diff SMART/Java/Python/GetReadDistribution.py @ 18:94ab73e8a190

Uploaded

--- a/SMART/Java/Python/GetReadDistribution.py	Mon Apr 22 11:11:10 2013 -0400
+++ b/SMART/Java/Python/GetReadDistribution.py	Mon Apr 29 03:20:15 2013 -0400
@@ -39,7 +39,7 @@
 
 LOG_DEPTH      = "smart"
 DEFAULT_REGION = "_all_"
-MULTIPLE_STR   = {1: "", 1000: " (in kpb)", 1000000: " (in Gbp)"}
+MULTIPLE_STR   = {1: "", 1000: " (in kbp)", 1000000: " (in Gbp)"}
 
 class GetReadDistribution(object):
 
@@ -56,8 +56,10 @@
 		self.tmpDatName   = None
 		self.tmpRName     = None
 		self.quorum       = 1
+		self.strands      = False
 		self.width        = 800
 		self.height       = 300
+		self.arial        = False
 
 	def setNames(self, names):
 		self.names = names
@@ -82,7 +84,10 @@
 		self.colors = colors
 
 	def setFactors(self, factors):
-		self.factors = dict(zip(self.names, factors))
+		if factors == None:
+			self.factors = dict([name, 1.0] for name in self.names)
+		else:
+			self.factors = dict(zip(self.names, factors))
 
 	def setMultiple(self, boolean):
 		self.multiple = boolean
@@ -100,6 +105,12 @@
 		if fileName != None:
 			self._loadRegions(fileName)
 
+	def setBothStrands(self, strands):
+		self.strands = strands
+
+	def setArial(self, arial):
+		self.arial = arial
+
 	def _checkOptions(self):
 		if not self.parsers:
 			self.logAndRaise("ERROR: Missing input file names")
@@ -156,14 +167,17 @@
 				chromosome  = transcript.getChromosome()
 				nbElements  = float(transcript.getTagValue("nbElements")) if "nbElements" in transcript.getTagNames() else 1
 				nbElements *= self.factors.get(name, 1)
+				strand      = transcript.getDirection() if self.strands else 1
 				if chromosome not in self.distribution[region][name]:
 					self.distribution[region][name][chromosome] = {}
+				if strand not in self.distribution[region][name][chromosome]:
+					self.distribution[region][name][chromosome][strand] = {}
 				previousBin = None
 				for exon in transcript.getExons():
 					for pos in range(exon.getStart(), exon.getEnd()+1):
 						bin = pos / self.binSize
 						if bin != previousBin:
-							self.distribution[region][name][chromosome][bin] = self.distribution[region][name][chromosome].get(bin, 0) + nbElements
+							self.distribution[region][name][chromosome][strand][bin] = self.distribution[region][name][chromosome][strand].get(bin, 0) + nbElements
 							previousBin = bin
 			progress.inc()
 		progress.done()
@@ -171,22 +185,23 @@
 	def _checkQuorum(self, region):
 		if self.quorum == None:
 			return True
-		return max([max([max(self.distribution[region][name][chromosome].values()) for chromosome in self.distribution[region][name]]) for name in self.distribution[region]]) >= self.quorum
+		return max([max([max([max(self.distribution[region][name][chromosome][strand].values()) for strand in self.distribution[region][name][chromosome]]) for chromosome in self.distribution[region][name]]) for name in self.distribution[region]])
 
 	def _writeData(self, region):
 		self.tmpDatName = "tmpFile%d.dat" % (self.number)
 		handle          = open(self.tmpDatName, "w")
-		handle.write("Chr\tPos\tCount\tSample\n")
+		handle.write("Chr\tPos\tStrand\tCount\tSample\n")
 		for name in self.distribution[region]:
 			for chromosome in sorted(self.distribution[region][name].keys()):
-				for pos in sorted(self.distribution[region][name][chromosome].keys()):
-					handle.write("%s\t%d\t%d\t\"%s\"\n" % (chromosome, pos * self.binSize, self.distribution[region][name][chromosome].get(pos, 0), name))
+				for strand in sorted(self.distribution[region][name][chromosome].keys()):
+					for pos in sorted(self.distribution[region][name][chromosome][strand].keys()):
+						handle.write("%s\t%d\t%d\t%d\t\"%s\"\n" % (chromosome, pos * self.binSize, strand, self.distribution[region][name][chromosome][strand].get(pos, 0) * strand, name))
 		handle.close()
 
 	def _findMultiple(self, region):
 		if not self.multiple:
 			return 1
-		maxPosition = max([self.distribution[region][name][chromosome].keys() for name in self.distribution[region] for chromosome in self.distribution[region][name]])
+		maxPosition = max([max([max([max(self.distribution[region][name][chromosome][strand].keys()) for strand in self.distribution[region][name][chromosome]]) for chromosome in self.distribution[region][name]]) for name in self.distribution[region]]) * self.binSize
 		if maxPosition > 2000000:
 			return 1000000
 		elif maxPosition > 2000:
@@ -197,23 +212,24 @@
 		self.tmpRName = "tmpFile%d.R" % (self.number)
 		fileName      = self.outputFileName if region == DEFAULT_REGION else "%s_%s.png" % (os.path.splitext(self.outputFileName)[0], region)
 		colors        = "scale_fill_brewer(palette=\"Set1\") + scale_color_brewer(palette=\"Set1\")" if self.colors == None else "scale_fill_manual(values = c(%s)) + scale_color_manual(values = c(%s))" % (", ".join(["\"%s\"" % (color) for color in self.colors]), ", ".join(["\"%s\"" % (color) for color in self.colors]))
-		title         = "" if region == DEFAULT_REGION else " of %s" % (region)
+		title         = "" if region == DEFAULT_REGION else " + labs(title = \"Distribution of %s\") " % (region)
 		facet         = "Sample ~ Chr" if region == DEFAULT_REGION else "Sample ~ ."
 		handle        = open(self.tmpRName, "w")
 		multiple      = self._findMultiple(region)
+		arial         = ", text = element_text(family=\"Arial\", size=20)" if self.arial else ""
+		if self.arial:
+			handle.write("library(extrafont)\nloadfonts()\n")
 		handle.write("library(ggplot2)\n")
 		handle.write("data <- read.table(\"%s\", header = T)\n" % (self.tmpDatName))
 		handle.write("data$Sample <- factor(data$Sample, levels=c(%s))\n" % (", ".join(["\"%s\"" % (name) for name in self.names])))
 		handle.write("png(\"%s\", width = %d, height = %d)\n" % (fileName, self.width, self.height))
-		handle.write("ggplot(data, aes(x = Pos/%d, y = Count, fill = Sample, color = Sample)) + opts(title = \"Distribution%s\") + geom_bar(stat = \"identity\") + facet_grid(%s, space=\"free\") + xlab(\"%s%s\") + ylab(\"%s\") + %s + opts(legend.position = \"none\", panel.grid.major = theme_blank(), panel.grid.minor = theme_blank(), panel.background = theme_blank())\n" % (multiple, title, facet, self.xLab, MULTIPLE_STR[multiple], self.yLab, colors))
+		handle.write("ggplot(data, aes(x = Pos/%d, y = Count, fill = Sample, color = Sample)) %s + geom_bar(stat = \"identity\") + facet_grid(%s, space=\"free\") + xlab(\"%s%s\") + ylab(\"%s\") + %s + theme(legend.position = \"none\", panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank()%s)\n" % (multiple, title, facet, self.xLab, MULTIPLE_STR[multiple], self.yLab, colors, arial))
 		handle.write("dev.off()\n")
 
 	def _runR(self):
-		rCommand = "R"
-		if "SMARTRPATH" in os.environ:
-			rCommand = os.environ["SMARTRPATH"]
-		command = "\"%s\" CMD BATCH %s" % (rCommand, self.tmpRName)
-		status = subprocess.call(command, shell=True)
+		rCommand = os.environ["SMARTRPATH"] if "SMARTRPATH" in os.environ else "R"
+		command  = "\"%s\" CMD BATCH %s" % (rCommand, self.tmpRName)
+		status   = subprocess.call(command, shell=True)
 		if status != 0:
 			raise Exception("Problem with the execution of script file %s, status is: %s" % (self.tmpRName, status))
 
@@ -261,10 +277,12 @@
 	parser.add_option("-c", "--colors",    dest="colors",          action="store",      default=None,      type="string", help="colors of the bars, separated by commas  [format: string]")
 	parser.add_option("-a", "--factors",   dest="factors",         action="store",      default=None,      type="string", help="normalization factors, separated by commas  [format: string]")
 	parser.add_option("-r", "--regions",   dest="regionsFileName", action="store",      default=None,      type="string", help="regions to plot [format: transcript file in GFF format]")
-	parser.add_option("-m", "--multiple",  dest="multiple",        action="store_true", default=False,                    help="print position using multiples (k, G) [format: boolean] [default: False]")
+	parser.add_option("-2", "--strands",   dest="strands",         action="store_true", default=False,                    help="plot negative strands on the negative x-axis [format: boolean] [default: False]")
+	parser.add_option("-m", "--multiple",  dest="multiple",        action="store_true", default=False,                    help="use human readable genomic positions (k, G) [format: boolean] [default: False]")
 	parser.add_option("-q", "--quorum",    dest="quorum",          action="store",      default=1,         type="int",    help="minimum number of intervals to plot a region [format: int] [default: 1]")
 	parser.add_option("-z", "--width",     dest="width",           action="store",      default=800,       type="int",    help="width of the image [format: int] [default: 800]")
 	parser.add_option("-Z", "--height",    dest="height",          action="store",      default=300,       type="int",    help="height of the image [format: int] [default: 300]")
+	parser.add_option("-A", "--arial",     dest="arial",           action="store_true", default=False,                    help="use Arial font [format: boolean] [default: false]")
 	parser.add_option("-v", "--verbosity", dest="verbosity",       action="store",      default=1,         type="int",    help="trace level [format: int]")
 	options = parser.parse_args()[0]
 	iGetReadDistribution = GetReadDistribution(options.verbosity)
@@ -279,5 +297,7 @@
 	iGetReadDistribution.setMultiple(options.multiple)
 	iGetReadDistribution.setQuorum(options.quorum)
 	iGetReadDistribution.setImageSize(options.width, options.height)
+	iGetReadDistribution.setBothStrands(options.strands)
+	iGetReadDistribution.setArial(options.arial)
 	iGetReadDistribution.run()