<jupyter_start><jupyter_text>Predictive performance comparisonThe idea of this notebook is to take a look at the predictive performance on cell lines for all the drugs. The idea is two-fold:     Assessing that the source top PVs can yield same predictive performance as a direct ridge on the source data. It would mean that the top PVs contain the relevant information for drug response prediction.     Taking a look at which drug gets predicted using both the PV duos and the consensus representation.We here use all the cell line data for the domain adaptation. Other settings can be imagined as well.  Parameters (to change)<jupyter_code># None for 'rnaseq', 'fpkm' for FPKM
type_data = 'rnaseq'
normalization = 'TMM'
transformation = 'log'
mean_center = True
std_unit = False
filter_mytochondrial = False
protein_coding_only = True
d_test = [40]
n_factors = 70
same_pv_pca = True
drug_file = 'input/drug_list_small.txt' # To change to drug_list.txt for full-scale analysis
n_jobs=5
import os, sys
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import scipy
from sklearn.model_selection import GroupKFold, GridSearchCV
from sklearn.linear_model import ElasticNet, Ridge
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.externals.joblib import Parallel, delayed
import pickle
plt.style.use('ggplot')

#Import src implementations
os.environ['OMP_NUM_THREADS'] = '1'
os.environ['KMP_DUPLICATE_LIB_OK']='True'
from data_reader.read_data import read_data
from data_reader.read_drug_response import read_drug_response
from data_reader.read_cna_tumors import read_cna_tumors
from normalization_methods.feature_engineering import feature_engineering
import precise
from precise import DrugResponsePredictor, ConsensusRepresentation<jupyter_output><empty_output><jupyter_text>Read all the drug from the file and load all the data<jupyter_code>with open(drug_file,'r') as drug_file_reader:
    drug_file_content = drug_file_reader.read()
    drug_file_content = drug_file_content.split('\n')
    drug_file_content = [e.split(',') for e in drug_file_content]
    
    # drug_IDs and tumor tissues are ordered in the same way
    drug_IDs = np.array(list(zip(*drug_file_content))[0]).astype(int)
    tumor_tissues = np.array(list(zip(*drug_file_content))[1])

unique_tumor_tissues = np.unique(tumor_tissues)
target_raw_data = dict()
source_raw_data = dict()
target_barcodes = dict()
source_names = dict()
target_data = dict()
source_data = dict()
source_data_filtered = dict()
source_response_data = dict()
source_names_filtered = dict()
drug_names = dict()
target_primary_site = dict()
# Load cell line data
# /!\ Due to some mismatch in the genes available in TCGA, cell line data has to be loaded all the time
for tissue_name in unique_tumor_tissues:
    print(tissue_name)
    
    if tissue_name in target_raw_data:
        continue
        
    X_target, X_source, _, s, target_names = read_data('cell_line',
                                                       'tumor',
                                                       'count',
                                                        None,
                                                        tissue_name,
                                                        filter_mytochondrial)
    
    target_raw_data[tissue_name] = X_target
    source_raw_data[tissue_name] = X_source
    target_barcodes[tissue_name] = target_names
    source_names[tissue_name] = s
# Normalize the data
for tissue_name in unique_tumor_tissues:
    print(tissue_name)
    if tissue_name in target_data:
        continue
    
    target_data[tissue_name] = feature_engineering(target_raw_data[tissue_name],
                                                  normalization,
                                                  transformation,
                                                  mean_center,
                                                  std_unit)
    
    # source data is not mean-centered as it will be done during cross-validation procedure.
    source_data[tissue_name] = feature_engineering(source_raw_data[tissue_name],
                                                  normalization,
                                                  transformation,
                                                  False,
                                                  False)
# Normalize for variance
for tissue_name in unique_tumor_tissues:
    print(tissue_name)
    if tissue_name in target_data:
        continue
    
    target_total_variance = np.sqrt(np.sum(np.var(target_data[tissue_name], 0)))
    target_data[tissue_name] = target_data[tissue_name] / target_total_variance * 10**3
    
    source_total_variance = np.sqrt(np.sum(np.var(source_data[tissue_name], 0)))
    source_data[tissue_name] = source_data[tissue_name] / source_total_variance * 10**3
# Read drug response
for i, (ID, tissue) in enumerate(zip(drug_IDs, tumor_tissues)):
    if (ID, tissue) in source_data_filtered:
        continue
        
    x, y, s, name = read_drug_response(ID,
                                       source_data[tissue],
                                       source_names[tissue],
                                      'count')
    
    source_data_filtered[(ID, tissue)] = x
    source_response_data[(ID, tissue)] = y
    drug_names[(ID, tissue)] = name
    source_names_filtered[(ID, tissue)] = s<jupyter_output><empty_output><jupyter_text>Principal vector testHere we compute the predictive performance for several different drugs using either the osurce, the target of both principal vector. The latter one is still biases towards the source.  Consensus representation<jupyter_code>l1_ratio  = 0

for ID, tissue in zip(drug_IDs, tumor_tissues):
    print(ID, tissue)
    
    X_source = source_data_filtered[ID, tissue]
    y_source = source_response_data[ID, tissue]
    X_target = target_data[tissue]
    
    pickle_file = 'consensus_drug_%s_tissue_%s_l1_ratio_%s_n_factors_%s.pkl'%(ID,
                                                                      tissue,
                                                                      l1_ratio,
                                                                      n_factors)
    if pickle_file in os.listdir('./output/pred_performance/'):
        print('%s, %s ALREADY COMPUTED'%(ID, tissue))
        continue
        
    with open('./output/pred_performance/%s'%(pickle_file), 'wb') as f:
        pickle.dump(dict(), f, pickle.HIGHEST_PROTOCOL)
    
    pred_performance = {}
    for d in d_test:
        print(d)

        predictor = DrugResponsePredictor(source_data=source_data[tissue][~np.isin(source_names[tissue], source_names_filtered[(ID, tissue)])],\
                                            method='consensus',\
                                            n_representations = 100,\
                                            target_data=X_target,\
                                            n_pv=d,\
                                            n_factors=n_factors,\
                                            n_jobs=n_jobs,\
                                            mean_center=mean_center,\
                                            std_unit=std_unit,\
                                            l1_ratio=l1_ratio)
        predictor.alpha_values = list(np.logspace(-2,10,17))
        predictor.verbose = 5
        predictor.fit(X_source, y_source, use_data=True)
        pred_performance[d] = predictor.compute_predictive_performance(X_source, y_source)
        plt.plot(predictor.alpha_values, predictor.regression_model_.cv_results_['mean_test_score'], '+-')
        plt.title(pred_performance[d])
        plt.xscale('log')
        plt.show()
    
    with open('./output/pred_performance/%s'%(pickle_file), 'wb') as f:
        pickle.dump(pred_performance, f, pickle.HIGHEST_PROTOCOL)<jupyter_output><empty_output><jupyter_text>ElasticNet/Ridge comparison<jupyter_code>from sklearn.model_selection import GroupKFold

l1_ratio  = 0.

pickle_file = 'elasticnet_drug_l1_ratio_%s_std.pkl'%(l1_ratio)
if pickle_file in os.listdir('./output/pred_performance/'):
    with open('./output/pred_performance/%s'%(pickle_file), 'rb') as f:
        elasticnet_perf = pickle.load(f)

for ID, tissue in zip(drug_IDs, tumor_tissues):
    print(ID, tissue)
    
    pickle_file = 'en_std_drug_%s_tissue_%s_l1_ratio_%s_n_factors_%s.pkl'%(ID,
                                                                          tissue,
                                                                          l1_ratio,
                                                                          n_factors)
    
    if pickle_file in os.listdir('./output/pred_performance/'):
        print('%s, %s ALREADY COMPUTED'%(ID, tissue))
        continue
    
    if (ID, tissue) in elasticnet_perf:
        continue
        
    with open('./output/pred_performance/%s'%(pickle_file), 'wb') as f:
        pickle.dump(dict(), f, pickle.HIGHEST_PROTOCOL)
    
    X_source = source_data_filtered[ID, tissue]
    y_source = source_response_data[ID, tissue]
    X_target = target_data[tissue]

    #Parameters for the grid search
    alpha_values = np.logspace(-5,10,16)
    param_grid ={
        'regression__alpha': alpha_values
    }

    #Grid search setup

    k_fold_split = GroupKFold(10)
    y_predicted = np.zeros(X_source.shape[0])
    
    for train_index, test_index in k_fold_split.split(X_source, y_source, y_source):
        grid_en = GridSearchCV(Pipeline([
                                ('normalization', StandardScaler(with_mean=mean_center, with_std=True)),
                                ('regression', ElasticNet(l1_ratio) if l1_ratio > 0 else Ridge())
                            ]),\
                            cv=10, n_jobs=30, param_grid=param_grid, verbose=1, scoring='neg_mean_squared_error')
        grid_en.fit(X_source[train_index], y_source[train_index])
        y_predicted[test_index] = grid_en.predict(X_source[test_index])
    
    #Fit grid search
    grid_en.fit(X_source, y_source)
    elasticnet_perf[ID, tissue] = scipy.stats.pearsonr(y_predicted, y_source)[0]
    print(elasticnet_perf[ID, tissue])
    
    with open('./output/pred_performance/%s'%(pickle_file), 'wb') as f:
        pickle.dump(elasticnet_perf[ID, tissue], f, pickle.HIGHEST_PROTOCOL)<jupyter_output><empty_output><jupyter_text>Load pickle and look at results<jupyter_code>l1_ratio = 0
l1_ratio_en = 0.

two_pv_results = dict()
consensus_pv_results = dict()
source_pv_results = dict()
target_pv_results = dict()
en_results_std = dict()

def sort_dictionary(d):
    return {e:d[e] for e in sorted(d)}

for ID, tissue in zip(drug_IDs, tumor_tissues):
    print(ID, tissue)
        
    # Read results of consensus PVs
    pickle_file = 'consensus_drug_%s_tissue_%s_l1_ratio_%s_n_factors_%s.pkl'%(ID,
                                                                      tissue,
                                                                      l1_ratio,
                                                                      n_factors)
    with open('./output/pred_performance/%s'%(pickle_file), 'rb') as f:
        consensus_pv_results[ID,tissue] = sort_dictionary(pickle.load(f))
        
    
    # Read results of EN
    pickle_file = 'en_std_drug_%s_tissue_%s_l1_ratio_%s_n_factors_%s.pkl'%(ID,
                                                                          tissue,
                                                                          '0.0',
                                                                          n_factors)
    
    with open('./output/pred_performance/%s'%(pickle_file), 'rb') as f:
        en_results_std[ID,tissue] = pickle.load(f)
        print(en_results[ID, tissue])
for ID, tissue in zip(drug_IDs, tumor_tissues):
    # Plot for a specific number of PV
    plt.plot([e[0] for e in consensus_pv_results[ID,tissue].items()],
             [e[1] for e in consensus_pv_results[ID,tissue].items()],
             label='consensus', linewidth=3, alpha=0.5, marker='+')
    plt.plot([e[0] for e in source_pv_results[ID,tissue].items()],
             [e[1] for e in source_pv_results[ID,tissue].items()],
             label='source', linewidth=3, alpha=0.5, marker='+')
    plt.plot([e[0] for e in target_pv_results[ID,tissue].items()],
             [e[1] for e in target_pv_results[ID,tissue].items()],
             label='target', linewidth=3, alpha=0.5, marker='+')
    plt.plot([e[0] for e in two_pv_results[ID,tissue].items()],
             [e[1] for e in two_pv_results[ID,tissue].items()],
             label='2 pv', linewidth=3, alpha=0.5, marker='+')
    plt.hlines(en_results[ID,tissue], xmin=0, xmax=plt.xlim()[1], label='Ridge', linewidth=3, alpha=0.7)

    plt.title(drug_names[ID, tissue] + ' '+ tissue)
    plt.xlabel('Number of Principal Vectors', fontsize=15)
    plt.ylabel('Predictive Performance', fontsize=15)
    plt.legend()    
    plt.show()
n_pv = 40

perf_scatter = []
for ID, tissue in zip(drug_IDs, tumor_tissues):
    #print(ID, tissue)
    if n_pv not in consensus_pv_results[ID,tissue]:
        print(ID, tissue)
        continue
    plt.scatter(en_results_std[ID,tissue],
                consensus_pv_results[ID,tissue][n_pv],
               color='blue', marker='x', alpha=0.7)
    perf_scatter.append([en_results_std[ID,tissue], consensus_pv_results[ID,tissue][n_pv]])
    
plt.xlabel('ElasticNet', fontsize=20)
plt.ylabel('Consensus \n representation', fontsize=20)
plt.xticks(fontsize=15, color='black')
plt.yticks(fontsize=15, color='black')
plt.tight_layout()
plt.xlim(0.1,0.8)
plt.ylim(0.1,0.8)
plt.plot(plt.xlim(), plt.xlim(), linewidth=3, alpha=0.5)
#plt.savefig('./figures/fig4_pred_perf_consensus_%s_en_%s.png'%(l1_ratio, l1_ratio_en), dpi=300)
plt.show()

perf_scatter = np.array(perf_scatter)
p = scipy.stats.pearsonr(perf_scatter[:,0], perf_scatter[:,1])
print('Pearson Correlation: %s, %s'%(p[0], p[1]))
plt.scatter(perf_scatter[:,1], (perf_scatter[:,0] - perf_scatter[:,1])/perf_scatter[:,0])

np.median((perf_scatter[:,0] - perf_scatter[:,1])/perf_scatter[:,0])

#for e in en_results:
#    print(e, en_results[e], consensus_pv_results[e])
for ID, tissue in zip(drug_IDs, tumor_tissues):
    #print(ID, tissue)
    if n_pv not in consensus_pv_results[ID,tissue]:
        print(ID, tissue)
        continue
    plt.scatter(en_results[ID,tissue],
                en_results_std[ID,tissue],
               color='blue', marker='x', alpha=0.7)
    #perf_scatter.append([en_results[ID,tissue], consensus_pv_results[ID,tissue][n_pv]])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10050525740816606, 'fraction_numerical': 0.009217533797623925, 'mean_word_length': 2.313800904977376, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 22, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 83, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '16549', 'license': 'Apache-2.0', 'n_tokens_mistral': 4583, 'n_tokens_neox': 4382, 'n_words': 920, 'path': 'fig4_predictive_performance_comparison.ipynb', 'repo_name': 'NKI-CCB/PRECISE_figures', 'size': '14646'}