- if doCache and flagRDS:
- hitRDS.saveCacheDB(hitfile)
-
- writeLog(logfilename, versionString, "returned %d region counts for %s (%.2f M reads)" % (len(labelList), hitfile, totalCount / 1000000.))
+ writeLog(logfilename, versionString, "returned %d region counts (%.2f M reads)" % (len(labelList), totalCount / 1000000.))
+
+
+def getRegionReadCounts(bamfile, chr, start, end, doUniqs=True, doMulti=False, doSplices=False):
+ uniques = 0
+ multis = 0.0
+ uniqueSplice = 0
+ multiSplice = 0.0
+ for alignedread in bamfile.fetch(chr, start, end):
+ readMultiplicity = alignedread.opt('NH')
+ if doSplices and isSpliceEntry(alignedread.cigar):
+ if readMultiplicity == 1 and doUniqs:
+ uniqueSplice += 1
+ elif doMulti:
+ multiSplice += 1.0/readMultiplicity
+ elif readMultiplicity == 1 and doUniqs:
+ uniques += 1
+ elif doMulti:
+ multis += 1.0/readMultiplicity
+
+ totalReads = uniques + multis + uniqueSplice + multiSplice
+
+ return totalReads