msi_mut_features_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#returns the name of the gene and left most position spanning the msi indel\n",
    "#uses re stuff\n",
    "def get_left_aligned_allele_name(hgvsNames):\n",
    "    \n",
    "    positions = []\n",
    "    geneName = hgvsNames[0].split('_p.')[0]\n",
    "    for entry in hgvsNames:\n",
    "        if len(entry.split('_p.')) == 2: #ignore weirdly formatted hgvs names\n",
    "            variantNotation = entry.split('_p.')[1]\n",
    "            number = variantNotation[1:]\n",
    "            refAA = variantNotation[0]\n",
    "            position = re.match('\\d*', number).group(0)\n",
    "            positions.append((position, refAA))\n",
    "    \n",
    "    if len(positions) == 0: return None #if all the hgvs names were ill formatted return None\n",
    "    \n",
    "    minEntry = sorted(positions)[0] #this is a sorted list of tuples the first thing is the position second is the reference aa\n",
    "    return geneName + '_p.' + str(minEntry[1]) + str(minEntry[0])\n",
    "\n",
    "#collapses all indels within 1 bp of each other for a start to be at the same location/name for matching\n",
    "def standardize_allele_names(msiLengthInfo, observedMuts):\n",
    "    \n",
    "    neverObservedSites = set([]) #all the names of sites from criags msi file we cant match with the real maf\n",
    "    msiSitesToNameMapping = {} #a dictionary mapping each msi site allele from craigs file to its corrected name\n",
    "    mafMsiSiteToNameMapping = {} #a dictionary mapping each msi site allele from the maf to its corrected name\n",
    "    \n",
    "    cntr = 0.0\n",
    "    for hgvs in set(msiLengthInfo['allele']):\n",
    "        \n",
    "        cntr += 1\n",
    "        if cntr%500 == 0: print 100*(cntr/len(set(msiLengthInfo['allele']))), 'percent done'\n",
    "        \n",
    "        startPos = msiLengthInfo[msiLengthInfo['allele'] == hgvs]['Start_Position']\n",
    "        \n",
    "        \n",
    "        if startPos.shape[0] == 1:\n",
    "            #we want all names given to indels near (within 1position) of the start position of the MSI site in Craig's file\n",
    "            putativeVariantNames = list(set(observedMuts[(abs(observedMuts['Start_Position'] - int(startPos)) < 2)\n",
    "                                                    & (observedMuts['Variant_Type'].isin(set(['INS', 'DEL'])))]['allele']))\n",
    "            \n",
    "            if len(putativeVariantNames) == 0:\n",
    "                neverObservedSites.add(hgvs) #if it cant be matched in the MAF we add it to never observed sites\n",
    "                #note some of there are likely to be actually matched but missed by my method\n",
    "            else:\n",
    "                trueVariantName = get_left_aligned_allele_name(putativeVariantNames)\n",
    "                \n",
    "                #NOW WE PROPERLY create the mappings\n",
    "                for putativeVariantName in putativeVariantNames:\n",
    "                    mafMsiSiteToNameMapping[putativeVariantName] = trueVariantName\n",
    "                msiSitesToNameMapping[hgvs] = trueVariantName\n",
    "        else:\n",
    "            pass #ignore variants with multiple start position in the msi info file\n",
    "        \n",
    "    return neverObservedSites, msiSitesToNameMapping, mafMsiSiteToNameMapping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda2/lib/python2.7/site-packages/IPython/core/interactiveshell.py:2714: DtypeWarning: Columns (0,76,97) have mixed types. Specify dtype option on import or set low_memory=False.\n",
      "  interactivity=interactivity, compiler=compiler, result=result)\n"
     ]
    }
   ],
   "source": [
    "msiSummary = pd.read_table('/Users/friedman/Desktop/impact_microsatellites.txt')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda2/lib/python2.7/site-packages/IPython/core/interactiveshell.py:2878: DtypeWarning: Columns (36) have mixed types. Specify dtype option on import or set low_memory=False.\n",
      "  exec(code_obj, self.user_global_ns, self.user_ns)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.64 percent done\n",
      "7.27 percent done\n",
      "10.91 percent done\n",
      "14.55 percent done\n",
      "18.18 percent done\n",
      "21.82 percent done\n",
      "25.45 percent done\n",
      "29.09 percent done\n",
      "32.73 percent done\n",
      "36.36 percent done\n",
      "40.0 percent done\n",
      "43.64 percent done\n",
      "47.27 percent done\n",
      "50.91 percent done\n",
      "54.55 percent done\n",
      "58.18 percent done\n",
      "61.82 percent done\n",
      "65.45 percent done\n",
      "69.09 percent done\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda2/lib/python2.7/site-packages/IPython/core/interactiveshell.py:2878: DtypeWarning: Columns (85) have mixed types. Specify dtype option on import or set low_memory=False.\n",
      "  exec(code_obj, self.user_global_ns, self.user_ns)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "72.73 percent done\n",
      "76.36 percent done\n",
      "80.0 percent done\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda2/lib/python2.7/site-packages/IPython/core/interactiveshell.py:2878: DtypeWarning: Columns (36,39,136) have mixed types. Specify dtype option on import or set low_memory=False.\n",
      "  exec(code_obj, self.user_global_ns, self.user_ns)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "83.64 percent done\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda2/lib/python2.7/site-packages/IPython/core/interactiveshell.py:2878: DtypeWarning: Columns (13,36) have mixed types. Specify dtype option on import or set low_memory=False.\n",
      "  exec(code_obj, self.user_global_ns, self.user_ns)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "87.27 percent done\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda2/lib/python2.7/site-packages/IPython/core/interactiveshell.py:2878: DtypeWarning: Columns (13) have mixed types. Specify dtype option on import or set low_memory=False.\n",
      "  exec(code_obj, self.user_global_ns, self.user_ns)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "90.91 percent done\n",
      "94.55 percent done\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda2/lib/python2.7/site-packages/IPython/core/interactiveshell.py:2878: DtypeWarning: Columns (39) have mixed types. Specify dtype option on import or set low_memory=False.\n",
      "  exec(code_obj, self.user_global_ns, self.user_ns)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "98.18 percent done\n",
      "101.82 percent done\n",
      "performing big concat then returning\n"
     ]
    }
   ],
   "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": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "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": [
    "neverObservedSites, msiSitesToNameMapping, mafMsiSiteToNameMapping =  standardize_allele_names(msiSummary, allImpactMuts)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "allImpactMuts['correctedAllele'] = allImpactMuts['allele'].apply(lambda x:\n",
    "                                                                 mafMsiSiteToNameMapping[x] if x in mafMsiSiteToNameMapping else None)\n",
    "msiSummary['correctedAllele'] = msiSummary['allele'].apply(lambda x:\n",
    "                                                           msiSitesToNameMapping[x] if x in msiSitesToNameMapping else None)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "def summarize_mutation_counts_and_msi_allele_relationship(impactMuts, msiInfo):\n",
    "    listOfDicts = []\n",
    "    for correctedAllele in set(msiInfo['correctedAllele']):\n",
    "        msiInfoAlleleDf = msiInfo[msiInfo['correctedAllele'] == correctedAllele]\n",
    "        if msiInfoAlleleDf.shape[0] > 0:\n",
    "            nRepeats = msiInfoAlleleDf['repeat_times'].iloc[0]\n",
    "            gene = msiInfoAlleleDf[msiInfoAlleleDf['correctedAllele'] == correctedAllele]['Hugo_Symbol'].iloc[0]\n",
    "            basePair = msiInfoAlleleDf[msiInfoAlleleDf['correctedAllele'] == correctedAllele]['Tumor_Seq_Allele2'].iloc[0]\n",
    "            nOccurences = impactMuts[impactMuts['correctedAllele'] == correctedAllele].shape[0]\n",
    "            listOfDicts.append({'Hugo_Symbol': gene, 'repeat_length': nRepeats,\n",
    "                                'nOccurences': nOccurences, 'allele': correctedAllele, 'basePair': basePair})\n",
    "    return pd.DataFrame(listOfDicts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "occurenceSummaryDf = summarize_mutation_counts_and_msi_allele_relationship(allImpactMuts, msiSummary)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [],
   "source": [
    "occurenceSummaryDf['repeatLengthsLabel'] = occurenceSummaryDf['repeat_length'].apply(lambda x: '_' + str(x) if x < 10 else '>=10')\n",
    "tumorSuppressors = 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'])\n",
    "occurenceSummaryDf['isTumorSuppresor'] = occurenceSummaryDf['Hugo_Symbol'].apply(lambda x: True if x in tumorSuppressors else False)\n",
    "occurenceSummaryDf['label'] = occurenceSummaryDf.apply(lambda row: row['allele'] if row['nOccurences'] >= 75 else None, axis=1)\n",
    "occurenceSummaryDf['basePair'] = occurenceSummaryDf['basePair'].apply(lambda x: x if x in set(['A', 'C', 'T', 'G']) else '_other')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [],
   "source": [
    "#orderingValDict = {'A': .1, 'T':.2, 'C':.3, 'G':.4, '_other': None}\n",
    "#occurenceSummaryDf['basePairAndLength'] = occurenceSummaryDf.apply(lambda row: row['repeatLengthsLabel'] + '_' + row['basePair'], axis=1)\n",
    "#occurenceSummaryDf['orderingVal'] = occurenceSummaryDf.apply(lambda row: row['repeat_length'] + orderingValDict[row['basePair']] if orderingValDict[row['basePair']] != None else None, axis=1)\n",
    "\n",
    "occurenceSummaryDf = occurenceSummaryDf[occurenceSummaryDf['nOccurences'].notnull()]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [],
   "source": [
    "occurenceSummaryDf.to_csv('/Users/friedman/Desktop/WORK/msiSiteOccurenceInfo.tsv', index=False, sep='\\t')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**SUMMARIZE DOUBLE HIT MUTATIONS** <br/>\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [],
   "source": [
    "def summarize_double_mut_msi_character(allelesAndCounts, msiInfo):\n",
    "    listOfDicts = []\n",
    "    for correctedAllele, count in allelesAndCounts.items():\n",
    "        msiInfoAlleleDf = msiInfo[msiInfo['correctedAllele'] == correctedAllele]\n",
    "        if msiInfoAlleleDf.shape[0] > 0:\n",
    "            nRepeats = msiInfoAlleleDf['repeat_times'].iloc[0]\n",
    "            gene = msiInfoAlleleDf[msiInfoAlleleDf['correctedAllele'] == correctedAllele]['Hugo_Symbol'].iloc[0]\n",
    "            basePair = msiInfoAlleleDf[msiInfoAlleleDf['correctedAllele'] == correctedAllele]['Tumor_Seq_Allele2'].iloc[0]\n",
    "        listOfDicts.append({'Hugo_Symbol': gene, 'repeat_length': nRepeats,\n",
    "                            'allele': correctedAllele, 'basePair': basePair, 'nOccurences': count})\n",
    "    df = pd.DataFrame(listOfDicts)\n",
    "    return df\n",
    "        \n",
    "        \n",
    "        \n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "msiAlleleNames = {}\n",
    "observedDoubleMutsEndoColoHyper = Counter({'RNF43_p.G659Vfs*41': 44, 'TCF7L2_p.K468Sfs*23': 14, 'MSH3_p.K383Rfs*32': 12, 'JAK1_p.K860Nfs*16': 10, 'PTEN_p.K267Rfs*9': 6, 'TGFBR2_p.K153Afs*3': 5, 'INPPL1_p.R1156Gfs*46': 4, 'ASXL1_p.G645Vfs*58': 3, 'KMT2C_p.X2481_splice': 2, 'TP53BP1_p.X1769_splice': 2, 'B2M_p.V69Wfs*34': 2, 'APC_p.T1556Nfs*3': 2, 'KMT2D_p.X2683_splice': 2, 'FLCN_p.H429Tfs*39': 2, 'ERF_p.G299Efs*12': 1, 'HLA-A_p.W157*': 1, 'EP300_p.H2324Tfs*29': 1, 'TGFBR2_p.K153Sfs*35': 1, 'APC_p.R856Nfs*6': 1, 'NF1_p.R1362*': 1, 'FUBP1_p.S11Lfs*43': 1, 'CYLD_p.N719Mfs*13': 1, 'ASXL2_p.X135_splice': 1, 'SOCS1_p.P36Vfs*75': 1, 'CIC_p.P1248Hfs*54': 1, 'NF1_p.N78Ifs*7': 1, 'ARID1A_p.P224Rfs*8': 1, 'ATM_p.K2811Sfs*46': 1, 'TSC1_p.R420Gfs*20': 1, 'KMT2C_p.K2797Rfs*26': 1, 'CREBBP_p.X608_splice': 1, 'MSH6_p.F1088Lfs*5': 1, 'APC_p.R2204*': 1, 'TP53_p.V73Wfs*50': 1, 'PIK3R1_p.R348*': 1, 'RAD50_p.K722Rfs*14': 1, 'PTEN_p.R130Q': 1, 'MSH6_p.E1322*': 1, 'MGA_p.E1249Rfs*42': 1, 'HNF1A_p.P291Qfs*51': 1, 'HLA-B_p.X25_splice': 1, 'NF1_p.I679Dfs*21': 1, 'ARID1A_p.D1850Tfs*33': 1, 'SMAD4_p.R361C': 1, 'CTNNB1_p.S45F': 1, 'B2M_p.L15Ffs*41': 1, 'PTEN_p.X267_splice': 1, 'DICER1_p.X504_splice': 1, 'PTEN_p.R173C': 1, 'DICER1_p.X301_splice': 1, 'JAK2_*505*': 1, 'B2M_p.T93Lfs*10': 1, 'NF1_p.Y628Tfs*3': 1, 'PTPRD_p.X1751_splice': 1, 'PTEN_p.D268Gfs*30': 1, 'TOP1_p.R140*': 1})\n",
    "for val in observedDoubleMutsEndoColoHyper.keys():\n",
    "    if val in mafMsiSiteToNameMapping:\n",
    "        msiAlleleNames[mafMsiSiteToNameMapping[val]] = observedDoubleMutsEndoColoHyper[val]\n",
    "        \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = summarize_double_mut_msi_character(msiAlleleNames, msiSummary)\n",
    "df['label'] = df.apply(lambda row: row['allele'] if row['nOccurences'] > 2 else None, axis=1)\n",
    "df['basePair'] = df['basePair'].apply(lambda x: x if x in set(['A', 'C', 'T', 'G']) else 'other')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.to_csv('/Users/friedman/Desktop/WORK/dataForLocalPlotting/doubleMutationMSIInfo.tsv', index=False, sep='\\t')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Look for MSI Signatures** <br/> <br/> <br/> <br/> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "metadata": {},
   "outputs": [],
   "source": [
    "def summarize_msi_repeat_info(impactMsiMuts, msiInfo, impactSigs, mmrGeneSummary):\n",
    "    listOfDicts = []\n",
    "    cntr = 0\n",
    "    for case in set(impactMsiMuts['Tumor_Sample_Barcode']):\n",
    "        if cntr%100 == 0: print cntr\n",
    "        cntr +=1\n",
    "        \n",
    "        caseMuts = impactMsiMuts[impactMsiMuts['Tumor_Sample_Barcode'] == case]\n",
    "        alleles = set(caseMuts[caseMuts['correctedAllele'].notnull()]['correctedAllele'])\n",
    "        msiInfoAllelesInCase = msiInfo[msiInfo['correctedAllele'].isin(alleles)] #get msi site info about the case\n",
    "        basePairCounter = Counter(msiInfoAllelesInCase['Tumor_Seq_Allele2'])\n",
    "        nCG = basePairCounter['C'] + basePairCounter['G']\n",
    "        nAT = basePairCounter['A'] + basePairCounter['T']\n",
    "        nOther = sum(basePairCounter.values()) - nCG - nAT\n",
    "        \n",
    "        #cg ratio and indel length\n",
    "        cg_atRATIO = None\n",
    "        if nAT > 0: cg_atRATIO = (1.0*nCG)/nAT\n",
    "        caseAverageIndelLength = np.nanmean(msiInfoAllelesInCase['repeat_times'])    \n",
    "        \n",
    "        #gene info\n",
    "        mmrGene = 'None'\n",
    "        pid = caseMuts['pid'].iloc[0]\n",
    "        if pid in set(mmrGeneSummary['Patient_ID']):\n",
    "            mmrGene = mmrGeneSummary[mmrGeneSummary['Patient_ID'] == pid]['Hugo_Symbol'].iloc[0]\n",
    "        \n",
    "        dominantSig = impactSigs[impactSigs['Tumor_Sample_Barcode'] == case]['dominantSiganture'].iloc[0]\n",
    "        \n",
    "        listOfDicts.append({'Tumor_Sample_Barcode': case,'dominantSig': dominantSig, 'mmrGene': mmrGene,\n",
    "                            'nCG': nCG, 'nAT': nAT,\n",
    "                            'nC': basePairCounter['C'], 'nG': basePairCounter['G'],\n",
    "                            'nA': basePairCounter['A'], 'nT': basePairCounter['T'],\n",
    "                            'cg_atRATIO': cg_atRATIO, 'averageIndelLen': caseAverageIndelLength, 'nOtherBP': nOther\n",
    "                           })\n",
    "    return pd.DataFrame(listOfDicts)\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda2/lib/python2.7/site-packages/IPython/core/interactiveshell.py:2818: DtypeWarning: Columns (5,14,18,24,25,26,28,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47) have mixed types. Specify dtype option on import or set low_memory=False.\n",
      "  if self.run_code(code, result):\n"
     ]
    }
   ],
   "source": [
    "allImpactMuts = analysis_utils.map_cases_to_msi_sensor_class(allImpactMuts, msiSensorInfo=pathPrefix + '/ifs/work/taylorlab/friedman/mskImpactAsOfMarch2019/dmp/mskimpact/data_clinical_sample.txt')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "metadata": {},
   "outputs": [],
   "source": [
    "allImpactMuts['pid'] = allImpactMuts['Tumor_Sample_Barcode'].apply(lambda x: x[:9])\n",
    "msiOnlyMuts = allImpactMuts[allImpactMuts['caseMsiClass'] == 'Instable']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "preethiMMRSummary = pd.read_table('/Users/friedman/Downloads/unique_MMR-4genes.txt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 176,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "100\n",
      "200\n",
      "300\n",
      "400\n",
      "500\n",
      "600\n",
      "700\n",
      "800\n"
     ]
    }
   ],
   "source": [
    "df = summarize_msi_repeat_info(msiOnlyMuts, msiSummary, impactSigs, preethiMMRSummary)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 177,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['dominantSig'] = df['dominantSig'].apply(lambda x: x if x in set(['mean_6', 'mean_1', 'mean_14', 'mean_15', 'mean_20', \n",
    "                                                                'mean_21', 'mean_26']) else 'other')\n",
    "df.to_csv('/Users/friedman/Desktop/WORK/dataForLocalPlotting/msiIndelFeatures.tsv', index=False, sep='\\t')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**New section find out info about msi mut co-occurence** <br/><br/><br/><br/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "msiMuts = allImpactMuts[allImpactMuts['correctedAllele'].notnull()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "tumorSuppressors = 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'])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda2/lib/python2.7/site-packages/ipykernel_launcher.py:1: 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",
      "  \"\"\"Entry point for launching an IPython kernel.\n",
      "/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": [
    "msiMuts['geneType'] = msiMuts['Hugo_Symbol'].apply(lambda x: 'Tumor Suppressor' if x in tumorSuppressors else 'Oncogene')\n",
    "msiMuts['isFrameshift'] = msiMuts['Consequence'].apply(lambda x: True if x == 'frameshift_variant' else False)\n",
    "msiSummary['geneType'] = msiSummary['Hugo_Symbol'].apply(lambda x: 'Tumor Suppressor' if x in tumorSuppressors else 'Oncogene')\n",
    "msiSummary['isFrameshift'] = msiSummary['Consequence'].apply(lambda x: True if x == 'frameshift_variant' else False)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.stats as stats\n",
    "\n",
    "def do_fisher_test_for_genes_indels(obsMuts, possibleMuts):\n",
    "    listOfDicts = []\n",
    "    for gene in set(obsMuts['Hugo_Symbol']):\n",
    "        obsGene = obsMuts[obsMuts['Hugo_Symbol'] == gene]\n",
    "        possgene = possibleMuts[possibleMuts['Hugo_Symbol'] == gene]\n",
    "        \n",
    "        nObsFrameshift = obsGene[obsGene['isFrameshift'] == True].shape[0]\n",
    "        nObsNotFrameshift = obsGene[obsGene['isFrameshift'] == False].shape[0]\n",
    "        nPossFrameshift = possibleMuts[possibleMuts['isFrameshift'] == True].shape[0]\n",
    "        nPossNotFrameshift = possibleMuts[possibleMuts['isFrameshift'] == False].shape[0]\n",
    "        \n",
    "        oddsRatio, pValue = stats.fisher_exact([[nPossFrameshift, nPossNotFrameshift],\n",
    "                                                                    [nObsFrameshift, nObsNotFrameshift]])\n",
    "        \n",
    "        listOfDicts.append({'Hugo_Symbol': gene, 'oddsRatio': oddsRatio, 'pValue': pValue})\n",
    "    return pd.DataFrame(listOfDicts)\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = do_fisher_test_for_genes_indels(msiMuts, msiSummary)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['displayLabel'] = df.apply(lambda row: row['Hugo_Symbol'] if row['pValue'] < 1e-8 else None, axis=1)\n",
    "df['isTumorS'] = df['Hugo_Symbol'].apply(lambda x: True if x in tumorSuppressors else False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.to_csv('/Users/friedman/Desktop/WORK/dataForLocalPlotting/frameShiftTestMsi.tsv', index=False, sep='\\t')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "allImpactMuts['quadNuc'] = allImpactMuts.apply(lambda row: mutationSigUtils.create_reference_four_nuc(row['Ref_Tri'], row['Reference_Allele'], row['Tumor_Seq_Allele2'], row['Variant_Type']), axis=1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "impactSigs = pd.read_table(pathPrefix + '/ifs/work/taylorlab/friedman/hypermutationAnalysisProj/projectDataAndConfigFiles/signatures_from_unfiltered_maf.txt')\n",
    "\n",
    "#print allImpactMuts[allImpactMuts['HGVSp_Short'] == 'p.Y27C']['quadNuc']"
   ]
  }
 ],
 "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
}