analyze_gtex_data.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 re\n",
    "import scipy.stats as stats\n",
    "\n",
    "from collections import Counter\n",
    "\n",
    "pathPrefix = '/Users/friedman/Desktop/mnt'\n",
    "\n",
    "sys.path.append(pathPrefix + '/juno/work/taylorlab/friedman/myUtils')\n",
    "import analysis_utils \n",
    "import mutationSigUtils\n",
    "import get_gene_and_cohort_list_utils"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 150,
   "metadata": {},
   "outputs": [],
   "source": [
    "def summarize_rates(maf, sizeInfo, expressedGenes, notExpressedGenes, eosGenes):\n",
    "    \n",
    "    listOfDicts = []\n",
    "    maf = maf[~maf['Hugo_Symbol'].isin(eosGenes)]\n",
    "    \n",
    "    expressedGeneMuts = maf[maf['Hugo_Symbol'].isin(expressedGenes)]\n",
    "    notExpressedGeneMuts = maf[maf['Hugo_Symbol'].isin(notExpressedGenes)]\n",
    "    \n",
    "    nExpressedGeneBases = sum([sizeInfo[gene] for gene in expressedGenes if gene in sizeInfo and gene not in eosGenes])\n",
    "    nNotExpressedGeneBases = sum([sizeInfo[gene] for gene in notExpressedGenes if gene in sizeInfo and gene not in eosGenes])\n",
    "    \n",
    "    expressedGeneCounts = expressedGeneMuts['SAMPLE_ID'].value_counts()\n",
    "    notExpressedGeneCounts = notExpressedGeneMuts['SAMPLE_ID'].value_counts()\n",
    "    \n",
    "    for case in set(maf['SAMPLE_ID']):\n",
    "        expressedRate = 0\n",
    "        if case in expressedGeneCounts:\n",
    "            expressedRate = 1.0*expressedGeneCounts[case]/nExpressedGeneBases\n",
    "        notExpressedRate = 0\n",
    "        if case in notExpressedGeneCounts:\n",
    "            notExpressedRate = 1.0*notExpressedGeneCounts[case]/nNotExpressedGeneBases\n",
    "        listOfDicts.append({'case': case, 'type': 'expressed', 'rate': expressedRate})\n",
    "        listOfDicts.append({'case': case, 'type': 'notExpressed', 'rate': notExpressedRate})\n",
    "        \n",
    "    df = pd.DataFrame(listOfDicts)\n",
    "    return df\n",
    "   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#Do the same analysis for TCGA\n",
    "mc3maf = analysis_utils.load_in_df_with_progress(pathPrefix + '/ifs/res/taylorlab/ang46/ext/mafs/mc3/mc3.v0.2.8.PUBLIC.LAML_PATCH_prepped_facets_oncokb.maf', nLinesFile= 2699035)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "gtexDf = pd.read_table(pathPrefix + '/juno/work/taylorlab/friedman/myAdjustedDataFiles/GTEx_Analysis_2017-06-05_v8_RNASeQCv1.1.9_gene_median_tpm.gct', skiprows=[0,1])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 151,
   "metadata": {},
   "outputs": [],
   "source": [
    "impact341Genes = set(['ABL1', 'AKT1', 'AKT2', 'AKT3', 'ALK', 'ALOX12B', 'APC', 'AR', 'ARAF', 'ARID1A', 'ARID1B', 'ARID2', 'ARID5B', 'ASXL1', 'ASXL2', 'ATM', 'ATR', 'ATRX', 'AURKA', 'AURKB', 'AXIN1', 'AXIN2', 'AXL', 'B2M', 'BAP1', 'BARD1', 'BBC3', 'BCL2', 'BCL2L1', 'BCL2L11', 'BCL6', 'BCOR', 'BLM', 'BMPR1A', 'BRAF', 'BRCA1', 'BRCA2', 'BRD4', 'BRIP1', 'BTK', 'CARD11', 'CASP8', 'CBFB', 'CBL', 'CCND1', 'CCND2', 'CCND3', 'CCNE1', 'CD274', 'CD276', 'CD79B', 'CDC73', 'CDH1', 'CDK12', 'CDK4', 'CDK6', 'CDK8', 'CDKN1A', 'CDKN1B', 'CDKN2A', 'CDKN2B', 'CDKN2C', 'CHEK1', 'CHEK2', 'CIC', 'CREBBP', 'CRKL', 'CRLF2', 'CSF1R', 'CTCF', 'CTLA4', 'CTNNB1', 'CUL3', 'DAXX', 'DCUN1D1', 'DDR2', 'DICER1', 'DIS3', 'DNMT1', 'DNMT3A', 'DNMT3B', 'DOT1L', 'E2F3', 'EED', 'EGFL7', 'EGFR', 'EIF1AX', 'EP300', 'EPCAM', 'EPHA3', 'EPHA5', 'EPHB1', 'ERBB2', 'ERBB3', 'ERBB4', 'ERCC2', 'ERCC3', 'ERCC4', 'ERCC5', 'ERG', 'ESR1', 'ETV1', 'ETV6', 'EZH2', 'FAM123B', 'FAM175A', 'FAM46C', 'FANCA', 'FANCC', 'FAT1', 'FBXW7', 'FGF19', 'FGF3', 'FGF4', 'FGFR1', 'FGFR2', 'FGFR3', 'FGFR4', 'FH', 'FLCN', 'FLT1', 'FLT3', 'FLT4', 'FOXA1', 'FOXL2', 'FOXP1', 'FUBP1', 'GATA1', 'GATA2', 'GATA3', 'GNA11', 'GNAQ', 'GNAS', 'GREM1', 'GRIN2A', 'GSK3B', 'H3F3C', 'HGF', 'HIST1H1C', 'HIST1H2BD', 'HIST1H3B', 'HNF1A', 'HRAS', 'ICOSLG', 'IDH1', 'IDH2', 'IFNGR1', 'IGF1', 'IGF1R', 'IGF2', 'IKBKE', 'IKZF1', 'IL10', 'IL7R', 'INPP4A', 'INPP4B', 'INSR', 'IRF4', 'IRS1', 'IRS2', 'JAK1', 'JAK2', 'JAK3', 'JUN', 'KDM5A', 'KDM5C', 'KDM6A', 'KDR', 'KEAP1', 'KIT', 'KLF4', 'KRAS', 'LATS1', 'LATS2', 'LMO1', 'MAP2K1', 'MAP2K2', 'MAP2K4', 'MAP3K1', 'MAP3K13', 'MAPK1', 'MAX', 'MCL1', 'MDC1', 'MDM2', 'MDM4', 'MED12', 'MEF2B', 'MEN1', 'MET', 'MITF', 'MLH1', 'MLL', 'MLL2', 'MLL3', 'MPL', 'MRE11A', 'MSH2', 'MSH6', 'MTOR', 'MUTYH', 'MYC', 'MYCL1', 'MYCN', 'MYD88', 'MYOD1', 'NBN', 'NCOR1', 'NF1', 'NF2', 'NFE2L2', 'NKX2-1', 'NKX3-1', 'NOTCH1', 'NOTCH2', 'NOTCH3', 'NOTCH4', 'NPM1', 'NRAS', 'NSD1', 'NTRK1', 'NTRK2', 'NTRK3', 'PAK1', 'PAK7', 'PALB2', 'PARK2', 'PARP1', 'PAX5', 'PBRM1', 'PDCD1', 'PDGFRA', 'PDGFRB', 'PDPK1', 'PHOX2B', 'PIK3C2G', 'PIK3C3', 'PIK3CA', 'PIK3CB', 'PIK3CD', 'PIK3CG', 'PIK3R1', 'PIK3R2', 'PIK3R3', 'PIM1', 'PLK2', 'PMAIP1', 'PMS1', 'PMS2', 'PNRC1', 'POLE', 'PPP2R1A', 'PRDM1', 'PRKAR1A', 'PTCH1', 'PTEN', 'PTPN11', 'PTPRD', 'PTPRS', 'PTPRT', 'RAC1', 'RAD50', 'RAD51', 'RAD51B', 'RAD51C', 'RAD51D', 'RAD52', 'RAD54L', 'RAF1', 'RARA', 'RASA1', 'RB1', 'RBM10', 'RECQL4', 'REL', 'RET', 'RFWD2', 'RHOA', 'RICTOR', 'RIT1', 'RNF43', 'ROS1', 'RPS6KA4', 'RPS6KB2', 'RPTOR', 'RUNX1', 'RYBP', 'SDHA', 'SDHAF2', 'SDHB', 'SDHC', 'SDHD', 'SETD2', 'SF3B1', 'SH2D1A', 'SHQ1', 'SMAD2', 'SMAD3', 'SMAD4', 'SMARCA4', 'SMARCB1', 'SMARCD1', 'SMO', 'SOCS1', 'SOX17', 'SOX2', 'SOX9', 'SPEN', 'SPOP', 'SRC', 'STAG2', 'STK11', 'STK40', 'SUFU', 'SUZ12', 'SYK', 'TBX3', 'TERT', 'TET1', 'TET2', 'TGFBR1', 'TGFBR2', 'TMEM127', 'TMPRSS2', 'TNFAIP3', 'TNFRSF14', 'TOP1', 'TP53', 'TP63', 'TRAF7', 'TSC1', 'TSC2', 'TSHR', 'U2AF1', 'VHL', 'VTCN1', 'WT1', 'XIAP', 'XPO1', 'YAP1', 'YES1'])\n",
    "impactGenes = ['ABL1', 'ACVR1', 'AGO2', 'AKT1', 'AKT2', 'AKT3', 'ALK', 'ALOX12B', 'AMER1', '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', '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', 'GTF2I', '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', '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', 'MLL', 'MLL2', 'MLL3', 'MLL4', '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', 'RAD51B', 'RAD51C', 'RAD51D', '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', 'SETD8', '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",
    "\n",
    "impactGtex = gtexDf[gtexDf['Description'].isin(impactGenes)]\n",
    "#impactGtex = gtexDf[gtexDf['Description'].isin(impact341Genes)]\n",
    "impactGeneLengths = analysis_utils.get_cds_size_targeted_by_impact(infoFilePath = pathPrefix + '/ifs/work/taylorlab/friedman/myAdjustedDataFiles/impact_gene_reference_signatures.tsv')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "metadata": {},
   "outputs": [],
   "source": [
    "#identify the genes that are not expressed in the cancer types I care about\n",
    "#todo identify them in a more sophisticated manner\n",
    "expressionThresh = 1\n",
    "colonNotExpressedGenes = set(impactGtex[impactGtex['Colon - Transverse'] < expressionThresh]['Description']) & set(impactGtex[impactGtex['Colon - Sigmoid'] < expressionThresh]['Description'])\n",
    "endometrialNotExpressedGenes = set(impactGtex[impactGtex['Uterus'] < expressionThresh]['Description'])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 187,
   "metadata": {},
   "outputs": [],
   "source": [
    "#GET TCGA IDs by signatures\n",
    "mmrIds, poleIds = get_gene_and_cohort_list_utils.get_tcga_pole_mmr_hypermutator_ids(tcgaSigsPath = pathPrefix + '/juno/work/taylorlab/friedman/myAdjustedDataFiles/tcgaSigsCombined.txt')\n",
    "uterineIds = get_gene_and_cohort_list_utils.get_tcga_cancer_type_info(tcgaInfoPath = pathPrefix + '/juno/work/taylorlab/friedman/myAdjustedDataFiles/tcgaCancerTypeInfo.txt', cancerTypes = ['TCGA-UCEC'])\n",
    "colorectalIds = get_gene_and_cohort_list_utils.get_tcga_cancer_type_info(tcgaInfoPath = pathPrefix + '/juno/work/taylorlab/friedman/myAdjustedDataFiles/tcgaCancerTypeInfo.txt', cancerTypes = ['TCGA-COAD', 'TCGA-READ'])\n",
    "\n",
    "#strip the numbers from the signature ids\n",
    "mmrIds = set([i[:-3] for i in mmrIds])\n",
    "poleIds = set([i[:-3] for i in poleIds])\n",
    "\n",
    "colorectalMMRIds = mmrIds & colorectalIds\n",
    "endometrialMMRIds = mmrIds & uterineIds\n",
    "colorectalPOLEIds = poleIds & colorectalIds\n",
    "endometrialPOLEIds = poleIds & uterineIds\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mc3maf['SAMPLE_ID_short'] = mc3maf['SAMPLE_ID'].apply(lambda x: x[:-3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 216,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "59 139\n",
      "50 97\n"
     ]
    }
   ],
   "source": [
    "#QC by removing cases with MMR signature but not enough indels\n",
    "mmrFrameshiftIndels = mc3maf[(mc3maf['Consequence'] == 'frameshift_variant') \n",
    "    & (mc3maf['SAMPLE_ID_short'].isin(colorectalMMRIds | endometrialMMRIds))]\n",
    "\n",
    "minNTCGAIndels = 75\n",
    "enoughMuts = mmrFrameshiftIndels['SAMPLE_ID_short'].value_counts() > minNTCGAIndels\n",
    "casesWithTooFewIndels = set(enoughMuts[enoughMuts == False].index)\n",
    "\n",
    "colorectalMMRIds = colorectalMMRIds - casesWithTooFewIndels\n",
    "endometrialMMRIds = endometrialMMRIds - casesWithTooFewIndels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 217,
   "metadata": {},
   "outputs": [],
   "source": [
    "#get the mutations we care about from the mc3 maf\n",
    "mmrFrameshiftIndelsColorectal = mc3maf[(mc3maf['Consequence'] == 'frameshift_variant') \n",
    "    & (mc3maf['SAMPLE_ID_short'].isin(colorectalMMRIds))]\n",
    "\n",
    "mmrFrameshiftIndelsEndometrial = mc3maf[(mc3maf['Consequence'] == 'frameshift_variant') \n",
    "    & (mc3maf['SAMPLE_ID_short'].isin(endometrialMMRIds))]\n",
    "\n",
    "poleStopGainMutsEndometrial = mc3maf[(mc3maf['Consequence'] == 'stop_gained')\n",
    "    & (mc3maf['SAMPLE_ID_short'].isin(endometrialPOLEIds))]\n",
    "\n",
    "poleStopGainMutsColorectal = mc3maf[(mc3maf['Consequence'] == 'stop_gained')\n",
    "    & (mc3maf['SAMPLE_ID_short'].isin(colorectalPOLEIds))]\n",
    "                        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [],
   "source": [
    "eosDict = {'Endometrial Cancer': ['CASP8', 'MSH6', 'JAK1', 'INPPL1', 'MSH2', 'MAP3K1', 'ARID5B', 'RNF43', 'APC', 'RB1',\n",
    "                      'CCND1', 'NFE2L2', 'RASA1', 'ARID1A', 'PTEN', 'PIK3R1', 'ATM', 'EPHA3', 'B2M', 'TP53',\n",
    "                      'BCOR', 'CTNNB1', 'KRAS', 'FGFR2', 'PPP2R1A', 'POLE', 'NF1', 'FBXW7', 'XPO1', 'ERBB2', 'ERBB3',\n",
    "                      'PTPRD', 'CTCF', 'PIK3CA'],\n",
    "\n",
    "            'Colorectal Cancer': ['CASP8', 'MSH2', 'HLA-A', 'HLA-B', 'PIK3R1', 'B2M', 'KRAS', 'PTEN', 'APC',\n",
    "                     'PIK3CA', 'FBXW7', 'RNF43', 'AMER1', 'TCF7L2', 'SMAD4', 'NRAS', 'BRAF', 'SOX9',\n",
    "                      'TP53', 'CDKN2A', 'SMAD2', 'ARID1A', 'POLE', 'MSH3', 'MLH1', 'MSH6'],\n",
    "              \n",
    "              'Glioma': ['IDH1', 'TP53', 'TERT', 'MSH6', 'CIC', 'PIK3CA', 'MSH6',\n",
    "                 'MSH2', 'NF1', 'TSC2', 'PTEN', 'ATRX', 'EGFR', 'PIK3R1', 'RB1']\n",
    "              }\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "essentialGenes = get_gene_and_cohort_list_utils.get_essential_genes(depMapPath = pathPrefix + '/juno/work/taylorlab/friedman/myAdjustedDataFiles/achillesDepMap.csv')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geneLengthInfo = pd.read_table(pathPrefix + '/juno/work/taylorlab/friedman/myAdjustedDataFiles/coding_sequences_all_exome_genes_generated_by_alexG.tsv')\n",
    "geneLengthDict = dict(zip(geneLengthInfo['hgnc_symbol'], geneLengthInfo['nt.length']))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 204,
   "metadata": {},
   "outputs": [],
   "source": [
    "#subset the gtex information by gene type\n",
    "gtexImpact = gtexDf[(gtexDf['Description'].isin(impactGenes))]\n",
    "\n",
    "gtexNotEssentialNotImpact = gtexDf[(~gtexDf['Description'].isin(essentialGenes)) &\n",
    "                                  (~gtexDf['Description'].isin(impactGenes))]\n",
    "\n",
    "allNotEssentialNotImpactGenes = set(gtexNotEssentialNotImpact['Description'])\n",
    "expressionThresh = 1\n",
    "\n",
    "colonNotExpressedGenesGtexIMPACT = set(gtexImpact[(gtexImpact['Colon - Transverse'] < expressionThresh) & (gtexImpact['Colon - Sigmoid'] < expressionThresh)]['Description'])\n",
    "endometrialNotExpressedGenesGtexIMPACT = set(gtexImpact[gtexImpact['Uterus'] < expressionThresh]['Description'])\n",
    "colonNotExpressedGenesGtexNOTIMPACT = set(gtexNotEssentialNotImpact[(gtexNotEssentialNotImpact['Colon - Transverse'] < expressionThresh) & (gtexNotEssentialNotImpact['Colon - Sigmoid'] < expressionThresh)]['Description'])\n",
    "endometrialNotExpressedGenesGtexNOTIMPACT = set(gtexNotEssentialNotImpact[gtexNotEssentialNotImpact['Uterus'] < expressionThresh]['Description'])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 218,
   "metadata": {},
   "outputs": [],
   "source": [
    "#BIG CALL: get the summary for all the cancer types and genes\n",
    "colorectalMMRSummaryIMGenes = summarize_rates(mmrFrameshiftIndelsColorectal[mmrFrameshiftIndelsColorectal['Hugo_Symbol'].isin(impactGenes)],\n",
    "    impactGeneLengths, set(impactGenes) - colonNotExpressedGenesGtexIMPACT, colonNotExpressedGenesGtexIMPACT, eosDict['Colorectal Cancer'])\n",
    "\n",
    "colorectalPOLESummaryIMGenes = summarize_rates(poleStopGainMutsColorectal[poleStopGainMutsColorectal['Hugo_Symbol'].isin(impactGenes)],\n",
    "    impactGeneLengths, set(impactGenes) - colonNotExpressedGenesGtexIMPACT, colonNotExpressedGenesGtexIMPACT, eosDict['Colorectal Cancer'])\n",
    "\n",
    "endometrialMMRSummaryIMGenes = summarize_rates(mmrFrameshiftIndelsEndometrial[mmrFrameshiftIndelsEndometrial['Hugo_Symbol'].isin(impactGenes)],\n",
    "    impactGeneLengths, set(impactGenes) - endometrialNotExpressedGenesGtexIMPACT, endometrialNotExpressedGenesGtexIMPACT, eosDict['Endometrial Cancer'])\n",
    "\n",
    "endometrialPOLESummaryIMGenes = summarize_rates(poleStopGainMutsEndometrial[poleStopGainMutsEndometrial['Hugo_Symbol'].isin(impactGenes)],\n",
    "    impactGeneLengths, set(impactGenes) - endometrialNotExpressedGenesGtexIMPACT, endometrialNotExpressedGenesGtexIMPACT, eosDict['Endometrial Cancer'])\n",
    "\n",
    "############################\n",
    "\n",
    "colorectalMMRSummaryNOTIMGenes = summarize_rates(mmrFrameshiftIndelsColorectal[~mmrFrameshiftIndelsColorectal['Hugo_Symbol'].isin(impactGenes)],\n",
    "    geneLengthDict, allNotEssentialNotImpactGenes - colonNotExpressedGenesGtexNOTIMPACT, colonNotExpressedGenesGtexNOTIMPACT, [])\n",
    "\n",
    "colorectalPOLESummaryNOTIMGenes = summarize_rates(poleStopGainMutsColorectal[~poleStopGainMutsColorectal['Hugo_Symbol'].isin(impactGenes)],\n",
    "    geneLengthDict, allNotEssentialNotImpactGenes - colonNotExpressedGenesGtexNOTIMPACT, colonNotExpressedGenesGtexNOTIMPACT, [])\n",
    "\n",
    "endometrialMMRSummaryNOTIMGenes = summarize_rates(mmrFrameshiftIndelsEndometrial[~mmrFrameshiftIndelsEndometrial['Hugo_Symbol'].isin(impactGenes)],\n",
    "    geneLengthDict, allNotEssentialNotImpactGenes - endometrialNotExpressedGenesGtexNOTIMPACT, endometrialNotExpressedGenesGtexNOTIMPACT, [])\n",
    "\n",
    "endometrialPOLESummaryNOTIMGenes = summarize_rates(poleStopGainMutsEndometrial[~poleStopGainMutsEndometrial['Hugo_Symbol'].isin(impactGenes)],\n",
    "    geneLengthDict, allNotEssentialNotImpactGenes - endometrialNotExpressedGenesGtexNOTIMPACT, endometrialNotExpressedGenesGtexNOTIMPACT, [])\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 219,
   "metadata": {},
   "outputs": [],
   "source": [
    "impactMMRSummary = pd.concat([colorectalMMRSummaryIMGenes, endometrialMMRSummaryIMGenes])\n",
    "impactMMRSummary['sigType'] = 'MMR'\n",
    "impactPOLESummary = pd.concat([colorectalPOLESummaryIMGenes, endometrialPOLESummaryIMGenes])\n",
    "impactPOLESummary['sigType'] = 'POLE'\n",
    "impactGenesSummary = pd.concat([impactMMRSummary, impactPOLESummary])\n",
    "impactGenesSummary['genes'] = 'IMPACT genes\\n(sans EOS genes)'\n",
    "\n",
    "notIMgenesMMRSummary = pd.concat([colorectalMMRSummaryNOTIMGenes, endometrialMMRSummaryNOTIMGenes])\n",
    "notIMgenesMMRSummary['sigType'] = 'MMR'\n",
    "notIMgenesPOLESummary = pd.concat([colorectalPOLESummaryNOTIMGenes, endometrialPOLESummaryNOTIMGenes])\n",
    "notIMgenesPOLESummary['sigType'] = 'POLE'\n",
    "notIMGenesSummary = pd.concat([notIMgenesPOLESummary, notIMgenesMMRSummary])\n",
    "notIMGenesSummary['genes'] = 'All non IMPACT &\\nnon essential genes'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 220,
   "metadata": {},
   "outputs": [],
   "source": [
    "finalCombinedDf = pd.concat([impactGenesSummary, notIMGenesSummary])\n",
    "finalCombinedDf['plottingLabel'] = finalCombinedDf.apply(lambda row: row['sigType'] + '_' + row['genes'], axis=1)\n",
    "finalCombinedDf.to_csv('/Users/friedman/Desktop/WORK/dataForLocalPlotting/gtexInfo.tsv', index=False, sep='\\t')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Now do gene based analysis in IMPACT**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 265,
   "metadata": {},
   "outputs": [],
   "source": [
    "def summarize_gene_mut_rates(mutsMaf, expressedTsgs, notExpressedTsgs, expressedOncogenes, notExpressedOncogenes, geneLengths):\n",
    "    listOfDicts = []\n",
    "    nCases = len(set(mutsMaf['Tumor_Sample_Barcode']))\n",
    "    \n",
    "    for gene in set(expressedTsgs) & set(geneLengths.keys()):\n",
    "        geneMutMaf = mutsMaf[mutsMaf['Hugo_Symbol'] == gene]\n",
    "        listOfDicts.append({'Gene': gene, 'mutRate': (1.0*geneMutMaf.shape[0])/(geneLengths[gene]*nCases), 'expression': 'expressed', 'geneType': 'TSG'})\n",
    "    \n",
    "    for gene in set(notExpressedTsgs) & set(geneLengths.keys()):\n",
    "        geneMutMaf = mutsMaf[mutsMaf['Hugo_Symbol'] == gene]\n",
    "        listOfDicts.append({'Gene': gene, 'mutRate': (1.0*geneMutMaf.shape[0])/(geneLengths[gene]*nCases), 'expression': 'notExpressed', 'geneType': 'TSG'})\n",
    "    \n",
    "    for gene in set(expressedOncogenes) & set(geneLengths.keys()):\n",
    "        geneMutMaf = mutsMaf[mutsMaf['Hugo_Symbol'] == gene]\n",
    "        listOfDicts.append({'Gene': gene, 'mutRate': (1.0*geneMutMaf.shape[0])/(geneLengths[gene]*nCases), 'expression': 'expressed', 'geneType': 'Oncogene'})\n",
    "    \n",
    "    for gene in set(notExpressedOncogenes) & set(geneLengths.keys()):\n",
    "        geneMutMaf = mutsMaf[mutsMaf['Hugo_Symbol'] == gene]\n",
    "        listOfDicts.append({'Gene': gene, 'mutRate': (1.0*geneMutMaf.shape[0])/(geneLengths[gene]*nCases), 'expression': 'notExpressed', 'geneType': 'Oncogene'})\n",
    "        \n",
    "    return pd.DataFrame(listOfDicts)   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 231,
   "metadata": {},
   "outputs": [],
   "source": [
    "poleCases, mmrCases, tmzCases = get_gene_and_cohort_list_utils.get_hypermutator_signature_cohorts(impactSigsPath = pathPrefix + '/juno/work/taylorlab/friedman/hypermutationAnalysisProj/projectDataAndConfigFiles/signatures_from_unfiltered_maf.txt', pathPrefix=pathPrefix)\n",
    "\n",
    "endometrialIds = analysis_utils.get_ids_by_hypermutant_status(hypermutantIdDir=pathPrefix + '/juno/work/taylorlab/friedman/hypermutationAnalysisProj/projectDataAndConfigFiles/hypermutationStatusIds', cancerType='Endometrial Cancer', hypermutantStatus = 'Hypermutated')\n",
    "colorectalIds = analysis_utils.get_ids_by_hypermutant_status(hypermutantIdDir=pathPrefix + '/juno/work/taylorlab/friedman/hypermutationAnalysisProj/projectDataAndConfigFiles/hypermutationStatusIds', cancerType='Colorectal Cancer', hypermutantStatus = 'Hypermutated')\n",
    "\n",
    "poleEndometrial = poleCases & endometrialIds\n",
    "mmrEndometrial = mmrCases & endometrialIds\n",
    "mmrColorectal = mmrCases & colorectalIds\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 235,
   "metadata": {},
   "outputs": [],
   "source": [
    "#just do it in POLE endometrial to start\n",
    "endometrialExpressedGenes = set(impactGenes) - endometrialNotExpressedGenesGtexIMPACT\n",
    "tsgs = 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",
    "oncogenes = set(impactGenes) - tsgs\n",
    "\n",
    "endometrialExpressedTsgs = endometrialExpressedGenes & tsgs\n",
    "endometrialNotExpressedTsgs = endometrialNotExpressedGenesGtexIMPACT & tsgs\n",
    "endometrialExpressedOncogenes = endometrialExpressedGenes & oncogenes\n",
    "endometrialNotExpressedOncogenes = endometrialNotExpressedGenesGtexIMPACT & oncogenes\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 266,
   "metadata": {},
   "outputs": [],
   "source": [
    "poleTruncatingMuts = allImpactMuts[(allImpactMuts['Tumor_Sample_Barcode'].isin(poleEndometrial)) &\n",
    "                                  (allImpactMuts['Consequence'] == 'stop_gained')]\n",
    "\n",
    "poleMutSummary = summarize_gene_mut_rates(poleTruncatingMuts, endometrialExpressedTsgs, endometrialNotExpressedTsgs, endometrialExpressedOncogenes, endometrialNotExpressedOncogenes, impactGeneLengths)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 267,
   "metadata": {},
   "outputs": [],
   "source": [
    "poleMutSummary['plottingLabel'] = poleMutSummary.apply(lambda row: row['geneType'] + '_' + row['expression'], axis=1)\n",
    "poleMutSummary.to_csv('/Users/friedman/Desktop/WORK/dataForLocalPlotting/impactGenesByExpressionPole.tsv', index=False, sep='\\t')\n"
   ]
  }
 ],
 "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
}