4 from commoncode import getConfigParser, getConfigBoolOption, writeLog, getHeaderComment
5 from checkrmask import checkrmask
6 from geneMrnaCounts import geneMrnaCounts
7 from geneMrnaCountsWeighted import geneMrnaCountsWeighted
8 from regionCounts import regionCounts
9 from rnafarPairs import rnaFarPairs
10 from normalizeFinalExonic import normalizeFinalExonic
11 from normalizeExpandedExonic import normalizeExpandedExonic
12 from findall import findall, RegionFinder
20 print "runRNAPairedAnalysis: version %s" % VERSION
21 usage = "usage: python runRNAPairedAnalysis.py genome rdsprefix repeatmaskdb [modelfile] [--replacemodels]"
23 parser = getParser(usage)
24 (options, args) = parser.parse_args(argv[1:])
32 repeatmaskdb = args[2]
38 runRNAPairedAnalysis(genome, rdsprefix, repeatmaskdb, modelfile=modelfile, replacemodels=options.replacemodels)
42 parser = optparse.OptionParser(usage=usage)
43 parser.add_option("--replacemodels", action="store_true", dest="replacemodels")
45 configParser = getConfigParser()
46 section = "RunRNAPairedAnalysis"
47 replacemodels = getConfigBoolOption(configParser, section, "replacemodels", False)
48 parser.set_defaults(replacemodels=replacemodels)
53 def runRNAPairedAnalysis(genome, fileprefix, repeatmaskdb, modelfile="", replacemodels=False):
54 """ based on original script runRNAPairedAnalysis.sh
55 usage:runRNAPairedAnalysis.sh genome fileprefix repeatmaskdb [modelfile] [--replacemodels]
56 where fileprefix is the name of the bam file without the .bam extension
57 use "none" for the repeatmaskdb if you do not have one
60 bamfilename = "%s.bam" % fileprefix
61 bamfile = pysam.Samfile(bamfilename, "rb")
65 message = "with parameters: %s %s %s" % (genome, fileprefix, repeatmaskdb)
66 writeLog(logfile, "runRNAPairedAnalysis.py", message)
68 # count the unique reads falling on the gene models ; the nomatch files are
69 # mappable reads that fell outside of the Cistematic gene models and not the
70 # unmappable of Eland (i.e, the "NM" reads)
72 # These will need to be marked by running the BAM preprocessor with the markGID option
73 uniquecountfilename = "%s.uniqs.count" % fileprefix
74 geneMrnaCounts(genome, bamfile, uniquecountfilename, extendGenome=modelfile, replaceModels=replacemodels)
76 # calculate a first-pass RPKM to re-weigh the unique reads,
77 # using 'none' for the splice count
78 initialrpkmfilename = "%s.firstpass.rpkm" % fileprefix
80 readCounts["uniq"] = int(getHeaderComment(bamfile.header, "Unique"))
81 readCounts["splice"] = int(getHeaderComment(bamfile.header, "UniqueSplices"))
82 readCounts["multi"] = int(getHeaderComment(bamfile.header, "Multis"))
83 normalizeExpandedExonic(genome, readCounts["uniq"], uniquecountfilename, "none", initialrpkmfilename, doCache=True, extendGenome=modelfile, replaceModels=replacemodels)
85 # recount the unique reads with weights calculated during the first pass
86 uniquerecountfilename = "%s.uniqs.recount" % fileprefix
87 geneMrnaCountsWeighted(genome, bamfile, initialrpkmfilename, uniquerecountfilename, withUniqs=True, cachePages=1, extendGenome=modelfile, replaceModels=replacemodels)
90 splicecountfilename = "%s.splices.count" % fileprefix
91 geneMrnaCounts(genome, bamfile, splicecountfilename, doSplices=True, doUniqs=False, extendGenome=modelfile, replaceModels=replacemodels)
93 # find new regions outside of gene models with reads piled up
94 newregionfilename = "%s.newregions.txt" % fileprefix
95 regionFinder = RegionFinder("RNAFAR", minHits=1, withFlag="NM")
96 findall(regionFinder, bamfilename, newregionfilename, logfilename=logfile, rnaSettings=True, useMulti=False)
98 # filter out new regions that overlap repeats more than a certain fraction
99 outFileName = "%s.newregions.repstatus" % fileprefix
100 regionfilename = "%s.newregions.checked" % fileprefix
101 checkrmask(repeatmaskdb, newregionfilename, outFileName, regionfilename, startField=1, cachePages=1, logfilename=logfile)
103 # calculate the read densities
104 goodFileName = "%s.newregions.good" % fileprefix
105 regionCounts(regionfilename, bamfile, goodFileName, flagRDS=True, cachePages=1, logfilename=logfile)
107 # map all candidate regions that have paired ends overlapping with known genes
108 candidatefilename = "%s.candidates.txt" % fileprefix
109 rnaFarPairs(genome, goodFileName, bamfile, candidatefilename, doCache=True)
111 expandedRPKMfilename = "%s.expanded.rpkm" % fileprefix
112 # calculate expanded exonic read density
113 acceptedfilename = "%s.accepted.rpkm" % fileprefix
114 normalizeExpandedExonic(genome, readCounts["uniq"], uniquerecountfilename, splicecountfilename, expandedRPKMfilename, candidatefilename=candidatefilename,
115 acceptedfilename=acceptedfilename, doCache=True, extendGenome=modelfile, replaceModels=replacemodels)
118 multicountfilename = "%s.multi.count" % fileprefix
119 acceptfile = "%s.accepted.rpkm" % fileprefix
120 geneMrnaCountsWeighted(genome, bamfile, expandedRPKMfilename, multicountfilename, withMulti=True, acceptfile=acceptfile, cachePages=1,
121 extendGenome=modelfile, replaceModels=replacemodels)
123 # calculate final exonic read density
124 outfilename = "%s.final.rpkm" % fileprefix
125 normalizeFinalExonic(readCounts, expandedRPKMfilename, multicountfilename, outfilename, reportFraction=True, doCache=True, writeGID=True)
128 if __name__ == "__main__":