bamPositionalCoverage.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pysam\n",
    "import multiprocessing as mp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {
    "code_folding": [
     2,
     40,
     56,
     63,
     77,
     113
    ]
   },
   "outputs": [],
   "source": [
    "### General Tool for querying bam files for positional basepair content from a bed ###\n",
    "\n",
    "class Position:\n",
    "    ''' This is a class that is used to calculate alignment statistics over a specified position in a bam file.\n",
    "        It is suppose to be initiated by reading a bedFile where the reference and alternate alleles are alread specified.\n",
    "        There are two methods assosicated with this class that will get the positional info over the specified position,\n",
    "        and another that will format a string for output to a tsv file.'''\n",
    "    def __init__(self,chromosome,position,ref,alt):\n",
    "        self.chromosome = chromosome\n",
    "        self.position = position\n",
    "        self.ref = ref\n",
    "        self.alt = alt\n",
    "    ##################################\n",
    "    def getInfo(self,samFile,mapQ=60):\n",
    "        self.mappingQuality = mapQ\n",
    "        self.nucDict = {\"A\":0,\"C\":0,\"T\":0,\"G\":0,\"N\":0}\n",
    "        self.coverage,self.indelReads = 0.0,0.0\n",
    "        self.properPairs,self.mapQ = 0.0,0.0\n",
    "        pysamAdjustedPos = self.position - 1\n",
    "        for pileupcolumn in samFile.pileup(reference=self.chromosome,max_depth=100000000,start=pysamAdjustedPos,end=self.position):\n",
    "            if pileupcolumn.pos == (pysamAdjustedPos):\n",
    "                for pileUpRead in pileupcolumn.pileups:\n",
    "                    ### Check if the read has a \"deletion over this position\" ###\n",
    "                    if pileUpRead.is_del == 1:\n",
    "                        self.indelReads += 1.0\n",
    "                    else:\n",
    "                        ### Deal with reads that actually have coverage over this position and not a deletion ###\n",
    "                        self.nucDict[pileUpRead.alignment.query_sequence[pileUpRead.query_position]] += 1.0\n",
    "                        self.coverage += 1\n",
    "                    if pileUpRead.alignment.is_proper_pair:\n",
    "                        self.properPairs += 1\n",
    "                    ### For bwa mem, this means a read is of the highest mapping quality and uniquely mapped \n",
    "                    if pileUpRead.alignment.mapping_quality >= mapQ:\n",
    "                        self.mapQ += 1.0\n",
    "    ########################\n",
    "    def getAnnotation(self):\n",
    "        self.ann = self.chromosome+'\\t'+str(self.position)+'\\t'+self.ref+'\\t'+'\\t'+\"Coverage:\"+str(self.coverage)+'\\t'\n",
    "        self.ann = self.ann+'A'+\":\"+str(self.nucDict[\"A\"])+'\\t'+'C'+\":\"+str(self.nucDict[\"C\"])+'\\t'+'T'+\":\"+str(self.nucDict[\"T\"])+'\\t'+'G'+\":\"+str(self.nucDict[\"G\"])+'\\t'\n",
    "        self.ann = self.ann+'N'+\":\"+str(self.nucDict[\"N\"])+'\\t'\"INDELReads:\"+str(self.indelReads)+'\\t'+\"MAPQ>=\"+str(self.mappinqQuality)+\":\"+str(self.mapQ)+'\\t'+\"PROPER_PAIRED:\"+str(self.properPairs)\n",
    "\n",
    "def reportPositionStats(ls,samMan,mapQ=60):\n",
    "    '''This function takes in a list of bed file entries and will calculate positional statistics over each entry.\n",
    "       It inititally opens up a bamFile via pysam and then loops through the bedFile entries calling the Position class.object and getInfo method from the object.\n",
    "       It then returns a list of Position objects where each one has had the getInfo method called on it'''\n",
    "    samFile = pysam.AlignmentFile(samMan,\"rb\")\n",
    "    returnArray = []\n",
    "    for line in ls:\n",
    "        cols = line.strip('\\r').strip('\\n').split('\\t')\n",
    "        chrom,start,ref,pAlt = cols[0],int(cols[1]),cols[2],cols[3]\n",
    "        posMan = Position(chrom,start,ref,pAlt)\n",
    "        posMan.getInfo(samFile,mapQ=mapQ)\n",
    "        returnArray.append(posMan)\n",
    "    ###############\n",
    "    samFile.close()\n",
    "    return returnArray\n",
    "\n",
    "def readBedFileToMem(bedFile):\n",
    "    '''Pretty straight foreward. Reads through a bed file and appends each line to a list'''\n",
    "    bedFile = open(bedFile,'r')\n",
    "    bedLines = [ line for line in bedFile ]\n",
    "    bedFile.close()\n",
    "    return bedLines\n",
    "\n",
    "def divideWorkLoad(items,proc=1):\n",
    "    '''Clever little script that I paritally stole from the internet, that divides up a list into equal parts.\n",
    "       Or as equal as possible which is ideal for \"parallel processing\" because each node will get an equal work load.\n",
    "       Other methods Ive found tend to have 1 node that gets an unequal distribution of work so it just sits there.'''\n",
    "    if proc == 1:\n",
    "        return items\n",
    "    ######################################\n",
    "    baskets = [ [] for i in xrange(proc) ]\n",
    "    for i,item in enumerate(items):\n",
    "        baskets[ i % proc ].append(item)\n",
    "    ##############################\n",
    "    baskets = filter(None,baskets)\n",
    "    return baskets\n",
    "    \n",
    "def grabCoverageSNPs(bedFile,samFileLocation,proc=1):\n",
    "    '''Master function that will report position statistics for each entry in a bed file.\n",
    "       Uses multiprocessing library to distribute the work load evenly among cores.\n",
    "       Runtime is essentially linear with respect to the number of processors which is ideal.\n",
    "       Will write out the results to a file.\n",
    "       bedFile = fullFilePathToBedFile. In format of chromosome\"\\t\"position\"\\t\"position\"\\t\"referenceAllele\"\\t\"alternateAllele\n",
    "       samFileLocation = fullFilePathToBamFile\n",
    "       I refer to the input bam file as samFileLocation because pysam opens the bam file as a sam essentially.'''\n",
    "    startTime = time.time()\n",
    "    bedLineLists = divideWorkLoad(readBedFileToMem(bedFile),proc=proc)\n",
    "    output = mp.Queue()\n",
    "    pool = mp.Pool(processes=proc)\n",
    "    print \"Workload divided\"\n",
    "    print \"Calculating position stats across \"+str(proc)+\" cores...\"\n",
    "    results = [ pool.apply_async(reportPositionStats,args=(bedLines,samFileLocation)) for bedLines in bedLineLists ]\n",
    "    output = [ p.get() for p in results ]\n",
    "    print \"Positional stats calculated\"\n",
    "    print \"Ordering output...\"\n",
    "    ############\n",
    "    ### Now we want to order our positionObjects by chromosome and position (smalles position first) ###\n",
    "    ### This just makes the output easier to understand/follow instead of it just being essentially random ###\n",
    "    chrDict = {}\n",
    "    for posObjs in output:\n",
    "        for posObj in posObjs:\n",
    "            if posObj.chromosome not in chrDict:\n",
    "                chrDict.update({ posObj.chromosome:[posObj] })\n",
    "            else:\n",
    "                chrDict[posObj.chromosome].append(posObj)\n",
    "    #################\n",
    "    for c in chrDict:\n",
    "        chrDict[c] = sorted(chrDict[c],key=lambda posObj: posObj.position)\n",
    "    ######################\n",
    "    stopTime = time.time()\n",
    "    print \"RunTime:\"+'\\t'+str(stopTime-startTime)\n",
    "    return chrDict\n",
    "\n",
    "def writePositionsToFile(outFile,chrDict,chromosomeOrderList=\"NA\"):\n",
    "    '''Also pretty self explanatory. Writes the results of the positional stats objects to a file.\n",
    "       outFile = fullFilePathToOutputFile. Should include the name you want to name the outFile\n",
    "       chromosomeOrderList = [ chrName, chrName, ... ] if you want to order the output file in a chromosome speficic manner you can'''\n",
    "    if chromosomeOrderList == \"NA\":\n",
    "        chromosomeOrderList = [ c for c in chrDict ]\n",
    "    ###########################\n",
    "    outFile = open(outFile,'w')\n",
    "    for c in chromosomeOrderList:\n",
    "        for posObj in chrDict[c]:\n",
    "            posObj.getAnnotation()\n",
    "            outFile.write(posObj.ann+'\\n')\n",
    "    ###############\n",
    "    outFile.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Example as a one liner ###\n",
    "writePositionsToFile(\"pathToOutputFile\",(grabCoverageSNPS(\"pathToBedFile\",\"pathToBamFile\",proc=10),chromosomeOrderList=\"NA\"))\n",
    "\n",
    "### Or a bit easier to read ###\n",
    "chromosomePositionDict = grabCoverageSNPS(\"pathToBedFile\",\"pathToBamFile\",proc=10)\n",
    "writePositionsToFile(\"pathToOutputFile\",chromosomePositionDict,chromosomeOrderList=\"NA\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}