Datasets
Models
Downloads: 1
Switch to side-by-side view
--- a +++ b/utils.py @@ -0,0 +1,892 @@ +# Base / Native +import math +import os +import pickle +import re +import warnings +warnings.filterwarnings('ignore') + +# Numerical / Array +import lifelines +from lifelines.utils import concordance_index +from lifelines import CoxPHFitter +from lifelines.datasets import load_regression_dataset +from lifelines.utils import k_fold_cross_validation +from lifelines.statistics import logrank_test +from imblearn.over_sampling import RandomOverSampler +import matplotlib as mpl +import matplotlib.pyplot as plt +import matplotlib.font_manager as font_manager +import numpy as np +import pandas as pd +from PIL import Image +import pylab +import scipy +import seaborn as sns +from sklearn import preprocessing +from sklearn.model_selection import train_test_split, KFold +from sklearn.metrics import average_precision_score, auc, f1_score, roc_curve, roc_auc_score +from sklearn.preprocessing import LabelBinarizer + +from scipy import interp +mpl.rcParams['axes.linewidth'] = 3 #set the value globally + +# Torch +import torch +import torch.nn as nn +from torch.nn import init, Parameter +from torch.utils.data._utils.collate import * +from torch.utils.data.dataloader import default_collate +import torch_geometric +from torch_geometric.data import Batch + + + +################ +# Regularization +################ +def regularize_weights(model, reg_type=None): + l1_reg = None + + for W in model.parameters(): + if l1_reg is None: + l1_reg = torch.abs(W).sum() + else: + l1_reg = l1_reg + torch.abs(W).sum() # torch.abs(W).sum() is equivalent to W.norm(1) + return l1_reg + + +def regularize_path_weights(model, reg_type=None): + l1_reg = None + + for W in model.module.classifier.parameters(): + if l1_reg is None: + l1_reg = torch.abs(W).sum() + else: + l1_reg = l1_reg + torch.abs(W).sum() # torch.abs(W).sum() is equivalent to W.norm(1) + + for W in model.module.linear.parameters(): + if l1_reg is None: + l1_reg = torch.abs(W).sum() + else: + l1_reg = l1_reg + torch.abs(W).sum() # torch.abs(W).sum() is equivalent to W.norm(1) + + return l1_reg + + +def regularize_MM_weights(model, reg_type=None): + l1_reg = None + + if model.module.__hasattr__('omic_net'): + for W in model.module.omic_net.parameters(): + if l1_reg is None: + l1_reg = torch.abs(W).sum() + else: + l1_reg = l1_reg + torch.abs(W).sum() # torch.abs(W).sum() is equivalent to W.norm(1) + + if model.module.__hasattr__('linear_h_path'): + for W in model.module.linear_h_path.parameters(): + if l1_reg is None: + l1_reg = torch.abs(W).sum() + else: + l1_reg = l1_reg + torch.abs(W).sum() # torch.abs(W).sum() is equivalent to W.norm(1) + + if model.module.__hasattr__('linear_h_omic'): + for W in model.module.linear_h_omic.parameters(): + if l1_reg is None: + l1_reg = torch.abs(W).sum() + else: + l1_reg = l1_reg + torch.abs(W).sum() # torch.abs(W).sum() is equivalent to W.norm(1) + + if model.module.__hasattr__('linear_h_grph'): + for W in model.module.linear_h_grph.parameters(): + if l1_reg is None: + l1_reg = torch.abs(W).sum() + else: + l1_reg = l1_reg + torch.abs(W).sum() # torch.abs(W).sum() is equivalent to W.norm(1) + + if model.module.__hasattr__('linear_z_path'): + for W in model.module.linear_z_path.parameters(): + if l1_reg is None: + l1_reg = torch.abs(W).sum() + else: + l1_reg = l1_reg + torch.abs(W).sum() # torch.abs(W).sum() is equivalent to W.norm(1) + + if model.module.__hasattr__('linear_z_omic'): + for W in model.module.linear_z_omic.parameters(): + if l1_reg is None: + l1_reg = torch.abs(W).sum() + else: + l1_reg = l1_reg + torch.abs(W).sum() # torch.abs(W).sum() is equivalent to W.norm(1) + + if model.module.__hasattr__('linear_z_grph'): + for W in model.module.linear_z_grph.parameters(): + if l1_reg is None: + l1_reg = torch.abs(W).sum() + else: + l1_reg = l1_reg + torch.abs(W).sum() # torch.abs(W).sum() is equivalent to W.norm(1) + + if model.module.__hasattr__('linear_o_path'): + for W in model.module.linear_o_path.parameters(): + if l1_reg is None: + l1_reg = torch.abs(W).sum() + else: + l1_reg = l1_reg + torch.abs(W).sum() # torch.abs(W).sum() is equivalent to W.norm(1) + + if model.module.__hasattr__('linear_o_omic'): + for W in model.module.linear_o_omic.parameters(): + if l1_reg is None: + l1_reg = torch.abs(W).sum() + else: + l1_reg = l1_reg + torch.abs(W).sum() # torch.abs(W).sum() is equivalent to W.norm(1) + + if model.module.__hasattr__('linear_o_grph'): + for W in model.module.linear_o_grph.parameters(): + if l1_reg is None: + l1_reg = torch.abs(W).sum() + else: + l1_reg = l1_reg + torch.abs(W).sum() # torch.abs(W).sum() is equivalent to W.norm(1) + + if model.module.__hasattr__('encoder1'): + for W in model.module.encoder1.parameters(): + if l1_reg is None: + l1_reg = torch.abs(W).sum() + else: + l1_reg = l1_reg + torch.abs(W).sum() # torch.abs(W).sum() is equivalent to W.norm(1) + + if model.module.__hasattr__('encoder2'): + for W in model.module.encoder2.parameters(): + if l1_reg is None: + l1_reg = torch.abs(W).sum() + else: + l1_reg = l1_reg + torch.abs(W).sum() # torch.abs(W).sum() is equivalent to W.norm(1) + + if model.module.__hasattr__('classifier'): + for W in model.module.classifier.parameters(): + if l1_reg is None: + l1_reg = torch.abs(W).sum() + else: + l1_reg = l1_reg + torch.abs(W).sum() # torch.abs(W).sum() is equivalent to W.norm(1) + + return l1_reg + + +def regularize_MM_omic(model, reg_type=None): + l1_reg = None + + if model.module.__hasattr__('omic_net'): + for W in model.module.omic_net.parameters(): + if l1_reg is None: + l1_reg = torch.abs(W).sum() + else: + l1_reg = l1_reg + torch.abs(W).sum() # torch.abs(W).sum() is equivalent to W.norm(1) + + return l1_reg + + + +################ +# Network Initialization +################ +def init_weights(net, init_type='orthogonal', init_gain=0.02): + """Initialize network weights. + + Parameters: + net (network) -- network to be initialized + init_type (str) -- the name of an initialization method: normal | xavier | kaiming | orthogonal + init_gain (float) -- scaling factor for normal, xavier and orthogonal. + + We use 'normal' in the original pix2pix and CycleGAN paper. But xavier and kaiming might + work better for some applications. Feel free to try yourself. + """ + def init_func(m): # define the initialization function + classname = m.__class__.__name__ + if hasattr(m, 'weight') and (classname.find('Conv') != -1 or classname.find('Linear') != -1): + if init_type == 'normal': + init.normal_(m.weight.data, 0.0, init_gain) + elif init_type == 'xavier': + init.xavier_normal_(m.weight.data, gain=init_gain) + elif init_type == 'kaiming': + init.kaiming_normal_(m.weight.data, a=0, mode='fan_in') + elif init_type == 'orthogonal': + init.orthogonal_(m.weight.data, gain=init_gain) + else: + raise NotImplementedError('initialization method [%s] is not implemented' % init_type) + if hasattr(m, 'bias') and m.bias is not None: + init.constant_(m.bias.data, 0.0) + elif classname.find('BatchNorm2d') != -1: # BatchNorm Layer's weight is not a matrix; only normal distribution applies. + init.normal_(m.weight.data, 1.0, init_gain) + init.constant_(m.bias.data, 0.0) + + print('initialize network with %s' % init_type) + net.apply(init_func) # apply the initialization function <init_func> + + +def init_max_weights(module): + for m in module.modules(): + if type(m) == nn.Linear: + stdv = 1. / math.sqrt(m.weight.size(1)) + m.weight.data.normal_(0, stdv) + m.bias.data.zero_() + + +def init_net(net, init_type='normal', init_gain=0.02, gpu_ids=[]): + """Initialize a network: 1. register CPU/GPU device (with multi-GPU support); 2. initialize the network weights + Parameters: + net (network) -- the network to be initialized + init_type (str) -- the name of an initialization method: normal | xavier | kaiming | orthogonal + gain (float) -- scaling factor for normal, xavier and orthogonal. + gpu_ids (int list) -- which GPUs the network runs on: e.g., 0,1,2 + + Return an initialized network. + """ + if len(gpu_ids) > 0: + assert(torch.cuda.is_available()) + net.to(gpu_ids[0]) + net = torch.nn.DataParallel(net, gpu_ids) # multi-GPUs + + if init_type != 'max' and init_type != 'none': + print("Init Type:", init_type) + init_weights(net, init_type, init_gain=init_gain) + elif init_type == 'none': + print("Init Type: Not initializing networks.") + elif init_type == 'max': + print("Init Type: Self-Normalizing Weights") + return net + + + +################ +# Freeze / Unfreeze +################ +def unfreeze_unimodal(opt, model, epoch): + if opt.mode == 'graphomic': + if epoch == 5: + dfs_unfreeze(model.module.omic_net) + print("Unfreezing Omic") + if epoch == 5: + dfs_unfreeze(model.module.grph_net) + print("Unfreezing Graph") + elif opt.mode == 'pathomic': + if epoch == 5: + dfs_unfreeze(model.module.omic_net) + print("Unfreezing Omic") + elif opt.mode == 'pathgraph': + if epoch == 5: + dfs_unfreeze(model.module.grph_net) + print("Unfreezing Graph") + elif opt.mode == "pathgraphomic": + if epoch == 5: + dfs_unfreeze(model.module.omic_net) + print("Unfreezing Omic") + if epoch == 5: + dfs_unfreeze(model.module.grph_net) + print("Unfreezing Graph") + elif opt.mode == "omicomic": + if epoch == 5: + dfs_unfreeze(model.module.omic_net) + print("Unfreezing Omic") + elif opt.mode == "graphgraph": + if epoch == 5: + dfs_unfreeze(model.module.grph_net) + print("Unfreezing Graph") + + +def dfs_freeze(model): + for name, child in model.named_children(): + for param in child.parameters(): + param.requires_grad = False + dfs_freeze(child) + + +def dfs_unfreeze(model): + for name, child in model.named_children(): + for param in child.parameters(): + param.requires_grad = True + dfs_unfreeze(child) + + +def print_if_frozen(module): + for idx, child in enumerate(module.children()): + for param in child.parameters(): + if param.requires_grad == True: + print("Learnable!!! %d:" % idx, child) + else: + print("Still Frozen %d:" % idx, child) + + +def unfreeze_vgg_features(model, epoch): + epoch_schedule = {30:45} + unfreeze_index = epoch_schedule[epoch] + for idx, child in enumerate(model.features.children()): + if idx > unfreeze_index: + print("Unfreezing %d:" %idx, child) + for param in child.parameters(): + param.requires_grad = True + else: + print("Still Frozen %d:" %idx, child) + continue + + + +################ +# Collate Utils +################ +def mixed_collate(batch): + elem = batch[0] + elem_type = type(elem) + transposed = zip(*batch) + return [Batch.from_data_list(samples, []) if type(samples[0]) is torch_geometric.data.data.Data else default_collate(samples) for samples in transposed] + + + +################ +# Survival Utils +################ +def CoxLoss(survtime, censor, hazard_pred, device): + # This calculation credit to Travers Ching https://github.com/traversc/cox-nnet + # Cox-nnet: An artificial neural network method for prognosis prediction of high-throughput omics data + current_batch_len = len(survtime) + R_mat = np.zeros([current_batch_len, current_batch_len], dtype=int) + for i in range(current_batch_len): + for j in range(current_batch_len): + R_mat[i,j] = survtime[j] >= survtime[i] + + R_mat = torch.FloatTensor(R_mat).to(device) + theta = hazard_pred.reshape(-1) + exp_theta = torch.exp(theta) + loss_cox = -torch.mean((theta - torch.log(torch.sum(exp_theta*R_mat, dim=1))) * censor) + return loss_cox + + +def accuracy(output, labels): + preds = output.max(1)[1].type_as(labels) + correct = preds.eq(labels).double() + correct = correct.sum() + return correct / len(labels) + + +def accuracy_cox(hazardsdata, labels): + # This accuracy is based on estimated survival events against true survival events + median = np.median(hazardsdata) + hazards_dichotomize = np.zeros([len(hazardsdata)], dtype=int) + hazards_dichotomize[hazardsdata > median] = 1 + correct = np.sum(hazards_dichotomize == labels) + return correct / len(labels) + + +def cox_log_rank(hazardsdata, labels, survtime_all): + median = np.median(hazardsdata) + hazards_dichotomize = np.zeros([len(hazardsdata)], dtype=int) + hazards_dichotomize[hazardsdata > median] = 1 + idx = hazards_dichotomize == 0 + T1 = survtime_all[idx] + T2 = survtime_all[~idx] + E1 = labels[idx] + E2 = labels[~idx] + results = logrank_test(T1, T2, event_observed_A=E1, event_observed_B=E2) + pvalue_pred = results.p_value + return(pvalue_pred) + + +def CIndex(hazards, labels, survtime_all): + concord = 0. + total = 0. + N_test = labels.shape[0] + for i in range(N_test): + if labels[i] == 1: + for j in range(N_test): + if survtime_all[j] > survtime_all[i]: + total += 1 + if hazards[j] < hazards[i]: concord += 1 + elif hazards[j] < hazards[i]: concord += 0.5 + + return(concord/total) + + +def CIndex_lifeline(hazards, labels, survtime_all): + return(concordance_index(survtime_all, -hazards, labels)) + + + +################ +# Data Utils +################ +def addHistomolecularSubtype(data): + """ + Molecular Subtype: IDHwt == 0, IDHmut-non-codel == 1, IDHmut-codel == 2 + Histology Subtype: astrocytoma == 0, oligoastrocytoma == 1, oligodendroglioma == 2, glioblastoma == 3 + """ + subtyped_data = data.copy() + subtyped_data.insert(loc=0, column='Histomolecular subtype', value=np.ones(len(data))) + idhwt_ATC = np.logical_and(data['Molecular subtype'] == 0, np.logical_or(data['Histology'] == 0, data['Histology'] == 3)) + subtyped_data.loc[idhwt_ATC, 'Histomolecular subtype'] = 'idhwt_ATC' + + idhmut_ATC = np.logical_and(data['Molecular subtype'] == 1, np.logical_or(data['Histology'] == 0, data['Histology'] == 3)) + subtyped_data.loc[idhmut_ATC, 'Histomolecular subtype'] = 'idhmut_ATC' + + ODG = np.logical_and(data['Molecular subtype'] == 2, data['Histology'] == 2) + subtyped_data.loc[ODG, 'Histomolecular subtype'] = 'ODG' + return subtyped_data + + +def changeHistomolecularSubtype(data): + """ + Molecular Subtype: IDHwt == 0, IDHmut-non-codel == 1, IDHmut-codel == 2 + Histology Subtype: astrocytoma == 0, oligoastrocytoma == 1, oligodendroglioma == 2, glioblastoma == 3 + """ + data = data.drop(['Histomolecular subtype'], axis=1) + subtyped_data = data.copy() + subtyped_data.insert(loc=0, column='Histomolecular subtype', value=np.ones(len(data))) + idhwt_ATC = np.logical_and(data['Molecular subtype'] == 0, np.logical_or(data['Histology'] == 0, data['Histology'] == 3)) + subtyped_data.loc[idhwt_ATC, 'Histomolecular subtype'] = 'idhwt_ATC' + + idhmut_ATC = np.logical_and(data['Molecular subtype'] == 1, np.logical_or(data['Histology'] == 0, data['Histology'] == 3)) + subtyped_data.loc[idhmut_ATC, 'Histomolecular subtype'] = 'idhmut_ATC' + + ODG = np.logical_and(data['Molecular subtype'] == 2, data['Histology'] == 2) + subtyped_data.loc[ODG, 'Histomolecular subtype'] = 'ODG' + return subtyped_data + + +def getCleanAllDataset(dataroot='./data/TCGA_GBMLGG/', ignore_missing_moltype=False, ignore_missing_histype=False, use_rnaseq=False): + ### 1. Joining all_datasets.csv with grade data. Looks at columns with misisng samples + metadata = ['Histology', 'Grade', 'Molecular subtype', 'TCGA ID', 'censored', 'Survival months'] + all_dataset = pd.read_csv(os.path.join(dataroot, 'all_dataset.csv')).drop('indexes', axis=1) + all_dataset.index = all_dataset['TCGA ID'] + + all_grade = pd.read_csv(os.path.join(dataroot, 'grade_data.csv')) + all_grade['Histology'] = all_grade['Histology'].str.replace('astrocytoma (glioblastoma)', 'glioblastoma', regex=False) + all_grade.index = all_grade['TCGA ID'] + assert pd.Series(all_dataset.index).equals(pd.Series(sorted(all_grade.index))) + + all_dataset = all_dataset.join(all_grade[['Histology', 'Grade', 'Molecular subtype']], how='inner') + cols = all_dataset.columns.tolist() + cols = cols[-3:] + cols[:-3] + all_dataset = all_dataset[cols] + + if use_rnaseq: + gbm = pd.read_csv(os.path.join(dataroot, 'mRNA_Expression_z-Scores_RNA_Seq_RSEM.txt'), sep='\t', skiprows=1, index_col=0) + lgg = pd.read_csv(os.path.join(dataroot, 'mRNA_Expression_Zscores_RSEM.txt'), sep='\t', skiprows=1, index_col=0) + gbm = gbm[gbm.columns[~gbm.isnull().all()]] + lgg = lgg[lgg.columns[~lgg.isnull().all()]] + glioma_RNAseq = gbm.join(lgg, how='inner').T + glioma_RNAseq = glioma_RNAseq.dropna(axis=1) + glioma_RNAseq.columns = [gene+'_rnaseq' for gene in glioma_RNAseq.columns] + glioma_RNAseq.index = [patname[:12] for patname in glioma_RNAseq.index] + glioma_RNAseq = glioma_RNAseq.iloc[~glioma_RNAseq.index.duplicated()] + glioma_RNAseq.index.name = 'TCGA ID' + all_dataset = all_dataset.join(glioma_RNAseq, how='inner') + + pat_missing_moltype = all_dataset[all_dataset['Molecular subtype'].isna()].index + pat_missing_idh = all_dataset[all_dataset['idh mutation'].isna()].index + pat_missing_1p19q = all_dataset[all_dataset['codeletion'].isna()].index + print("# Missing Molecular Subtype:", len(pat_missing_moltype)) + print("# Missing IDH Mutation:", len(pat_missing_idh)) + print("# Missing 1p19q Codeletion:", len(pat_missing_1p19q)) + assert pat_missing_moltype.equals(pat_missing_idh) + assert pat_missing_moltype.equals(pat_missing_1p19q) + pat_missing_grade = all_dataset[all_dataset['Grade'].isna()].index + pat_missing_histype = all_dataset[all_dataset['Histology'].isna()].index + print("# Missing Histological Subtype:", len(pat_missing_histype)) + print("# Missing Grade:", len(pat_missing_grade)) + assert pat_missing_histype.equals(pat_missing_grade) + + ### 2. Impute Missing Genomic Data: Removes patients with missing molecular subtype / idh mutation / 1p19q. Else imputes with median value of each column. Fills missing Molecular subtype with "Missing" + if ignore_missing_moltype: + all_dataset = all_dataset[all_dataset['Molecular subtype'].isna() == False] + for col in all_dataset.drop(metadata, axis=1).columns: + all_dataset['Molecular subtype'] = all_dataset['Molecular subtype'].fillna('Missing') + all_dataset[col] = all_dataset[col].fillna(all_dataset[col].median()) + + ### 3. Impute Missing Histological Data: Removes patients with missing histological subtype / grade. Else imputes with "missing" / grade -1 + if ignore_missing_histype: + all_dataset = all_dataset[all_dataset['Histology'].isna() == False] + else: + all_dataset['Grade'] = all_dataset['Grade'].fillna(1) + all_dataset['Histology'] = all_dataset['Histology'].fillna('Missing') + all_dataset['Grade'] = all_dataset['Grade'] - 2 + + ### 4. Adds Histomolecular subtype + ms2int = {'Missing':-1, 'IDHwt':0, 'IDHmut-non-codel':1, 'IDHmut-codel':2} + all_dataset[['Molecular subtype']] = all_dataset[['Molecular subtype']].applymap(lambda s: ms2int.get(s) if s in ms2int else s) + hs2int = {'Missing':-1, 'astrocytoma':0, 'oligoastrocytoma':1, 'oligodendroglioma':2, 'glioblastoma':3} + all_dataset[['Histology']] = all_dataset[['Histology']].applymap(lambda s: hs2int.get(s) if s in hs2int else s) + all_dataset = addHistomolecularSubtype(all_dataset) + metadata.extend(['Histomolecular subtype']) + all_dataset['censored'] = 1 - all_dataset['censored'] + return metadata, all_dataset + + + +################ +# Analysis Utils +################ +def count_parameters(model): + return sum(p.numel() for p in model.parameters() if p.requires_grad) + + +def hazard2grade(hazard, p): + if hazard < p[0]: + return 0 + elif hazard < p[1]: + return 1 + return 2 + + +def p(n): + def percentile_(x): + return np.percentile(x, n) + percentile_.__name__ = 'p%s' % n + return percentile_ + + +def natural_sort(l): + convert = lambda text: int(text) if text.isdigit() else text.lower() + alphanum_key = lambda key: [ convert(c) for c in re.split('([0-9]+)', key) ] + return sorted(l, key = alphanum_key) + + +def CI_pm(data, confidence=0.95): + a = 1.0 * np.array(data) + n = len(a) + m, se = np.mean(a), scipy.stats.sem(a) + h = se * scipy.stats.t.ppf((1 + confidence) / 2., n-1) + return str("{0:.4f} ± ".format(m) + "{0:.3f}".format(h)) + + +def CI_interval(data, confidence=0.95): + a = 1.0 * np.array(data) + n = len(a) + m, se = np.mean(a), scipy.stats.sem(a) + h = se * scipy.stats.t.ppf((1 + confidence) / 2., n-1) + return str("{0:.3f}, ".format(m-h) + "{0:.3f}".format(m+h)) + + +def poolSurvTestPD(ckpt_name='./checkpoints/TCGA_GBMLGG/surv_15_rnaseq/', model='pathgraphomic_fusion', split='test', zscore=False, agg_type='Hazard_mean'): + all_dataset_regstrd_pooled = [] + ignore_missing_moltype = 1 if 'omic' in model else 0 + ignore_missing_histype = 1 if 'grad' in ckpt_name else 0 + use_patch, roi_dir, use_vgg_features = ('_patch_', 'all_st_patches_512', 1) if ((('path' in model) or ('graph' in model)) and ('cox' not in model)) else ('_', 'all_st', 0) + use_rnaseq = '_rnaseq' if ('rnaseq' in ckpt_name and 'path' != model and 'pathpath' not in model and 'graph' != model and 'graphgraph' not in model) else '' + + for k in range(1,16): + pred = pickle.load(open(ckpt_name+'/%s/%s_%d%spred_%s.pkl' % (model, model, k, use_patch, split), 'rb')) + + if 'cox' not in model: + surv_all = pd.DataFrame(np.stack(np.delete(np.array(pred), 3))).T + surv_all.columns = ['Hazard', 'Survival months', 'censored', 'Grade'] + data_cv = pickle.load(open('./data/TCGA_GBMLGG/splits/gbmlgg15cv_%s_%d_%d_%d%s.pkl' % (roi_dir, ignore_missing_moltype, ignore_missing_histype, use_vgg_features, use_rnaseq), 'rb')) + data_cv_splits = data_cv['cv_splits'] + data_cv_split_k = data_cv_splits[k] + assert np.all(data_cv_split_k[split]['t'] == pred[1]) # Data is correctly registered + all_dataset = data_cv['data_pd'].drop('TCGA ID', axis=1) + all_dataset_regstrd = all_dataset.loc[data_cv_split_k[split]['x_patname']] # Subset of "all_datasets" (metadata) that is registered with "pred" (predictions) + assert np.all(np.array(all_dataset_regstrd['Survival months']) == pred[1]) + assert np.all(np.array(all_dataset_regstrd['censored']) == pred[2]) + assert np.all(np.array(all_dataset_regstrd['Grade']) == pred[4]) + all_dataset_regstrd.insert(loc=0, column='Hazard', value = np.array(surv_all['Hazard'])) + all_dataset_regstrd.index.name = 'TCGA ID' + hazard_agg = all_dataset_regstrd.groupby('TCGA ID').agg({'Hazard': ['mean', 'median', max, p(0.25), p(0.75)]}) + hazard_agg.columns = ["_".join(x) for x in hazard_agg.columns.ravel()] + hazard_agg = hazard_agg[[agg_type]] + hazard_agg.columns = ['Hazard'] + pred = hazard_agg.join(all_dataset, how='inner') + + if zscore: pred['Hazard'] = scipy.stats.zscore(np.array(pred['Hazard'])) + all_dataset_regstrd_pooled.append(pred) + + all_dataset_regstrd_pooled = pd.concat(all_dataset_regstrd_pooled) + all_dataset_regstrd_pooled = changeHistomolecularSubtype(all_dataset_regstrd_pooled) + return all_dataset_regstrd_pooled + + +def getAggHazardCV(ckpt_name='./checkpoints/TCGA_GBMLGG/surv_15_rnaseq/', model='pathgraphomic_fusion', split='test', agg_type='Hazard_mean'): + result = [] + + ignore_missing_moltype = 1 if 'omic' in model else 0 + ignore_missing_histype = 1 if 'grad' in ckpt_name else 0 + use_patch, roi_dir, use_vgg_features = ('_patch_', 'all_st_patches_512', 1) if (('path' in model) or ('graph' in model)) else ('_', 'all_st', 0) + use_rnaseq = '_rnaseq' if ('rnaseq' in ckpt_name and 'path' != model and 'pathpath' not in model and 'graph' != model and 'graphgraph' not in model) else '' + + for k in range(1,16): + pred = pickle.load(open(ckpt_name+'/%s/%s_%d%spred_%s.pkl' % (model, model, k, use_patch, split), 'rb')) + surv_all = pd.DataFrame(np.stack(np.delete(np.array(pred), 3))).T + surv_all.columns = ['Hazard', 'Survival months', 'censored', 'Grade'] + data_cv = pickle.load(open('./data/TCGA_GBMLGG/splits/gbmlgg15cv_%s_%d_%d_%d%s.pkl' % (roi_dir, ignore_missing_moltype, ignore_missing_histype, use_vgg_features, use_rnaseq), 'rb')) + data_cv_splits = data_cv['cv_splits'] + data_cv_split_k = data_cv_splits[k] + assert np.all(data_cv_split_k[split]['t'] == pred[1]) # Data is correctly registered + all_dataset = data_cv['data_pd'].drop('TCGA ID', axis=1) + all_dataset_regstrd = all_dataset.loc[data_cv_split_k[split]['x_patname']] # Subset of "all_datasets" (metadata) that is registered with "pred" (predictions) + assert np.all(np.array(all_dataset_regstrd['Survival months']) == pred[1]) + assert np.all(np.array(all_dataset_regstrd['censored']) == pred[2]) + assert np.all(np.array(all_dataset_regstrd['Grade']) == pred[4]) + all_dataset_regstrd.insert(loc=0, column='Hazard', value = np.array(surv_all['Hazard'])) + all_dataset_regstrd.index.name = 'TCGA ID' + hazard_agg = all_dataset_regstrd.groupby('TCGA ID').agg({'Hazard': ['mean', max, p(0.75)]}) + hazard_agg.columns = ["_".join(x) for x in hazard_agg.columns.ravel()] + hazard_agg = hazard_agg[[agg_type]] + hazard_agg.columns = ['Hazard'] + all_dataset_hazard = hazard_agg.join(all_dataset, how='inner') + cin = CIndex_lifeline(all_dataset_hazard['Hazard'], all_dataset_hazard['censored'], all_dataset_hazard['Survival months']) + result.append(cin) + + return result + + +def calcGradMetrics(ckpt_name='./checkpoints/grad_15/', model='pathgraphomic_fusion', split='test', avg='micro'): + auc_all = [] + ap_all = [] + f1_all = [] + f1_gradeIV_all = [] + + ignore_missing_moltype = 1 if 'omic' in model else 0 + ignore_missing_histype = 1 if 'grad' in ckpt_name else 0 + use_patch, roi_dir, use_vgg_features = ('_patch_', 'all_st_patches_512', 1) if (('path' in model) or ('graph' in model)) else ('_', 'all_st', 0) + + for k in range(1,16): + pred = pickle.load(open(ckpt_name+'/%s/%s_%d%spred_%s.pkl' % (model, model, k, use_patch, split), 'rb')) + grade_pred, grade = np.array(pred[3]), np.array(pred[4]) + enc = LabelBinarizer() + enc.fit(grade) + grade_oh = enc.transform(grade) + rocauc = roc_auc_score(grade_oh, grade_pred, avg) + ap = average_precision_score(grade_oh, grade_pred, average=avg) + f1 = f1_score(grade_pred.argmax(axis=1), grade, average=avg) + f1_gradeIV = f1_score(grade_pred.argmax(axis=1), grade, average=None)[2] + + auc_all.append(rocauc) + ap_all.append(ap) + f1_all.append(f1) + f1_gradeIV_all.append(f1_gradeIV) + + return np.array([CI_pm(auc_all), CI_pm(ap_all), CI_pm(f1_all), CI_pm(f1_gradeIV_all)]) + + + +################ +# Plot Utils +################ +def makeKaplanMeierPlot(ckpt_name='./checkpoints/surv_15_rnaseq/', model='omic', split='test', zscore=False, agg_type='Hazard_mean'): + def hazard2KMCurve(data, subtype): + p = np.percentile(data['Hazard'], [33, 66]) + if p[0] == p[1]: p[0] = 2.99997 + data.insert(0, 'grade_pred', [hazard2grade(hazard, p) for hazard in data['Hazard']]) + kmf_pred = lifelines.KaplanMeierFitter() + kmf_gt = lifelines.KaplanMeierFitter() + + def get_name(model): + mode2name = {'pathgraphomic':'Pathomic F.', 'pathomic':'Pathomic F.', 'graphomic':'Pathomic F.', 'path':'Histology CNN', 'graph':'Histology GCN', 'omic':'Genomic SNN'} + for mode in mode2name.keys(): + if mode in model: return mode2name[mode] + return 'N/A' + + fig = plt.figure(figsize=(10, 10), dpi=600) + ax = plt.subplot() + censor_style = {'ms': 20, 'marker': '+'} + + temp = data[data['Grade']==0] + kmf_gt.fit(temp['Survival months']/365, temp['censored'], label="Grade II") + kmf_gt.plot(ax=ax, show_censors=True, ci_show=False, c='g', linewidth=3, ls='--', markerfacecolor='black', censor_styles=censor_style) + temp = data[data['grade_pred']==0] + kmf_pred.fit(temp['Survival months']/365, temp['censored'], label="%s (Low)" % get_name(model)) + kmf_pred.plot(ax=ax, show_censors=True, ci_show=False, c='g', linewidth=4, ls='-', markerfacecolor='black', censor_styles=censor_style) + + temp = data[data['Grade']==1] + kmf_gt.fit(temp['Survival months']/365, temp['censored'], label="Grade III") + kmf_gt.plot(ax=ax, show_censors=True, ci_show=False, c='b', linewidth=3, ls='--', censor_styles=censor_style) + temp = data[data['grade_pred']==1] + kmf_pred.fit(temp['Survival months']/365, temp['censored'], label="%s (Mid)" % get_name(model)) + kmf_pred.plot(ax=ax, show_censors=True, ci_show=False, c='b', linewidth=4, ls='-', censor_styles=censor_style) + + if subtype != 'ODG': + temp = data[data['Grade']==2] + kmf_gt.fit(temp['Survival months']/365, temp['censored'], label="Grade IV") + kmf_gt.plot(ax=ax, show_censors=True, ci_show=False, c='r', linewidth=3, ls='--', censor_styles=censor_style) + temp = data[data['grade_pred']==2] + kmf_pred.fit(temp['Survival months']/365, temp['censored'], label="%s (High)" % get_name(model)) + kmf_pred.plot(ax=ax, show_censors=True, ci_show=False, c='r', linewidth=4, ls='-', censor_styles=censor_style) + + ax.set_xlabel('') + ax.set_ylim(0, 1) + ax.set_yticks(np.arange(0, 1.001, 0.5)) + + ax.tick_params(axis='both', which='major', labelsize=40) + plt.legend(fontsize=32, prop=font_manager.FontProperties(family='Arial', style='normal', size=32)) + if subtype != 'idhwt_ATC': ax.get_legend().remove() + return fig + + data = poolSurvTestPD(ckpt_name, model, split, zscore, agg_type) + for subtype in ['idhwt_ATC', 'idhmut_ATC', 'ODG']: + fig = hazard2KMCurve(data[data['Histomolecular subtype'] == subtype], subtype) + fig.savefig(ckpt_name+'/%s_KM_%s.png' % (model, subtype)) + + fig = hazard2KMCurve(data, 'all') + fig.savefig(ckpt_name+'/%s_KM_%s.png' % (model, 'all')) + + +def makeHazardSwarmPlot(ckpt_name='./checkpoints/surv_15_rnaseq/', model='path', split='test', zscore=True, agg_type='Hazard_mean'): + mpl.rcParams['font.family'] = "arial" + data = poolSurvTestPD(ckpt_name=ckpt_name, model=model, split=split, zscore=zscore, agg_type=agg_type) + data = data[data['Grade'] != -1] + data = data[data['Histomolecular subtype'] != -1] + data['Grade'] = data['Grade'].astype(int).astype(str) + data['Grade'] = data['Grade'].str.replace('0', 'Grade II', regex=False) + data['Grade'] = data['Grade'].str.replace('1', 'Grade III', regex=False) + data['Grade'] = data['Grade'].str.replace('2', 'Grade IV', regex=False) + data['Histomolecular subtype'] = data['Histomolecular subtype'].str.replace('idhwt_ATC', 'IDH-wt \n astryocytoma', regex=False) + data['Histomolecular subtype'] = data['Histomolecular subtype'].str.replace('idhmut_ATC', 'IDH-mut \n astrocytoma', regex=False) + data['Histomolecular subtype'] = data['Histomolecular subtype'].str.replace('ODG', 'Oligodendroglioma', regex=False) + + fig, ax = plt.subplots(dpi=600) + ax.set_ylim([-2, 2.5]) # plt.ylim(-2, 2) + ax.spines['right'].set_visible(False) + ax.spines['top'].set_visible(False) + ax.set_yticks(np.arange(-2, 2.001, 1)) + + sns.swarmplot(x = 'Histomolecular subtype', y='Hazard', data=data, hue='Grade', + palette={"Grade II":"#AFD275" , "Grade III":"#7395AE", "Grade IV":"#E7717D"}, + size = 4, alpha = 0.9, ax=ax) + + ax.set_xlabel('') # ax.set_xlabel('Histomolecular subtype', size=16) + ax.set_ylabel('') # ax.set_ylabel('Hazard (Z-Score)', size=16) + ax.tick_params(axis='y', which='both', labelsize=20) + ax.tick_params(axis='x', which='both', labelsize=15) + ax.tick_params(axis='x', which='both', labelbottom='off') # doesn't work?? + ax.legend(prop={'size': 8}) + fig.savefig(ckpt_name+'/%s_HSP.png' % (model)) + + +def makeHazardBoxPlot(ckpt_name='./checkpoints/surv_15_rnaseq/', model='omic', split='test', zscore=True, agg_type='Hazard_mean'): + mpl.rcParams['font.family'] = "arial" + data = poolSurvTestPD(ckpt_name, model, split, zscore, 'Hazard_mean') + data['Grade'] = data['Grade'].astype(int).astype(str) + data['Grade'] = data['Grade'].str.replace('0', 'II', regex=False) + data['Grade'] = data['Grade'].str.replace('1', 'III', regex=False) + data['Grade'] = data['Grade'].str.replace('2', 'IV', regex=False) + + fig, axes = plt.subplots(nrows=1, ncols=3, gridspec_kw={'width_ratios': [3, 3, 2]}, dpi=600) + plt.subplots_adjust(wspace=0, hspace=0) + plt.ylim(-2, 2) + plt.yticks(np.arange(-2, 2.001, 1)) + #color_dict = {0: '#CF9498', 1: '#8CC7C8', 2: '#AAA0C6'} + #color_dict = {0: '#F76C6C', 1: '#A8D0E6', 2: '#F8E9A1'} + color_dict = ['#F76C6C', '#A8D0E6', '#F8E9A1'] + subtypes = ['idhwt_ATC', 'idhmut_ATC', 'ODG'] + + for i in range(len(subtypes)): + axes[i].spines["top"].set_visible(False) + axes[i].spines["right"].set_visible(False) + axes[i].xaxis.grid(False) + axes[i].yaxis.grid(False) + + if i > 0: + axes[i].get_yaxis().set_visible(False) + axes[i].spines["left"].set_visible(False) + + order = ["II","III","IV"] if subtypes[i] != 'ODG' else ["II", "III"] + + axes[i].xaxis.label.set_visible(False) + axes[i].yaxis.label.set_visible(False) + axes[i].tick_params(axis='y', which='both', labelsize=20) + axes[i].tick_params(axis='x', which='both', labelsize=15) + datapoints = data[data['Histomolecular subtype'] == subtypes[i]] + sns.boxplot(y='Hazard', x="Grade", data=datapoints, ax = axes[i], color=color_dict[i], order=order) + sns.stripplot(y='Hazard', x='Grade', data=datapoints, alpha=0.2, jitter=0.2, color='k', ax = axes[i], order=order) + axes[i].set_ylim(-2.5, 2.5) + axes[i].set_yticks(np.arange(-2.0, 2.1, 1)) + + #axes[2].legend(prop={'size': 10}) + fig.savefig(ckpt_name+'/%s_HBP.png' % (model)) + + +def makeAUROCPlot(ckpt_name='./checkpoints/grad_15/', model_list=['path', 'omic', 'pathgraphomic_fusion'], split='test', avg='micro', use_zoom=False): + mpl.rcParams['font.family'] = "arial" + colors = {'path':'dodgerblue', 'graph':'orange', 'omic':'green', 'pathgraphomic_fusion':'crimson'} + names = {'path':'Histology CNN', 'graph':'Histology GCN', 'omic':'Genomic SNN', 'pathgraphomic_fusion':'Pathomic F.'} + zoom_params = {0:([0.2, 0.4], [0.8, 1.0]), + 1:([0.25, 0.45], [0.75, 0.95]), + 2:([0.0, 0.2], [0.8, 1.0]), + 'micro':([0.15, 0.35], [0.8, 1.0])} + mean_fpr = np.linspace(0, 1, 100) + classes = [0, 1, 2, avg] + ### 1. Looping over classes + for i in classes: + print("Class: " + str(i)) + fi = pylab.figure(figsize=(10,10), dpi=600, linewidth=0.2) + axi = plt.subplot() + + ### 2. Looping over models + for m, model in enumerate(model_list): + ignore_missing_moltype = 1 if 'omic' in model else 0 + ignore_missing_histype = 1 if 'grad' in ckpt_name else 0 + use_patch, roi_dir, use_vgg_features = ('_patch_', 'all_st_patches_512', 1) if (('path' in model) or ('graph' in model)) else ('_', 'all_st', 0) + + ###. 3. Looping over all splits + tprs, pres, aucrocs, rocaucs, = [], [], [], [] + for k in range(1,16): + pred = pickle.load(open(ckpt_name+'/%s/%s_%d%spred_%s.pkl' % (model, model, k, use_patch, split), 'rb')) + grade_pred, grade = np.array(pred[3]), np.array(pred[4]) + enc = LabelBinarizer() + enc.fit(grade) + grade_oh = enc.transform(grade) + + if i != avg: + pres.append(average_precision_score(grade_oh[:, i], grade_pred[:, i])) # from https://scikit-learn.org/stable/auto_examples/model_selection/plot_precision_recall.html + fpr, tpr, thresh = roc_curve(grade_oh[:,i], grade_pred[:,i], drop_intermediate=False) + aucrocs.append(auc(fpr, tpr)) # https://scikit-learn.org/stable/auto_examples/model_selection/plot_roc.html + rocaucs.append(roc_auc_score(grade_oh[:,i], grade_pred[:,i])) # https://scikit-learn.org/stable/modules/generated/sklearn.metrics.roc_auc_score.html#sklearn.metrics.roc_auc_score + tprs.append(interp(mean_fpr, fpr, tpr)) + tprs[-1][0] = 0.0 + else: + # A "micro-average": quantifying score on all classes jointly + pres.append(average_precision_score(grade_oh, grade_pred, average=avg)) + fpr, tpr, thresh = roc_curve(grade_oh.ravel(), grade_pred.ravel()) + aucrocs.append(auc(fpr, tpr)) + rocaucs.append(roc_auc_score(grade_oh, grade_pred, avg)) + tprs.append(interp(mean_fpr, fpr, tpr)) + tprs[-1][0] = 0.0 + + mean_tpr = np.mean(tprs, axis=0) + mean_tpr[-1] = 1.0 + #mean_auc = auc(mean_fpr, mean_tpr) + mean_auc = np.mean(aucrocs) + std_auc = np.std(aucrocs) + print('\t'+'%s - AUC: %0.3f ± %0.3f' % (model, mean_auc, std_auc)) + + if use_zoom: + alpha, lw = (0.8, 6) if model =='pathgraphomic_fusion' else (0.5, 6) + plt.plot(mean_fpr, mean_tpr, color=colors[model], + label=r'%s (AUC = %0.3f $\pm$ %0.3f)' % (names[model], mean_auc, std_auc), lw=lw, alpha=alpha) + std_tpr = np.std(tprs, axis=0) + tprs_upper = np.minimum(mean_tpr + std_tpr, 1) + tprs_lower = np.maximum(mean_tpr - std_tpr, 0) + plt.fill_between(mean_fpr, tprs_lower, tprs_upper, color=colors[model], alpha=0.1) + plt.xlim([zoom_params[i][0][0]-0.005, zoom_params[i][0][1]+0.005]) + plt.ylim([zoom_params[i][1][0]-0.005, zoom_params[i][1][1]+0.005]) + axi.set_xticks(np.arange(zoom_params[i][0][0], zoom_params[i][0][1]+0.001, 0.05)) + axi.set_yticks(np.arange(zoom_params[i][1][0], zoom_params[i][1][1]+0.001, 0.05)) + axi.tick_params(axis='both', which='major', labelsize=26) + else: + alpha, lw = (0.8, 4) if model =='pathgraphomic_fusion' else (0.5, 3) + plt.plot(mean_fpr, mean_tpr, color=colors[model], + label=r'%s (AUC = %0.3f $\pm$ %0.3f)' % (names[model], mean_auc, std_auc), lw=lw, alpha=alpha) + std_tpr = np.std(tprs, axis=0) + tprs_upper = np.minimum(mean_tpr + std_tpr, 1) + tprs_lower = np.maximum(mean_tpr - std_tpr, 0) + plt.fill_between(mean_fpr, tprs_lower, tprs_upper, color=colors[model], alpha=0.1) + plt.xlim([-0.05, 1.05]) + plt.ylim([-0.05, 1.05]) + axi.set_xticks(np.arange(0, 1.001, 0.2)) + axi.set_yticks(np.arange(0, 1.001, 0.2)) + axi.legend(loc="lower right", prop={'size': 20}) + axi.tick_params(axis='both', which='major', labelsize=30) + #plt.plot([0, 1], [0, 1], linestyle='--', lw=2, color='navy', alpha=.8) + + figures = [manager.canvas.figure + for manager in mpl._pylab_helpers.Gcf.get_all_fig_managers()] + + zoom = '_zoom' if use_zoom else '' + for i, fig in enumerate(figures): + fig.savefig(ckpt_name+'/AUC_%s%s.png' % (classes[i], zoom)) \ No newline at end of file