Analysis_Source_Code.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import libraries\n",
    "import os\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import random as rand"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "setlen=lambda x:len(set(x))\n",
    "def chgdir(dirname):\n",
    "    \"\"\"This wrapper creates a directory, if not already present, and then shifts code to that working directory.\"\"\"\n",
    "    if os.path.isdir(dirname):\n",
    "        os.chdir(dirname)\n",
    "        print('Changing directory to '+dirname)\n",
    "    else:\n",
    "        print('Directory '+dirname+' not found.\\nCreating new directory.\\nChanging to new directory.')\n",
    "        os.mkdir(dirname)\n",
    "        os.chdir(dirname)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "cwd=os.getcwd()\n",
    "cwd_input_matrices=os.path.join(cwd,'input_matrices')\n",
    "cwd_interim=os.path.join(cwd,'interim_tmp_files')\n",
    "cwd_genomics=os.path.join(cwd,'Genomics_Analysis')\n",
    "cwd_SEER=os.path.join(cwd,'SEER_Analysis')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import data from input_matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfC=pd.read_csv(os.path.join(cwd_genomics,'Genomics_Output_Processed.txt'),sep='\\t',index_col='Hugo_Symbol')\n",
    "dfCred=dfC.drop(['All'],axis='columns').drop('Total') # genomic matrix without the total rows\n",
    "dfCp=dfC.div(dfC.loc['Total':,:].values[0]) # P(gene|hist): Calculate probability of having a gene mutated given a histology.\n",
    "dfCp=(dfCp.drop(['All'], axis=1)).drop(['Total'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate P(hist) from epidemiological studies\n",
    "dfS0=pd.read_csv(os.path.join(cwd_SEER,'Output_SEER.txt'),sep='\\t',index_col='HISTO CODE')\n",
    "if (dfS0[-1:].values)[0,0] is np.nan:\n",
    "    dfS0=dfS0[:-1] # remove the last (NAN) line.\n",
    "dfS=dfS0.drop(index='CUSTOM SITES') # remove custom site names so histologies are only known by the codes (as in case of dfC)\n",
    "dfS=dfS.astype('float') # convert all numbers to float\n",
    "\n",
    "#P(hist AND loc): Calculate probability of incidence of a malignant cancer in a given location and histology in US population.\n",
    "sumStot=dfS.sum().sum()\n",
    "dfSp=((dfS.transpose()).div(sumStot)).fillna(0.) # The only way the last row would become unity is if the division occured accurately\n",
    "dfSp=dfSp.astype('float')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0000000000000002\n"
     ]
    }
   ],
   "source": [
    "# P(hist AND loc) including only ROSETTA histologies also incorporated in genomic data. For histologies with no sequenced cases, the genomic data has no information to offer.\n",
    "dfSC=dfS[dfCp.columns]\n",
    "sumSCtot=dfSC.sum().sum()\n",
    "dfSCp=((dfSC.transpose()).div(sumSCtot)).fillna(0.) # The only way the last row would become unity is if the division occured accurately\n",
    "dfSCp=dfSCp.astype('float')\n",
    "print(dfSCp.sum().sum()) # sum over histologies and locations for P(hist and tiss) =1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "dfp=dfCp@dfSCp # P(gene|loc)= sum over histologies in (p(gene|hist)*p(hist AND location))\n",
    "\n",
    "dfNhist=pd.DataFrame({'Incidence':dfSC.sum(axis='index')}) # P(hist)= sum over locations for P(hist and loc) in epidemiological data (SEER)\n",
    "dfSCNhist=dfCp.mul(dfNhist['Incidence'])# p(gene AND hist)=p(gene|hist)*p(hist) # probability reweighted to epidemiological data.\n",
    "dfSCNhist.columns=dfS0.loc['CUSTOM SITES'][dfSCNhist.columns]\n",
    "\n",
    "dfIncd=pd.DataFrame({'Incidence':dfp.sum(axis='columns')}) # p(gene)=sum over locations on dfp \n",
    "dfIncdsort=dfIncd.sort_values(by='Incidence',ascending=False)*100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create In-Silico samples for genomics studies and find statistical measures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of genes above cut-off incidence 0.01%: 20451 \n",
      "of total number of non-zero genes: 21271\n"
     ]
    }
   ],
   "source": [
    "# Only genes above cut off percent incidence are considered for statistical analysis - This is both due to significance of these genes and statistical analysis is less reliable for low incidence genes due to systematic effects.\n",
    "CutOffPercent=0.01 # select genes with incidence above a cut-off. Any percent level value cuts off the long tail of nearly zero incidence genes.\n",
    "NSilicoStudies=2000 # Select Number of silico studies ~ 2000 is enough to obtain stable central values for common, 1% level genes.\n",
    "dfIncdsort_sub=dfIncdsort[dfIncdsort>CutOffPercent].dropna()\n",
    "print('Number of genes above cut-off incidence '+str(CutOffPercent)+'%:',len(dfIncdsort_sub),'\\nof total number of non-zero genes:',(dfIncdsort>0).sum()[0])"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# Uncomment and run this line to skip de-novo generation of in-silico data below and speed processing along.\n",
    "df_Incdstats=pd.read_csv(os.path.join(cwd_interim,'df_Incdstats.csv'),sep='\\t',index_col='Hugo_Symbol')\n",
    "df_Incdstats=df_Incdstats.astype(float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "dfinp=(dfCred.copy(deep=True).loc[dfIncdsort_sub.index]).astype(int) # select genes with incidence above a cut-off. Any percent level value cuts off the long tail of nearly zero incidence genes.\n",
    "\n",
    "df_SilicoStudies=dfinp.applymap(lambda x:np.random.poisson(x,NSilicoStudies) if x>0 else np.array([0]))\n",
    "print('Silico Studies Generated')\n",
    "df_Silico_cases_ALL_Total=df_SilicoStudies+(dfC.drop(columns='All').loc['Total'].astype(int)-dfinp)\n",
    "print('Denominator defined')\n",
    "# For those rare cases where there are very few sequencing cases available.\n",
    "dfAllboolean=df_Silico_cases_ALL_Total.applymap(lambda x:x<1).applymap(lambda x:any(x))\n",
    "icollist=df_Silico_cases_ALL_Total.columns[dfAllboolean.any()]\n",
    "indexlist=df_Silico_cases_ALL_Total.index[dfAllboolean.transpose().any()]\n",
    "df_Silico_cases_ALL_Total.loc[indexlist,icollist]=df_Silico_cases_ALL_Total.loc[indexlist,icollist].applymap(lambda ilist:[x if x>0 else 1 for x in ilist])\n",
    "print('Denominator Updated')\n",
    "# Calculate conditional probability for silico samples\n",
    "df_Silico_Cp=df_SilicoStudies.div(df_Silico_cases_ALL_Total) #divide the each case count by total respective sample sizes for each histology.\n",
    "print('Conditional Probabilities Generated')\n",
    "del dfAllboolean\n",
    "del df_SilicoStudies\n",
    "\n",
    "df_Silico_p=df_Silico_Cp@dfSCp # Calcualte (p(gene|hist)*p(hist AND tissue)) for each in-silico study. Note that SEER epidimiological data is obtained from datasets orders of magnitude larger so errors are ignored.\n",
    "print('Total Probabilty Samples Generated')\n",
    "del df_Silico_Cp\n",
    "df_Silico_Incdsort=df_Silico_p.sum(axis='columns')*100\n",
    "del df_Silico_p\n",
    "print('Done.')\n",
    "df_Incdstats=pd.DataFrame(df_Silico_Incdsort,columns=['Incidence_Silico'])\n",
    "df_Incdstats['Incidence']=dfIncdsort.loc[df_Incdstats.index,'Incidence']\n",
    "df_Incdstats['CL97.5%']=df_Incdstats['Incidence_Silico'].apply(lambda arr:np.percentile(arr,97.5))\n",
    "df_Incdstats['CL2.5%']=df_Incdstats['Incidence_Silico'].apply(lambda arr:np.percentile(arr,2.5))\n",
    "df_Incdstats['Mean']=df_Incdstats['Incidence_Silico'].apply(np.mean)\n",
    "df_Incdstats=df_Incdstats[['Incidence','CL97.5%','CL2.5%','Mean']]\n",
    "df_Incdstats.to_csv(os.path.join(cwd_interim,'df_Incdstats.csv'),sep='\\t')\n",
    "del df_Silico_Incdsort"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import list of genes in RAS pathway, Kinome, Consensus 1 and 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "dfGA=pd.read_csv(os.path.join(cwd_input_matrices,'List_Ras_Kinome.txt'),sep='\\t') # annotated lists\n",
    "dfCT=pd.read_csv(os.path.join(cwd_input_matrices,'Lists_CT1and2.txt'),sep='\\t') # Updated 2020 CT1,2 lists.\n",
    "# remove NAN values and isolate the proteins for which length is available.\n",
    "kinomelist=set([elem for elem in dfGA['Kinome'] if elem is not np.nan])\n",
    "CT1list=set([elem for elem in dfCT['Consensus_Tier1'] if elem is not np.nan])\n",
    "CT2list=set([elem for elem in dfCT['Consensus_Tier2'] if elem is not np.nan])\n",
    "CT12list=CT1list.union(CT2list)\n",
    "Raslist=set([elem for elem in dfGA['RAS_Pathway'] if elem is not np.nan])\n",
    "\n",
    "def checkindices(l1):\n",
    "    \"\"\"Returns TRUE if all items present in dfIncdsort else returns a list of indices absent\"\"\"\n",
    "    flag=True\n",
    "    listfalse=[]\n",
    "    for item in l1:\n",
    "        flag=flag and (item in dfIncdsort.index)\n",
    "        if not (item in dfIncdsort.index):\n",
    "            listfalse=listfalse+[item]\n",
    "    return flag if flag else listfalse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Import validation data from Kandoth et al and Bailey et al."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfTCGA13=pd.read_csv(os.path.join(cwd_input_matrices,'TCGA2013_Kandoth_Supplementary_Table_4.txt'),sep='\\t',index_col='Gene')\n",
    "dfTCGA13.index=dfTCGA13.index.map(lambda x:x.upper())\n",
    "dfTCGA18=pd.read_csv(os.path.join(cwd_input_matrices,'TCGA_2018_Bailey.txt'),sep='\\t',index_col='Gene')\n",
    "dfTCGA18.index=dfTCGA18.index.map(lambda x:x.upper())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import alternate gene nomenclature file from cbioportal: https://docs.cbioportal.org/3.-cbioportal-maintenance/updating-gene-and-gene_alias-tables\n",
    "#Homo_sapien.gene_info.gz ftp://ftp.ncbi.nih.gov/gene/DATA/GENE_INFO/Mammalia/Homo_sapiens.gene_info.gz\n",
    "dfGeneNames=pd.read_csv(os.path.join(cwd_input_matrices,'Homo_sapiens_gene_info_GM.txt'),sep='\\t')\n",
    "dfGeneNames=dfGeneNames.astype(str).applymap(lambda x:x.upper())\n",
    "dfGeneNames.Synonyms=[str(row).split('|') for row in dfGeneNames.Synonyms.values]\n",
    "listSynonyms=[elem for row in dfGeneNames.Synonyms.values for elem in row]\n",
    "def FindGeneName(igene):\n",
    "    # This function has already been applied on the input gene list. In that application, care was taken to match chromosomes before deciding that two genes are same. This function therefore does not filter any names already present in the input genome matrix dfC.\n",
    "    igene=igene.upper()\n",
    "    retgene=igene\n",
    "    if (igene in set(dfGeneNames.Symbol)) or (igene not in set(listSynonyms) or (igene in dfC.index.values)):\n",
    "        return retgene\n",
    "    else:\n",
    "        retgene=dfGeneNames[[igene in row for row in dfGeneNames.Synonyms]].Symbol.values[0]\n",
    "        return retgene"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Changing directory to Results\n"
     ]
    }
   ],
   "source": [
    "def chgdir(dirname):\n",
    "    if os.path.isdir(dirname):\n",
    "        os.chdir(dirname)\n",
    "        print('Changing directory to '+dirname)\n",
    "    else:\n",
    "        print('Directory '+dirname+' not found.\\nCreating new directory.\\nChanging to new directory.')\n",
    "        os.mkdir(dirname)\n",
    "        os.chdir(dirname)\n",
    "chgdir('Results')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figure 2A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "('Rows and columns in Genomics Matrix:', 21271, 73)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"Rows and columns in Genomics Matrix:\",len(dfC.index)-1, len(dfC.columns)-1 # One extra \"All\" column for genomic matrix is removed and Total row for index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "('The number of samples sequenced in Genomics studies included:', 19181)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "'The number of samples sequenced in Genomics studies included:', dfC.loc['Total','All']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "('Total number of cancer diagnoses included (SEER):', 7167808.0)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "'Total number of cancer diagnoses included (SEER):',round(dfS.sum().sum(),2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total number of histologies included in ICD03 <or> ROSETTA: 716 \n",
      "Number of histologies in ROSETTA with at least one case included: 370\n"
     ]
    }
   ],
   "source": [
    "print('Total number of histologies included in ICD03 <or> ROSETTA:', len(dfS.sum()),'\\nNumber of histologies in ROSETTA with at least one case included:',len(dfS.sum()[dfS.sum()>0])) # There were 670 histologies in ICD-03 hist/behav codes in SEER data.+12 represents histologies with zero cases but included in ROSETTA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "92.95"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "round(100*sumSCtot/sumStot,2) # Percentage of cases from SEER database which are included in our analyses."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figure 2B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "127 127\n"
     ]
    }
   ],
   "source": [
    "# Replace Gene Names From Redundant gene names list\n",
    "print(len(dfTCGA13.index),len(set([FindGeneName(igene) for igene in dfTCGA13.index])))\n",
    "if len(dfTCGA13.index)==len(set([FindGeneName(igene) for igene in dfTCGA13.index])):\n",
    "    dfTCGA13.index=[FindGeneName(igene) for igene in dfTCGA13.index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['MALAT1', 'MIR142']"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "checkindices(dfTCGA13.index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "if checkindices(dfTCGA13.index)!=True: # check that all indices are present in our prevalance vector\n",
    "    dfTCGA13=dfTCGA13.drop(checkindices(dfTCGA13.index))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "checkindices(dfTCGA13.index) # check that all indices are present in our prevalance vector"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figure 2D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfTCGA13.drop('TP53',inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "font = {'family' : 'Arial',\n",
    "        'weight' : 'normal',\n",
    "        'size'   : 16}\n",
    "plt.rc('font', **font)\n",
    "plt.figure(figsize=(4,4))\n",
    "dfIncdsortTCGA13genes=df_Incdstats.loc[dfTCGA13.index] # extract the TCGA genes from our prevalance vector\n",
    "xvals=dfTCGA13['Pancan12 Freq%'].values\n",
    "yvals=dfIncdsortTCGA13genes['Incidence'].values\n",
    "yerr=[abs(dfIncdsortTCGA13genes['CL2.5%']-dfIncdsortTCGA13genes['Incidence']).values,(dfIncdsortTCGA13genes['CL97.5%']-dfIncdsortTCGA13genes['Incidence']).values]\n",
    "dfdiff=pd.concat([dfTCGA13[['Pancan12 Freq%']],dfIncdsortTCGA13genes[['Incidence','CL97.5%','CL2.5%']],abs((dfTCGA13.rename(columns={'Pancan12 Freq%':'Incidence'}))[['Incidence']]-dfIncdsortTCGA13genes[['Incidence']]).rename(columns={'Incidence':'abs(PrevDiff)'})],axis=1)\n",
    "xyline=np.arange(0,20)\n",
    "plt.plot(xyline,xyline,color='gray')\n",
    "plt.errorbar(x=xvals,y=yvals,yerr=yerr,color='k',fmt='o',ms=6,elinewidth=3)\n",
    "plt.xlabel('Proportion of TCGA Samples (%)\\n (Kandoth 2013)',size=20)\n",
    "plt.ylabel('Estimated Mutation Frequency\\n in U.S. population (%)',size=20)\n",
    "dfdiffsort=dfdiff.loc[dfdiff['abs(PrevDiff)']>3]\n",
    "# dictextloc={'PIK3CA':1.9,'PTEN':1.5, 'VHL':1, 'KRAS':1.55, 'BRAF':1.45, 'KMT2D':1.85}\n",
    "for id1 in dfdiffsort.index:\n",
    "    plt.text(dfdiffsort['Pancan12 Freq%'][id1]+0.2,dfdiffsort['Incidence'][id1]+0.2,id1,fontsize=16,color='k')\n",
    "plt.xlim([0,20])\n",
    "plt.ylim([0,20])\n",
    "plt.gca().spines['right'].set_linewidth(3)\n",
    "plt.gca().spines['top'].set_linewidth(3)\n",
    "plt.gca().spines['bottom'].set_linewidth(3)\n",
    "plt.gca().spines['left'].set_linewidth(3)\n",
    "plt.gca().xaxis.set_tick_params(width=3,size=10,labelsize=20)\n",
    "plt.gca().yaxis.set_tick_params(width=3,size=10,labelsize=20)\n",
    "\n",
    "fname='Figure2D_TCGApancanKandoth_scatter'\n",
    "dfdiff.to_excel(fname+'.xlsx')\n",
    "plt.savefig(fname+'.pdf',bbox_inches='tight')\n",
    "plt.close() # replace with plt.show() to display chart"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figure 2E"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "299 299\n"
     ]
    }
   ],
   "source": [
    "# Replace Gene Names From Redundant gene names list\n",
    "print(len(dfTCGA18.index),len(set([FindGeneName(igene) for igene in dfTCGA18.index])))\n",
    "if len(dfTCGA18.index)==len(set([FindGeneName(igene) for igene in dfTCGA18.index])):\n",
    "    dfTCGA18.index=[FindGeneName(igene) for igene in dfTCGA18.index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "checkindices(dfTCGA18.index) # check that all indices are present in our prevalance vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfTCGA18=dfTCGA18[dfTCGA18.sum(axis='columns')>0.]*100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfTCGA18.drop('TP53',inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "font = {'family' : 'Arial',\n",
    "        'weight' : 'normal',\n",
    "        'size'   : 16}\n",
    "plt.rc('font', **font)\n",
    "plt.figure(figsize=(4,4))\n",
    "plt.rcParams['pdf.fonttype'] = 42 # To make output Adobe Illustrator friendly for formatting the output pie charts\n",
    "dfIncdsortTCGA18genes=df_Incdstats.loc[dfTCGA18.index] # extract the TCGA genes from our prevalance vector\n",
    "xvals=dfTCGA18['Pancan Freq%'].values\n",
    "yvals=dfIncdsortTCGA18genes['Incidence'].values\n",
    "yerr=[abs(dfIncdsortTCGA18genes['CL2.5%']-dfIncdsortTCGA18genes['Incidence']).values,(dfIncdsortTCGA18genes['CL97.5%']-dfIncdsortTCGA18genes['Incidence']).values]\n",
    "dfdiff=pd.concat([dfTCGA18[['Pancan Freq%']],dfIncdsortTCGA18genes[['Incidence','CL97.5%','CL2.5%']],abs((dfTCGA18.rename(columns={'Pancan Freq%':'Incidence'}))[['Incidence']]-dfIncdsortTCGA18genes[['Incidence']]).rename(columns={'Incidence':'abs(PrevDiff)'})],axis=1)\n",
    "xyline=np.arange(0,20)\n",
    "plt.plot(xyline,xyline,color='gray')\n",
    "plt.errorbar(x=xvals,y=yvals,yerr=yerr,color='k',fmt='o',ms=6,elinewidth=3)\n",
    "plt.xlabel('Proportion of TCGA Samples (%)\\n (Bailey 2018)',size=20)\n",
    "plt.ylabel('Estimated Mutation Proportion\\n in U.S. population (%)',size=20)\n",
    "dfdiffsort=dfdiff.loc[dfdiff['abs(PrevDiff)']>3]\n",
    "dictextloc={'PIK3CA':1.9,'PTEN':1.5, 'VHL':1, 'KRAS':0.2, 'BRAF':1.45, 'KMT2D':1.85,'APC':2.,'IDH1':0.9}\n",
    "for id1 in dfdiffsort.index:\n",
    "    plt.text(dfdiffsort['Pancan Freq%'][id1]-dictextloc[id1],dfdiffsort['Incidence'][id1]+0.5,id1,fontsize=16,color='k')\n",
    "plt.xlim([0,20])\n",
    "plt.ylim([0,20])\n",
    "plt.gca().spines['right'].set_linewidth(3)\n",
    "plt.gca().spines['top'].set_linewidth(3)\n",
    "plt.gca().spines['bottom'].set_linewidth(3)\n",
    "plt.gca().spines['left'].set_linewidth(3)\n",
    "plt.gca().xaxis.set_tick_params(width=3,size=10,labelsize=20)\n",
    "plt.gca().yaxis.set_tick_params(width=3,size=10,labelsize=20)\n",
    "\n",
    "fname='Figure2E_TCGApancanBailey_scatter'\n",
    "dfdiff.to_excel(fname+'.xlsx')\n",
    "plt.savefig(fname+'.pdf',bbox_inches='tight')\n",
    "plt.close() # replace with plt.show() to display chart"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Compare impact of epidemiological reweighting for our genomic dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of genes above cut-off incidence 0.01%: 20451 \n",
      "of total number of non-zero genes: 21271\n",
      "Wall time: 8.63 s\n"
     ]
    }
   ],
   "source": [
    "%%time \n",
    "# In-Silico Studies for Statistical Analysis\n",
    "# Only genes above cut off percent incidence are considered for statistical analysis - This is both due to significance of these genes and statistical analysis is less reliable for low incidence genes due to systematic effects.\n",
    "print('Number of genes above cut-off incidence '+str(CutOffPercent)+'%:',len(dfIncdsort_sub),'\\nof total number of non-zero genes:',(dfIncdsort>0).sum()[0])\n",
    "\n",
    "def fun_Silico_Unwtd(dflist):\n",
    "    \"\"\"Inputs a dataframe of numbers and creates Poisson in-silico studies\"\"\"\n",
    "    import numpy as np\n",
    "    inpdf=dflist[0]\n",
    "    dfC_AllTot=dflist[1]\n",
    "    NumSilicoStudies=dflist[2] # NUMBER of IN-SILICO SAMPLES generated for each gene and ROSETTA histology combination\n",
    "    df_SilicoStudies=inpdf.apply(lambda x:np.random.poisson(x,NumSilicoStudies))\n",
    "    df_Silicodiff=df_SilicoStudies-inpdf\n",
    "    df_Silico_cases_ALL_Total=df_Silicodiff+dfC_AllTot\n",
    "    df_Silico_Cp=df_SilicoStudies.div(df_Silico_cases_ALL_Total)*100 #divide the each case count by total respective sample sizes for each histology.\n",
    "    del df_SilicoStudies\n",
    "    del df_Silicodiff\n",
    "    del df_Silico_cases_ALL_Total\n",
    "    return df_Silico_Cp\n",
    "\n",
    "#Process in-silico studies in parallel\n",
    "dfinp=dfCred.sum(axis='columns').loc[dfIncdsort_sub.index] # select genes with incidence above a cut-off. Any percent level value cuts off the long tail of nearly zero incidence genes.\n",
    "datalist=[dfinp,dfC.All.Total,NSilicoStudies] \n",
    "\n",
    "df_Silico_Cp=fun_Silico_Unwtd(datalist)\n",
    "\n",
    "# stats for unweighted genomic Incdalence rates\n",
    "df_unwtd_Incdstats=pd.DataFrame(df_Silico_Cp,columns=['Incidence_Silico'])\n",
    "df_unwtd_Incdstats['Mean']=df_unwtd_Incdstats['Incidence_Silico'].apply(np.mean)\n",
    "df_unwtd_Incdstats['CL97.5%']=df_unwtd_Incdstats['Incidence_Silico'].apply(lambda arr:np.percentile(arr,97.5))\n",
    "df_unwtd_Incdstats['CL2.5%']=df_unwtd_Incdstats['Incidence_Silico'].apply(lambda arr:np.percentile(arr,2.5))\n",
    "df_unwtd_Incdstats[['Mean','CL97.5%','CL2.5%']].to_csv(os.path.join(cwd_interim,'df_unwtd_Incdstats.csv'),sep='\\t')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "font = {'family' : 'Arial',\n",
    "        'weight' : 'normal',\n",
    "        'size'   : 20}\n",
    "plt.rc('font', **font)\n",
    "plt.rcParams['pdf.fonttype'] = 42 # To make output Adobe Illustrator friendly for formatting the output pie charts\n",
    "\n",
    "plt.figure(figsize=(4,4))\n",
    "genelist_G=dfIncdsort_sub.index.values\n",
    "df_xplt=df_unwtd_Incdstats.loc[genelist_G].copy(deep=True)\n",
    "df_yplt=df_Incdstats.loc[genelist_G].copy(deep=True)\n",
    "xvals=df_xplt.Mean.loc[genelist_G].values\n",
    "yvals=df_yplt.Mean.loc[genelist_G].values\n",
    "plt.scatter(x=xvals,y=yvals,marker='o',color='k',s=10)\n",
    "xyline=np.arange(0,40)\n",
    "plt.plot(xyline,xyline,color='gray')\n",
    "\n",
    "genesublist_G=[row for row in genelist_G if df_xplt.Mean[row]>10 or df_yplt.Mean[row]>10]\n",
    "xvals=df_xplt.Mean.loc[genesublist_G].values\n",
    "yvals=df_yplt.Mean.loc[genesublist_G].values\n",
    "xerror=[abs(df_xplt.loc[genesublist_G,'CL2.5%']-df_xplt.loc[genesublist_G,'Mean']).values,(df_xplt.loc[genesublist_G,'CL97.5%']-df_xplt.loc[genesublist_G,'Mean']).values]\n",
    "yerror=[abs(df_yplt.loc[genesublist_G,'CL2.5%']-df_yplt.loc[genesublist_G,'Mean']).values,(df_yplt.loc[genesublist_G,'CL97.5%']-df_yplt.loc[genesublist_G,'Mean']).values]\n",
    "plt.errorbar(x=xvals,y=yvals,xerr=xerror,yerr=yerror,color='k',fmt='o',ms=1,elinewidth=3)\n",
    "\n",
    "labelgenes=[igene for igene in genelist_G if abs(df_xplt.Mean.loc[igene]-df_yplt.Mean.loc[igene])>3.5]\n",
    "# dicloctexty={ 'TTN', 'MUC16', 'CSMD3', 'PIK3CA', 'LRP1B', 'RYR2', 'USH2A', 'ZFHX4', 'APC','TP53':0.4,'APC':0.3,'NOTCH1':0.2,'NFE2L2':0.2}\n",
    "# dicloctextx={'FAT1':-4.5, 'KMT2D':0.2, 'LRP1B':1., 'KRAS':1, 'CDKN2A':0.4, 'TP53':0.7,'APC':0.2,'NOTCH1':0.2,'NFE2L2':0.2}\n",
    "for id1 in labelgenes:\n",
    "    plt.text(x=df_xplt.loc[id1,'Mean']-7,y=df_yplt.loc[id1,'Mean']+1.5,s=id1,fontsize=16,color='k')\n",
    "plt.text(x=36,y=33,s='x = y',fontsize=20,color='gray',rotation=30)\n",
    "plt.xlim([0,40])\n",
    "plt.ylim([0,40])\n",
    "plt.gca().spines['right'].set_linewidth(2)\n",
    "plt.gca().spines['top'].set_linewidth(2)\n",
    "plt.gca().spines['bottom'].set_linewidth(2)\n",
    "plt.gca().spines['left'].set_linewidth(2)\n",
    "plt.gca().xaxis.set_tick_params(width=3,size=10,labelsize=20)\n",
    "plt.gca().yaxis.set_tick_params(width=3,size=10,labelsize=20)\n",
    "\n",
    "plt.xlabel('Unweighted Genomics Data Mutation Proportion (%)',size=20)\n",
    "plt.ylabel('Weighted Mutation Proportion\\n in U.S. Population (%)',size=20)\n",
    "fname='Figure_GenomicsCompare_scatter'\n",
    "plt.savefig(fname+'.pdf',bbox_inches='tight')\n",
    "plt.close() # replace with plt.show() to display chart"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_xplt.columns=[icol+'_Unwtd' for icol in df_xplt.columns]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfdiff=pd.concat([df_xplt[['Mean_Unwtd','CL97.5%_Unwtd','CL2.5%_Unwtd']],df_yplt[['Mean','CL97.5%','CL2.5%']],abs(df_xplt[['Mean_Unwtd']].rename(columns={'Mean_Unwtd':'Mean'})-df_yplt[['Mean']]).rename(columns={'Mean':'abs(MeanDiff)'}),(abs(df_xplt[['Mean_Unwtd']].rename(columns={'Mean_Unwtd':'Mean'})-df_yplt[['Mean']]).div(df_xplt[['Mean_Unwtd']].rename(columns={'Mean_Unwtd':'Mean'}))).rename(columns={'Mean':'MeanDiffRatio(%)'})*100],axis=1)\n",
    "dfdiff.loc[dfdiff['abs(MeanDiff)']>0.5].loc[dfdiff['MeanDiffRatio(%)']>40].sort_values(by='Mean',ascending=False).to_excel('dfdiff_Reweighting_Impact40.xlsx')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figure 3A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "CT1list=[FindGeneName(igene) for igene in CT1list]\n",
    "CT2list=[FindGeneName(igene) for igene in CT2list]\n",
    "CT12list=[FindGeneName(igene) for igene in CT12list]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True ['HMGN2P46', 'MALAT1'] ['MALAT1', 'HMGN2P46']\n"
     ]
    }
   ],
   "source": [
    "# identify Incidence rates by CT1, CT2 and CT12 lists , SORT and store as percentages\n",
    "print(checkindices(CT1list),checkindices(CT2list),checkindices(CT12list))\n",
    "CT2list=[item for item in CT2list if item not in checkindices(CT2list)]\n",
    "CT12list=[item for item in CT12list if item not in checkindices(CT12list)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfCT1prev=(dfIncd.loc[CT1list].sort_values(by=['Incidence'],ascending=False)*100)\n",
    "dfCT1prev['CL97.5%']=[df_Incdstats.loc[igene,'CL97.5%'] if igene in df_Incdstats.index else np.nan for igene in dfCT1prev.index]\n",
    "dfCT1prev['CL2.5%']=[df_Incdstats.loc[igene,'CL2.5%'] if igene in df_Incdstats.index else np.nan for igene in dfCT1prev.index]\n",
    "dfCT1prev['Mean_Silico']=[df_Incdstats.loc[igene,'Mean'] if igene in df_Incdstats.index else np.nan for igene in dfCT1prev.index]\n",
    "dfCT1prev.round(5).to_excel('Figure3_Incidences_CT1.xlsx')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfCT2prev=(dfIncd.loc[CT2list].sort_values(by=['Incidence'],ascending=False)*100)\n",
    "dfCT2prev['CL97.5%']=[df_Incdstats.loc[igene,'CL97.5%'] if igene in df_Incdstats.index else np.nan for igene in dfCT2prev.index]\n",
    "dfCT2prev['CL2.5%']=[df_Incdstats.loc[igene,'CL2.5%'] if igene in df_Incdstats.index else np.nan for igene in dfCT2prev.index]\n",
    "dfCT2prev['Mean_Silico']=[df_Incdstats.loc[igene,'Mean'] if igene in df_Incdstats.index else np.nan for igene in dfCT2prev.index]\n",
    "dfCT2prev.round(5).to_excel('Figure3_Incidences_CT2.xlsx')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfCT12prev=(dfIncd.loc[CT12list].sort_values(by=['Incidence'],ascending=False)*100)\n",
    "dfCT12prev['CL97.5%']=[df_Incdstats.loc[igene,'CL97.5%'] if igene in df_Incdstats.index else np.nan for igene in dfCT12prev.index]\n",
    "dfCT12prev['CL2.5%']=[df_Incdstats.loc[igene,'CL2.5%'] if igene in df_Incdstats.index else np.nan for igene in dfCT12prev.index]\n",
    "dfCT12prev['Mean_Silico']=[df_Incdstats.loc[igene,'Mean'] if igene in df_Incdstats.index else np.nan for igene in dfCT12prev.index]\n",
    "dfCT12prev.round(5).to_excel('Figure3_Incidences_CT12.xlsx')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Figure3A. Bar Chart of top 50 CT1 genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(714, 714)"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(CT12list),setlen(CT12list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "axis_font = {'size':'10'}\n",
    "dfbarplt=dfCT1prev[:50]\n",
    "barlabels=dfbarplt.index.values\n",
    "barpos=np.arange(len(barlabels))\n",
    "bardata=dfbarplt['Incidence'].values\n",
    "barerr=[abs(dfbarplt['CL2.5%']-dfbarplt['Incidence']).values,(dfbarplt['CL97.5%']-dfbarplt['Incidence']).values]\n",
    "figobj,axobj=plt.subplots(figsize=(15,10))\n",
    "axobj.bar(barpos,bardata,yerr=barerr,align='center',alpha=0.5,edgecolor='k')\n",
    "plt.xticks(barpos,barlabels,rotation=90,**axis_font)\n",
    "plt.ylabel('% Incidence in US population')\n",
    "plt.savefig('Figure3A_Top50_CT1.png',dpi=300,bbox_inches='tight')\n",
    "plt.close() # replace with plt.show() to display chart"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figure 3B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(514, 514)"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "setlen(kinomelist),len(kinomelist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['CAMK1B', 'MAP3K14', 'CCRK', 'PRKY']\n"
     ]
    }
   ],
   "source": [
    "# identify Incidence rates by kinomelist, SORT and store as percentages and stats\n",
    "kinomelist=[FindGeneName(igene) for igene in kinomelist]\n",
    "print(checkindices(kinomelist))\n",
    "kinomelist=[item for item in kinomelist if item not in checkindices(kinomelist)]\n",
    "dfkinomeprev=(dfIncd.loc[kinomelist].sort_values(by=['Incidence'],ascending=False)*100).round(5)\n",
    "dfkinomeprev['CL97.5%']=[df_Incdstats.loc[igene,'CL97.5%'] if igene in df_Incdstats.index else np.nan for igene in dfkinomeprev.index]\n",
    "dfkinomeprev['CL2.5%']=[df_Incdstats.loc[igene,'CL2.5%'] if igene in df_Incdstats.index else np.nan for igene in dfkinomeprev.index]\n",
    "dfkinomeprev['Mean_Silico']=[df_Incdstats.loc[igene,'Mean'] if igene in df_Incdstats.index else np.nan for igene in dfkinomeprev.index]\n",
    "dfkinomeprev.to_excel('Figure3B_Incidences_kinome.xlsx')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "axis_font = {'size':'10'}\n",
    "dfbarplt=dfkinomeprev[:50]\n",
    "barlabels=dfbarplt.index.values\n",
    "barpos=np.arange(len(barlabels))\n",
    "bardata=dfbarplt['Incidence'].values\n",
    "barerr=[abs(dfbarplt['CL2.5%']-dfbarplt['Incidence']).values,(dfbarplt['CL97.5%']-dfbarplt['Incidence']).values]\n",
    "figobj,axobj=plt.subplots(figsize=(15,10))\n",
    "axobj.bar(barpos,bardata,yerr=barerr,align='center',alpha=0.5,color='g',edgecolor='k')\n",
    "plt.xticks(barpos,barlabels,rotation=90,**axis_font)\n",
    "plt.ylabel('% Incidence in US population')\n",
    "plt.savefig('Figure3B_Top50_kinome.png',dpi=150,bbox_inches='tight')\n",
    "plt.close() # replace with plt.show() to display chart"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figure 3C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['RASSF10']\n"
     ]
    }
   ],
   "source": [
    "# identify Incidence rates by raslist, SORT and store as percentages and stats\n",
    "Raslist=[FindGeneName(igene) for igene in Raslist]\n",
    "print(checkindices(Raslist))\n",
    "Raslist=[item for item in Raslist if item not in checkindices(Raslist)]\n",
    "dfrasprev=(dfIncd.loc[Raslist].sort_values(by=['Incidence'],ascending=False)*100).round(5)\n",
    "dfrasprev['CL97.5%']=[df_Incdstats.loc[igene,'CL97.5%'] if igene in df_Incdstats.index else np.nan for igene in dfrasprev.index]\n",
    "dfrasprev['CL2.5%']=[df_Incdstats.loc[igene,'CL2.5%'] if igene in df_Incdstats.index else np.nan for igene in dfrasprev.index]\n",
    "dfrasprev['Mean_Silico']=[df_Incdstats.loc[igene,'Mean'] if igene in df_Incdstats.index else np.nan for igene in dfrasprev.index]\n",
    "dfrasprev.to_excel('Figure3C_Incidences_RasPathway.xlsx')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(226, 226)"
      ]
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "setlen(dfrasprev.index.values),len(dfrasprev.index.values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "axis_font = {'size':'10'}\n",
    "dfbarplt=dfrasprev[:50]\n",
    "barlabels=dfbarplt.index.values\n",
    "barpos=np.arange(len(barlabels))\n",
    "bardata=dfbarplt['Incidence'].values\n",
    "barerr=[abs(dfbarplt['CL2.5%']-dfbarplt['Incidence']).values,(dfbarplt['CL97.5%']-dfbarplt['Incidence']).values]\n",
    "figobj,axobj=plt.subplots(figsize=(15,10))\n",
    "axobj.bar(barpos,bardata,yerr=barerr,align='center',alpha=0.5,color='r',edgecolor='k')\n",
    "plt.xticks(barpos,barlabels,rotation=90,**axis_font)\n",
    "plt.ylabel('% Incidence in US population')\n",
    "plt.savefig('Figure3C_Top50_ras.png',dpi=150,bbox_inches='tight')\n",
    "plt.close() # replace with plt.show() to display chart"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of adenocarcinomas in SEER epidemiological data: 4676878.54917131 \n",
      "Percent of SEER histologies which are adenocarcinoma: 65.2483792698034\n",
      "Adenocarcinoma ROSETTA histologies NOT in current genomic data: [['70407: Other Adenocarcinoma']]\n"
     ]
    }
   ],
   "source": [
    "# Define Adenocarcinoma histologies by code\n",
    "adenocarclist=['81603','83173','83703','70117','70267','70397','70187','70197','70277','70287','70297','70307','70317','70327','70337','70347','70357','70367','70387','70407']\n",
    "print('Number of adenocarcinomas in SEER epidemiological data:',dfS[adenocarclist].sum().sum(),'\\nPercent of SEER histologies which are adenocarcinoma:', dfS[adenocarclist].sum().sum()/sumStot*100)\n",
    "aclist=[val for val in adenocarclist if val in dfCred.columns]\n",
    "print('Adenocarcinoma ROSETTA histologies NOT in current genomic data:',[[val+': '+dfS0.loc['CUSTOM SITES',val]] for val in adenocarclist if val not in dfCred.columns])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0\n"
     ]
    }
   ],
   "source": [
    "histsredAC=aclist\n",
    "dfCp_ac=dfCp[histsredAC] # genomic conditional probablity matrix without All genes and total histology vectors\n",
    "# P(hist AND loc) including only histologies also incorporated in genomic data. For histologies not incorporated, the genomic data has no additional information to offer.\n",
    "sumSCtot_ac=dfSC[histsredAC].sum().sum()\n",
    "dfSCp_ac=dfSCp.loc[histsredAC]\n",
    "dfSCp_ac=(dfSCp_ac.div(sumSCtot_ac)).mul(sumSCtot) # Multiply by total cases over all histologies and divide by total cases in included histologies\n",
    "print(dfSCp_ac.sum().sum()) #Cross-check: sum over histologies and tissues for P(hist and tiss) =1\n",
    "dfp_ac=dfCp_ac@dfSCp_ac # sum over histologies for (p(gene|hist)*p(hist AND tissue))\n",
    "\n",
    "# dfpG=dfC.drop()\n",
    "\n",
    "# N(hist)\n",
    "dfNloc_ac=pd.DataFrame({'Incidence':dfSC[histsredAC].sum(axis='columns')})\n",
    "dfIncd_ac=pd.DataFrame({'Incidence':dfp_ac.sum(axis='columns')}) # p(mut)=sum over locations on dfp \n",
    "dfIncdsort_ac=dfIncd_ac.sort_values(by='Incidence',ascending=False)*100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 441 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# Uncomment and run when creating Table 1. (In jupyter, the key 'y' uncomments the whole block.)\n",
    "dfSCredNorm=dfSCNhist.loc[CT1list].div(dfSCNhist.loc[CT1list].sum(axis=1),axis='index')*100 # Histological case distribution (SEER reweighted) for each gene\n",
    "commonhistdic={}\n",
    "for id1 in dfSCredNorm.index:\n",
    "    tophistset=dfSCredNorm.loc[id1].sort_values(ascending=False)\n",
    "    tophistset=tophistset.sort_values(ascending=False)[:3]\n",
    "    tophiststr=[id2+'('+str(tophistset[id2].round(1))+'%) ' for id2 in tophistset.index] # replace histology ROSETTA number with name\n",
    "    commonhistdic[id1]=''.join(tophiststr)\n",
    "df_commonhist=pd.DataFrame.from_dict(commonhistdic, orient='index',columns=['Common Cancers'])\n",
    "# Run to calculate Table 1\n",
    "dfIncdsort_ac.columns=['Incidence_AdenoCarcinomas']\n",
    "df_Incid_GlobalvsAdeno=pd.concat([dfIncdsort.loc[CT1list].round(1),dfIncdsort_ac.loc[CT1list].round(1),df_commonhist.loc[CT1list]],axis='columns').sort_values(by='Incidence',ascending=False)\n",
    "df_Incid_GlobalvsAdeno.to_excel('Table1_ConsensusGenesFreq.xlsx',index_label='Gene')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figure 4A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [],
   "source": [
    "squamouslist=dfS0.columns[dfS0.columns>='70207']\n",
    "squamouslist=squamouslist[squamouslist<='70257']\n",
    "sqlist=[elem for elem in squamouslist if elem in dfCred.columns] # extract those squamous ROSETTA histologies which are present in genomic data.\n",
    "\n",
    "Melanomalist=['70097','70107','93633','87443']\n",
    "Melanomalist=[elem for elem in Melanomalist if elem in dfCred.columns]# extract those melanoma ROSETTA histologies which are present in genomic data.\n",
    "\n",
    "TCClist=['70167'] # To analyze any other single histology - this is an easy trick. enter histology here and output TCC data.\n",
    "TCClist=[elem for elem in TCClist if elem in dfCred.columns]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Percent of all cancers in epidemiological studies that are within major sub-classifications of adenocarcinoma, squamous cell carcinoma, melanoma, or transitional cell carcinoma:\n",
      " [64.0, 9.0, 5.0, 5.0]\n"
     ]
    }
   ],
   "source": [
    "TypeList=[aclist,squamouslist,Melanomalist,TCClist]\n",
    "print(\"Percent of all cancers in epidemiological studies that are within major sub-classifications of adenocarcinoma, squamous cell carcinoma, melanoma, or transitional cell carcinoma:\\n\",[round(100*dfS[histType].sum().sum()/sumStot,0) for histType in TypeList ])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figure 4B,C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [],
   "source": [
    "# A few categories are manually combined in the manuscript figure. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfbar=dfNloc_ac.sort_values(by='Incidence',ascending=False)\n",
    "dfbar=100*dfbar/dfbar.sum()\n",
    "axis_font = {'size':'10'}\n",
    "barlabels=dfbar[:20].index.values\n",
    "barpos=np.arange(len(barlabels))\n",
    "bardata=dfbar[:20]['Incidence'].values\n",
    "figobj,axobj=plt.subplots(figsize=(10,10))\n",
    "axobj.bar(barpos,bardata,align='center',alpha=0.5)\n",
    "plt.xticks(barpos,barlabels,rotation=90,**axis_font)\n",
    "plt.ylabel('Percent Incidence')\n",
    "dfbar.to_excel('Figure4B_Adenocarcinoma_Types.xlsx')\n",
    "plt.close() # replace with plt.show() to display chart"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9999999999999998\n"
     ]
    }
   ],
   "source": [
    "histsredSQ=sqlist\n",
    "dfCp_sq=dfCp[histsredSQ] # genomic conditional probablity matrix without All genes and total histology vectors\n",
    "# P(hist AND loc) including only histologies also incorporated in genomic data. For histologies not incorporated, the genomic data has no additional information to offer.\n",
    "sumSCtot_sq=dfSC[histsredSQ].sum().sum()\n",
    "dfSCp_sq=dfSCp.loc[histsredSQ]\n",
    "dfSCp_sq=(dfSCp_sq.div(sumSCtot_sq)).mul(sumSCtot) # Multiply by total cases over all histologies and divide by total cases in included histologies\n",
    "print(dfSCp_sq.sum().sum()) #Cross-check: sum over histologies and tissues for P(hist and tiss) =1\n",
    "dfp_sq=dfCp_sq@dfSCp_sq # sum over histologies for (p(gene|hist)*p(hist AND tissue))\n",
    "\n",
    "# dfpG=dfC.drop()\n",
    "\n",
    "# N(hist)\n",
    "dfNloc_sq=pd.DataFrame({'Incidence':dfSC[histsredSQ].sum(axis='columns')})\n",
    "dfIncd_sq=pd.DataFrame({'Incidence':dfp_sq.sum(axis='columns')}) # p(mut)=sum over locations on dfp \n",
    "dfIncdsort_sq=dfIncd_sq.sort_values(by='Incidence',ascending=False)*100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfbar=dfNloc_sq.sort_values(by='Incidence',ascending=False)\n",
    "dfbar=100*dfbar/dfbar.sum()\n",
    "axis_font = {'size':'10'}\n",
    "barlabels=dfbar[:10].index.values\n",
    "barpos=np.arange(len(barlabels))\n",
    "bardata=dfbar[:10]['Incidence'].values\n",
    "figobj,axobj=plt.subplots(figsize=(10,10))\n",
    "axobj.bar(barpos,bardata,align='center',alpha=0.5)\n",
    "plt.xticks(barpos,barlabels,rotation=90,**axis_font)\n",
    "plt.ylabel('Percent Incidence')\n",
    "dfbar.to_excel('Figure4C_SCC_Types.xlsx')\n",
    "plt.close() # replace with plt.show() to display chart"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figure 4D,E,G,H"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [],
   "source": [
    "def genIncdsrtdat(histinplist,igenelist):\n",
    "    \"\"\"Output: Dataframe which includes Silico_Incidence_mean,  CL97.5, CL2.5\\n Input:List of sub-histologies of interest.\"\"\"\n",
    "    histsred=histinplist # use aclist for adenocarcinomas and sqlist for squamous cell carcinoma, TCClist for transitional cell carcinoma and Melanomalist for melanoma. use dfSC.columns for ALL HISTOLOGIES\n",
    "    dfCp1=dfCp.loc[igenelist,histsred] # genomic conditional probablity matrix without All genes and total histology vectors\n",
    "    # P(hist AND loc) including only histologies also incorporated in genomic data. For histologies not incorporated, the genomic data has no additional information to offer.\n",
    "    sumSCtot1=dfSC[histsred].sum().sum()\n",
    "    dfSCp1=dfSCp.loc[histsred]\n",
    "    dfSCp1=(dfSCp1/sumSCtot1)*sumSCtot # Multiply by total cases over all histologies and divide by total cases in included histologies\n",
    "    dfp1=dfCp1@dfSCp1 # sum over histologies for (p(gene|hist)*p(hist AND tissue))\n",
    "\n",
    "    dfIncd1=pd.DataFrame({'Incidence':dfp1.sum(axis='columns')}) # p(mut)=sum over locations on dfp \n",
    "    dfIncdsort1=dfIncd1.sort_values(by='Incidence',ascending=False)*100\n",
    "    \n",
    "    dfinp=dfCred.loc[igenelist,histsred].astype(int)#.loc[dfIncdsort_sub1.index,histsred].astype(int)\n",
    "    \n",
    "    # Create silico study samples\n",
    "    df_SilicoStudies1=dfinp.applymap(lambda x:np.random.poisson(x,NSilicoStudies))\n",
    "    df_Silico_cases_ALL_Total=df_SilicoStudies1+(dfC.loc['Total',histsred].astype(int)-dfinp)\n",
    "\n",
    "    # For those rare cases where there are very few sequencing cases available.\n",
    "    dfAllboolean=df_Silico_cases_ALL_Total.applymap(lambda x:x<1).applymap(lambda x:any(x))\n",
    "    icollist=df_Silico_cases_ALL_Total.columns[dfAllboolean.any()]\n",
    "    indexlist=df_Silico_cases_ALL_Total.index[dfAllboolean.transpose().any()]\n",
    "    df_Silico_cases_ALL_Total.loc[indexlist,icollist]=df_Silico_cases_ALL_Total.loc[indexlist,icollist].applymap(lambda ilist:[x if x>0 else 1 for x in ilist])\n",
    "\n",
    "    # Calculate conditional probability for silico samples\n",
    "    df_Silico_Cp1=df_SilicoStudies1.div(df_Silico_cases_ALL_Total) #divide the each case count by total respective sample sizes for each histology.\n",
    "    del df_Silico_cases_ALL_Total\n",
    "    del df_SilicoStudies1\n",
    "    \n",
    "    df_Silico_p1=df_Silico_Cp1@dfSCp1 # Calcualte (p(gene|hist)*p(hist AND tissue)) for each in-silico study. Note that SEER epidimiological data is obtained from datasets orders of magnitude larger so errors are ignored.\n",
    "    del df_Silico_Cp1\n",
    "    \n",
    "    df_Silico_Incdsort1=df_Silico_p1.sum(axis='columns')*100\n",
    "    del df_Silico_p1\n",
    "    \n",
    "    df_Incdstats1=pd.DataFrame(df_Silico_Incdsort1,columns=['Incidence_Silico'])\n",
    "    df_Incdstats1['Incidence']=dfIncdsort1.loc[df_Incdstats1.index,'Incidence']\n",
    "    df_Incdstats1['CL97.5%']=df_Incdstats1['Incidence_Silico'].apply(lambda arr:np.percentile(arr,97.5))\n",
    "    df_Incdstats1['CL2.5%']=df_Incdstats1['Incidence_Silico'].apply(lambda arr:np.percentile(arr,2.5))\n",
    "    df_Incdstats1['Mean_Silico']=df_Incdstats1['Incidence_Silico'].apply(np.mean)\n",
    "    df_Incdstats1=df_Incdstats1.drop(columns='Incidence_Silico')\n",
    "    df_Incdstats1=df_Incdstats1.sort_values(by='Incidence',ascending=False)\n",
    "    del df_Silico_Incdsort1\n",
    "    \n",
    "    return df_Incdstats1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "Fig4histdic={'D':['Adenocarcinoma',aclist],'E':['Squamous Cell Carcinoma',sqlist],'G':['Melanoma',Melanomalist],'H':['Transitional Cell Carcinoma',TCClist]}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Done:  Adenocarcinoma\n",
      "Done:  Squamous Cell Carcinoma\n",
      "Done:  Melanoma\n",
      "Done:  Transitional Cell Carcinoma\n",
      "Wall time: 10.7 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "for figlabel in ['D','E','G','H']: # Change this label to D,E,G,H to generate respective plots\n",
    "    dfIncd_hist=genIncdsrtdat(Fig4histdic[figlabel][1],CT1list)\n",
    "#     dfIncd_hist=dfIncd_hist.loc[[igene for igene in CT1list if igene in dfIncd_hist.index]]\n",
    "    print('Done: ',Fig4histdic[figlabel][0])\n",
    "    dfIncd_hist=dfIncd_hist.sort_values(by='Incidence',ascending=False)\n",
    "    dfbarplt=dfIncd_hist[:25]\n",
    "    axis_font = {'size':'10'}\n",
    "    barlabels=dfbarplt.index.values\n",
    "    barpos=np.arange(len(barlabels))\n",
    "    bardata=dfbarplt['Incidence'].values\n",
    "    barerr=[abs(dfbarplt['CL2.5%']-dfbarplt['Incidence']).values,(dfbarplt['CL97.5%']-dfbarplt['Incidence']).values]\n",
    "    figobj,axobj=plt.subplots(figsize=(10,10))\n",
    "    axobj.bar(barpos,bardata,yerr=barerr,align='center',alpha=0.5)\n",
    "    plt.xticks(barpos,barlabels,rotation=90,**axis_font)\n",
    "    plt.ylabel('Mutation Frequency in '+Fig4histdic[figlabel][0]+' (%)')\n",
    "    figname='Figure4'+figlabel+'_Top25_'+Fig4histdic[figlabel][0]\n",
    "    dfbarplt.to_excel(figname+'.xlsx')\n",
    "    ST4='SuppTable4_'+Fig4histdic[figlabel][0]\n",
    "    dfIncd_hist.to_excel(ST4+'.xlsx')\n",
    "    plt.savefig(figname+'.png',dpi=300,bbox_inches='tight')\n",
    "    plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figure 4F"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Difference between two poisson distribution is a skellam distribution. However, each gene is distributed by a rewighted combination of poisson distributions which means that we need to explicitly calcualte the distribution of the difference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wall time: 133 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "df_Incdstats_ac=pd.read_excel('SuppTable4_Adenocarcinoma.xlsx',index_col='Hugo_Symbol')\n",
    "df_Incdstats_sq=pd.read_excel('SuppTable4_Squamous Cell Carcinoma.xlsx',index_col='Hugo_Symbol')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "genelist_CT1acsq=[igene for igene in CT1list if ((igene in df_Incdstats_ac.index) and (igene in df_Incdstats_sq.index))]\n",
    "genelist_acsq=list(df_Incdstats_ac.index[df_Incdstats_ac.index.isin(df_Incdstats_sq.index)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "def nonoverlapping(igene):\n",
    "    \"\"\"returns True if adeno and squamous 95% cl are non-overlapping for input gene\"\"\"\n",
    "    acmin=df_Incdstats_ac.loc[igene,'CL2.5%']\n",
    "    sqmin=df_Incdstats_sq.loc[igene,'CL2.5%']\n",
    "    acmax=df_Incdstats_ac.loc[igene,'CL97.5%']\n",
    "    sqmax=df_Incdstats_sq.loc[igene,'CL97.5%']\n",
    "    if (acmin>sqmax) or (sqmin>acmax):\n",
    "        return True\n",
    "    else:\n",
    "        return False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(87, 570)"
      ]
     },
     "execution_count": 83,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SigGenes_AcvsSq_CT1=[igene for igene in genelist_CT1acsq if nonoverlapping(igene)]\n",
    "len(SigGenes_AcvsSq_CT1),len(genelist_CT1acsq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [],
   "source": [
    "font = {'family' : 'Arial',\n",
    "        'weight' : 'normal',\n",
    "        'size'   : 20}\n",
    "plt.rc('font', **font)\n",
    "plt.rcParams['pdf.fonttype'] = 42 # To make output Adobe Illustrator friendly for formatting the output pie charts\n",
    "\n",
    "plt.figure(figsize=(4,4))\n",
    "df_xplt=df_Incdstats_ac.copy()\n",
    "df_yplt=df_Incdstats_sq.copy()\n",
    "\n",
    "xvals=df_xplt.Incidence.loc[genelist_CT1acsq].values\n",
    "yvals=df_yplt.Incidence.loc[genelist_CT1acsq].values\n",
    "plt.scatter(x=xvals,y=yvals,marker='o',color='k',s=50)\n",
    "xyline=np.arange(0,40)\n",
    "plt.plot(xyline,xyline,color='gray')\n",
    "\n",
    "genesublist_CT1acsq=[row for row in SigGenes_AcvsSq_CT1 if df_Incdstats_ac.Incidence[row]>6 or df_Incdstats_sq.Incidence[row]>5 ]\n",
    "xvals=df_xplt.Incidence.loc[genesublist_CT1acsq].values\n",
    "yvals=df_yplt.Incidence.loc[genesublist_CT1acsq].values\n",
    "xerror=[abs(df_xplt.loc[genesublist_CT1acsq,'CL2.5%']-df_xplt.loc[genesublist_CT1acsq,'Incidence']).values,(df_xplt.loc[genesublist_CT1acsq,'CL97.5%']-df_xplt.loc[genesublist_CT1acsq,'Incidence']).values]\n",
    "yerror=[abs(df_yplt.loc[genesublist_CT1acsq,'CL2.5%']-df_yplt.loc[genesublist_CT1acsq,'Incidence']).values,(df_yplt.loc[genesublist_CT1acsq,'CL97.5%']-df_yplt.loc[genesublist_CT1acsq,'Incidence']).values]\n",
    "plt.errorbar(x=xvals,y=yvals,xerr=xerror,yerr=yerror,color='k',fmt='o',ms=8,elinewidth=2)\n",
    "\n",
    "labelgenes=[igene for igene in SigGenes_AcvsSq_CT1 if abs(df_xplt.Incidence.loc[igene]-df_yplt.Incidence.loc[igene])>7]\n",
    "dicloctexty={'FAT1':0.4, 'KMT2D':0.5, 'LRP1B':0.3, 'KRAS':0.3, 'CDKN2A':-0.4, 'TP53':0.4,'APC':0.3,'NOTCH1':0.2,'NFE2L2':0.2}\n",
    "dicloctextx={'FAT1':-4.5, 'KMT2D':0.2, 'LRP1B':1., 'KRAS':1, 'CDKN2A':0.4, 'TP53':0.7,'APC':0.2,'NOTCH1':0.2,'NFE2L2':0.2}\n",
    "for id1 in labelgenes:\n",
    "    if id1 not in ['NFE2L2','CDKN2A','FAT1']:\n",
    "        plt.text(x=df_xplt.loc[id1,'Incidence']+dicloctextx[id1],y=df_yplt.loc[id1,'Incidence']+dicloctexty[id1],s=id1,fontsize=16,color='k')\n",
    "# plt.text(x=36,y=33,s='x = y',fontsize=20,color='gray',rotation=30)\n",
    "plt.xlim([0,40])\n",
    "plt.ylim([0,70])\n",
    "plt.gca().spines['right'].set_linewidth(3)\n",
    "plt.gca().spines['top'].set_linewidth(3)\n",
    "plt.gca().spines['bottom'].set_linewidth(3)\n",
    "plt.gca().spines['left'].set_linewidth(3)\n",
    "plt.gca().xaxis.set_tick_params(width=3,size=10,labelsize=20)\n",
    "plt.gca().yaxis.set_tick_params(width=3,size=10,labelsize=20)\n",
    "\n",
    "plt.xlabel('Adenocarcinoma (%)',size=20)\n",
    "plt.ylabel('Squamous Cell Carcinoma (%)',size=20)\n",
    "fname='Figure4F_SCCvsAdeno_scatter_CT1'\n",
    "plt.savefig(fname+'.pdf',bbox_inches='tight')\n",
    "plt.close() # replace with plt.show() to display chart"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[['NOTCH1', 2.9, 10.43],\n",
       " ['LRP1B', 10.46, 23.12],\n",
       " ['CDKN2A', 1.71, 12.49],\n",
       " ['KMT2D', 6.86, 16.47],\n",
       " ['FAT1', 5.35, 14.71],\n",
       " ['NFE2L2', 1.2, 9.85],\n",
       " ['TP53', 33.7, 62.94]]"
      ]
     },
     "execution_count": 85,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# How many genes have more incidence in SCC vs Adenocarcinoma?\n",
    "[[igene,round(df_xplt.Incidence.loc[igene],2),round(df_yplt.Incidence.loc[igene],2)] for igene in SigGenes_AcvsSq_CT1 if df_xplt.Incidence.loc[igene]-df_yplt.Incidence.loc[igene]<-5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_Incdstats_ac.loc[genelist_CT1acsq].to_excel(fname+'AC.xlsx')\n",
    "df_Incdstats_sq.loc[genelist_CT1acsq].to_excel(fname+'SQ.xlsx')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figure 4I"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_Incdstats_mel=pd.read_excel('SuppTable4_Melanoma.xlsx',index_col='Hugo_Symbol')\n",
    "df_Incdstats_tcc=pd.read_excel('SuppTable4_Transitional Cell Carcinoma.xlsx',index_col='Hugo_Symbol')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>KRAS</th>\n",
       "      <th>NRAS</th>\n",
       "      <th>HRAS</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>TCC</th>\n",
       "      <td>[4.252, 6.01, 2.763]</td>\n",
       "      <td>[2.041, 3.193, 1.031]</td>\n",
       "      <td>[6.122, 8.0, 4.333]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M</th>\n",
       "      <td>[1.394, 2.317, 0.589]</td>\n",
       "      <td>[24.272, 26.908, 21.523]</td>\n",
       "      <td>[1.889, 2.936, 0.954]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SCC</th>\n",
       "      <td>[1.101, 1.708, 0.564]</td>\n",
       "      <td>[0.753, 1.266, 0.316]</td>\n",
       "      <td>[3.212, 4.145, 2.353]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>AC</th>\n",
       "      <td>[13.978, 14.486, 13.44]</td>\n",
       "      <td>[1.488, 1.732, 1.251]</td>\n",
       "      <td>[0.595, 0.76, 0.438]</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                        KRAS                      NRAS                   HRAS\n",
       "TCC     [4.252, 6.01, 2.763]     [2.041, 3.193, 1.031]    [6.122, 8.0, 4.333]\n",
       "M      [1.394, 2.317, 0.589]  [24.272, 26.908, 21.523]  [1.889, 2.936, 0.954]\n",
       "SCC    [1.101, 1.708, 0.564]     [0.753, 1.266, 0.316]  [3.212, 4.145, 2.353]\n",
       "AC   [13.978, 14.486, 13.44]     [1.488, 1.732, 1.251]   [0.595, 0.76, 0.438]"
      ]
     },
     "execution_count": 88,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_RASstats=round(pd.concat([df_Incdstats_ac.loc[['KRAS','NRAS','HRAS']].rename(columns={'Incidence_AdenoCarcinomas':'Incidence'}),df_Incdstats_sq.loc[['KRAS','NRAS','HRAS']],df_Incdstats_mel.loc[['KRAS','NRAS','HRAS']],df_Incdstats_tcc.loc[['KRAS','NRAS','HRAS']]]),3)\n",
    "df_RASstats['Hist_Type']=['AC']*3+['SCC']*3+['M']*3+['TCC']*3\n",
    "df_RAS=pd.DataFrame(index=['AC','SCC','M','TCC'],columns=['KRAS','NRAS','HRAS'])\n",
    "for icol in df_RAS.columns:\n",
    "    for idx in df_RAS.index:\n",
    "        df_RAS.loc[idx,icol]=df_RASstats[df_RASstats.Hist_Type==idx].loc[icol,['Incidence','CL97.5%','CL2.5%']].values\n",
    "df_RAS.loc[['TCC','M','SCC','AC']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Supplementary Tables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Directory Supp_Tables/ not found.\n",
      "Creating new directory.\n",
      "Changing to new directory.\n"
     ]
    }
   ],
   "source": [
    "chgdir('Supp_Tables/')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(20451, 20451)"
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "setlen(df_Incdstats.index),len(df_Incdstats.index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfHistIncd=dfSp.sum(axis='columns') #p(hist) from SEER data\n",
    "dfHistIncd=dfHistIncd.sort_values(ascending=False)*100\n",
    "dfHPsorted=pd.DataFrame(columns=['CUSTOM SITES','Incidence'],index=dfHistIncd.index)\n",
    "dfHPsorted['CUSTOM SITES']=dfS0.loc['CUSTOM SITES']\n",
    "dfHPsorted['Incidence']=dfHistIncd\n",
    "dfHPsorted.to_excel('SuppTable1_ROSETTA_Abundance.xlsx')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_Incdstats.sort_values(by='Incidence',ascending=False).round(5).to_excel('SuppTableIII_AllGenes_Incidences.xlsx')"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# Supp Table III is a combination of Figure3_ CT1, RAS, Kinome and AllGenes txt files already output to the Results folder.\n",
    "# Supp Table IV is a combination of four histological incidence txt files already output to the Results folder with names starting with- SuppTable4_."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}