Mercurial > repos > yufei-luo > s_mart

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

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author m-zytnicki
date Mon, 29 Apr 2013 03:20:15 -0400
parents 769e306b7933
children
comparison
equal deleted inserted replaced
17:b0e8584489e6 18:94ab73e8a190
37 from commons.core.LoggerFactory import LoggerFactory 37 from commons.core.LoggerFactory import LoggerFactory
38 from commons.core.utils.RepetOptionParser import RepetOptionParser 38 from commons.core.utils.RepetOptionParser import RepetOptionParser
39 39
40 LOG_DEPTH = "smart" 40 LOG_DEPTH = "smart"
41 DEFAULT_REGION = "_all_" 41 DEFAULT_REGION = "_all_"
42 MULTIPLE_STR = {1: "", 1000: " (in kpb)", 1000000: " (in Gbp)"} 42 MULTIPLE_STR = {1: "", 1000: " (in kbp)", 1000000: " (in Gbp)"}
43 43
44 class GetReadDistribution(object): 44 class GetReadDistribution(object):
45 45
46 def __init__(self, verbosity = 0): 46 def __init__(self, verbosity = 0):
47 self.xLab = "" 47 self.xLab = ""
54 self.factors = {} 54 self.factors = {}
55 self.regions = None 55 self.regions = None
56 self.tmpDatName = None 56 self.tmpDatName = None
57 self.tmpRName = None 57 self.tmpRName = None
58 self.quorum = 1 58 self.quorum = 1
59 self.strands = False
59 self.width = 800 60 self.width = 800
60 self.height = 300 61 self.height = 300
62 self.arial = False
61 63
62 def setNames(self, names): 64 def setNames(self, names):
63 self.names = names 65 self.names = names
64 66
65 def setInputFiles(self, fileNames, format): 67 def setInputFiles(self, fileNames, format):
80 82
81 def setColors(self, colors): 83 def setColors(self, colors):
82 self.colors = colors 84 self.colors = colors
83 85
84 def setFactors(self, factors): 86 def setFactors(self, factors):
85 self.factors = dict(zip(self.names, factors)) 87 if factors == None:
88 self.factors = dict([name, 1.0] for name in self.names)
89 else:
90 self.factors = dict(zip(self.names, factors))
86 91
87 def setMultiple(self, boolean): 92 def setMultiple(self, boolean):
88 self.multiple = boolean 93 self.multiple = boolean
89 94
90 def setImageSize(self, width, height): 95 def setImageSize(self, width, height):
97 self.quorum = quorum 102 self.quorum = quorum
98 103
99 def setRegionsFile(self, fileName): 104 def setRegionsFile(self, fileName):
100 if fileName != None: 105 if fileName != None:
101 self._loadRegions(fileName) 106 self._loadRegions(fileName)
107
108 def setBothStrands(self, strands):
109 self.strands = strands
110
111 def setArial(self, arial):
112 self.arial = arial
102 113
103 def _checkOptions(self): 114 def _checkOptions(self):
104 if not self.parsers: 115 if not self.parsers:
105 self.logAndRaise("ERROR: Missing input file names") 116 self.logAndRaise("ERROR: Missing input file names")
106 117
154 if name not in self.distribution[region]: 165 if name not in self.distribution[region]:
155 self.distribution[region][name] = {} 166 self.distribution[region][name] = {}
156 chromosome = transcript.getChromosome() 167 chromosome = transcript.getChromosome()
157 nbElements = float(transcript.getTagValue("nbElements")) if "nbElements" in transcript.getTagNames() else 1 168 nbElements = float(transcript.getTagValue("nbElements")) if "nbElements" in transcript.getTagNames() else 1
158 nbElements *= self.factors.get(name, 1) 169 nbElements *= self.factors.get(name, 1)
170 strand = transcript.getDirection() if self.strands else 1
159 if chromosome not in self.distribution[region][name]: 171 if chromosome not in self.distribution[region][name]:
160 self.distribution[region][name][chromosome] = {} 172 self.distribution[region][name][chromosome] = {}
173 if strand not in self.distribution[region][name][chromosome]:
174 self.distribution[region][name][chromosome][strand] = {}
161 previousBin = None 175 previousBin = None
162 for exon in transcript.getExons(): 176 for exon in transcript.getExons():
163 for pos in range(exon.getStart(), exon.getEnd()+1): 177 for pos in range(exon.getStart(), exon.getEnd()+1):
164 bin = pos / self.binSize 178 bin = pos / self.binSize
165 if bin != previousBin: 179 if bin != previousBin:
166 self.distribution[region][name][chromosome][bin] = self.distribution[region][name][chromosome].get(bin, 0) + nbElements 180 self.distribution[region][name][chromosome][strand][bin] = self.distribution[region][name][chromosome][strand].get(bin, 0) + nbElements
167 previousBin = bin 181 previousBin = bin
168 progress.inc() 182 progress.inc()
169 progress.done() 183 progress.done()
170 184
171 def _checkQuorum(self, region): 185 def _checkQuorum(self, region):
172 if self.quorum == None: 186 if self.quorum == None:
173 return True 187 return True
174 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 188 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]])
175 189
176 def _writeData(self, region): 190 def _writeData(self, region):
177 self.tmpDatName = "tmpFile%d.dat" % (self.number) 191 self.tmpDatName = "tmpFile%d.dat" % (self.number)
178 handle = open(self.tmpDatName, "w") 192 handle = open(self.tmpDatName, "w")
179 handle.write("Chr\tPos\tCount\tSample\n") 193 handle.write("Chr\tPos\tStrand\tCount\tSample\n")
180 for name in self.distribution[region]: 194 for name in self.distribution[region]:
181 for chromosome in sorted(self.distribution[region][name].keys()): 195 for chromosome in sorted(self.distribution[region][name].keys()):
182 for pos in sorted(self.distribution[region][name][chromosome].keys()): 196 for strand in sorted(self.distribution[region][name][chromosome].keys()):
183 handle.write("%s\t%d\t%d\t\"%s\"\n" % (chromosome, pos * self.binSize, self.distribution[region][name][chromosome].get(pos, 0), name)) 197 for pos in sorted(self.distribution[region][name][chromosome][strand].keys()):
198 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))
184 handle.close() 199 handle.close()
185 200
186 def _findMultiple(self, region): 201 def _findMultiple(self, region):
187 if not self.multiple: 202 if not self.multiple:
188 return 1 203 return 1
189 maxPosition = max([self.distribution[region][name][chromosome].keys() for name in self.distribution[region] for chromosome in self.distribution[region][name]]) 204 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
190 if maxPosition > 2000000: 205 if maxPosition > 2000000:
191 return 1000000 206 return 1000000
192 elif maxPosition > 2000: 207 elif maxPosition > 2000:
193 return 1000 208 return 1000
194 return 1 209 return 1
195 210
196 def _writeScript(self, region): 211 def _writeScript(self, region):
197 self.tmpRName = "tmpFile%d.R" % (self.number) 212 self.tmpRName = "tmpFile%d.R" % (self.number)
198 fileName = self.outputFileName if region == DEFAULT_REGION else "%s_%s.png" % (os.path.splitext(self.outputFileName)[0], region) 213 fileName = self.outputFileName if region == DEFAULT_REGION else "%s_%s.png" % (os.path.splitext(self.outputFileName)[0], region)
199 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])) 214 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]))
200 title = "" if region == DEFAULT_REGION else " of %s" % (region) 215 title = "" if region == DEFAULT_REGION else " + labs(title = \"Distribution of %s\") " % (region)
201 facet = "Sample ~ Chr" if region == DEFAULT_REGION else "Sample ~ ." 216 facet = "Sample ~ Chr" if region == DEFAULT_REGION else "Sample ~ ."
202 handle = open(self.tmpRName, "w") 217 handle = open(self.tmpRName, "w")
203 multiple = self._findMultiple(region) 218 multiple = self._findMultiple(region)
219 arial = ", text = element_text(family=\"Arial\", size=20)" if self.arial else ""
220 if self.arial:
221 handle.write("library(extrafont)\nloadfonts()\n")
204 handle.write("library(ggplot2)\n") 222 handle.write("library(ggplot2)\n")
205 handle.write("data <- read.table(\"%s\", header = T)\n" % (self.tmpDatName)) 223 handle.write("data <- read.table(\"%s\", header = T)\n" % (self.tmpDatName))
206 handle.write("data$Sample <- factor(data$Sample, levels=c(%s))\n" % (", ".join(["\"%s\"" % (name) for name in self.names]))) 224 handle.write("data$Sample <- factor(data$Sample, levels=c(%s))\n" % (", ".join(["\"%s\"" % (name) for name in self.names])))
207 handle.write("png(\"%s\", width = %d, height = %d)\n" % (fileName, self.width, self.height)) 225 handle.write("png(\"%s\", width = %d, height = %d)\n" % (fileName, self.width, self.height))
208 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)) 226 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))
209 handle.write("dev.off()\n") 227 handle.write("dev.off()\n")
210 228
211 def _runR(self): 229 def _runR(self):
212 rCommand = "R" 230 rCommand = os.environ["SMARTRPATH"] if "SMARTRPATH" in os.environ else "R"
213 if "SMARTRPATH" in os.environ: 231 command = "\"%s\" CMD BATCH %s" % (rCommand, self.tmpRName)
214 rCommand = os.environ["SMARTRPATH"] 232 status = subprocess.call(command, shell=True)
215 command = "\"%s\" CMD BATCH %s" % (rCommand, self.tmpRName)
216 status = subprocess.call(command, shell=True)
217 if status != 0: 233 if status != 0:
218 raise Exception("Problem with the execution of script file %s, status is: %s" % (self.tmpRName, status)) 234 raise Exception("Problem with the execution of script file %s, status is: %s" % (self.tmpRName, status))
219 235
220 def _plot(self): 236 def _plot(self):
221 progress = Progress(len(self.distribution), "Plotting data", self.verbosity) 237 progress = Progress(len(self.distribution), "Plotting data", self.verbosity)
259 parser.add_option("-l", "--xLabel", dest="xLab", action="store", default="", type="string", help="x-axis label name [format: string]") 275 parser.add_option("-l", "--xLabel", dest="xLab", action="store", default="", type="string", help="x-axis label name [format: string]")
260 parser.add_option("-L", "--yLabel", dest="yLab", action="store", default="# reads", type="string", help="y-axis label name [format: string] [default: Reads]") 276 parser.add_option("-L", "--yLabel", dest="yLab", action="store", default="# reads", type="string", help="y-axis label name [format: string] [default: Reads]")
261 parser.add_option("-c", "--colors", dest="colors", action="store", default=None, type="string", help="colors of the bars, separated by commas [format: string]") 277 parser.add_option("-c", "--colors", dest="colors", action="store", default=None, type="string", help="colors of the bars, separated by commas [format: string]")
262 parser.add_option("-a", "--factors", dest="factors", action="store", default=None, type="string", help="normalization factors, separated by commas [format: string]") 278 parser.add_option("-a", "--factors", dest="factors", action="store", default=None, type="string", help="normalization factors, separated by commas [format: string]")
263 parser.add_option("-r", "--regions", dest="regionsFileName", action="store", default=None, type="string", help="regions to plot [format: transcript file in GFF format]") 279 parser.add_option("-r", "--regions", dest="regionsFileName", action="store", default=None, type="string", help="regions to plot [format: transcript file in GFF format]")
264 parser.add_option("-m", "--multiple", dest="multiple", action="store_true", default=False, help="print position using multiples (k, G) [format: boolean] [default: False]") 280 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]")
281 parser.add_option("-m", "--multiple", dest="multiple", action="store_true", default=False, help="use human readable genomic positions (k, G) [format: boolean] [default: False]")
265 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]") 282 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]")
266 parser.add_option("-z", "--width", dest="width", action="store", default=800, type="int", help="width of the image [format: int] [default: 800]") 283 parser.add_option("-z", "--width", dest="width", action="store", default=800, type="int", help="width of the image [format: int] [default: 800]")
267 parser.add_option("-Z", "--height", dest="height", action="store", default=300, type="int", help="height of the image [format: int] [default: 300]") 284 parser.add_option("-Z", "--height", dest="height", action="store", default=300, type="int", help="height of the image [format: int] [default: 300]")
285 parser.add_option("-A", "--arial", dest="arial", action="store_true", default=False, help="use Arial font [format: boolean] [default: false]")
268 parser.add_option("-v", "--verbosity", dest="verbosity", action="store", default=1, type="int", help="trace level [format: int]") 286 parser.add_option("-v", "--verbosity", dest="verbosity", action="store", default=1, type="int", help="trace level [format: int]")
269 options = parser.parse_args()[0] 287 options = parser.parse_args()[0]
270 iGetReadDistribution = GetReadDistribution(options.verbosity) 288 iGetReadDistribution = GetReadDistribution(options.verbosity)
271 iGetReadDistribution.setNames(options.names.split(",")) 289 iGetReadDistribution.setNames(options.names.split(","))
272 iGetReadDistribution.setInputFiles(options.inputFileNames.split(","), options.format) 290 iGetReadDistribution.setInputFiles(options.inputFileNames.split(","), options.format)
277 iGetReadDistribution.setFactors(None if options.factors == None else map(float, options.factors.split(","))) 295 iGetReadDistribution.setFactors(None if options.factors == None else map(float, options.factors.split(",")))
278 iGetReadDistribution.setRegionsFile(options.regionsFileName) 296 iGetReadDistribution.setRegionsFile(options.regionsFileName)
279 iGetReadDistribution.setMultiple(options.multiple) 297 iGetReadDistribution.setMultiple(options.multiple)
280 iGetReadDistribution.setQuorum(options.quorum) 298 iGetReadDistribution.setQuorum(options.quorum)
281 iGetReadDistribution.setImageSize(options.width, options.height) 299 iGetReadDistribution.setImageSize(options.width, options.height)
300 iGetReadDistribution.setBothStrands(options.strands)
301 iGetReadDistribution.setArial(options.arial)
282 iGetReadDistribution.run() 302 iGetReadDistribution.run()
283 303