runStandardAnalysis.py

   1 import sys
   2 import optparse
   3 import ReadDataset
   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
  12
  13 VERSION = "1.0"
  14
  15 def main(argv=None):
  16     if not argv:
  17         argv = sys.argv
  18
  19     print "runStandardAnalysis: version %s" % VERSION
  20     usage = "usage: python runStandardAnalysis.py genome rdsprefix repeatmaskdb bpradius [modelfile] [--replacemodels]"
  21
  22     parser = getParser(usage)
  23     (options, args) = parser.parse_args(argv[1:])
  24
  25     if len(args) < 4:
  26         print usage
  27         sys.exit(1)
  28
  29     genome = args[0]
  30     rdsprefix = args[1]
  31     repeatmaskdb = args[2]
  32     bpradius = args[3]
  33     try:
  34         modelfile = args[4]
  35     except IndexError:
  36         modelfile = ""
  37
  38     runStandardAnalysis(genome, rdsprefix, repeatmaskdb, bpradius, modelfile=modelfile, replacemodels=options.replacemodels)
  39
  40
  41 def getParser(usage):
  42     parser = optparse.OptionParser(usage=usage)
  43     parser.add_option("--replacemodels", action="store_true", dest="replacemodels")
  44
  45     configParser = getConfigParser()
  46     section = "RunStandardAnalysis"
  47     replacemodels = getConfigBoolOption(configParser, section, "replacemodels", False)
  48     parser.set_defaults(replacemodels=replacemodels)
  49
  50     return parser
  51
  52
  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
  58     """
  59
  60     rdsfile = "%s.rds" % rdsprefix
  61     logfile = "rna.log"
  62
  63     # log the parameters
  64     message = "with parameters: %s %s %s %s" % (genome, rdsprefix, repeatmaskdb, bpradius)
  65     writeLog(logfile, "runStandardAnalysis.py", message)
  66
  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)
  72
  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)
  79     readCounts = {}
  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)
  84
  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)
  88
  89     # count splice reads
  90     splicecountfilename = "%s.splices.count" % rdsprefix
  91     geneMrnaCounts(genome, rdsfile, splicecountfilename, doSplices=True, doUniqs=False, extendGenome=modelfile, replaceModels=replacemodels, cachePages=1)
  92
  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)
  97
  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)
 102
 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)
 106
 107     expandedRPKMfilename = "%s.expanded.rpkm" % rdsprefix
 108     # calculate expanded exonic read density
 109     acceptedfilename = "%s.accepted.rpkm" % rdsprefix
 110     try:
 111         candidatefile = open(candidatefilename)
 112         candidateLines = candidatefile.readlines()
 113         candidatefile.close()
 114     except IOError:
 115         candidateLines = []
 116
 117     normalizeExpandedExonic(genome, readCounts["uniq"], uniquerecountfilename, splicecountfilename, expandedRPKMfilename, candidateLines=candidateLines, acceptedfilename=acceptedfilename,
 118                             doCache=True, extendGenome=modelfile, replaceModels=replacemodels)
 119
 120     # weigh multi-reads
 121     multicountfilename = "%s.multi.count" % rdsprefix
 122     geneMrnaCountsWeighted(genome, rdsfile, expandedRPKMfilename, multicountfilename, withMulti=True, acceptfile=acceptedfilename, cachePages=1, extendGenome=modelfile,
 123                            replaceModels=replacemodels)
 124
 125     # calculate final exonic read density
 126     outfilename = "%s.final.rpkm" % rdsprefix
 127     normalizeFinalExonic(readCounts, expandedRPKMfilename, multicountfilename, outfilename, reportFraction=True, doCache=True, writeGID=True)
 128
 129
 130 if __name__ == "__main__":
 131     main(sys.argv)