<jupyter_start><jupyter_text>Cancer Genomics Analysis<NAME>2020-12-02Objective: Explore the cancer genomic data (tumor/normal whole exome sequencing) from 50 patients that received the same type of treatment, half of whom responded.The data for this analysis from the Van Allen lab at the Dana-Farber Cancer Institute: Google Drive.**Overall question: Can we identify mutations that are associated with the treatment response?**<jupyter_code>#!/usr/bin/env python

""" Basic Utils """
import warnings
import glob
import os
from collections import Counter

""" Data Analysis and Visualization """
import pandas as pd
import numpy as np
import scipy.stats as stats
import seaborn as sns
import matplotlib.pyplot as plt
import math

INPUT_PATH = "mafs/"<jupyter_output><empty_output><jupyter_text>Step 1: Download the Dataset* import and store each MAF (Mutation Annotation Format) file into a pandas dataframe* create a dictionary to hold the dataframes; the dictionary will be indexed by the patient ID<jupyter_code>patient_maf_df_collection = {}

for filename in glob.glob(INPUT_PATH + '*.maf'):
    base = os.path.basename(filename)
    patient_id = base.split(".")[0]
    patient_maf_df_collection[patient_id] = pd.read_csv(filename, sep="\t")<jupyter_output><empty_output><jupyter_text>Step 2: Subset For Nonsynonymous Mutations* subset for mutations that are not of the Variant Classification “Silent”    * get rid of silent mutations<jupyter_code>subset_patient_maf_df_collection = {}

for key in patient_maf_df_collection.keys():
    df = patient_maf_df_collection[key]
    subset_df = df.loc[df['Variant_Classification'] != "Silent"]
    subset_patient_maf_df_collection[key] = subset_df.applymap(str)<jupyter_output><empty_output><jupyter_text>We can check the dataframe before and after to verify that the solution is correct.Notice that the rows corresponding to Silent mutations are absent after subsetting.<jupyter_code>print("=== Original Dataset ===")
print(patient_maf_df_collection['Patient-0']['Variant_Classification'].value_counts())
print("")
print("=== Subset Dataset ===")
print(subset_patient_maf_df_collection['Patient-0']['Variant_Classification'].value_counts())<jupyter_output>=== Original Dataset ===
Missense_Mutation    208
Silent                94
Nonsense_Mutation     20
Splice_Site           12
Name: Variant_Classification, dtype: int64

=== Subset Dataset ===
Missense_Mutation    208
Nonsense_Mutation     20
Splice_Site           12
Name: Variant_Classification, dtype: int64<jupyter_text>Step 3: Find the 15 Most Common Genes names are included in the column Hugo_Symbol and protein changes are stored in the column Protein_Change.<jupyter_code>def most_common(column_name):
    value_counter = []
    for k,v in subset_patient_maf_df_collection.items():
        a = v[column_name].values
        unique, counts = np.unique(a, return_counts=True)
        value_counter.append(Counter(dict(zip(unique, counts))))

    counter = sum(value_counter, Counter())
    return(counter)<jupyter_output><empty_output><jupyter_text>Most Common Gene Mutations<jupyter_code>gene_mutations = most_common('Hugo_Symbol')
most_common_gene_mutations = gene_mutations.most_common()[:15]
pd.DataFrame(np.array(most_common_gene_mutations).reshape(15,2), columns=('Gene Mutation', 'Frequency'))<jupyter_output><empty_output><jupyter_text>Most Common Protein Changes<jupyter_code>protein_changes = most_common('Protein_Change')
most_common_protein_changes = protein_changes.most_common()[0:16]
pd.DataFrame(np.array(most_common_protein_changes).reshape(16,2), columns=('Protein Change', 'Frequency'))<jupyter_output><empty_output><jupyter_text>Step 4: Perform Statistical Test (Gene Enrichment Analysis)* Perform a statistical test to explore if any mutated genes are enriched in patients who either responded or not.I will perform the Fisher's Exact Test by constructing following contingency table:                  Responder      Non-Responder                  Gene A mutant                          Gene A wild-type                  There are 50 patients in the study therefore the cells of the contingency table should add up to 50.**Overall Question: Are there any mutations that are associated with treatment response?**$H_{o}$: For a given gene, an observed imbalance in response vs non-response is not statistically significant.$H_{A}$: For a given gene, an observed imbalance in response vs non-response is significantly significant.<jupyter_code>mutated_gene_list = gene_mutations.keys()
patient_sample_list = subset_patient_maf_df_collection.keys()

sample_information = pd.read_csv('sample-information.tsv', sep="\t")
def get_contingency_table(patient_id,gene):
    gene_count = np.sum(subset_patient_maf_df_collection[patient_id]['Hugo_Symbol'] == gene)
    response = sample_information[sample_information['Patient_ID'] == patient_id]['Response'].values[0]

    if gene_count == 0 and response == 'Non-Responder':
        return(np.array([[0,0],[0,1]]))
    if gene_count > 0 and response == 'Non-Responder':
        return(np.array([[0,1],[0,0]]))
    if gene_count == 0 and response == 'Responder':
        return(np.array([[0,0],[1,0]]))
    if gene_count > 0 and response == 'Responder':
        return(np.array([[1,0],[0,0]]))
fisher_exact_pval = {}             # the fisher's exact test p-values
fisher_exact_oddsratio = {}        # the fisher's exact test odds radios
mutated_patients_per_gene = {}     # the number of mutated patients per gene

for gene in mutated_gene_list:
    contingency_table = np.array([[0,0],[0,0]])
    for patient_id in patient_sample_list:
        contingency_table += get_contingency_table(patient_id,gene)
    
    # sanity check -- if the cell values of the contingency table sum greater than 50 then display an error message
    if contingency_table.sum() != 50:
        print("ERROR!! gene name = " + gene)
        break
        
    number_of_mutated_patients = sum(contingency_table[0])
    oddsratio, pvalue = stats.fisher_exact(contingency_table)
    fisher_exact_pval[gene] = pvalue
    fisher_exact_oddsratio[gene] = oddsratio
    mutated_patients_per_gene[gene] = number_of_mutated_patients<jupyter_output><empty_output><jupyter_text>Bonferroni Multiple Test CorrectionWe're performing a Fisher's Exact Test for every individual per gene. There are 50 individual patients in the study. Therefore, there are 50 hypothesis tests per gene. We need to perform a multiple test correction.One option is the Bonferroni Multiple Test Correction. It is calculated by dividing the specified Type I error (alpha = 0.05) by the number of hypothesis tests. Bonferroni correction $=α/N$Genes having a p-value lower than the Bonferroni correction (the new alpha) will be determined as significant.Note: The Bonferroni correction tends to be a bit too conservative. An alternative multiple test correction is the False Discovery Rate.<jupyter_code>bonferroni_correction = 0.05 / 50
bonferroni_correction<jupyter_output><empty_output><jupyter_text>Step 5: Scatterplot of the Results (Gene Enrichment Analysis)Create a scatter plot of genes with the number of mutated patients on the x-axis and your results from question 4 on the y-axis.Can the axes be scaled or transformed in any way to improve readability?<jupyter_code>""" create a dataframe of the Fisher's Exact Test Results and Number of Mutated Patients"""
x_num = [0] * len(mutated_gene_list)
y_pval = [0] * len(mutated_gene_list)
gene_name = [0] * len(mutated_gene_list)

for i, gene in enumerate(fisher_exact_pval.keys(), 0):
    x_num[i] = mutated_patients_per_gene[gene]
    y_pval[i] = fisher_exact_pval[gene]
    gene_name[i] = gene

gene_enrichment_results = df = pd.DataFrame(list(zip(x_num, y_pval, gene_name)), columns =['Mutated_Patients', 'p-value', 'Gene_Name'])
print("There are {0} nonsynonymous mutations in the study.".format(str(gene_enrichment_results.shape[0])))
gene_enrichment_results.head()
""" Create the Scatter Plot """ 
fig, ax = plt.subplots(figsize=(16,10))
ax.tick_params(axis='both', which='major', labelsize=15)
ax = sns.scatterplot(x="Mutated_Patients", y="p-value", data=gene_enrichment_results)
ax.set_title('Fisher\'s Exact Test: P-Value vs Number of Mutated Patients',y=1.02, fontsize=20)
ax.set_xlabel('Number of Mutated Patients', fontsize=15, y=1.08)
ax.set_ylabel('p-value', fontsize=15)
ax.axhline(0.05, color = 'blue', ls='--')
ax.axhline(bonferroni_correction, color = 'red', ls='--')
bbox_props = dict(boxstyle="round", fc="w", ec="0.5", alpha=0.9)

ax.text(14, -0.2, "Note: There are over 6000 nonsynonymous mutations in the study but many points overlap due to \
having the same p-value and number of mutated patients", ha="center", va="center", size=14, bbox=bbox_props)
plt.legend(("0.05 alpha significance level", "Bonferroni Correction", "genes"), prop={'size': 14})<jupyter_output><empty_output><jupyter_text>Note: We can improve the visualization by increasing the radius of the data points to reflect the number of mutations that cluster. Here is an Example. Applying Log Transformation to p-valuesWe're interested in genes which have p-values < 0.05. These genes are currently positioned towards the bottom of the plot.We can improve readibility by transforming the y-axis onto a logarithmic scale: -log10(p-value). This will help us better visualize significant genes by effectively moving them to top of the plot.<jupyter_code>""" Apply Log Transformation to the p-values """
gene_enrichment_results['negative_log10_pval'] = -np.log10(gene_enrichment_results['p-value']) + 0

""" Display results """
gene_enrichment_results.sort_values(by=['negative_log10_pval'],  ascending=False)[0:9]<jupyter_output><empty_output><jupyter_text>Most Significantly Enriched GeneThe most significantly enriched gene is ERCC2.<jupyter_code>""" Create the Scatter Plot """ 
fig, ax2 = plt.subplots(figsize=(16,10))
ax2.tick_params(axis='both', which='major', labelsize=15)
ax2 = sns.scatterplot(x="Mutated_Patients", y="negative_log10_pval", data=gene_enrichment_results)
ax2.set_title('Fisher\'s Exact Test: P-Value vs Number of Mutated Patients',y=1.02, fontsize=18)
ax2.set_xlabel('Number of Mutated Patients', fontsize=14, y=1.08)
ax2.set_ylabel('-log10(p-value)', fontsize=14)
ax2.axhline(-math.log10(0.05), color = 'blue', ls='--')
ax2.axhline(-math.log10(0.01), color = 'red', ls='--')

ax2.text(14, -0.5, "Note: There are over 6000 nonsynonymous mutations in the study but many points overlap due to \
having the same p-value and the same number of mutated patients", ha="center", va="center", size=14, bbox=bbox_props)
plt.legend(("0.05 alpha significance level", "0.01 alpha significance level", "genes"), prop={'size': 14})

""" gene annotation """
ax2.annotate('ERCC2', xy = (9, -math.log10(0.001631)), xytext = (11, 2.6), arrowprops = dict(facecolor = 'black', shrink = 0.05, width=3), fontsize = 14 )
ax2.annotate('HERC1', xy = (6, -math.log10(0.022290)), xytext = (8, 1.7), arrowprops = dict(facecolor = 'black', shrink = 0.05, width=3), fontsize = 14 )
ax2.annotate('AKAP9', xy = (6, -math.log10(0.022290)), xytext = (3, 1.7), arrowprops = dict(facecolor = 'black', shrink = 0.05, width=3), fontsize = 14 )
ax2.annotate('HECTD1', xy = (6, -math.log10(0.022290)), xytext = (3, 2.1), arrowprops = dict(facecolor = 'black', shrink = 0.05, width=3), fontsize = 14 )
ax2.annotate('MROH2B', xy = (6, -math.log10(0.022290)), xytext = (7, 2.1), arrowprops = dict(facecolor = 'black', shrink = 0.05, width=3), fontsize = 14 )
ax2.annotate('MACF1', xy = (6, -math.log10(0.022290)), xytext = (7.5, 1.8), arrowprops = dict(facecolor = 'black', shrink = 0.05, width=3), fontsize = 14 )
ax2.annotate('KMT2C', xy = (8, -math.log10(0.048797)), xytext = (10, 1.5), arrowprops = dict(facecolor = 'black', shrink = 0.05, width=3), fontsize = 14 )
bbox_props = dict(boxstyle="round", fc="w", ec="0.5", alpha=0.9)<jupyter_output><empty_output><jupyter_text>Step 6: Nonsynonymous Mutations Per MegabaseHow many samples are wild-type versus mutant with respect to the most significantly enriched gene from Question 4?<jupyter_code>gene = 'ERCC2'
gene_phenotype = []
contingency_table = np.array([[0,0],[0,0]])
for patient_id in patient_sample_list:
    contingency_table += get_contingency_table(patient_id,gene)
    status = 'wild-type' if np.sum(subset_patient_maf_df_collection[patient_id]['Hugo_Symbol'] == gene) == 0 else 'mutant'
    gene_phenotype.append(status)

contingency_table
gene_phenotype.count('mutant')<jupyter_output><empty_output><jupyter_text>Responder      Non-Responder                  Gene A mutant      9      0              Gene A wild-type      16      25        41 samples are wild-type and 9 samples are mutant with respect to the most significantly enriched gene ERCC2  Examine Nonsynonymous Mutatations Per Megabase in Mutant vs Wild-Type SamplesPlot the number of nonsynonymous mutations per megabase in the mutant vs. wild-type samples. Y-axis is nonsynonymous mutations per megabase. X-axis is categorical: mutant vs. wild-type samples (two different columns)<jupyter_code>df = sample_information[['Patient_ID', 'Response', 'Nonsynonymous_mutations_per_Mb']]
df.insert(loc=3, column='gene_phenotype', value=gene_phenotype)
df.head()
fig, ax3 = plt.subplots(figsize=(16,10))

""" create categorical scatter plot """ 
ax3 = sns.swarmplot(x="gene_phenotype", y="Nonsynonymous_mutations_per_Mb", hue="Response",
              palette=["r", "c"], data=df)

ax3.set_title('Nonsynonymous Mutations in Patient Samples with Respect to the ERCC2 Gene',y=1.02, fontsize=18)
ax3.set_xlabel('Sample Phenotype with respect to ERCC2', fontsize=14, y=1.08)
ax3.set_ylabel('Nonsynonymous Mutations Per Mb', fontsize=14)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11533974919801691, 'fraction_numerical': 0.030183727034120734, 'mean_word_length': 5.433864915572233, 'pattern_counts': {'":': 0, '<': 58, '<?xml version=': 0, '>': 59, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '633757', 'license': 'MIT', 'n_tokens_mistral': 4737, 'n_tokens_neox': 4302, 'n_words': 1246, 'path': 'Dana_Farber_LiuLab_SMoijueh.ipynb', 'repo_name': 'smoijueh/Cancer-Genomic-Analysis', 'size': '13724'}