annotationFinalForStrVCTVRE.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import pybedtools\n",
    "import pyBigWig\n",
    "import numpy as np\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Develop a function that takes the top N intervals from an SV, and averages those, instead of just doing max.\n",
    "\n",
    "def topUsage(df,n):\n",
    "    df['topUsage'] = df.sort_values('avgUsage',ascending=False).iloc[0:n,:]['avgUsage'].mean()\n",
    "    df['topExp'] = df.sort_values('avgExp',ascending=False).iloc[0:n,:]['avgExp'].mean()\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Define two functions that determine how far into the amino acids the SV begins\n",
    "\n",
    "# assumes exon rank ordered.\n",
    "def cdsRank(df):\n",
    "    l = df.shape[0]\n",
    "    dff = df.iloc[0,:] # get the first entry\n",
    "    # first separate by strand\n",
    "    if dff['strand'] == 1:\n",
    "        if dff['start'] - dff['cStart'] < 0: # if sv begins before the first exon\n",
    "            df['early'] = dff['cdsCount']\n",
    "            return df\n",
    "        else:\n",
    "            df['early'] = dff['start'] - dff['cStart'] + dff['cdsCount']\n",
    "            return df\n",
    "    # reverse strand\n",
    "    if dff['strand'] == -1:\n",
    "        if dff['cStop'] - dff['stop'] < 0: # if sv begins before the first exon\n",
    "            df['early'] = dff['cdsCount']\n",
    "            return df\n",
    "        else:\n",
    "            df['early'] = dff['cStop'] - dff['stop'] + dff['cdsCount']\n",
    "            return df\n",
    "        \n",
    "# assumes reverse exon rank ordered.\n",
    "def cdsEnd(df):\n",
    "    l = df.shape[0]\n",
    "    dff = df.iloc[0,:]\n",
    "    # first separate by strand\n",
    "    if dff['strand'] == 1:\n",
    "        if dff['stop'] > dff['cStop']: # if sv begins before the first exon\n",
    "            df['late'] = dff['cdsCount'] + dff['size']\n",
    "            return df\n",
    "        else:\n",
    "            df['late'] = dff['stop'] - dff['cStart'] + dff['cdsCount']\n",
    "            return df\n",
    "    # reverse strand\n",
    "    if dff['strand'] == -1:\n",
    "        if dff['start'] < dff['cStart']: # if sv ends after the last overlapped exon\n",
    "            df['late'] = dff['cdsCount'] + dff['size']\n",
    "            return df\n",
    "        else:\n",
    "            df['late'] = dff['cStop'] - dff['start'] + dff['cdsCount']\n",
    "            return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for def:\n",
    "def annotateSVs(inpath, outpath, phylopPath, tempdir):\n",
    "\n",
    "    # read csv file into dataframe\n",
    "    \n",
    "    df = pd.read_csv(inpath)\n",
    "    \n",
    "    # Do all exon-level and gene-level features\n",
    "\n",
    "    exons = pybedtools.BedTool('data/exons_Appris_featurized_transcript_Chr1-Y_loeuf.sorted.bed')\n",
    "    df['ID'] = 'sv' + pd.Series(df.index.values).apply(str)\n",
    "    df[['chrom','start','end','ID']].to_csv(os.path.join(tempdir,'df.bed'),sep='\\t', index=False,header=False)\n",
    "    a = pybedtools.BedTool(os.path.join(tempdir,'df.bed'))\n",
    "    b = a.intersect(exons, wa=True, wb=True).saveas(os.path.join(tempdir,'dfExonOverlap.bed'))\n",
    "    del a\n",
    "    del exons\n",
    "    del b\n",
    "    exonOverlap = pd.read_csv(os.path.join(tempdir,'dfExonOverlap.bed'), sep='\\t', header=None, \n",
    "                              names=['chrom', 'start', 'stop', 'ID', 'eChrom', 'eStart', 'eStop', 'gene', 'exonRank', 'skippable', 'exonsInGene', 'const','pLI','loeuf'])\n",
    "\n",
    "    exonOverlap['numExonsFinal'] = exonOverlap.groupby(by='ID').eStart.transform('size')\n",
    "    exonOverlap['allSkippable'] = exonOverlap.groupby(by='ID').skippable.transform(lambda x: all(x))\n",
    "    exonOverlap['lowestExonRank'] = exonOverlap.groupby(by='ID').exonRank.transform('min')\n",
    "    exonOverlap['lowestExonsInGene'] = exonOverlap.groupby(by='ID').exonsInGene.transform('min')\n",
    "    exonOverlap['anyConstExon'] = exonOverlap.groupby(by='ID').const.transform('max')\n",
    "    exonOverlap['pLIMax'] = exonOverlap.groupby(by='ID').pLI.transform('max')\n",
    "    exonOverlap['loeufMin'] = exonOverlap.groupby(by='ID').loeuf.transform('min')\n",
    "\n",
    "    exonOverlap.drop_duplicates(subset='ID', inplace=True)\n",
    "    df = df.merge(exonOverlap[['ID', 'numExonsFinal', 'allSkippable', 'lowestExonRank', 'lowestExonsInGene', 'anyConstExon','pLIMax','loeufMin']], how='left', on='ID')\n",
    "    del exonOverlap\n",
    "    # numExons = the total number of exons that an SV overlaps, across all genes\n",
    "\n",
    "    # allSkippable = 1 if all exons overlapped start and end in same phase, 0 otherwise\n",
    "\n",
    "    # lowestExonRank = the minimum rank of all exons overlapped\n",
    "\n",
    "    # lowestExonsInGene = the number of exons in the gene overlapped, minimum if multiple genes\n",
    "\n",
    "    # anyConstExon = 1 if any exon overlapped is constitutive, 0 otherwise\n",
    "\n",
    "    # pLIMax = the maximum pLI of all overlapped genes (high pLI is more intolerant)\n",
    "\n",
    "    # loeufMin = the minimum LOEUF of all overlapped genes (low loeuf is more intolerant)\n",
    "\n",
    "    # Calculate conservation feature using phyloP as average of top 400 most conserved position \n",
    "\n",
    "    size = 400\n",
    "\n",
    "    with open(\"data/hg38chromsizes.tsv\") as f:\n",
    "        chrms = dict((k, v) for k,v in (line.split() for line in f))\n",
    "\n",
    "    consBW = pyBigWig.open(phylopPath)\n",
    "    # get phyloP value for each position in the SV\n",
    "    x = []\n",
    "    for i in range(df.shape[0]):\n",
    "        if int(df.loc[i,'end']) - int(df.loc[i,'start']) > 1000000:\n",
    "            x.append(np.array([15.0]))\n",
    "        else:\n",
    "            try:\n",
    "                x.append(np.nan_to_num(np.array(consBW.values(df.loc[i,'chrom'], int(df.loc[i,'start']), int(df.loc[i,'end'])))))\n",
    "            except:\n",
    "                print(df.loc[i,'start'])\n",
    "                x.append(np.array([0.5]))\n",
    "    del consBW\n",
    "    x = np.asarray(x)\n",
    "    # get the mean of the top 100 most conserved positions\n",
    "    cons = [np.mean(y[np.argsort(y)[-size:]]) for y in x]\n",
    "    del x\n",
    "    df['phyloP'] = pd.Series(cons)\n",
    "    del cons\n",
    "    # Add TAD features\n",
    "\n",
    "    tads = pybedtools.BedTool('data/rep12tadsMergedhg38.bed')\n",
    "    df[['chrom','start','end','ID']].to_csv(os.path.join(tempdir,'df.bed'),sep='\\t', index=False,header=False)\n",
    "    a = pybedtools.BedTool(os.path.join(tempdir,'df.bed'))\n",
    "    b = a.intersect(tads, wa=True, wb=True).saveas(os.path.join(tempdir,'dfTadOverlap.bed'))\n",
    "    tadOverlap = pd.read_csv(os.path.join(tempdir,'dfTadOverlap.bed'), sep='\\t', header=None, \n",
    "                              names=['chrom', 'start', 'stop', 'ID', 'tChrom', 'tStart', 'tStop', 'strength'])\n",
    "\n",
    "    tadOverlap['maxStrength'] = tadOverlap.groupby(by='ID').strength.transform('max')\n",
    "    tadOverlap.drop_duplicates(subset='ID', inplace=True)\n",
    "    df = df.merge(tadOverlap[['ID', 'maxStrength']], how='left', on='ID')\n",
    "    df['maxStrength'].fillna(value=0, inplace=True)\n",
    "    del tads\n",
    "    del tadOverlap\n",
    "    ## Add amino acid features\n",
    "\n",
    "    cds = pybedtools.BedTool('data/exons_CDS_Chr1-Y.sorted.bed')\n",
    "    df[['chrom','start','end','ID']].to_csv(os.path.join(tempdir,'df.bed'),sep='\\t', index=False,header=False)\n",
    "    a = pybedtools.BedTool(os.path.join(tempdir,'df.bed'))\n",
    "    b = a.intersect(cds, wa=True, wb=True).saveas(os.path.join(tempdir,'dfCDSOverlap.bed'))\n",
    "    cdsOverlap = pd.read_csv(os.path.join(tempdir,'dfCDSOverlap.bed'), sep='\\t', header=None, \n",
    "                              names=['chrom', 'start', 'stop', 'ID', 'cChrom', 'cStart', 'cStop', 'CDSLength', 'size', 'exonRank', 'strand','gene', 'cdsCount', 'pLI','loeuf'])\n",
    "    del cds\n",
    "    # use if statement to address possible scenario in which all given variants are in UTR and don't overlap a CDS\n",
    "    if cdsOverlap.shape[0] != 0:\n",
    "        # apply above functions to the SVs that were previously intersected with coding exons\n",
    "        out = cdsOverlap.sort_values('exonRank').groupby(['ID','gene']).apply(cdsRank)\n",
    "        out = out.sort_values('exonRank', ascending=False).groupby(['ID','gene']).apply(cdsEnd)\n",
    "\n",
    "        # get shape, but don't drop duplicates (not in place)\n",
    "        out.drop_duplicates(subset='ID').shape[0]\n",
    "\n",
    "        # Featurize above information into features normalized by cds length\n",
    "\n",
    "        out['cdsFracStart'] = out['early']/out['CDSLength']\n",
    "        out['cdsFracEnd'] = out['late']/out['CDSLength']\n",
    "        out['cdsFrac'] = (out['late'] - out['early'])/out['CDSLength'] \n",
    "\n",
    "        # This is an experimental feature, which gives the max pLI and loeuf of the genes which are signficantly disrupted by the SV.\n",
    "\n",
    "        out['pLI_max25'] = out[(out['cdsFracStart'] == 0) | (out['cdsFrac'] > 0.25)].groupby('ID')['pLI'].transform('max')\n",
    "        #out['pLI_max25'].fillna(value=0, inplace=True)\n",
    "        out['loeuf_min25']= out[(out['cdsFracStart'] == 0) | (out['cdsFrac'] > 0.25)].groupby('ID')['loeuf'].transform('min')\n",
    "        #out['loeuf_min25'].fillna(value=0, inplace=True)\n",
    "\n",
    "        # but we now need to fill in all the cells with the max loeuf_max25 in their ID\n",
    "        out['pLI_max25_ID'] = out.groupby('ID')['pLI_max25'].transform('max')\n",
    "        out['loeuf_min25_ID'] = out.groupby('ID')['loeuf_min25'].transform('max')\n",
    "\n",
    "        out['cdsFracMax'] = out.groupby('ID')['cdsFrac'].transform('max')\n",
    "        out['cdsFracStartMin'] = out.groupby('ID')['cdsFracStart'].transform('min')\n",
    "        out['cdsFracEndMax'] = out.groupby('ID')['cdsFracEnd'].transform('max')\n",
    "\n",
    "        out.drop_duplicates(subset='ID', inplace=True)\n",
    "        \n",
    "        final = df.merge(out[['ID', 'cdsFracStartMin', 'cdsFracEndMax', 'cdsFracMax', 'pLI_max25_ID', 'loeuf_min25_ID']], how='left')\n",
    "        del out\n",
    "    else:\n",
    "        final = df.copy()\n",
    "        del df\n",
    "        final['cdsFracStartMin'] = float('NaN')\n",
    "        final['cdsFracEndMax'] = float('NaN')\n",
    "        final['cdsFracMax'] = float('NaN')\n",
    "        final['pLI_max25_ID'] = float('NaN')\n",
    "        final['loeuf_min25_ID'] = float('NaN')\n",
    "\n",
    "    del cdsOverlap\n",
    "    final['cdsFracStartMin'].fillna(value=2, inplace=True)\n",
    "    final['cdsFracEndMax'].fillna(value=-1, inplace=True)\n",
    "    final['cdsFracMax'].fillna(value=-1, inplace=True)\n",
    "    final['pLI_max25_ID'].fillna(value=-1, inplace=True)\n",
    "    final['loeuf_min25_ID'].fillna(value=3, inplace=True)\n",
    "\n",
    "    # Add exon inclusion features\n",
    "\n",
    "    usage = pybedtools.BedTool('data/summary_exon_usage_hg38.sorted.bed')\n",
    "    final[['chrom','start','end','ID']].to_csv(os.path.join(tempdir,'df.bed'),sep='\\t', index=False,header=False)\n",
    "    a = pybedtools.BedTool(os.path.join(tempdir,'df.bed'))\n",
    "    b = a.intersect(usage, wa=True, wb=True).saveas(os.path.join(tempdir,'dfUsageOverlap.bed'))\n",
    "    usageOverlap = pd.read_csv(os.path.join(tempdir,'dfUsageOverlap.bed'), sep='\\t', header=None, \n",
    "                              names=['chrom', 'start', 'stop', 'ID', 'uChrom', 'uStart', 'uStop', 'avgUsage', 'avgExp'])\n",
    "    out = usageOverlap.groupby('ID').apply(topUsage,n=size)\n",
    "    out.drop_duplicates(subset='ID', inplace=True)\n",
    "    del a\n",
    "    del b\n",
    "    \n",
    "    if out.shape[0] == 0: \n",
    "        # try looking for features that are nearby, using bedtools window, and window size of 1000\n",
    "        b = a.window(usage, w=1000).saveas(os.path.join(tempdir,'dfUsageOverlap.bed'))\n",
    "        usageOverlap = pd.read_csv(os.path.join(tempdir,'dfUsageOverlap.bed'), sep='\\t', header=None, \n",
    "                              names=['chrom', 'start', 'stop', 'ID', 'uChrom', 'uStart', 'uStop', 'avgUsage', 'avgExp'])\n",
    "        pd.read_csv('data/summary_exon_usage_hg38.sorted.bed')\n",
    "        out = usageOverlap.groupby('ID').apply(topUsage,n=size)\n",
    "        out.drop_duplicates(subset='ID', inplace=True)\n",
    "        \n",
    "        if out.shape[0] == 0:\n",
    "            # if still zero, assume too far from exon to be found, give the median values of all exons \n",
    "            exp = pd.read_csv('data/summary_exon_usage_hg38.sorted.bed', names=['uChrom', 'uStart', 'uStop','Usage','Exp'],sep='\\t')\n",
    "            usageMed = exp.Usage.median()\n",
    "            expMed = exp.Exp.median()\n",
    "            out = final[['ID']].copy()\n",
    "            out['topUsage'] = usageMed\n",
    "            out['topExp'] = expMed\n",
    "            \n",
    "    final2 = final.merge(out[['ID', 'topUsage', 'topExp']], how='left', on='ID')\n",
    "    del final\n",
    "    del out\n",
    "    final2['topExp'].fillna(value=final2['topExp'].median(), inplace=True)\n",
    "    final2['topUsage'].fillna(value=final2['topUsage'].median(), inplace=True)\n",
    "    \n",
    "    final2.to_csv(outpath,index=False)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "StrVCTVRE_py3user",
   "language": "python",
   "name": "strvctvre_py3user"
  },
  "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.7.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}