getSNPs.py

   1 #
   2 #  getSNPs.py
   3 #  ENRAGE
   4 #
   5 # Originally written by: Wendy Lee
   6 # Last modified: May 11th, 2009 by Ali Mortazavi
   7
   8 """
   9     Get the matches and mismatches from the RDS file, and calculate the SNP thresholds uniqStartMin (Sl * readlength) and and totalRatio (Cl).
  10     For each mismatch, choose the base change that occur most frequently (ie: has the highest number
  11     of independent reads)
  12     Threshold of Sl and Cl are from user input
  13     Sl = # of independent reads supporting a base change at position S
  14     Cl = total # of all reads supporting a base change at position S / # of all # reads that pass through position S
  15
  16     usage: python getSNPs.py samplerdsfile uniqStartMin totalRatioMin outfile [--nosplices] [--enforceChr] [--cache pages] where
  17
  18     uniqStartMin = # of independent reads supporting a base change at position S
  19     totalRatioMin = total # of reads supporting a base change at position S / total # reads that pass through position S
  20 """
  21
  22 import sys
  23 import optparse
  24 from commoncode import writeLog, getConfigParser, getConfigBoolOption, getConfigIntOption
  25 import ReadDataset
  26
  27 print "getSNPs: version 3.6"
  28
  29 try:
  30     import psyco
  31     psyco.full()
  32 except:
  33     print "psyco is not running"
  34     pass
  35
  36 def usage():
  37     print __doc__
  38
  39
  40 def main(argv=None):
  41     if not argv:
  42         argv = sys.argv
  43
  44     usage = __doc__
  45
  46     parser = makeParser(usage)
  47     (options, args) = parser.parse_args(argv[1:])
  48
  49     if len(args) < 4:
  50         usage()
  51         sys.exit(2)
  52
  53     hitfile = args[0]
  54     uniqStartMin = float(args[1])
  55     totalRatioMin = float(args[2])
  56     outfilename = args[3]
  57
  58     if options.cachePages > 0:
  59         doCache = True
  60     else:
  61         doCache = False
  62
  63     writeSNPsToFile(hitfile, uniqStartMin, totalRatioMin, outfilename, doCache, options.cachePages, options.doSplices, options.forceChr)
  64
  65
  66 def makeParser(usage=""):
  67     parser = optparse.OptionParser(usage=usage)
  68     parser.add_option("--nosplices", action="store_false", dest="doSplices")
  69     parser.add_option("--enforceChr", action="store_true", dest="forceChr")
  70     parser.add_option("--cache", type="int", dest="cachePages")
  71
  72     configParser = getConfigParser()
  73     section = "getSNPs"
  74     doSplices = getConfigBoolOption(configParser, section, "doSplices", True)
  75     forceChr = getConfigBoolOption(configParser, section, "forceChr", False)
  76     cachePages = getConfigIntOption(configParser, section, "cachePages", 0)
  77
  78     parser.set_defaults(doSplices=True, forceChr=False, cachePages=0)
  79
  80     return parser
  81
  82
  83 def writeSNPsToFile(hitfile, uniqStartMin, totalRatioMin, outfilename, doCache, cachePages=0, doSplices=True, forceChr=False):
  84     writeLog("snp.log", sys.argv[0], "rdsfile: %s uniqStartMin: %1.2f totalRatioMin: %1.2f" % (hitfile, uniqStartMin, totalRatioMin))
  85
  86     outfile  = open(outfilename, "w")
  87     header = "#Sl\tCl\tchrom\tpos\tmatch\tuniqMis\t\ttotalMis\tchange"
  88     outfile.write(header + "\n")
  89
  90     snpPropertiesList = getSNPs(hitfile, uniqStartMin, totalRatioMin, doCache, cachePages, doSplices, forceChr)
  91     for snpEntry in snpPropertiesList:
  92         outline = "%1.2f\t%1.2f\t%s\t%d\t%d\t%d\t\t%d\t%s\n" % snpEntry
  93         print outline
  94         outfile.write(outline + "\n")
  95         outfile.flush()
  96
  97     outfile.close()
  98
  99     writeLog("snp.log", sys.argv[0], "%d candidate SNPs\n" % len(snpPropertiesList))
 100
 101
 102 def getSNPs(hitfile, uniqStartMin, totalRatioMin, doCache, cachePages=0, doSplices=True, forceChr=False):
 103
 104     hitRDS = ReadDataset.ReadDataset(hitfile, verbose=True, cache=doCache)
 105     if cachePages > 20000:
 106         hitRDS.setDBcache(cachePages)
 107
 108     snpPropertiesList = []
 109     readLength = hitRDS.getReadSize()
 110     chromList = hitRDS.getChromosomes()
 111
 112     for chrom in chromList:
 113         if doNotProcessChromosome(forceChr, chrom):
 114             continue
 115
 116         matchDict = getMatchDict(hitRDS, chrom, doSplices)
 117         print "got match dict for %s " % chrom
 118         mismatchDict = getMismatchDict(hitRDS, chrom, doSplices)
 119         print "got mismatch dict for %s " % chrom
 120         mismatchPositions = mismatchDict.keys()
 121         mismatchPositions.sort()
 122         for position in mismatchPositions:
 123             totalCount = mismatchDict[position]["totalCount"]
 124             uniqBaseDict = mismatchDict[position]["uniqBaseDict"]
 125             totalBaseDict = mismatchDict[position]["totalBaseDict"]
 126             highestCount = 0
 127             highestBaseChange = "N-N"
 128             highestTotalCount = 0
 129             for baseChange in uniqBaseDict:
 130                 if totalBaseDict[baseChange] > highestTotalCount:
 131                     highestBaseChange = baseChange
 132                     highestCount = uniqBaseDict[baseChange]
 133                     highestTotalCount = totalBaseDict[baseChange]
 134
 135             Cl = 0.
 136             matchCount = 0
 137             if highestCount >= uniqStartMin:
 138                 for matchpos in xrange(position - readLength + 1, position + 1):
 139                     try:
 140                         matchCount += len([mstop for mstop in matchDict[matchpos] if position <= mstop])
 141                     except:
 142                         pass
 143
 144                 matchCount -= totalCount
 145                 if matchCount < 0:
 146                     matchCount = 0
 147
 148                 Sl = highestCount/float(readLength)
 149                 Cl = highestTotalCount/float(highestTotalCount + matchCount)
 150                 if Cl >= totalRatioMin:
 151                     snpProperties = (Sl, Cl, chrom, position, matchCount, highestCount, highestTotalCount, highestBaseChange)
 152                     snpPropertiesList.append(snpProperties)
 153
 154     return snpPropertiesList
 155
 156
 157 def doNotProcessChromosome(forceChr, chromosome):
 158     if forceChr:
 159         if chromosome[:3] != "chr":
 160             return True
 161     else:
 162         return False
 163
 164
 165 def getMatchDict(rds, chrom, withSplices=True):
 166     spliceDict = {}
 167     readDict = {}
 168     finalDict = {}
 169
 170     try:
 171         readDict = rds.getReadsDict(fullChrom=True, bothEnds=True, noSense=True, chrom=chrom)
 172     except:
 173         readDict[chrom] = []
 174
 175     for read in readDict[chrom]:
 176         start = read["start"]
 177         stop = read["stop"]
 178         if finalDict.has_key(start):
 179             finalDict[start].append(stop)
 180         else:
 181             finalDict[start] = [stop]
 182
 183     if withSplices:
 184         try:
 185             spliceDict = rds.getSplicesDict(noSense=True, fullChrom=True, chrom=chrom, splitRead=True)
 186         except:
 187             spliceDict[chrom] = []
 188
 189         for read in spliceDict[chrom]:
 190             try:
 191                 start = read["startL"]
 192                 stop = read["stopL"]
 193             except KeyError:
 194                 start = read["startR"]
 195                 stop = read["stopR"]
 196
 197             if finalDict.has_key(start):
 198                 finalDict[start].append(stop)
 199             else:
 200                 finalDict[start] = [stop]
 201
 202     return finalDict
 203
 204
 205 def getMismatchDict(rds, chrom, withSplices=True):
 206     mismatchDict = {}
 207     spliceDict = rds.getMismatches(mischrom=chrom, useSplices=withSplices)
 208     for (start, change_at, change_base, change_from) in spliceDict[chrom]:
 209         change = "%s-%s" % (change_base, change_from)
 210         uniqueReadCount = 1
 211         totalCount = 1
 212         back = "%s:%s" % (str(start), change)
 213         uniqBaseDict = {change: 1}
 214         totalBaseDict = {change: 1}
 215         if mismatchDict.has_key(change_at):
 216             (uniqueReadCount, totalCount, back, uniqBaseDict, totalBaseDict) = mismatchDict[change_at]
 217             pos = "%s:%s" % (str(start), change)
 218             totalCount += 1
 219             if totalBaseDict.has_key(change):
 220                 totalBaseDict[change] += 1
 221
 222             if pos not in back:
 223                 uniqueReadCount += 1
 224                 if uniqBaseDict.has_key(change):
 225                     uniqBaseDict[change] += 1 # dict contains total unique read counts
 226
 227                 back = "%s,%s" % (back, pos)
 228
 229         mismatchDict[change_at] = {"uniqueReadCount": uniqueReadCount,
 230                                    "totalCount": totalCount,
 231                                    "back": back,
 232                                    "uniqBaseDict": uniqBaseDict,
 233                                    "totalBaseDict": totalBaseDict
 234         }
 235
 236     return mismatchDict
 237
 238
 239 if __name__ == "__main__":
 240     main(sys.argv)