summarize_nmut_and_noncogenic_info_for_mutation_simulation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import argparse\n",
    "import os\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from tqdm import tqdm\n",
    "tqdm.pandas()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "pathPrefix = '/Users/friedman/Desktop/mnt'\n",
    "sys.path.append(pathPrefix + '/ifs/work/taylorlab/friedman/myUtils')\n",
    "import analysis_utils\n",
    "import mutationSigUtils\n",
    "import maf_analysis_utils"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#reload modules as needed in interactive dev\n",
    "reload(analysis_utils)\n",
    "reload(maf_analysis_utils)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**QC SCRIPTS**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def summarize_qc_metrics_for_simulated_mafs(simMafSummaryPath = pathPrefix + '/ifs/work/taylorlab/friedman/hypermutationAnalysisProj/projectDataAndConfigFiles/simulatedMutationSummary.tsv',\n",
    "                                            observedMafSummaryPath = pathPrefix + '/ifs/work/taylorlab/friedman/hypermutationAnalysisProj/projectDataAndConfigFiles/observedMutationCounts.tsv'):\n",
    "    print 'loading in simulated maf summary'\n",
    "    dfSimulation = pd.read_table(simMafSummaryPath)\n",
    "    print 'loading in observed summary'\n",
    "    dfObserved = pd.read_table(observedMafSummaryPath)\n",
    "    print 'doing mapping'\n",
    "    chromosomeDict = analysis_utils.map_impact_genes_to_chromosome(impactRegionsFile = pathPrefix + '/ifs/work/taylorlab/friedman/myAdjustedDataFiles/genelist.with_aa.interval_list')\n",
    "    \n",
    "    simulatedOncogenicDict = dict(zip(dfSimulation['Hugo_Symbol'], dfSimulation['totalNOncogenic']))\n",
    "    observedOncogenicDict = dict(zip(dfObserved['gene'], dfObserved['nDistinctOncogenicMutations']))\n",
    "    \n",
    "    for gene in simulatedOncogenicDict.keys():\n",
    "        nSimOnc = simulatedOncogenicDict[gene]\n",
    "        if gene in observedOncogenicDict:\n",
    "            nDistinctObservedOnc = observedOncogenicDict[gene]\n",
    "            if nDistinctObservedOnc > nSimOnc and nSimOnc == 0:\n",
    "                if gene in chromosomeDict:\n",
    "                    print 'problem with: ', gene, chromosomeDict[gene], '****** ',  nDistinctObservedOnc, '   vs', nSimOnc\n",
    "                else:\n",
    "                    print 'problem with: ', gene, ' unknown chr'\n",
    "        else:\n",
    "            print 'never observed oncogenic mut in ', gene\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "loading in simulated maf summary\n",
      "loading in observed summary\n",
      "doing mapping\n",
      "problem with:  GATA1 X ******  6    vs 0\n",
      "problem with:  INPP4A 2 ******  2    vs 0\n",
      "problem with:  CUL3 2 ******  4    vs 0\n",
      "problem with:  MDC1 6 ******  2    vs 0\n",
      "problem with:  CENPA 2 ******  1    vs 0\n",
      "problem with:  FLT1 13 ******  3    vs 0\n",
      "problem with:  FANCA 16 ******  79    vs 0\n",
      "problem with:  CDK8 13 ******  7    vs 0\n",
      "problem with:  HIST1H3G 6 ******  1    vs 0\n",
      "problem with:  HIST1H3E 6 ******  3    vs 0\n",
      "problem with:  PARK2 6 ******  64    vs 0\n",
      "problem with:  INHBA 7 ******  6    vs 0\n",
      "problem with:  SYK 9 ******  4    vs 0\n",
      "problem with:  NTRK2 9 ******  3    vs 0\n",
      "problem with:  HIST1H2BD 6 ******  10    vs 0\n",
      "problem with:  AGO2  unknown chr\n",
      "problem with:  TRAF2 9 ******  5    vs 0\n",
      "problem with:  TMPRSS2 21 ******  2    vs 0\n"
     ]
    }
   ],
   "source": [
    "summarize_qc_metrics_for_simulated_mafs(simMafSummaryPath = pathPrefix + '/ifs/work/taylorlab/friedman/hypermutationAnalysisProj/projectDataAndConfigFiles/simulatedMutationSummary.tsv')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "chromosomeDict = analysis_utils.map_impact_genes_to_chromosome(impactRegionsFile = pathPrefix + '/ifs/work/taylorlab/friedman/myAdjustedDataFiles/genelist.with_aa.interval_list')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "for key, value in chromosomeDict.items():\n",
    "    if value == '2': print key"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Count totals for nmut, and nmut trinuc across the big possible mutations maf**\n",
    "<br/>\n",
    "Loading in the maf takes forever--as long as 45 minutes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#wc -l /ifs/work/taylorlab/friedman/myAdjustedDataFiles/simulatedMafs/all_possible_impact_muts_limited_columns.maf\n",
    "nPossibleImpactMuts = 5008844"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#enumerate n muts per each quadnuc\n",
    "#this is currently about 1 hour to load and process using a mounted filesystem\n",
    "\n",
    "#df = pd.progress_read_table(pathPrefix + '/ifs/work/taylorlab/friedman/myAdjustedDataFiles/simulatedMafs/all_possible_impact_muts_limited_columns.maf')\n",
    "#filepath = pathPrefix + '/ifs/work/taylorlab/friedman/myAdjustedDataFiles/simulatedMafs/all_possible_impact_muts_limited_columns.maf'\n",
    "filepath = pathPrefix + '/ifs/work/taylorlab/friedman/myAdjustedDataFiles/simulatedMafs/allmut_impact_simulated_muts_mafAnno_trinuc_v2.maf'\n",
    "#get a realistic time estimate for read table\n",
    "\n",
    "nPossibleImpactMuts = 5008844\n",
    "allPossibleMutMaf = analysis_utils.load_in_df_with_progress(filepath, nPossibleImpactMuts)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "allPossibleMutMaf['quadNuc'] = allPossibleMutMaf.progress_apply(lambda row: mutationSigUtils.create_reference_four_nuc(row['Ref_Tri'], row['Reference_Allele'], row['Tumor_Seq_Allele2'], row['Variant_Type']), axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Iterates over genes and writes a per gene mutation info file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "impactGenes = set(['ABL1', 'ACVR1', 'AGO2', 'AKT1', 'AKT2', 'AKT3', 'ALK', 'ALOX12B', 'ANKRD11', 'APC', 'AR', 'ARAF', 'ARID1A', 'ARID1B', 'ARID2', 'ARID5B', 'ASXL1', 'ASXL2', 'ATM', 'ATR', 'ATRX', 'AURKA', 'AURKB', 'AXIN1', 'AXIN2', 'AXL', 'B2M', 'BABAM1', 'BAP1', 'BARD1', 'BBC3', 'BCL10', 'BCL2', 'BCL2L1', 'BCL2L11', 'BCL6', 'BCOR', 'BIRC3', 'BLM', 'BMPR1A', 'BRAF', 'BRCA1', 'BRCA2', 'BRD4', 'BRIP1', 'BTK', 'CALR', 'CARD11', 'CARM1', 'CASP8', 'CBFB', 'CBL', 'CCND1', 'CCND2', 'CCND3', 'CCNE1', 'CD274', 'CD276', 'CD79A', 'CD79B', 'CDC42', 'CDC73', 'CDH1', 'CDK12', 'CDK4', 'CDK6', 'CDK8', 'CDKN1A', 'CDKN1B', 'CDKN2A', 'CDKN2B', 'CDKN2C', 'CEBPA', 'CENPA', 'CHEK1', 'CHEK2', 'CIC', 'CREBBP', 'CRKL', 'CRLF2', 'CSDE1', 'CSF1R', 'CSF3R', 'CTCF', 'CTLA4', 'CTNNB1', 'CUL3', 'CXCR4', 'CYLD', 'CYSLTR2', 'DAXX', 'DCUN1D1', 'DDR2', 'DICER1', 'DIS3', 'DNAJB1', 'DNMT1', 'DNMT3A', 'DNMT3B', 'DOT1L', 'DROSHA', 'DUSP4', 'E2F3', 'EED', 'EGFL7', 'EGFR', 'EIF1AX', 'EIF4A2', 'EIF4E', 'ELF3', 'EP300', 'EPAS1', 'EPCAM', 'EPHA3', 'EPHA5', 'EPHA7', 'EPHB1', 'ERBB2', 'ERBB3', 'ERBB4', 'ERCC2', 'ERCC3', 'ERCC4', 'ERCC5', 'ERF', 'ERG', 'ERRFI1', 'ESR1', 'ETV1', 'ETV6', 'EZH1', 'EZH2', 'FAM123B', 'FAM175A', 'FAM46C', 'FAM58A', 'FANCA', 'FANCC', 'FAT1', 'FBXW7', 'FGF19', 'FGF3', 'FGF4', 'FGFR1', 'FGFR2', 'FGFR3', 'FGFR4', 'FH', 'FLCN', 'FLT1', 'FLT3', 'FLT4', 'FOXA1', 'FOXL2', 'FOXO1', 'FOXP1', 'FUBP1', 'FYN', 'GATA1', 'GATA2', 'GATA3', 'GLI1', 'GNA11', 'GNAQ', 'GNAS', 'GPS2', 'GREM1', 'GRIN2A', 'GSK3B', 'H3F3A', 'H3F3B', 'H3F3C', 'HGF', 'HIST1H1C', 'HIST1H2BD', 'HIST1H3A', 'HIST1H3B', 'HIST1H3C', 'HIST1H3D', 'HIST1H3E', 'HIST1H3F', 'HIST1H3G', 'HIST1H3H', 'HIST1H3I', 'HIST1H3J', 'HIST2H3C', 'HIST2H3D', 'HIST3H3', 'HLA-A', 'HLA-B', 'HNF1A', 'HOXB13', 'HRAS', 'ICOSLG', 'ID3', 'IDH1', 'IDH2', 'IFNGR1', 'IGF1', 'IGF1R', 'IGF2', 'IKBKE', 'IKZF1', 'IL10', 'IL7R', 'INHA', 'INHBA', 'INPP4A', 'INPP4B', 'INPPL1', 'INSR', 'IRF4', 'IRS1', 'IRS2', 'JAK1', 'JAK2', 'JAK3', 'JUN', 'KDM5A', 'KDM5C', 'KDM6A', 'KDR', 'KEAP1', 'KIT', 'KLF4', 'KMT2B', 'KMT5A', 'KNSTRN', 'KRAS', 'LATS1', 'LATS2', 'LMO1', 'LYN', 'MALT1', 'MAP2K1', 'MAP2K2', 'MAP2K4', 'MAP3K1', 'MAP3K13', 'MAP3K14', 'MAPK1', 'MAPK3', 'MAPKAP1', 'MAX', 'MCL1', 'MDC1', 'MDM2', 'MDM4', 'MED12', 'MEF2B', 'MEN1', 'MET', 'MGA', 'MITF', 'MLH1', 'KMT2A', 'KMT2B', 'KMT2C', 'MPL', 'MRE11A', 'MSH2', 'MSH3', 'MSH6', 'MSI1', 'MSI2', 'MST1', 'MST1R', 'MTOR', 'MUTYH', 'MYC', 'MYCL1', 'MYCN', 'MYD88', 'MYOD1', 'NBN', 'NCOA3', 'NCOR1', 'NEGR1', 'NF1', 'NF2', 'NFE2L2', 'NFKBIA', 'NKX2-1', 'NKX3-1', 'NOTCH1', 'NOTCH2', 'NOTCH3', 'NOTCH4', 'NPM1', 'NRAS', 'NSD1', 'NTHL1', 'NTRK1', 'NTRK2', 'NTRK3', 'NUF2', 'NUP93', 'PAK1', 'PAK7', 'PALB2', 'PARK2', 'PARP1', 'PAX5', 'PBRM1', 'PDCD1', 'PDCD1LG2', 'PDGFRA', 'PDGFRB', 'PDPK1', 'PGR', 'PHOX2B', 'PIK3C2G', 'PIK3C3', 'PIK3CA', 'PIK3CB', 'PIK3CD', 'PIK3CG', 'PIK3R1', 'PIK3R2', 'PIK3R3', 'PIM1', 'PLCG2', 'PLK2', 'PMAIP1', 'PMS1', 'PMS2', 'PNRC1', 'POLD1', 'POLE', 'PPARG', 'PPM1D', 'PPP2R1A', 'PPP4R2', 'PPP6C', 'PRDM1', 'PRDM14', 'PREX2', 'PRKAR1A', 'PRKCI', 'PRKD1', 'PTCH1', 'PTEN', 'PTP4A1', 'PTPN11', 'PTPRD', 'PTPRS', 'PTPRT', 'RAB35', 'RAC1', 'RAC2', 'RAD21', 'RAD50', 'RAD51', 'RAD51C', 'RAD51L1', 'RAD51L3', 'RAD52', 'RAD54L', 'RAF1', 'RARA', 'RASA1', 'RB1', 'RBM10', 'RECQL', 'RECQL4', 'REL', 'RET', 'RFWD2', 'RHEB', 'RHOA', 'RICTOR', 'RIT1', 'RNF43', 'ROS1', 'RPS6KA4', 'RPS6KB2', 'RPTOR', 'RRAGC', 'RRAS', 'RRAS2', 'RTEL1', 'RUNX1', 'RXRA', 'RYBP', 'SDHA', 'SDHAF2', 'SDHB', 'SDHC', 'SDHD', 'SESN1', 'SESN2', 'SESN3', 'SETD2', 'SF3B1', 'SH2B3', 'SH2D1A', 'SHOC2', 'SHQ1', 'SLX4', 'SMAD2', 'SMAD3', 'SMAD4', 'SMARCA4', 'SMARCB1', 'SMARCD1', 'SMO', 'SMYD3', 'SOCS1', 'SOS1', 'SOX17', 'SOX2', 'SOX9', 'SPEN', 'SPOP', 'SPRED1', 'SRC', 'SRSF2', 'STAG2', 'STAT3', 'STAT5A', 'STAT5B', 'STK11', 'STK19', 'STK40', 'SUFU', 'SUZ12', 'SYK', 'TAP1', 'TAP2', 'TBX3', 'TCEB1', 'TCF3', 'TCF7L2', 'TEK', 'TERT', 'TET1', 'TET2', 'TGFBR1', 'TGFBR2', 'TMEM127', 'TMPRSS2', 'TNFAIP3', 'TNFRSF14', 'TOP1', 'TP53', 'TP53BP1', 'TP63', 'TRAF2', 'TRAF7', 'TSC1', 'TSC2', 'TSHR', 'U2AF1', 'UPF1', 'VEGFA', 'VHL', 'VTCN1', 'WHSC1', 'WHSC1L1', 'WT1', 'WWTR1', 'XIAP', 'XPO1', 'XRCC2', 'YAP1', 'YES1', 'ZFHX3', 'ZRSR2'])\n",
    "writeDir = pathPrefix + '/ifs/work/taylorlab/friedman/myAdjustedDataFiles/simulatedMafs/geneMutMafs'\n",
    "cntr = 0\n",
    "for gene in impactGenes:\n",
    "    cntr += 1\n",
    "    geneMuts = allPossibleMutMaf[allPossibleMutMaf['Hugo_Symbol'] == gene]\n",
    "    \n",
    "    print gene, cntr, 'nmut:', geneMuts.shape[0]\n",
    "    \n",
    "    writePath = os.path.join(writeDir, gene + '_all_possible_snps_v2.maf')\n",
    "    geneMuts.to_csv(writePath, index=False, sep='\\t')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Now use functions to look at the files we wrote and assess number of oncogenic muts etc**<br/>\n",
    "note--currently we are using the files with the suffix v2 (be aware of this)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "cntr = 0\n",
    "nOncogenicMutDict = {}\n",
    "for gene in impactGenes:\n",
    "    if cntr %10 == 0: print cntr\n",
    "    cntr += 1\n",
    "    fileName = gene + '_all_possible_snps_v2.maf'\n",
    "    nmut = maf_analysis_utils.count_n_oncogenic_muts_in_maf(\n",
    "        pathPrefix + '/ifs/work/taylorlab/friedman/myAdjustedDataFiles/simulatedMafs/geneMutMafs/' + fileName)\n",
    "    nOncogenicMutDict[gene] = nmut\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "scrolled": true
   },
   "source": [
    "**PART 3--use simulated data to write out a summary file**<br/>\n",
    "Do python calculate_summary_information_from_simulated_mafs.py (hypermutation project directory)\n",
    "Its done this way because it is by far the fastest way to do it (too slow to deal with IO times on this mounted stuff) Takes about 1hr to runa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br> <br/>"
   ]
  }
 ],
 "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
}