runRNAPairedAnalysis.py

   1 import sys
   2 import optparse
   3 import pysam
   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
  13
  14 VERSION = "1.0"
  15
  16 def main(argv=None):
  17     if not argv:
  18         argv = sys.argv
  19
  20     print "runRNAPairedAnalysis: version %s" % VERSION
  21     usage = "usage: python runRNAPairedAnalysis.py genome rdsprefix repeatmaskdb [modelfile] [--replacemodels]"
  22
  23     parser = getParser(usage)
  24     (options, args) = parser.parse_args(argv[1:])
  25
  26     if len(args) < 3:
  27         print usage
  28         sys.exit(1)
  29
  30     genome = args[0]
  31     rdsprefix = args[1]
  32     repeatmaskdb = args[2]
  33     try:
  34         modelfile = args[3]
  35     except IndexError:
  36         modelfile = ""
  37
  38     runRNAPairedAnalysis(genome, rdsprefix, repeatmaskdb, 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 = "RunRNAPairedAnalysis"
  47     replacemodels = getConfigBoolOption(configParser, section, "replacemodels", False)
  48     parser.set_defaults(replacemodels=replacemodels)
  49
  50     return parser
  51
  52
  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
  58     """
  59
  60     bamfilename = "%s.bam" % fileprefix
  61     bamfile = pysam.Samfile(bamfilename, "rb")
  62     logfile = "rna.log"
  63
  64     # log the parameters
  65     message = "with parameters: %s %s %s" % (genome, fileprefix, repeatmaskdb)
  66     writeLog(logfile, "runRNAPairedAnalysis.py", message)
  67
  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)
  71     #
  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)
  75
  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
  79     readCounts = {}
  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)
  84
  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)
  88
  89     # count splice reads
  90     splicecountfilename = "%s.splices.count" % fileprefix
  91     geneMrnaCounts(genome, bamfile, splicecountfilename, doSplices=True, doUniqs=False, extendGenome=modelfile, replaceModels=replacemodels)
  92
  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)
  97
  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)
 102
 103     # calculate the read densities
 104     goodFileName = "%s.newregions.good" % fileprefix
 105     regionCounts(regionfilename, bamfile, goodFileName, flagRDS=True, cachePages=1, logfilename=logfile)
 106
 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)
 110
 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)
 116
 117     # weigh multi-reads
 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)
 122
 123     # calculate final exonic read density
 124     outfilename = "%s.final.rpkm" % fileprefix
 125     normalizeFinalExonic(readCounts, expandedRPKMfilename, multicountfilename, outfilename, reportFraction=True, doCache=True, writeGID=True)
 126
 127
 128 if __name__ == "__main__":
 129     main(sys.argv)