<jupyter_start><jupyter_text>Sample Prep<jupyter_code>HGNC = pd.read_csv("../../../keiser/data/hgnc_mapping_ids/hgnc.txt", sep="\t")<jupyter_output><empty_output><jupyter_text>Define sample types<jupyter_code>from collections import OrderedDict 
hierarchy = OrderedDict({'Project ID':['TCGA-BRCA', 'TCGA-LUAD', 'TCGA-KIRC', 'TCGA-THCA', 'TCGA-PRAD', 'TCGA-SKCM'],
                         'Sample Type':['Primary Tumor', 'Solid Tissue Normal', 'Metastatic']})
from dutils import Experiment, train_test_split_case
exp = Experiment(samples_dir='../../data/TCGA/rna-seq_pan/meta/gdc_sample_sheet.2019-12-12.tsv',
                 hierarchy=hierarchy,
                 cases='Case ID',
                 min_samples=20)
exp.holdout(holdout=['TCGA-SKCM:Metastatic'])
exp.meta_dict<jupyter_output><empty_output><jupyter_text>Dataset Prep Need to ensure test set has health/diseased samples matched for each case<jupyter_code>exp.samples['meta'].value_counts()
exp.holdout_samples['meta'].value_counts()
from sklearn.preprocessing import StandardScaler

target = 'meta'

train, test = train_test_split_case(exp.samples, cases='Case ID')
# stratification is not quite perfect but close
# in order to preserve matched samples for each case together
# in train or test set
case_counts = exp.samples[target].value_counts()
train[target].value_counts()[case_counts.index.to_numpy()] / case_counts
test[target].value_counts()[case_counts.index.to_numpy()] / case_counts
import torch
from torch.optim import lr_scheduler
import torch.optim as optim
from torch.autograd import Variable
#torch.manual_seed(123)

from trainer import fit
import visualization as vis
import numpy as np
cuda = torch.cuda.is_available()
print("Cuda is available: {}".format(cuda))

from tcga_datasets import TCGA, SiameseTCGA
root_dir = "../../data/TCGA/rna-seq_pan/"

train_dataset = TCGA(root_dir, samples=train, meta_dict=exp.meta_dict, train=True, target=target, log=True)
train_dataset.featurize(HGNC=HGNC, features='all', proteins_only=True)

test_dataset = TCGA(root_dir, samples=test, meta_dict=exp.meta_dict, train=False, target=target, log=True)
test_dataset.featurize(HGNC=HGNC, features='all', proteins_only=True)
np.array_equal(train_dataset.data.columns, test_dataset.data.columns)
scaler = StandardScaler()
train_dataset.data = pd.DataFrame(scaler.fit_transform(train_dataset.data),
                                  index=train_dataset.data.index,
                                  columns=train_dataset.data.columns)
test_dataset.data = pd.DataFrame(scaler.transform(test_dataset.data),
                                 index=test_dataset.data.index,
                                 columns=test_dataset.data.columns)
train_dataset.data.shape
test_dataset.data.shape

# Still required in order to feed single samples for embedding extraction
batch_size = 1
kwargs = {'num_workers': 10, 'pin_memory': True} if cuda else {'num_workers': 10}
train_loader = torch.utils.data.DataLoader(train_dataset, batch_size=batch_size, shuffle=True, **kwargs)
test_loader = torch.utils.data.DataLoader(test_dataset, batch_size=batch_size, shuffle=False, **kwargs)<jupyter_output><empty_output><jupyter_text>Siamese Network<jupyter_code># Step 1 set up dataloader
siamese_train_dataset = SiameseTCGA(train_dataset) # Returns pairs of images and target same/different
siamese_test_dataset = SiameseTCGA(test_dataset)
batch_size = 8
kwargs = {'num_workers': 10, 'pin_memory': True} if cuda else {'num_workers': 10}
siamese_train_loader = torch.utils.data.DataLoader(siamese_train_dataset, batch_size=batch_size, shuffle=True, **kwargs)
siamese_test_loader = torch.utils.data.DataLoader(siamese_test_dataset, batch_size=batch_size, shuffle=False, **kwargs)

# Set up the network and training parameters
from tcga_networks import EmbeddingNet, SiameseNet
from losses import ContrastiveLoss
from metrics import AccumulatedAccuracyMetric

# Step 2
n_samples, n_features = siamese_train_dataset.data.shape
embedding_net = EmbeddingNet(n_features)
# Step 3
model = SiameseNet(embedding_net)
if cuda:
    model.cuda()
    
# Step 4
margin = 1.
loss_fn = ContrastiveLoss(margin)
lr = 1e-3
optimizer = optim.Adam(model.parameters(), lr=lr)
scheduler = lr_scheduler.StepLR(optimizer, 8, gamma=0.1, last_epoch=-1)
n_epochs = 1
# print training metrics every log_interval * batch_size
log_interval = round(len(siamese_train_dataset)/4/batch_size)
train_loss, val_loss = fit(siamese_train_loader, siamese_test_loader, model, loss_fn, optimizer, scheduler, 
    n_epochs, cuda, log_interval)
plt.plot(range(0, n_epochs), train_loss, 'rx-', label='train')
plt.plot(range(0, n_epochs), val_loss, 'bx-', label='validation')
plt.legend()
train_embeddings_cl, train_labels_cl = vis.extract_embeddings(train_loader, model)
vis.sns_plot_embeddings(train_embeddings_cl, train_labels_cl, siamese_train_dataset, hue='meta', style='Sample Type')
plt.title('PanCancer Train: Siamese')
val_embeddings_baseline, val_labels_baseline = vis.extract_embeddings(test_loader, model)
vis.sns_plot_embeddings(val_embeddings_baseline, val_labels_baseline, siamese_test_dataset, hue='meta', style='Sample Type')
plt.title('PanCancer Test: Siamese')<jupyter_output><empty_output><jupyter_text>Save Model<jupyter_code>save = False
if save:
    torch.save(model.state_dict(), "../../data/TCGA/models/220_seminar_results.pt")<jupyter_output><empty_output><jupyter_text>Kmeans Cluster and NMI Evaluation<jupyter_code>from sklearn.cluster import KMeans
from sklearn.metrics import normalized_mutual_info_score as NMI
from sklearn.metrics import adjusted_mutual_info_score as aNMI
n_clusters = len(siamese_test_dataset.meta_dict) - len(exp.holdout)
kmeans = KMeans(n_clusters=n_clusters, random_state=0)
siamese_clusters = kmeans.fit_predict(val_embeddings_baseline)
aNMI(siamese_clusters, val_labels_baseline)
# class Experiment():
#     """
#     Defines an experimental class hierarchy object.
#     """
#     def __init__(self, samples_dir, hierarchy, cases, min_samples):
#         assert isinstance(hierarchy, OrderedDict), "Argument of wrong type."
#         self.samples_dir = samples_dir
#         self.hierarchy = hierarchy
#         self.min_samples = min_samples
#         self.data = self.load_data(self.samples_dir)
#         self.cases = self.data[cases].unique()
        
#     def load_data(self, samples_dir):
#         samples = pd.read_csv(samples_dir, sep="\t")
#         samples['data'] = [val[1] for i,val in samples['File Name'].str.split(".").items()]
#         return samples
    
#     def process_data(self):
#         # downsample data
#         for key,val in self.hierarchy.items():
#             self.samples = self.data[self.data[key].isin(val)].copy()
#         # unique meta classes
#         self.samples['meta'] = self.samples[list(hierarchy.keys())].apply(lambda row: ':'.join(row.values.astype(str)), axis=1)
#         # filter meta classes
#         meta_counts = self.samples['meta'].value_counts()
#         meta_keep = meta_counts[meta_counts > self.min_samples].index
#         self.samples = self.samples[self.samples['meta'].isin(meta_keep)]
#         # generate class categories
#         self.samples['meta'] = self.samples['meta'].astype('category')
#         self.meta = {key:val for key,val in enumerate(self.samples['meta'].cat.categories.values)}<jupyter_output><empty_output><jupyter_text>Project Metastatic SKCM<jupyter_code># proj_hierarchy = OrderedDict({'Project ID':['TCGA-SKCM'],
#                               'Sample Type':['Metastatic']})
# Need to define Experiment method for use in adding classes for projection
# because categorical classes are predefined, need to add an additional class based
# on these numerical labels
# this will probably be helpful with plotting also...
# proj = Experiment(samples_dir='../../data/TCGA/rna-seq_pan/meta/gdc_sample_sheet.2019-12-12.tsv',
#                  hierarchy=proj_hierarchy,
#                  cases='Case ID',
#                  min_samples=20)
proj_dataset = TCGA(root_dir, samples=exp.holdout_samples, meta_dict=exp.meta_dict, train=False, target=target, log=True)
proj_dataset.featurize(HGNC=HGNC, features='all', proteins_only=True)
proj_dataset.data = pd.DataFrame(scaler.transform(proj_dataset.data),
                                 index=proj_dataset.data.index,
                                 columns=proj_dataset.data.columns)
# Still required in order to feed single samples for embedding extraction
batch_size = 1
kwargs = {'num_workers': 10, 'pin_memory': True} if cuda else {'num_workers': 10}
proj_loader = torch.utils.data.DataLoader(proj_dataset, batch_size=batch_size, shuffle=False, **kwargs)
proj_embeddings_baseline, proj_labels_baseline = vis.extract_embeddings(proj_loader, model)
vis.sns_plot_embeddings(proj_embeddings_baseline, proj_labels_baseline, proj_dataset, hue='meta', style='Sample Type')
plt.title('PanCancer Test: Siamese')
val_labels_baseline
proj_labels_baseline
def sns_plot_embeddings(embeddings, labels, classes, hue, style):
    sns.set_style("ticks")
    id_map = {key:val for key,val in enumerate(classes)}
    data = np.concatenate([embeddings, labels.reshape(-1,1)], axis=1)
    data = pd.DataFrame(data, columns=[0,1,'label'])
    data['label'] = data['label'].astype(int)
    data['meta'] = data['label'].map(id_map).astype('category')
    data[style] = [i[1] for i in data['meta'].str.split(":")]
    # plot
    sns.scatterplot(x=0, y=1, hue=hue, style=style, palette="Paired", data=data,
                s=100, alpha=0.8)<jupyter_output><empty_output><jupyter_text>SKCM Clinical Info<jupyter_code>skcm_clinical = pd.read_csv("../../data/TCGA/clinical/skcm/https:/api.gdc.cancer.gov/data/58cbbc07-5ec4-47c7-9295-11ccbf7693f4/nationwidechildrens.org_clinical_patient_skcm.txt", 
                            sep="\t")
# remove extra header
skcm_clinical = skcm_clinical.loc[2:, :].reset_index(drop=True)
skcm_clinical.shape
skcm_metas_dist = skcm_clinical[skcm_clinical['primary_location'] == 'Distant Metastasis']
skcm_metas_dist = skcm_clinical[skcm_clinical['metastatic_tumor_site'] == 'Lung']
skcm_metas_dist.shape
bool_skcm_metas_dist = projection_dataset.samples['Case ID'].isin(skcm_metas_dist['bcr_patient_barcode'])
bool_skcm_metas_dist.sum()
import seaborn as sns
sns.set(rc={'figure.figsize':(18.7,14.27)})
vis.sns_plot_embeddings(proj_embeddings_baseline, projection_dataset, hue='Meta ID', style='Sample Type')
plt.scatter(proj_embeddings_baseline[bool_skcm_metas_dist,0], 
            proj_embeddings_baseline[bool_skcm_metas_dist,1],
            marker="+", color='k', s=200)
plt.title('PanCancer Metastatic: Siamese')<jupyter_output><empty_output><jupyter_text>Project CCLE Cells- no normalization or scaling done here...<jupyter_code>ccle_meta = pd.read_csv("../../data/DepMap/sample_info_v2.csv", index_col=0)
ccle_meta.head()
msys = ['Lung Cancer', 'Breast Cancer', 'Kidney Cancer', 'Prostate Cancer', 'Skin Cancer', 'Thyroid Cancer']
ccle_meta_sub = ccle_meta[ccle_meta['disease'].isin(msys)]
ccle_meta_sub.shape
# get shared protein subset to training data
ccle_overlap_idx = ccle.columns.isin(train_dataset.data.columns)
ccle = ccle.loc[:, ccle_overlap_idx]
# reorder columns to match training columns
ccle = ccle[train_dataset.data.columns]
assert np.array_equal(ccle.columns.values, train_dataset.data.columns.values)
# reorder rows to match data and meta data
ccle_sub = ccle.loc[ccle_meta_sub.index,:]
assert np.array_equal(ccle_meta_sub.index, ccle_sub.index)
from tcga_datasets import CCLE
ccle_dataset = CCLE(ccle_sub, ccle_meta_sub, train=False, log=True)
ccle_dataset.data.head()
ccle_scaler = StandardScaler()
ccle_dataset.data = pd.DataFrame(ccle_scaler.fit_transform(ccle_dataset.data),
                                  index=ccle_dataset.data.index,
                                  columns=ccle_dataset.data.columns)
ccle_loader = torch.utils.data.DataLoader(ccle_dataset, batch_size=batch_size, shuffle=False, **kwargs)
ccle_embeddings_baseline, ccle_labels_baseline = vis.extract_embeddings(ccle_loader, model)
vis.sns_plot_embeddings(ccle_embeddings_baseline, ccle_dataset, hue="disease", style="primary_or_metastasis")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1212718757702736, 'fraction_numerical': 0.008791389368170241, 'mean_word_length': 5.073852295409182, 'pattern_counts': {'":': 1, '<': 37, '<?xml version=': 0, '>': 38, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '321171', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 4250, 'n_tokens_neox': 4049, 'n_words': 743, 'path': 'notebook/.ipynb_checkpoints/2020.03.03_clean_evalMetrics-proteinsSubset-noMatch-panCancer-Contrastive-checkpoint.ipynb', 'repo_name': 'wconnell/metrx', 'size': '12171'}