4 from commoncode import getConfigParser, getConfigBoolOption, writeLog
5 from checkrmask import checkrmask
6 from geneMrnaCounts import geneMrnaCounts
7 from geneMrnaCountsWeighted import geneMrnaCountsWeighted
8 from getallgenes import getallgenes
9 from normalizeFinalExonic import normalizeFinalExonic
10 from normalizeExpandedExonic import normalizeExpandedExonic
11 from findall import findall, RegionFinder
19 print "runStandardAnalysis: version %s" % VERSION
20 usage = "usage: python runStandardAnalysis.py genome rdsprefix repeatmaskdb bpradius [modelfile] [--replacemodels]"
22 parser = getParser(usage)
23 (options, args) = parser.parse_args(argv[1:])
31 repeatmaskdb = args[2]
38 runStandardAnalysis(genome, rdsprefix, repeatmaskdb, bpradius, 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 = "RunStandardAnalysis"
47 replacemodels = getConfigBoolOption(configParser, section, "replacemodels", False)
48 parser.set_defaults(replacemodels=replacemodels)
53 def runStandardAnalysis(genome, rdsprefix, repeatmaskdb, bpradius, modelfile="", replacemodels=False):
54 """ based on original script runStandardAnalysis.sh
55 usage: runStandardAnalysis.sh genome rdsprefix repeatmaskdb bpradius [modelfile] [--replacemodels]
56 where rdsprefix is the name of the rds file without the .rds extension
57 use "none" for the repeatmaskdb if you do not have one
60 rdsfile = "%s.rds" % rdsprefix
64 message = "with parameters: %s %s %s %s" % (genome, rdsprefix, repeatmaskdb, bpradius)
65 writeLog(logfile, "runStandardAnalysis.py", message)
67 # count the unique reads falling on the gene models ; the nomatch files are
68 # mappable reads that fell outside of the Cistematic gene models and not the
69 # unmappable of Eland (i.e, the "NM" reads)
70 uniquecountfilename = "%s.uniqs.count" % rdsprefix
71 geneMrnaCounts(genome, rdsfile, uniquecountfilename, extendGenome=modelfile, replaceModels=replacemodels, cachePages=1, markGID=True)
73 # calculate a first-pass RPKM to re-weigh the unique reads,
74 # using 'none' for the splice count
75 splicecountfilename = "none"
76 initialrpkmfilename = "%s.firstpass.rpkm" % rdsprefix
77 RDS = ReadDataset.ReadDataset(rdsfile, verbose=True, cache=True, reportCount=False)
78 (ucount, mcount, scount) = RDS.getCounts(multi=True, splices=True, reportCombined=False)
80 readCounts["uniq"] = ucount
81 readCounts["splice"] = mcount
82 readCounts["multi"] = scount
83 normalizeExpandedExonic(genome, readCounts["uniq"], uniquecountfilename, splicecountfilename, initialrpkmfilename, doCache=True, extendGenome=modelfile, replaceModels=replacemodels)
85 # recount the unique reads with weights calculated during the first pass
86 uniquerecountfilename = "%s.uniqs.recount" % rdsprefix
87 geneMrnaCountsWeighted(genome, rdsfile, initialrpkmfilename, uniquerecountfilename, withUniqs=True, cachePages=1, extendGenome=modelfile, replaceModels=replacemodels)
90 splicecountfilename = "%s.splices.count" % rdsprefix
91 geneMrnaCounts(genome, rdsfile, splicecountfilename, doSplices=True, doUniqs=False, extendGenome=modelfile, replaceModels=replacemodels, cachePages=1)
93 # Alternative 1: find new regions outside of gene models with reads piled up
94 newregionfilename = "%s.newregions.txt" % rdsprefix
95 regionFinder = RegionFinder("RNAFAR", minHits=1, withFlag="NM")
96 findall(regionFinder, rdsfile, newregionfilename, logfilename=logfile, rnaSettings=True, cachePages=1, useMulti=False)
98 # Alternative 1: filter out new regions that overlap repeats more than a certain fraction
99 outFileName = "%s.newregions.repstatus" % rdsprefix
100 goodFileName = "%s.newregions.good" % rdsprefix
101 checkrmask(repeatmaskdb, newregionfilename, outFileName, goodFileName, startField=1, cachePages=1, logfilename=logfile)
103 # map all candidate regions that are within a given radius of a gene in bp
104 candidatefilename = "%s.candidates.txt" % rdsprefix
105 getallgenes(genome, goodFileName, candidatefilename, maxRadius=bpradius, trackFar=True, doCache=True, extendGenome=modelfile, replaceModels=replacemodels)
107 expandedRPKMfilename = "%s.expanded.rpkm" % rdsprefix
108 # calculate expanded exonic read density
109 acceptedfilename = "%s.accepted.rpkm" % rdsprefix
111 candidatefile = open(candidatefilename)
112 candidateLines = candidatefile.readlines()
113 candidatefile.close()
117 normalizeExpandedExonic(genome, readCounts["uniq"], uniquerecountfilename, splicecountfilename, expandedRPKMfilename, candidateLines=candidateLines, acceptedfilename=acceptedfilename,
118 doCache=True, extendGenome=modelfile, replaceModels=replacemodels)
121 multicountfilename = "%s.multi.count" % rdsprefix
122 geneMrnaCountsWeighted(genome, rdsfile, expandedRPKMfilename, multicountfilename, withMulti=True, acceptfile=acceptedfilename, cachePages=1, extendGenome=modelfile,
123 replaceModels=replacemodels)
125 # calculate final exonic read density
126 outfilename = "%s.final.rpkm" % rdsprefix
127 normalizeFinalExonic(readCounts, expandedRPKMfilename, multicountfilename, outfilename, reportFraction=True, doCache=True, writeGID=True)
130 if __name__ == "__main__":