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
24 from commoncode import writeLog, getConfigParser, getConfigBoolOption, getConfigIntOption
27 print "getSNPs: version 3.6"
33 print "psyco is not running"
46 parser = makeParser(usage)
47 (options, args) = parser.parse_args(argv[1:])
54 uniqStartMin = float(args[1])
55 totalRatioMin = float(args[2])
58 if options.cachePages > 0:
63 writeSNPsToFile(hitfile, uniqStartMin, totalRatioMin, outfilename, doCache, options.cachePages, options.doSplices, options.forceChr)
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")
72 configParser = getConfigParser()
74 doSplices = getConfigBoolOption(configParser, section, "doSplices", True)
75 forceChr = getConfigBoolOption(configParser, section, "forceChr", False)
76 cachePages = getConfigIntOption(configParser, section, "cachePages", 0)
78 parser.set_defaults(doSplices=True, forceChr=False, cachePages=0)
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))
86 outfile = open(outfilename, "w")
87 header = "#Sl\tCl\tchrom\tpos\tmatch\tuniqMis\t\ttotalMis\tchange"
88 outfile.write(header + "\n")
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
94 outfile.write(outline + "\n")
99 writeLog("snp.log", sys.argv[0], "%d candidate SNPs\n" % len(snpPropertiesList))
102 def getSNPs(hitfile, uniqStartMin, totalRatioMin, doCache, cachePages=0, doSplices=True, forceChr=False):
104 hitRDS = ReadDataset.ReadDataset(hitfile, verbose=True, cache=doCache)
105 if cachePages > 20000:
106 hitRDS.setDBcache(cachePages)
108 snpPropertiesList = []
109 readLength = hitRDS.getReadSize()
110 chromList = hitRDS.getChromosomes()
112 for chrom in chromList:
113 if doNotProcessChromosome(forceChr, chrom):
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"]
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]
137 if highestCount >= uniqStartMin:
138 for matchpos in xrange(position - readLength + 1, position + 1):
140 matchCount += len([mstop for mstop in matchDict[matchpos] if position <= mstop])
144 matchCount -= totalCount
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)
154 return snpPropertiesList
157 def doNotProcessChromosome(forceChr, chromosome):
159 if chromosome[:3] != "chr":
165 def getMatchDict(rds, chrom, withSplices=True):
171 readDict = rds.getReadsDict(fullChrom=True, bothEnds=True, noSense=True, chrom=chrom)
175 for read in readDict[chrom]:
176 start = read["start"]
178 if finalDict.has_key(start):
179 finalDict[start].append(stop)
181 finalDict[start] = [stop]
185 spliceDict = rds.getSplicesDict(noSense=True, fullChrom=True, chrom=chrom, splitRead=True)
187 spliceDict[chrom] = []
189 for read in spliceDict[chrom]:
191 start = read["startL"]
194 start = read["startR"]
197 if finalDict.has_key(start):
198 finalDict[start].append(stop)
200 finalDict[start] = [stop]
205 def getMismatchDict(rds, chrom, withSplices=True):
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)
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)
219 if totalBaseDict.has_key(change):
220 totalBaseDict[change] += 1
224 if uniqBaseDict.has_key(change):
225 uniqBaseDict[change] += 1 # dict contains total unique read counts
227 back = "%s,%s" % (back, pos)
229 mismatchDict[change_at] = {"uniqueReadCount": uniqueReadCount,
230 "totalCount": totalCount,
232 "uniqBaseDict": uniqBaseDict,
233 "totalBaseDict": totalBaseDict
239 if __name__ == "__main__":