diff --git a/ras_analysis.sh b/ras_analysis.sh
index 5632baa..bb615b3 100755
--- a/ras_analysis.sh
+++ b/ras_analysis.sh
@@ -14,12 +14,16 @@ ras_diseases='BLCA,CESC,COAD,ESCA,HNSC,LUAD,LUSC,OV,PAAD,PCPG,READ,SKCM,STAD,TGC
# 1. PanCancer NF1 Classification
python scripts/pancancer_classifier.py --genes 'NF1' --drop --copy_number \
--diseases $nf1_diseases --alphas $alphas --l1_ratios $l1_mixing \
- --remove_hyper --alt_folder 'classifiers/NF1'
+ --remove_hyper --shuffled --alt_folder 'classifiers/NF1_shuffled' \
+ --keep_intermediate
# 2. PanCancer RAS Classification and predict NF1 using RAS classifier
python scripts/pancancer_classifier.py --genes 'KRAS,HRAS,NRAS' --drop \
- --remove_hyper --copy_number --alphas $alphas --l1_ratios $l1_mixing \
- --alt_genes 'NF1' --alt_diseases $nf1_diseases --alt_folder 'classifiers/RAS'
+ --diseases $ras_diseases --copy_number --remove_hyper \
+ --alphas $alphas --l1_ratios $l1_mixing \
+ --shuffled --alt_genes 'NF1' --alt_diseases $nf1_diseases \
+ --alt_folder 'classifiers/RAS_shuffled' \
+ --keep_intermediate
# 3. Within cancer-type NF1 Classification
python scripts/within_tissue_analysis.py --genes 'NF1' \
@@ -57,7 +61,8 @@ ras_no_thca_skcm=${ras_no_thca_skcm/THCA,}
python scripts/pancancer_classifier.py --genes 'KRAS,HRAS,NRAS' --drop \
--remove_hyper --copy_number --alphas $alphas --l1_ratio $l1_mixing \
- --diseases $ras_no_thca_skcm --alt_folder 'classifiers/RAS_noTHCASKCM'
+ --diseases $ras_no_thca_skcm --shuffled --alt_folder 'classifiers/RAS_noTHCASKCM_shuffled' \
+ --keep_intermediate
python scripts/apply_weights.py --classifier 'classifiers/RAS_noTHCASKCM' --copy_number
python scripts/map_mutation_class.py --scores 'classifiers/RAS_noTHCASKCM' \
@@ -69,3 +74,51 @@ Rscript --vanilla scripts/viz/ras_summary_figures.R
Rscript --vanilla scripts/viz/nf1_summary_figures.R
Rscript --vanilla scripts/viz/braf_summary_figures.R
+# 10. Additional Benchmarking Analysis
+# Randomly shuffle input RNAseq features and build a classifier
+python scripts/pancancer_classifier.py --genes 'KRAS,HRAS,NRAS' \
+ --diseases $ras_diseases --copy_number --remove_hyper \
+ --alphas $alphas --l1_ratios $l1_mixing \
+ --shuffled_before_training \
+ --keep_intermediate --alt_folder 'classifiers/RAS_shuffled_before_training'
+
+# Do not include copy number in the classifier construction
+# Note that the shuffled flag here makes classifier predictions on shuffled RNAseq data
+python scripts/pancancer_classifier.py --genes 'KRAS,HRAS,NRAS' --drop \
+ --diseases $ras_diseases \
+ --alphas $alphas --l1_ratios $l1_mixing --remove_hyper \
+ --shuffled --alt_folder 'classifiers/RAS_nocopy' --keep_intermediate
+
+# Do not include mutation in the classifier construction
+python scripts/pancancer_classifier.py --genes 'KRAS,HRAS,NRAS' --drop \
+ --diseases $ras_diseases --copy_number --no_mutation \
+ --alphas $alphas --l1_ratios $l1_mixing --remove_hyper \
+ --shuffled --alt_folder 'classifiers/RAS_nomutation' --keep_intermediate
+
+# Drop all Rasopathy genes
+python scripts/pancancer_classifier.py --genes 'KRAS,HRAS,NRAS' --drop --drop_rasopathy \
+ --diseases $ras_diseases --copy_number --remove_hyper \
+ --alphas $alphas --l1_ratios $l1_mixing --shuffled \
+ --alt_folder 'classifiers/RAS_droprasopathy'
+
+# Use only covariate information
+python scripts/pancancer_classifier.py --genes 'KRAS,HRAS,NRAS' --drop \
+ --diseases $ras_diseases --copy_number --remove_hyper \
+ --alphas $alphas --l1_ratios $l1_mixing \
+ --shuffled --keep_intermediate --drop_expression \
+ --alt_folder 'classifiers/RAS_onlycovariate'
+
+# Use only gene expression information
+python scripts/pancancer_classifier.py --genes 'KRAS,HRAS,NRAS' --drop \
+ --diseases $ras_diseases --copy_number --remove_hyper \
+ --alphas $alphas --l1_ratios $l1_mixing \
+ --shuffled --keep_intermediate --drop_covariate \
+ --alt_folder 'classifiers/RAS_onlyexpression'
+
+# Drop no genes
+python scripts/pancancer_classifier.py --genes 'KRAS,HRAS,NRAS' \
+ --diseases $ras_diseases --copy_number --remove_hyper \
+ --alphas $alphas --l1_ratios $l1_mixing \
+ --shuffled --keep_intermediate \
+ --alt_folder 'classifiers/RAS_nodrop' \
+
diff --git a/scripts/pancancer_classifier.py b/scripts/pancancer_classifier.py
index 75f15ab..514c4b1 100644
--- a/scripts/pancancer_classifier.py
+++ b/scripts/pancancer_classifier.py
@@ -53,6 +53,7 @@
import os
import sys
import warnings
+from random import sample
import pandas as pd
import csv
@@ -69,6 +70,7 @@
sys.path.insert(0, os.path.join('scripts', 'util'))
from tcga_util import get_threshold_metrics, integrate_copy_number
+from tcga_util import shuffle_columns
parser = argparse.ArgumentParser()
parser.add_argument('-g', '--genes',
@@ -111,6 +113,18 @@
help='Filename of features to use in model')
parser.add_argument('-y', '--y_matrix', default='default',
help='Either `default` or `xena` based on underlying data')
+parser.add_argument('-e', '--shuffled', action='store_true',
+ help='Shuffle the input gene expression matrix alongside')
+parser.add_argument('--shuffled_before_training', action='store_true',
+ help='Shuffle the gene expression matrix before training')
+parser.add_argument('-m', '--no_mutation', action='store_false',
+ help='Remove mutation data from y matrix')
+parser.add_argument('-z', '--drop_rasopathy', action='store_true',
+ help='Decision to drop all rasopathy genes from X matrix')
+parser.add_argument('-q', '--drop_expression', action='store_true',
+ help='Decision to drop all gene expression values from X')
+parser.add_argument('-j', '--drop_covariates', action='store_true',
+ help='Decision to drop all covariate information from X')
args = parser.parse_args()
# Load command arguments
@@ -118,6 +132,7 @@
diseases = args.diseases.split(',')
folds = int(args.folds)
drop = args.drop
+drop_rasopathy = args.drop_rasopathy
copy_number = args.copy_number
filter_count = int(args.filter_count)
filter_prop = float(args.filter_prop)
@@ -133,6 +148,11 @@
keep_inter = args.keep_intermediate
x_matrix = args.x_matrix
y_matrix = args.y_matrix
+shuffled = args.shuffled
+shuffled_before_training = args.shuffled_before_training
+no_mutation = args.no_mutation
+drop_expression = args.drop_expression
+drop_covariates = args.drop_covariates
warnings.filterwarnings('ignore',
message='Changing the shape of non-C contiguous array')
@@ -155,13 +175,13 @@
os.makedirs(disease_folder)
count_table_file = os.path.join(base_folder, 'summary_counts.csv')
-cv_heatmap_file = os.path.join(base_folder, 'cv_heatmap.svg')
-full_roc_file = os.path.join(base_folder, 'all_disease_roc.svg')
-full_pr_file = os.path.join(base_folder, 'all_disease_pr.svg')
+cv_heatmap_file = os.path.join(base_folder, 'cv_heatmap.pdf')
+full_roc_file = os.path.join(base_folder, 'all_disease_roc.pdf')
+full_pr_file = os.path.join(base_folder, 'all_disease_pr.pdf')
disease_roc_file = os.path.join(base_folder, 'disease', 'classifier_roc_')
disease_pr_file = os.path.join(base_folder, 'disease', 'classifier_pr_')
-dis_summary_auroc_file = os.path.join(base_folder, 'disease_auroc.svg')
-dis_summary_aupr_file = os.path.join(base_folder, 'disease_aupr.svg')
+dis_summary_auroc_file = os.path.join(base_folder, 'disease_auroc.pdf')
+dis_summary_aupr_file = os.path.join(base_folder, 'disease_aupr.pdf')
classifier_file = os.path.join(base_folder, 'classifier_coefficients.tsv')
roc_results_file = os.path.join(base_folder, 'pancan_roc_results.tsv')
@@ -170,9 +190,9 @@
args.alt_diseases.replace(',', '_'))
alt_count_table_file = os.path.join(base_folder, 'alt_summary_counts.csv')
alt_gene_auroc_file = os.path.join(base_folder,
- '{}_auroc_bar.svg'.format(alt_gene_base))
+ '{}_auroc_bar.pdf'.format(alt_gene_base))
alt_gene_aupr_file = os.path.join(base_folder,
- '{}_aupr_bar.svg'.format(alt_gene_base))
+ '{}_aupr_bar.pdf'.format(alt_gene_base))
alt_gene_summary_file = os.path.join(base_folder,
'{}_summary.tsv'.format(alt_gene_base))
@@ -220,6 +240,13 @@
if x_matrix == 'raw':
rnaseq_full_df.drop(common_genes, axis=1, inplace=True)
+if drop_rasopathy:
+ rasopathy_genes = set(['BRAF', 'CBL', 'HRAS', 'KRAS', 'MAP2K1', 'MAP2K2',
+ 'NF1', 'NRAS', 'PTPN11', 'RAF1', 'SHOC2', 'SOS1',
+ 'SPRED1', 'RIT1'])
+ rasopathy_drop = list(rasopathy_genes.intersection(rnaseq_full_df.columns))
+ rnaseq_full_df.drop(rasopathy_drop, axis=1, inplace=True)
+
# Incorporate copy number for gene activation/inactivation
if copy_number:
# Load copy number matrices
@@ -236,7 +263,8 @@
y = integrate_copy_number(y=y, cancer_genes_df=cancer_genes,
genes=common_genes, loss_df=copy_loss_df,
- gain_df=copy_gain_df)
+ gain_df=copy_gain_df,
+ include_mutation=no_mutation)
# Process y matrix
y = y.assign(total_status=y.max(axis=1))
@@ -295,6 +323,21 @@
x_df_update.index = x_df.index
x_df = x_df_update.merge(covar, left_index=True, right_index=True)
+# Remove information from the X matrix given input arguments
+if drop_expression:
+ x_df = x_df.iloc[:, num_features_kept:]
+elif drop_covariates:
+ x_df = x_df.iloc[:, 0:num_features_kept]
+
+# Shuffle expression matrix _before_ training - this can be used as NULL model
+if shuffled_before_training:
+ # Shuffle genes
+ x_train_genes = x_df.iloc[:, 0:num_features_kept]
+ rnaseq_shuffled_df = x_train_genes.apply(shuffle_columns, axis=1)
+
+ x_train_cov = x_df.iloc[:, num_features_kept:]
+ x_df = pd.concat([rnaseq_shuffled_df, x_train_cov], axis=1)
+
# Build classifier pipeline
x_train, x_test, y_train, y_test = train_test_split(x_df, y_df, test_size=0.1,
random_state=0,
@@ -326,7 +369,7 @@
ax.set_xlabel('Regularization strength multiplier (alpha)')
ax.set_ylabel('Elastic net mixing parameter (l1_ratio)')
plt.tight_layout()
-plt.savefig(cv_heatmap_file, dpi=600, format='svg', bbox_inches='tight')
+plt.savefig(cv_heatmap_file, dpi=600, bbox_inches='tight')
plt.close()
# Get predictions
@@ -345,6 +388,21 @@
metrics_cv = get_threshold_metrics(y_train, y_cv,
drop_intermediate=keep_inter)
+# Determine shuffled predictive ability of shuffled gene expression matrix
+# representing a test of inflation of ROC metrics. Be sure to only shuffle
+# gene names, retain covariate information (tissue type and log10 mutations)
+if shuffled:
+ # Shuffle genes
+ x_train_genes = x_train.iloc[:, 0:num_features_kept]
+ rnaseq_shuffled_df = x_train_genes.apply(shuffle_columns, axis=1)
+
+ x_train_cov = x_train.iloc[:, num_features_kept:]
+ rnaseq_shuffled_df = pd.concat([rnaseq_shuffled_df, x_train_cov], axis=1)
+
+ y_predict_shuffled = cv_pipeline.decision_function(rnaseq_shuffled_df)
+ metrics_shuffled = get_threshold_metrics(y_train, y_predict_shuffled,
+ drop_intermediate=keep_inter)
+
# Decide to save ROC results to file
if keep_inter:
train_roc = metrics_train['roc_df']
@@ -354,16 +412,25 @@
cv_roc = metrics_cv['roc_df']
cv_roc = cv_roc.assign(train_type='cv')
full_roc_df = pd.concat([train_roc, test_roc, cv_roc])
+ if shuffled:
+ shuffled_roc = metrics_shuffled['roc_df']
+ shuffled_roc = shuffled_roc.assign(train_type='shuffled')
+ full_roc_df = pd.concat([full_roc_df, shuffled_roc])
full_roc_df = full_roc_df.assign(disease='PanCan')
# Plot ROC
sns.set_style("whitegrid")
plt.figure(figsize=(2.7, 2.4))
total_auroc = {}
-colors = ['blue', 'green', 'orange']
+colors = ['blue', 'green', 'orange', 'grey']
idx = 0
-for label, metrics in [('Training', metrics_train), ('Testing', metrics_test),
- ('CV', metrics_cv)]:
+
+metrics_list = [('Training', metrics_train), ('Testing', metrics_test),
+ ('CV', metrics_cv)]
+if shuffled:
+ metrics_list += [('Random', metrics_shuffled)]
+
+for label, metrics in metrics_list:
roc_df = metrics['roc_df']
plt.plot(roc_df.fpr, roc_df.tpr,
@@ -371,6 +438,8 @@
linewidth=2, c=colors[idx])
total_auroc[label] = metrics['auroc']
idx += 1
+
+plt.plot([0, 1], [0, 1], color='navy', linewidth=2, linestyle='--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('False Positive Rate', fontsize=8)
@@ -380,24 +449,29 @@
plt.legend(bbox_to_anchor=(0.2, -0.45, 0.7, .202), loc=0, borderaxespad=0.,
fontsize=7.5)
plt.tight_layout()
-plt.savefig(full_roc_file, dpi=600, format='svg', bbox_inches='tight')
+plt.savefig(full_roc_file, dpi=600, bbox_inches='tight')
plt.close()
# Plot PR
sns.set_style("whitegrid")
plt.figure(figsize=(2.7, 2.4))
total_aupr = {}
-colors = ['blue', 'green', 'orange']
+colors = ['blue', 'green', 'orange', 'grey']
idx = 0
-for label, metrics in [('Training', metrics_train), ('Testing', metrics_test),
- ('CV', metrics_cv)]:
+metrics_list = [('Training', metrics_train), ('Testing', metrics_test),
+ ('CV', metrics_cv)]
+if shuffled:
+ metrics_list += [('Random', metrics_shuffled)]
+
+for label, metrics in metrics_list:
pr_df = metrics['pr_df']
plt.plot(pr_df.recall, pr_df.precision,
label='{} (AUPR = {:.1%})'.format(label, metrics['aupr']),
linewidth=2, c=colors[idx])
total_aupr[label] = metrics['aupr']
idx += 1
+
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('Recall', fontsize=8)
@@ -407,7 +481,7 @@
plt.legend(bbox_to_anchor=(0.2, -0.45, 0.7, .202), loc=0, borderaxespad=0.,
fontsize=7.5)
plt.tight_layout()
-plt.savefig(full_pr_file, dpi=600, format='svg', bbox_inches='tight')
+plt.savefig(full_pr_file, dpi=600, bbox_inches='tight')
plt.close()
# disease specific performance
@@ -418,10 +492,14 @@
# Get true and predicted training labels
y_disease_train = y_train[y_train.index.isin(sample_sub)]
+ if y_disease_train.sum() < 1:
+ continue
y_disease_predict_train = y_predict_train[y_train.index.isin(sample_sub)]
# Get true and predicted testing labels
y_disease_test = y_test[y_test.index.isin(sample_sub)]
+ if y_disease_test.sum() < 1:
+ continue
y_disease_predict_test = y_predict_test[y_test.index.isin(sample_sub)]
# Get predicted labels for samples when they were in cross validation set
@@ -442,6 +520,14 @@
disease=disease,
drop_intermediate=keep_inter)
+ # Get predictions and metrics with shuffled gene expression
+ if shuffled:
+ y_dis_predict_shuf = y_predict_shuffled[y_train.index.isin(sample_sub)]
+ met_shuff_dis = get_threshold_metrics(y_disease_train,
+ y_dis_predict_shuf,
+ disease=disease,
+ drop_intermediate=keep_inter)
+
if keep_inter:
train_roc = met_train_dis['roc_df']
train_roc = train_roc.assign(train_type='train')
@@ -450,25 +536,41 @@
cv_roc = met_cv_dis['roc_df']
cv_roc = cv_roc.assign(train_type='cv')
full_dis_roc_df = train_roc.append(test_roc).append(cv_roc)
+
+ if shuffled:
+ shuffled_roc = met_shuff_dis['roc_df']
+ shuffled_roc = shuffled_roc.assign(train_type='shuffled')
+ full_dis_roc_df = full_dis_roc_df.append(shuffled_roc)
+
full_dis_roc_df = full_dis_roc_df.assign(disease=disease)
full_roc_df = full_roc_df.append(full_dis_roc_df)
# Store results in disease indexed dictionary
disease_metrics[disease] = [met_train_dis, met_test_dis, met_cv_dis]
+ if shuffled:
+ disease_metrics[disease] += [met_shuff_dis]
+
disease_auroc = {}
disease_aupr = {}
for disease, metrics_val in disease_metrics.items():
- met_train, met_test, met_cv = metrics_val
- disease_pr_sub_file = '{}_pred_{}.svg'.format(disease_pr_file, disease)
- disease_roc_sub_file = '{}_pred_{}.svg'.format(disease_roc_file, disease)
+
+ labels = ['Training', 'Testing', 'CV', 'Random']
+ met_list = []
+ idx = 0
+ for met in metrics_val:
+ lab = labels[idx]
+ met_list.append((lab, met))
+ idx += 1
+
+ disease_pr_sub_file = '{}_pred_{}.pdf'.format(disease_pr_file, disease)
+ disease_roc_sub_file = '{}_pred_{}.pdf'.format(disease_roc_file, disease)
# Plot disease specific PR
plt.figure(figsize=(2.7, 2.4))
aupr = []
idx = 0
- for label, metrics in [('Training', met_train), ('Testing', met_test),
- ('CV', met_cv)]:
+ for label, metrics in met_list:
pr_df = metrics['pr_df']
plt.plot(pr_df.recall, pr_df.precision,
label='{} (AUPR = {:.1%})'.format(label, metrics['aupr']),
@@ -486,16 +588,14 @@
plt.legend(bbox_to_anchor=(0.2, -0.45, 0.7, .202), loc=0, borderaxespad=0.,
fontsize=7.5)
plt.tight_layout()
- plt.savefig(disease_pr_sub_file, dpi=600, format='svg',
- bbox_inches='tight')
+ plt.savefig(disease_pr_sub_file, dpi=600, bbox_inches='tight')
plt.close()
# Plot disease specific ROC
plt.figure(figsize=(2.7, 2.4))
auroc = []
idx = 0
- for label, metrics in [('Training', met_train), ('Testing', met_test),
- ('CV', met_cv)]:
+ for label, metrics in met_list:
roc_df = metrics['roc_df']
plt.plot(roc_df.fpr, roc_df.tpr,
label='{} (AUROC = {:.1%})'.format(label, metrics['auroc']),
@@ -504,6 +604,7 @@
idx += 1
disease_auroc[disease] = auroc
+ plt.plot([0, 1], [0, 1], color='navy', linewidth=2, linestyle='--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('False Positive Rate', fontsize=8)
@@ -513,28 +614,30 @@
plt.legend(bbox_to_anchor=(0.2, -0.45, 0.7, .202), loc=0, borderaxespad=0.,
fontsize=7.5)
plt.tight_layout()
- plt.savefig(disease_roc_sub_file, dpi=600, format='svg',
- bbox_inches='tight')
+ plt.savefig(disease_roc_sub_file, dpi=600, bbox_inches='tight')
plt.close()
-disease_auroc_df = pd.DataFrame(disease_auroc, index=['Train', 'Test',
- 'Cross Validation']).T
+index_lab = ['Train', 'Test', 'Cross Validation']
+
+if shuffled:
+ index_lab += ['Random']
+
+disease_auroc_df = pd.DataFrame(disease_auroc, index=index_lab).T
disease_auroc_df = disease_auroc_df.sort_values('Cross Validation',
ascending=False)
ax = disease_auroc_df.plot(kind='bar', title='Disease Specific Performance')
ax.set_ylabel('AUROC')
plt.tight_layout()
-plt.savefig(dis_summary_auroc_file, dpi=600, format='svg', bbox_inches='tight')
+plt.savefig(dis_summary_auroc_file, dpi=600, bbox_inches='tight')
plt.close()
-disease_aupr_df = pd.DataFrame(disease_aupr, index=['Train', 'Test',
- 'Cross Validation']).T
+disease_aupr_df = pd.DataFrame(disease_aupr, index=index_lab).T
disease_aupr_df = disease_aupr_df.sort_values('Cross Validation',
ascending=False)
ax = disease_aupr_df.plot(kind='bar', title='Disease Specific Performance')
ax.set_ylabel('AUPR')
plt.tight_layout()
-plt.savefig(dis_summary_aupr_file, dpi=600, format='svg', bbox_inches='tight')
+plt.savefig(dis_summary_aupr_file, dpi=600, bbox_inches='tight')
plt.close()
# Save classifier coefficients
@@ -636,6 +739,10 @@
y_sub = y_test[y_test.index.isin(sample_dis)]
category = 'Holdout'
+ # If there are not enough classes do not proceed to plot
+ if y_sub.sum() < 1:
+ continue
+
neg, pos = y_sub.value_counts()
val_x_type[disease] = [category, neg, pos]
y_pred_alt = cv_pipeline.decision_function(x_sub)
@@ -670,8 +777,7 @@
ax.set_ylim([0, 1])
ax.set_ylabel('AUROC')
plt.tight_layout()
- plt.savefig(alt_gene_auroc_file, dpi=600, format='svg',
- bbox_inches='tight')
+ plt.savefig(alt_gene_auroc_file, dpi=600, bbox_inches='tight')
plt.close()
# Plot alternative gene cancer-type specific AUPR plots
@@ -683,8 +789,7 @@
ax.set_ylim([0, 1])
ax.set_ylabel('AUPR')
plt.tight_layout()
- plt.savefig(alt_gene_aupr_file, dpi=600, format='svg',
- bbox_inches='tight')
+ plt.savefig(alt_gene_aupr_file, dpi=600, bbox_inches='tight')
plt.close()
# Write a summary for the inputs and outputs of the classifier
diff --git a/scripts/ras_differential_expression.ipynb b/scripts/ras_differential_expression.ipynb
new file mode 100644
index 0000000..d98ae0f
--- /dev/null
+++ b/scripts/ras_differential_expression.ipynb
@@ -0,0 +1,416 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Ras Differential Expression Analysis\n",
+ "\n",
+ "Perform a t-test across all genes to determine differentially expressed genes between Ras wildtype and Ras mutant genes. Compare with what the model learns and output the data for visualization with R/ggplot2."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "/home/gway/anaconda3/lib/python3.5/site-packages/statsmodels/compat/pandas.py:56: FutureWarning: The pandas.core.datetools module is deprecated and will be removed in a future version. Please use the pandas.tseries module instead.\n",
+ " from pandas.core import datetools\n"
+ ]
+ }
+ ],
+ "source": [
+ "import os\n",
+ "import sys\n",
+ "import pandas as pd\n",
+ "import scipy\n",
+ "from sklearn.preprocessing import StandardScaler\n",
+ "\n",
+ "import matplotlib.pyplot as plt\n",
+ "import seaborn as sns\n",
+ "import plotnine as gg\n",
+ "\n",
+ "sys.path.insert(0, os.path.join('..', 'scripts', 'util'))\n",
+ "from tcga_util import integrate_copy_number"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "%matplotlib inline\n",
+ "plt.style.use('seaborn-notebook')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Load Data\n",
+ "\n",
+ "Loading RNAseq (X matrix), copy number and mutation (Y matrix), and coefficients data (W matrix from machine learning model)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "# Load RNAseq matrix\n",
+ "expr_file = os.path.join('..', 'data', 'pancan_rnaseq_freeze.tsv')\n",
+ "rnaseq_full_df = pd.read_table(expr_file, index_col=0)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "# Load Mutation matrix\n",
+ "mut_file = os.path.join('..', 'data', 'pancan_mutation_freeze.tsv')\n",
+ "mutation_df = pd.read_table(mut_file, index_col=0)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "# Load sample freeze data and cancer genes\n",
+ "sample_freeze_file = os.path.join('..', 'data', 'sample_freeze.tsv')\n",
+ "sample_freeze = pd.read_table(sample_freeze_file, index_col=0)\n",
+ "\n",
+ "cancer_gene_file = os.path.join('..', 'data', 'vogelstein_cancergenes.tsv')\n",
+ "cancer_genes = pd.read_table(cancer_gene_file)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "# Load copy number data to determine final status matrix\n",
+ "copy_loss_file = os.path.join('..', 'data', 'copy_number_loss_status.tsv')\n",
+ "copy_gain_file = os.path.join('..', 'data', 'copy_number_gain_status.tsv')\n",
+ "\n",
+ "copy_loss_df = pd.read_table(copy_loss_file, index_col=0)\n",
+ "copy_gain_df = pd.read_table(copy_gain_file, index_col=0)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "# Load Coefficients File\n",
+ "coef_file = os.path.join('..', 'classifiers', 'RAS', 'classifier_coefficients.tsv')\n",
+ "coef_df = pd.read_table(coef_file, index_col=0)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Process and Prep Data for Differential Expression Analysis"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "# Process y matrix\n",
+ "genes = ['KRAS', 'HRAS', 'NRAS']\n",
+ "y = mutation_df[genes]\n",
+ "\n",
+ "y_df = integrate_copy_number(y=y, cancer_genes_df=cancer_genes,\n",
+ " genes=genes, loss_df=copy_loss_df,\n",
+ " gain_df=copy_gain_df)\n",
+ "\n",
+ "y_df = y_df.assign(total_status=y.max(axis=1))\n",
+ "y_df = y_df.reset_index().merge(sample_freeze,\n",
+ " how='left').set_index('SAMPLE_BARCODE')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "0 8008\n",
+ "1 1066\n",
+ "Name: total_status, dtype: int64"
+ ]
+ },
+ "execution_count": 9,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "y_df['total_status'].value_counts()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "# Subset X matrix to coefficients\n",
+ "x_df = rnaseq_full_df.loc[:, coef_df['feature']]\n",
+ "x_df = x_df.dropna(axis=1)\n",
+ "\n",
+ "x_df_update = StandardScaler().fit_transform(x_df)\n",
+ "x_df_update = pd.DataFrame(x_df_update, columns=x_df.columns, index=x_df.index)\n",
+ "x_df = x_df_update"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "# Get two arrays to compare\n",
+ "ras_wt_samples = y_df[y_df['total_status'] == 0].index\n",
+ "ras_mut_samples = y_df[y_df['total_status'] == 1].index\n",
+ "\n",
+ "x_wt_df = x_df.loc[ras_wt_samples, :]\n",
+ "x_mut_df = x_df.loc[ras_mut_samples, :]"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Perform the analysis and output results"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "ttest_results = scipy.stats.ttest_ind(x_mut_df, x_wt_df)\n",
+ "t_stat = ttest_results.statistic\n",
+ "p_val = ttest_results.pvalue"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "ttest_df = pd.DataFrame(t_stat, columns=['stat'])\n",
+ "ttest_df = ttest_df.assign(pval = p_val)\n",
+ "ttest_df = ttest_df.assign(gene = x_wt_df.columns)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | \n",
+ " stat | \n",
+ " pval | \n",
+ " gene | \n",
+ " feature | \n",
+ " weight | \n",
+ " abs | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | 0 | \n",
+ " 6.812033 | \n",
+ " 1.023496e-11 | \n",
+ " PBX3 | \n",
+ " PBX3 | \n",
+ " 0.139761 | \n",
+ " 0.139761 | \n",
+ "
\n",
+ " \n",
+ " | 1 | \n",
+ " 21.001785 | \n",
+ " 1.164957e-95 | \n",
+ " SPRY2 | \n",
+ " SPRY2 | \n",
+ " 0.123713 | \n",
+ " 0.123713 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "