5 # Originally written by: Wendy Lee
6 # Last modified: May 11th, 2009 by Ali Mortazavi
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
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
16 usage: python getSNPs.py samplerdsfile uniqStartMin totalRatioMin outfile [--nosplices] [--enforceChr] [--cache pages] where
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
23 from commoncode import readDataset, writeLog
25 print "%prog: version 3.5"
31 print "psyco is not running"
44 parser = optparse.OptionParser(usage=usage)
45 parser.add_option("--nosplices", action="store_false", dest="doSplices")
46 parser.add_option("--enforceChr", action="store_true", dest="forceChr")
47 parser.add_option("--cache", type="int", dest="cachePages")
48 parser.set_defaults(doSplices=True, forceChr=False, cachePages=0)
49 (options, args) = parser.parse_args(argv[1:])
56 uniqStartMin = float(args[1])
57 totalRatioMin = float(args[2])
60 if options.cachePages > 0:
65 writeSNPsToFile(hitfile, uniqStartMin, totalRatioMin, outfilename, doCache, options.cachePages, options.doSplices, options.forceChr)
68 def writeSNPsToFile(hitfile, uniqStartMin, totalRatioMin, outfilename, doCache, cachePages=0, doSplices=True, forceChr=False):
69 writeLog("snp.log", sys.argv[0], "rdsfile: %s uniqStartMin: %1.2f totalRatioMin: %1.2f" % (hitfile, uniqStartMin, totalRatioMin))
71 outfile = open(outfilename, "w")
72 header = "#Sl\tCl\tchrom\tpos\tmatch\tuniqMis\t\ttotalMis\tchange"
73 outfile.write(header + "\n")
75 snpPropertiesList = getSNPs(hitfile, uniqStartMin, totalRatioMin, doCache, cachePages, doSplices, forceChr)
76 for snpEntry in snpPropertiesList:
77 outline = "%1.2f\t%1.2f\t%s\t%d\t%d\t%d\t\t%d\t%s\n" % snpEntry
79 outfile.write(outline + "\n")
84 writeLog("snp.log", sys.argv[0], "%d candidate SNPs\n" % len(snpPropertiesList))
87 def getSNPs(hitfile, uniqStartMin, totalRatioMin, doCache, cachePages=0, doSplices=True, forceChr=False):
89 hitRDS = readDataset(hitfile, verbose=True, cache=doCache)
90 if cachePages > 20000:
91 hitRDS.setDBcache(cachePages)
93 snpPropertiesList = []
94 readLength = hitRDS.getReadSize()
95 chromList = hitRDS.getChromosomes()
97 for chrom in chromList:
98 if doNotProcessChromosome(forceChr, chrom):
101 matchDict = getMatchDict(hitRDS, chrom, doSplices)
102 print "got match dict for %s " % chrom
103 mismatchDict = getMismatchDict(hitRDS, chrom, doSplices)
104 print "got mismatch dict for %s " % chrom
105 mismatchPositions = mismatchDict.keys()
106 mismatchPositions.sort()
107 for position in mismatchPositions:
108 totalCount = mismatchDict[position]["totalCount"]
109 uniqBaseDict = mismatchDict[position]["uniqBaseDict"]
110 totalBaseDict = mismatchDict[position]["totalBaseDict"]
112 highestBaseChange = "N-N"
113 highestTotalCount = 0
114 for baseChange in uniqBaseDict:
115 if totalBaseDict[baseChange] > highestTotalCount:
116 highestBaseChange = baseChange
117 highestCount = uniqBaseDict[baseChange]
118 highestTotalCount = totalBaseDict[baseChange]
122 if highestCount >= uniqStartMin:
123 for matchpos in xrange(position - readLength + 1, position + 1):
125 matchCount += len([mstop for mstop in matchDict[matchpos] if position <= mstop])
129 matchCount -= totalCount
133 Sl = highestCount/float(readLength)
134 Cl = highestTotalCount/float(highestTotalCount + matchCount)
135 if Cl >= totalRatioMin:
136 snpProperties = (Sl, Cl, chrom, position, matchCount, highestCount, highestTotalCount, highestBaseChange)
137 snpPropertiesList.append(snpProperties)
139 return snpPropertiesList
142 def doNotProcessChromosome(forceChr, chromosome):
144 if chromosome[:3] != "chr":
150 def getMatchDict(rds, chrom, withSplices=True):
156 readDict = rds.getReadsDict(fullChrom=True, bothEnds=True, noSense=True, chrom=chrom)
160 for (start, stop) in readDict[chrom]:
161 if finalDict.has_key(start):
162 finalDict[start].append(stop)
164 finalDict[start] = [stop]
168 spliceDict = rds.getSplicesDict(noSense=True, fullChrom=True, chrom=chrom, splitRead=True)
170 spliceDict[chrom] = []
172 for (start, stop) in spliceDict[chrom]:
173 if finalDict.has_key(start):
174 finalDict[start].append(stop)
176 finalDict[start] = [stop]
181 def getMismatchDict(rds, chrom, withSplices=True):
183 spliceDict = rds.getMismatches(mischrom=chrom, useSplices=withSplices)
184 for (start, change_at, change_base, change_from) in spliceDict[chrom]:
185 change = "%s-%s" % (change_base, change_from)
188 back = "%s:%s" % (str(start), change)
189 uniqBaseDict = {change: 1}
190 totalBaseDict = {change: 1}
191 if mismatchDict.has_key(change_at):
192 (uniqueReadCount, totalCount, back, uniqBaseDict, totalBaseDict) = mismatchDict[change_at]
193 pos = "%s:%s" % (str(start), change)
195 if totalBaseDict.has_key(change):
196 totalBaseDict[change] += 1
200 if uniqBaseDict.has_key(change):
201 uniqBaseDict[change] += 1 # dict contains total unique read counts
203 back = "%s,%s" % (back, pos)
205 mismatchDict[change_at] = {"uniqueReadCount": uniqueReadCount,
206 "totalCount": totalCount,
208 "uniqBaseDict": uniqBaseDict,
209 "totalBaseDict": totalBaseDict
215 if __name__ == "__main__":