test_StrVCTVRE_GRCh37.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# This notebook tests if StrVCTVRE is annotating correctly\n",
    "For a vcf entry to be annotated, it must have an END tag and a SVTYPE tag. Only exonic deletions and duplications will be annotated. Must be in GRCh38. Only annotates autosomes, X, and Y."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# may need to put each of these in it's own statement, throw an error if fails\n",
    "import sys\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pybedtools\n",
    "from cyvcf2 import VCF,Writer\n",
    "import annotationFinalForStrVCTVRE\n",
    "import liftover_hg19_to_hg38_public\n",
    "from joblib import dump, load\n",
    "import argparse\n",
    "import tempfile\n",
    "import shutil\n",
    "import os\n",
    "import filecmp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "parser = argparse.ArgumentParser(description='Test if StrVCTVRE is annotating correctly.')\n",
    "#parser.add_argument('-i','--input',help='Input file path',required=True,metavar = '/path/to/input/file',dest='pathIn')\n",
    "#parser.add_argument('-o','--output',help='Output file path',required=True,metavar = '/path/to/output/file',dest='pathOut')\n",
    "#parser.add_argument('-f','--format',help='Input file format, either vcf or bed, defaults to vcf when not provided',choices=['vcf','bed'],dest='formatIn',default='vcf')\n",
    "parser.add_argument('-p','--phyloP',help='phyloP file path, defaults to \\'data/hg38.phyloP100way.bw\\' when not provided',default='data/hg38.phyloP100way.bw',\n",
    "                    metavar = 'path/to/hg38.phyloP100way.bw',dest='phylopPath')\n",
    "#parser.add_argument('-a','--assembly',help='Genome assembly, either GRCh38 or GRCh37',choices=['GRCh37','GRCh38'])\n",
    "parser.add_argument('-l','--liftover',help='Liftover executable path, required if assembly is GRCh37',required=True,metavar='/path/to/liftover',dest='pathLiftover')\n",
    "\n",
    "args = parser.parse_args()\n",
    "args.formatIn = 'vcf'\n",
    "args.pathIn = 'data/test.GRCh37.vcf.gz'\n",
    "args.pathOut = 'data/test.GRCh37.annotated.vcf.gz'\n",
    "args.assembly = 'GRCh37'\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create temporary directory to store files created, deleted after finished running"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "td = tempfile.mkdtemp(prefix='StrVCTVRE.',suffix='.tmp')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "read VCF or BED into one large csv file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "reading VCF...\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# if VCF\n",
    "if args.formatIn == 'vcf':\n",
    "    print('\\nreading VCF...\\n')\n",
    "    toDf = []\n",
    "    for var in VCF(args.pathIn,gts012=True):\n",
    "        if var.INFO.get('END') and var.INFO.get('SVTYPE'):\n",
    "            entry = np.array([var.CHROM, var.POS, var.INFO['END'], var.INFO['SVTYPE']])\n",
    "            toDf.append(entry)\n",
    "    df = pd.DataFrame(toDf,columns=['chrom','start','end','svtype'])\n",
    "\n",
    "#if BED    \n",
    "else:\n",
    "    print('\\nreading BED...\\n')\n",
    "    toDf = []\n",
    "    df = pd.read_csv(args.pathIn,sep='\\t',names=['chrom','start','end','svtype'],header=None)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check bed file input has SVTYPE, an easy thing to forget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "if df['svtype'].isnull().all() & (args.formatIn == 'bed'):\n",
    "    sys.exit('ERROR: likely missing SVTYPE column from bed file')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print('\\nformatting VCF data...\\n')\n",
    "\n",
    "# make old index so we can annotate SVs rapidly at the end\n",
    "df['OldID'] = pd.Series(df.index.values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Confirm correct chromosomes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# check that the chroms all have chr in front\n",
    "if sum(df['chrom'].astype(str).str.startswith('chr',na=False))/df.shape[0] < 0.5:\n",
    "    df['chrom'] = 'chr' + df['chrom'].astype(str)\n",
    "\n",
    "acceptedChroms = ['chr1','chr2','chr3','chr4','chr5','chr6','chr7','chr8','chr9','chr10','chr11','chr12','chr13',\n",
    "                 'chr14','chr15','chr16','chr17','chr18','chr19','chr20','chr21','chr22','chrX','chrY']\n",
    "# keep only autosomes, X, and Y \n",
    "df = df[df['chrom'].isin(acceptedChroms)].copy()\n",
    "validChrom = df.copy()\n",
    "validChrom['validChrom'] = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Liftover data to hg38 if needed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "if args.assembly == 'GRCh37':\n",
    "    df['currentIndex'] = list(df.index.values)\n",
    "    df[['chrom','start','end','svtype','OldID','currentIndex']].to_csv(os.path.join(td,'svsForLiftover.csv'))\n",
    "    liftover_hg19_to_hg38_public.liftover(os.path.join(td,'svsForLiftover.csv'), td, True,os.path.join(td,'svsFromLiftover.csv'), args.pathLiftover)\n",
    "    df = pd.read_csv(os.path.join(td,'svsFromLiftover.csv'))\n",
    "    df['chrom'] = df['chrom_38']\n",
    "    df['start'] = df['start_38']\n",
    "    df['end'] = df['end_38']\n",
    "    df.index = list(df['currentIndex'])\n",
    "    df = df[['chrom','start','end','svtype','OldID']].copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Change formatting, keep only dels and dups"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "formatting VCF data...\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# remove all start and end values that are not numeric\n",
    "df = df[pd.to_numeric(df['start'], errors='coerce').notnull()].copy()\n",
    "df = df[pd.to_numeric(df['end'], errors='coerce').notnull()].copy()\n",
    "# convert from string to float (relevant to vcf only)\n",
    "df['start'] = df['start'].astype(float)\n",
    "df['end'] = df['end'].astype(float)\n",
    "# check all start and end values are integers\n",
    "df = df[df['start'] == df['start'] // 1]\n",
    "df = df[df['end'] == df['end'] // 1]\n",
    "df['start'] = df['start'].astype(int)\n",
    "df['end'] = df['end'].astype(int)\n",
    "validStartEnd = df.copy()\n",
    "\n",
    "# keep only SVs 50bp or longer\n",
    "df['length'] = df['end'].astype(int) - df['start'].astype(int)\n",
    "df = df[df['length'] > 49].copy()\n",
    "validLength = df.copy()\n",
    "validLength['validLength'] = True\n",
    "\n",
    "# keep only deletions and duplications\n",
    "df = df[((df['svtype'] == 'DEL') | (df['svtype'] == 'DUP'))].copy()\n",
    "validSVType = df.copy()\n",
    "validSVType['validSVType'] = True\n",
    "df['DEL'] = df['svtype'] == 'DEL'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Determine how many exons overlap each variant"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "identifying exonic deletions and duplications...\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print('\\nidentifying exonic deletions and duplications...\\n')\n",
    "\n",
    "exons = pybedtools.BedTool('data/exons_Appris_featurized_transcript_Chr1-Y_loeuf.sorted.bed')\n",
    "df[['chrom','start','end','OldID']].to_csv(os.path.join(td,'svs.bed'),sep='\\t', index=False,header=False)\n",
    "a = pybedtools.BedTool(os.path.join(td,'svs.bed'))\n",
    "b = a.intersect(exons, wa=True, wb=True).saveas(os.path.join(td,'svsExonOverlap.bed'))\n",
    "exonOverlap = pd.read_csv(os.path.join(td,'svsExonOverlap.bed'), sep='\\t', header=None, usecols=[0,1,2,3],\n",
    "                          names=['chrom', 'start', 'stop', 'OldID'])\n",
    "exonOverlap['numExons'] = exonOverlap.groupby(by='OldID').chrom.transform('size') # choice of chrom column here is arbitrary\n",
    "exonOverlap.drop_duplicates(subset='OldID', inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Drop variants that overlap no exons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "out = df.merge(exonOverlap[['numExons','OldID']],how='left',on='OldID')\n",
    "out = out[out['numExons'] > 0]\n",
    "validExon = out.set_index('OldID').copy()\n",
    "validExon['validExon'] = True\n",
    "# only annotate vars less than 3Mb\n",
    "out = out[out['length'] < 3000000]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "out[['chrom','start','end','OldID','DEL']].to_csv(os.path.join(td,'svsForAnnotation.csv'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Score each variant"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "scoring exonic deletions and duplications...\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print('\\nscoring exonic deletions and duplications...\\n')\n",
    "annotationFinalForStrVCTVRE.annotateSVs(os.path.join(td,'svsForAnnotation.csv'), os.path.join(td,'svsAnnotated.csv'), args.phylopPath, td)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "an = pd.read_csv(os.path.join(td,'svsAnnotated.csv'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# annotate SVs on each chromosome, using random forest trained on all other chroms, to avoid overfitting\n",
    "an['path'] = 0\n",
    "presentChroms = an['chrom'].value_counts().index.values\n",
    "rf = load('data/rfTrainedAllChromsPy3.joblib')\n",
    "X = an[['DEL','numExonsFinal','phyloP', 'lowestExonRank', 'allSkippable','lowestExonsInGene', 'anyConstExon','pLIMax','loeufMin', 'cdsFracStartMin', 'cdsFracEndMax', 'cdsFracMax', 'pLI_max25_ID', 'loeuf_min25_ID','topExp','topUsage','maxStrength']].copy()\n",
    "an['path'] = rf.predict_proba(X)[:,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "an.set_index('OldID', inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Annotate vcf with StrVCTVRE pathogenicity scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "writing annotated VCF...\n",
      "\n"
     ]
    }
   ],
   "source": [
    "if args.formatIn == 'vcf':\n",
    "    print('\\nwriting annotated VCF...\\n')\n",
    "\n",
    "    vcf = VCF(args.pathIn)\n",
    "    vcf.add_info_to_header({'ID':'StrVCTVRE','Description':'pathogenicity score for structural variants','Type':'Float','Number':'1'})\n",
    "\n",
    "    w = Writer(args.pathOut,vcf)\n",
    "\n",
    "    count = 0\n",
    "    for var in vcf:\n",
    "        if var.INFO.get('END') and var.INFO.get('SVTYPE'):\n",
    "            if count in an.index:\n",
    "                var.INFO['StrVCTVRE'] = str(round(an.loc[count,'path'],3))\n",
    "            elif count in validExon.index:\n",
    "                var.INFO['StrVCTVRE'] = '1.0'\n",
    "            elif count in validSVType.index:\n",
    "                var.INFO['StrVCTVRE'] = 'not_exonic'\n",
    "            elif count in validLength.index:\n",
    "                var.INFO['StrVCTVRE'] = 'not_dup_or_del'    \n",
    "            elif count in validStartEnd.index:\n",
    "                var.INFO['StrVCTVRE'] = 'less_than_50bp'\n",
    "            elif count in validChrom.index:\n",
    "                if args.assembly == 'GRCh37':\n",
    "                    var.INFO['StrVCTVRE'] = 'liftOver_to_GRCh38_failed'\n",
    "                else:\n",
    "                    var.INFO['StrVCTVRE'] = 'invalid_start_or_end'\n",
    "            else:\n",
    "                var.INFO['StrVCTVRE'] = 'not_valid_chrom'\n",
    "            count += 1\n",
    "        else:\n",
    "            var.INFO['StrVCTVRE'] = 'missing_END_or_SVTYPE'\n",
    "        w.write_record(var)\n",
    "    w.close();\n",
    "    vcf.close()\n",
    "    \n",
    "else:\n",
    "    print('\\nwriting annotated BED...\\n')\n",
    "\n",
    "    #bed = pd.read_csv(args.pathIn, sep='\\t', names=['chrom','start','end','svtype'], header=None,dtype = {'chrom':str,'start':int,'end':int,'svtype':str})\n",
    "    f = open(args.pathIn)\n",
    "    outf = open(args.pathOut,'w')\n",
    "    idx=0\n",
    "    for row in [x.strip() for x in f.readlines()]:\n",
    "        if idx in an.index:\n",
    "            score = str(round(an.loc[idx,'path'],3))\n",
    "        elif idx in validExon.index:\n",
    "            score = '1.0'\n",
    "        elif idx in validSVType.index:\n",
    "            score = 'not_exonic'\n",
    "        elif idx in validLength.index:\n",
    "            score = 'not_dup_or_del'   \n",
    "        elif idx in validStartEnd.index:\n",
    "            score = 'less_than_50bp'\n",
    "        elif idx in validChrom.index:\n",
    "            if args.assembly == 'GRCh37':\n",
    "                score = 'liftOver_to_GRCh38_failed'\n",
    "            else:\n",
    "                score = 'invalid_start_or_end'\n",
    "        else:\n",
    "            score = 'not_valid_chrom'\n",
    "        outf.write(row + '\\t' + str(score) + '\\n')\n",
    "        idx += 1\n",
    "    f.close()\n",
    "    outf.close();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "delete temporary files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "shutil.rmtree(td)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SUCCESS: StrVCTVRE is running correctly\n"
     ]
    }
   ],
   "source": [
    "if filecmp.cmp(args.pathOut,'data/test.GRCh37.correctAnnotation.vcf.gz',shallow=False):\n",
    "    print(\"SUCCESS: StrVCTVRE is running correctly\\n\")\n",
    "else:\n",
    "    print(\"ERROR: StrVCTVRE is not running correctly\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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
}