multiple_mutation_analyses.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import argparse\n",
    "import os\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import numpy\n",
    "\n",
    "from collections import Counter\n",
    "\n",
    "pathPrefix = '/Users/friedman/Desktop/mnt'\n",
    "\n",
    "sys.path.append(pathPrefix + '/ifs/work/taylorlab/friedman/myUtils')\n",
    "import analysis_utils \n",
    "import mutationSigUtils \n",
    "import maf_analysis_utils\n",
    "import clonality_analysis_util\n",
    "import re\n",
    "from scipy import stats\n",
    "import math\n",
    "import msi_analysis_utils"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def count_number_of_times_allele_is_in_multiplet(maf, allele):\n",
    "    casesWithAllele = set(maf[maf['allele'] == allele]['Tumor_Sample_Barcode'])\n",
    "    if len(casesWithAllele) == 0:\n",
    "        return None\n",
    "    alleleMaf = maf[maf['Tumor_Sample_Barcode'].isin(casesWithAllele)]\n",
    "    gene = allele.split('_')[0]\n",
    "    otherOncogenicGeneAlleles = alleleMaf[(alleleMaf['oncogenic'].notnull()) & (alleleMaf['Hugo_Symbol'] == gene) &(alleleMaf['allele'] != allele)]\n",
    "    return (1.0*len(set(otherOncogenicGeneAlleles['Tumor_Sample_Barcode'])))/len(casesWithAllele)\n",
    "\n",
    "#lambda function that marks whether an allele is part of a multiplet\n",
    "def is_allele_multiplet(maf, allele, tsb):\n",
    "    gene = allele.split('_')[0]\n",
    "    if maf[(maf['allele'] != allele) & (maf['oncogenic'].notnull()) & (maf['Hugo_Symbol'] == gene) & (maf['Tumor_Sample_Barcode'] == tsb)].shape[0] > 0:\n",
    "        return 1\n",
    "    else:\n",
    "        return 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "allImpactMuts = analysis_utils.load_in_df_with_progress(filePath = pathPrefix + '/ifs/work/taylorlab/friedman/hypermutationAnalysisProj/projectDataAndConfigFiles/all_impact_mutations_annotated_cohort.maf', nLinesFile = 275000)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "allImpactMuts = maf_analysis_utils.fix_mll_genes(allImpactMuts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "poleStopGainRecurrentAlleles = set(['PIK3R1_p.R348*', 'TP53_p.R213*', 'ASXL2_p.R312*', 'FBXW7_p.R658*', 'SMAD4_p.E520*', 'CASP8_p.E239*', 'ASXL2_p.E45*', 'APC_p.R2237*', 'ARID1A_p.R1989*', 'APC_p.R1114*', 'MSH2_p.E580*', 'ATM_p.R250*', 'ARID1A_p.E896*', 'RASA1_p.E429*', 'NF1_p.R2450*'])\n",
    "msiFrameshiftRecurrentAlleles = set(['ARID1A_p.Q758', 'SMO_p.P694', 'PTPRT_p.P1094', 'SMAD4_p.S32', 'HLA-B_p.P209', 'SOX9_p.P346', 'ARID1A_p.D1850', 'ZFHX3_p.E1892', 'TCF7L2_p.K468', 'MSH3_p.K383', 'CTCF_p.T204', 'AXIN2_p.E405', 'NSD1_p.M1531', 'MRE11A_p.N511', 'ERF_p.G299', 'KMT2B_p.G1879', 'MLL2_p.G1235', 'RNF43_p.R117', 'RNF43_p.G659', 'JAK1_p.K860', 'GLI1_p.G274', 'ARID1A_p.G276', 'ARID1A_p.F2141', 'B2M_p.L13', 'PDCD1_p.T36', 'CIC_p.P1248', 'KDM5A_p.G1200', 'PTEN_p.T321', 'ARID1A_p.Q372', 'IRS1_p.P368', 'PPM1D_p.N512', 'EP300_p.M1470', 'SOX9_p.V306', 'KMT2D_p.P2354', 'CD79A_p.R131', 'PTCH1_p.P1202', 'MSH6_p.F1088', 'APC_p.E1464', 'CREBBP_p.I1084', 'ARID5B_p.K967', 'KMT2D_p.V3089', 'ARID1A_p.T294', 'ARID1A_p.K1072', 'B2M_p.T93', 'SPRED1_p.D79', 'TGFBR2_p.K153', 'MLL2_p.P647', 'KMT2C_p.F4496', 'TCF7L2_p.L200', 'ARID1A_p.P224', 'FBXW7_p.S668', 'AXIN2_p.H661', 'ARAF_p.R255', 'ARID1A_p.P1326', 'PTEN_p.L265', 'B2M_p.V69', 'INPPL1_p.R1156', 'RPS6KA4_p.S772', 'SOX17_p.P328', 'PPM1D_p.L450', 'CARD11_p.R555', 'APC_p.E1554', 'AMER1_p.F173', 'NF1_p.Y628'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5.74382538771 percent done\n",
      "11.4876507754 percent done\n",
      "17.2314761631 percent done\n",
      "22.9753015508 percent done\n",
      "28.7191269385 percent done\n",
      "34.4629523262 percent done\n",
      "40.206777714 percent done\n",
      "45.9506031017 percent done\n",
      "51.6944284894 percent done\n",
      "57.4382538771 percent done\n",
      "63.1820792648 percent done\n",
      "68.9259046525 percent done\n",
      "74.6697300402 percent done\n",
      "80.4135554279 percent done\n",
      "86.1573808156 percent done\n",
      "91.9012062033 percent done\n",
      "97.645031591 percent done\n"
     ]
    }
   ],
   "source": [
    "#FIX MSI ALLELES\n",
    "reload(msi_analysis_utils)\n",
    "msiSummary = pd.read_table('/Users/friedman/Desktop/impact_microsatellites.txt')\n",
    "allImpactMuts['allele'] = allImpactMuts.apply(lambda row: str(row['Hugo_Symbol']) + '_' + str(row['HGVSp_Short']), axis=1)\n",
    "msiSummary['allele'] = msiSummary.apply(lambda row: str(row['Hugo_Symbol']) + '_' + str(row['HGVSp_Short']), axis=1)\n",
    "neverObservedSites, msiSitesToNameMapping, mafMsiSiteToNameMapping =  msi_analysis_utils.standardize_allele_names(msiSummary, allImpactMuts)\n",
    "\n",
    "allImpactMuts['correctedAllele'] = allImpactMuts['allele'].apply(lambda x:\n",
    "                                                                 mafMsiSiteToNameMapping[x] if x in mafMsiSiteToNameMapping else None)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "impactSigs = pd.read_table(pathPrefix + '/ifs/work/taylorlab/friedman/hypermutationAnalysisProj/projectDataAndConfigFiles/signatures_from_unfiltered_maf.txt')\n",
    "impactSigs = mutationSigUtils.merge_signature_columns(impactSigs)\n",
    "impactSigs['dominantSignature'] = impactSigs.apply(lambda row: \n",
    "        mutationSigUtils.get_dominant_signature(row.to_dict(), cols=None, prefix='mean', notEnoughMuts= True), axis=1)\n",
    "\n",
    "casesWithMSISignature = set(impactSigs[(impactSigs['dominantSignature'] == 'mean_MMR')\n",
    "                                      |((impactSigs['Nmut_Mb'] >= 30) & (impactSigs['dominantSignature'].isin(set(['mean_1']))))\n",
    "                                        ]['Tumor_Sample_Barcode'])\n",
    "\n",
    "casesWithPOLESignature = set(impactSigs[(impactSigs['dominantSignature'] == 'mean_10') & (impactSigs['Nmut_Mb'] > 30)]['Tumor_Sample_Barcode'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "msiPoleMuts = allImpactMuts[allImpactMuts['Tumor_Sample_Barcode'].isin(casesWithMSISignature | casesWithPOLESignature)]\n",
    "#alleleCounts = Counter(msiPoleMuts['allele'])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda2/lib/python2.7/site-packages/ipykernel_launcher.py:2: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
      "  \n"
     ]
    }
   ],
   "source": [
    "#THIS LINE TAKES A REALLY LONG TIME\n",
    "msiPoleMuts['multipletPresentInCase'] = msiPoleMuts.apply(lambda row:\n",
    "                is_allele_multiplet(msiPoleMuts, row['allele'], row['Tumor_Sample_Barcode']), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "tumorSuppresors = set(['ERRFI1', 'ASXL2', 'PMAIP1', 'ACTG1', 'SUFU', 'FBXO11', 'MEN1', 'FAM58A', 'B2M', 'RB1', 'DUSP22', 'SESN1', 'GPS2', 'RAD51D', 'SMG1', 'CDC73', 'MAP3K1', 'SMARCB1', 'INPP4B', 'PARK2', 'SMAD4', 'CBFB', 'CDH1', 'PPP6C', 'SETDB1', 'SETDB2', 'NF2', 'CDKN2B', 'CDKN2C', 'CDKN2A', 'DDX3X', 'PIK3R1', 'BARD1', 'PDS5B', 'KLF4', 'SPRED1', 'VHL', 'SMAD2', 'PMS1', 'PMS2', 'SETD2', 'GATA3', 'TBL1XR1', 'MUTYH', 'SOCS1', 'FAM175A', 'ROBO1', 'ARID1B', 'ARID1A', 'TCF7L2', 'STK11', 'FOXA1', 'PTEN', 'FAT1', 'FAS', 'CYLD', 'MAX', 'SH2D1A', 'APC', 'NTHL1', 'CTCF', 'KDM5C', 'KMT2C', 'ZFHX3', 'FOXP1', 'PIGA', 'CDKN1B', 'CDKN1A', 'FUBP1', 'MSH2', 'ID3', 'TNFRSF14', 'TRAF3', 'EP400', 'BRIP1', 'ARID4A', 'ARID4B', 'XRCC2', 'DAXX', 'SDHAF2', 'ASXL1', 'AMER1', 'RASA1', 'EGR1', 'MST1', 'SOX17', 'RUNX1', 'PIK3R3', 'NCOR1', 'NF1', 'JAK1', 'PTPRD', 'CHEK2', 'CHEK1', 'SMC1A', 'TMEM127', 'STAG1', 'RAD51', 'TCF3', 'STAG2', 'ARID2', 'RAD50', 'RNF43', 'PARP1', 'BLM', 'CUX1', 'RECQL', 'RAD21', 'PTPN2', 'PTPN1', 'SLX4', 'INHA', 'PAX5', 'IRF1', 'TP53', 'HLA-A', 'IRF8', 'CBL', 'TOP1', 'SHQ1', 'PRDM1', 'NSD1', 'ATXN2', 'CREBBP', 'HDAC4', 'SESN2', 'PPP2R1A', 'EPHA7', 'ATM', 'EPHA3', 'POT1', 'SMAD3', 'MOB3B', 'TBX3', 'POLE', 'ATR', 'FANCD2', 'FH', 'BCORL1', 'SOX9', 'IKZF3', 'TSC1', 'TP63', 'MRE11A', 'SDHC', 'BTG1', 'POLD1', 'CIITA', 'SMC3', 'SAMHD1', 'RTEL1', 'ECT2L', 'PIK3R2', 'CRBN', 'FANCC', 'NBN', 'FANCA', 'HLA-B', 'RECQL4', 'DUSP4', 'ERCC2', 'FBXW7', 'TGFBR2', 'TGFBR1', 'MSH3', 'RBM15', 'TET1', 'TET3', 'SESN3', 'MGA', 'LTB', 'FOXL2', 'SH2B3', 'BCOR', 'HIST1H1D', 'ATRX', 'EP300', 'RAD51C', 'RAD51B', 'HIST1H1B', 'TNFAIP3', 'DICER1', 'ARID5B', 'LATS2', 'FOXO1', 'KEAP1', 'EZH2', 'SP140', 'NKX3-1', 'PBRM1', 'PALB2', 'CIC', 'BRCA1', 'DTX1', 'FLCN', 'SPEN', 'CD58', 'ERCC3', 'ERCC4', 'MSH6', 'BCL11B', 'BMPR1A', 'ERF', 'BRCA2', 'NOTCH2', 'EED', 'MITF', 'ELF3', 'SMARCA4', 'BBC3', 'ANKRD11', 'CEBPA', 'BCL2L11', 'AXIN2', 'AXIN1', 'CDK12', 'ESCO2', 'MLH1', 'SDHB', 'MED12', 'HNF1A', 'RYBP', 'ATP6V1B2', 'DNMT3B', 'KMT2B', 'KMT2A', 'DNMT3A', 'NFKBIA', 'TRAF5', 'KMT2D', 'SPOP', 'RBM10', 'P2RY8', 'TP53BP1', 'TSC2', 'KDM6A', 'EPCAM', 'PHOX2B', 'NPM1', 'BCL10', 'LATS1', 'HOXB13', 'ARID3A', 'PTPRT', 'PTPRS', 'INPPL1', 'NOTCH4', 'TET2', 'NOTCH1', 'CASP8', 'NOTCH3', 'GRIN2A', 'MAP2K4', 'WT1', 'BACH2', 'SDHA', 'BAP1', 'PTCH1', 'SDHD'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda2/lib/python2.7/site-packages/ipykernel_launcher.py:3: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
      "  This is separate from the ipykernel package so we can avoid doing imports until\n"
     ]
    }
   ],
   "source": [
    "msiPoleMuts['alleleClass'] = msiPoleMuts.apply(lambda row: 'recurrentPoleStopGain' if row['allele'] in poleStopGainRecurrentAlleles\n",
    "                        else 'recurrentMsiIndel' if row['correctedAllele'] in msiFrameshiftRecurrentAlleles\n",
    "                        else 'notRecurrentTSG' if row['Hugo_Symbol'] in tumorSuppresors else 'Oncogene', axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "500 49865\n",
      "1000 49865\n",
      "1500 49865\n",
      "2000 49865\n",
      "2500 49865\n",
      "3000 49865\n",
      "3500 49865\n",
      "4000 49865\n",
      "4500 49865\n",
      "5000 49865\n",
      "5500 49865\n",
      "6000 49865\n",
      "6500 49865\n",
      "7000 49865\n",
      "7500 49865\n",
      "8000 49865\n",
      "8500 49865\n",
      "9000 49865\n",
      "9500 49865\n",
      "10000 49865\n",
      "10500 49865\n",
      "11000 49865\n",
      "11500 49865\n",
      "12000 49865\n",
      "12500 49865\n",
      "13000 49865\n",
      "13500 49865\n",
      "14000 49865\n",
      "14500 49865\n",
      "15000 49865\n",
      "15500 49865\n",
      "16000 49865\n",
      "16500 49865\n",
      "17000 49865\n",
      "17500 49865\n",
      "18000 49865\n",
      "18500 49865\n",
      "19000 49865\n",
      "19500 49865\n",
      "20000 49865\n",
      "20500 49865\n",
      "21000 49865\n",
      "21500 49865\n",
      "22000 49865\n",
      "22500 49865\n",
      "23000 49865\n",
      "23500 49865\n",
      "24000 49865\n",
      "24500 49865\n",
      "25000 49865\n",
      "25500 49865\n",
      "26000 49865\n",
      "26500 49865\n",
      "27000 49865\n",
      "27500 49865\n",
      "28000 49865\n",
      "28500 49865\n",
      "29000 49865\n",
      "29500 49865\n",
      "30000 49865\n",
      "30500 49865\n",
      "31000 49865\n",
      "31500 49865\n",
      "32000 49865\n",
      "32500 49865\n",
      "33000 49865\n",
      "33500 49865\n",
      "34000 49865\n",
      "34500 49865\n",
      "35000 49865\n",
      "35500 49865\n",
      "36000 49865\n",
      "36500 49865\n",
      "37000 49865\n",
      "37500 49865\n",
      "38000 49865\n",
      "38500 49865\n",
      "39000 49865\n",
      "39500 49865\n",
      "40000 49865\n",
      "40500 49865\n",
      "41000 49865\n",
      "41500 49865\n",
      "42000 49865\n",
      "42500 49865\n",
      "43000 49865\n",
      "43500 49865\n",
      "44000 49865\n",
      "44500 49865\n",
      "45000 49865\n",
      "45500 49865\n",
      "46000 49865\n",
      "46500 49865\n",
      "47000 49865\n",
      "47500 49865\n",
      "48000 49865\n",
      "48500 49865\n",
      "49000 49865\n",
      "49500 49865\n"
     ]
    }
   ],
   "source": [
    "#msiPoleMuts['multipletPresentInCase']\n",
    "\n",
    "alleleCounts = dict(msiPoleMuts['allele'].value_counts())\n",
    "multipletAlleleCounts = dict(msiPoleMuts[msiPoleMuts['multipletPresentInCase'] == 1]['allele'].value_counts())\n",
    "alleleClassMapping = dict(zip(msiPoleMuts['allele'], msiPoleMuts['alleleClass']))\n",
    "\n",
    "listOfDicts = []\n",
    "cntr = 0\n",
    "for allele in set(msiPoleMuts['allele']):\n",
    "    cntr += 1\n",
    "    if cntr%500 == 0:\n",
    "        print cntr, len(set(msiPoleMuts['allele']))\n",
    "    \n",
    "    nAlleles = alleleCounts[allele]\n",
    "    nMultiplet = 0\n",
    "    if allele in multipletAlleleCounts: nMultiplet = multipletAlleleCounts[allele]\n",
    "    fracMultiplet = 1.0*nMultiplet/nAlleles\n",
    "    #alleleMaf = msiPoleMuts[msiPoleMuts['allele'] == allele]\n",
    "    listOfDicts.append({'nOccurencesAllele': nAlleles, 'fracOccurencesMultiplet': fracMultiplet,\n",
    "                       'nMultiplet':nMultiplet, 'allele':allele, 'alleleClass': alleleClassMapping[allele]})\n",
    "    \n",
    "df = pd.DataFrame(listOfDicts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['orderingVal'] = df['alleleClass'].apply(lambda x: 1 if x == 'recurrentPoleStopGain' else 2 if x == 'recurrentMsiIndel'\n",
    "                        else 3 if x == 'notRecurrentTSG' else 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.to_csv('~/Desktop/WORK/dataForLocalPlotting/alleleRecurrences.tsv', index=False, sep='\\t')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Counter({'+': 275171})"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Counter(allImpactMuts['Strand'])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}