annotate SMART/Java/Python/plotCoverage.py @ 36:44d5973c188c

Uploaded
author m-zytnicki
date Tue, 30 Apr 2013 15:02:29 -0400
parents
children 169d364ddd91
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
36
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
1 #! /usr/bin/env python
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
2 #
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
3 # Copyright INRA-URGI 2009-2010
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
4 #
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
5 # This software is governed by the CeCILL license under French law and
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
6 # abiding by the rules of distribution of free software. You can use,
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
7 # modify and/ or redistribute the software under the terms of the CeCILL
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
8 # license as circulated by CEA, CNRS and INRIA at the following URL
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
9 # "http://www.cecill.info".
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
10 #
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
11 # As a counterpart to the access to the source code and rights to copy,
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
12 # modify and redistribute granted by the license, users are provided only
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
13 # with a limited warranty and the software's author, the holder of the
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
14 # economic rights, and the successive licensors have only limited
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
15 # liability.
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
16 #
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
17 # In this respect, the user's attention is drawn to the risks associated
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
18 # with loading, using, modifying and/or developing or reproducing the
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
19 # software by the user in light of its specific status of free software,
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
20 # that may mean that it is complicated to manipulate, and that also
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
21 # therefore means that it is reserved for developers and experienced
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
22 # professionals having in-depth computer knowledge. Users are therefore
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
23 # encouraged to load and test the software's suitability as regards their
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
24 # requirements in conditions enabling the security of their systems and/or
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
25 # data to be ensured and, more generally, to use and operate it in the
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
26 # same conditions as regards security.
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
27 #
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
28 # The fact that you are presently reading this means that you have had
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
29 # knowledge of the CeCILL license and that you accept its terms.
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
30 #
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
31 import os, os.path, subprocess, glob, random
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
32 from optparse import OptionParser
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
33 from SMART.Java.Python.structure.Interval import Interval
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
34 from SMART.Java.Python.structure.Transcript import Transcript
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
35 from commons.core.parsing.ParserChooser import ParserChooser
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
36 from SMART.Java.Python.misc.RPlotter import RPlotter
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
37 from SMART.Java.Python.misc.Progress import Progress
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
38 from commons.core.parsing.FastaParser import FastaParser
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
39
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
40 strands = [-1, 1]
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
41 colors = {-1: "blue", 1: "red", 0: "black"}
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
42 colorLine = "black"
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
43
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
44 def parseTargetField(field):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
45 strand = "+"
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
46 splittedFieldSpace = field.split()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
47 splittedFieldPlus = field.split("+", 4)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
48 if len(splittedFieldSpace) == 3:
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
49 id, start, end = splittedFieldSpace
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
50 elif len(splittedFieldSpace) == 4:
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
51 id, start, end, strand = splittedFieldSpace
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
52 elif len(splittedFieldPlus) == 3:
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
53 id, start, end = splittedFieldPlus
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
54 elif len(splittedFieldPlus) == 4:
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
55 id, start, end, strand = splittedFieldPlus
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
56 else:
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
57 raise Exception("Cannot parse Target field '%s'." % (field))
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
58 return (id, int(start), int(end), strand)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
59
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
60
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
61 class SimpleTranscript(object):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
62 def __init__(self, transcript1, transcript2, color = None):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
63 self.start = max(0, transcript1.getStart() - transcript2.getStart())
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
64 self.end = min(transcript2.getEnd() - transcript2.getStart(), transcript1.getEnd() - transcript2.getStart())
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
65 self.strand = transcript1.getDirection() * transcript2.getDirection()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
66 self.exons = []
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
67 for exon in transcript1.getExons():
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
68 if exon.getEnd() >= transcript2.getStart() and exon.getStart() <= transcript2.getEnd():
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
69 start = max(0, exon.getStart() - transcript2.getStart())
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
70 end = min(transcript2.getEnd() - transcript2.getStart(), exon.getEnd() - transcript2.getStart())
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
71 self.addExon(start, end, self.strand, color)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
72
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
73 def addExon(self, start, end, strand, color):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
74 exon = SimpleExon(start, end, strand, color)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
75 self.exons.append(exon)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
76
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
77 def getRScript(self, yOffset, height):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
78 rString = ""
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
79 previousEnd = None
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
80 for exon in sorted(self.exons, key=lambda exon: exon.start):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
81 if previousEnd != None:
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
82 rString += "segments(%.1f, %.1f, %.1f, %.1f, col = \"%s\")\n" % (previousEnd, yOffset + height / 4.0, exon.start, yOffset + height / 4.0, colorLine)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
83 rString += exon.getRScript(yOffset, height)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
84 previousEnd = exon.end
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
85 return rString
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
86
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
87
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
88 class SimpleExon(object):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
89 def __init__(self, start, end, strand, color = None):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
90 self.start = start
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
91 self.end = end
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
92 self.strand = strand
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
93 self.color = color
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
94
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
95 def getRScript(self, yOffset, height):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
96 color = self.color if self.color != None else colors[self.strand]
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
97 return "rect(%.1f, %.1f, %.1f, %.1f, col=\"%s\", border = \"%s\")\n" % (self.start, yOffset, self.end, yOffset + height / 2.0, color, colorLine)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
98
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
99
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
100 class Plotter(object):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
101
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
102 def __init__(self, seed, index, verbosity):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
103 self.seed = seed
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
104 self.index = index
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
105 self.verbosity = verbosity
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
106 self.maxCoverage = 0
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
107 self.maxOverlap = 0
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
108 self.log = ""
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
109 self.merge = False
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
110 self.width = 1500
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
111 self.heigth = 1000
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
112 self.xLabel = ""
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
113 self.yLabel = ""
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
114 self.title = None
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
115 self.absPath = os.getcwd()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
116 self.coverageDataFileName = "tmpFile_%d_%s.dat" % (seed, index)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
117 self.coverageScript = ""
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
118 self.overlapScript = ""
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
119 self.outputFileName = None
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
120
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
121 def setOutputFileName(self, fileName):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
122 self.outputFileName = fileName
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
123
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
124 def setTranscript(self, transcript):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
125 self.transcript = transcript
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
126 self.name = transcript.getName()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
127 self.size = transcript.getEnd() - transcript.getStart() + 1
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
128 if self.title == None:
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
129 self.title = self.name
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
130 else:
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
131 self.title += " " + self.name
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
132
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
133 def setTitle(self, title):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
134 self.title = title + " " + self.name
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
135
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
136 def setPlotSize(self, width, height):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
137 self.width = width
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
138 self.height = height
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
139
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
140 def setLabels(self, xLabel, yLabel):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
141 self.xLabel = xLabel
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
142 self.yLabel = yLabel
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
143
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
144 def setMerge(self, merge):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
145 self.merge = merge
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
146
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
147 def setCoverageData(self, coverage):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
148 outputCoveragePerStrand = dict([strand, 0] for strand in strands)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
149 outputCoverage = 0
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
150 dataFile = open(os.path.abspath(self.coverageDataFileName), "w")
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
151 for position in range(self.size+1):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
152 sumValue = 0
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
153 found = False
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
154 dataFile.write("%d\t" % (position))
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
155 for strand in strands:
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
156 value = coverage[strand].get(position, 0)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
157 sumValue += value
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
158 dataFile.write("%d\t" % (value))
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
159 if value > 0:
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
160 found = True
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
161 outputCoveragePerStrand[strand] += 1
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
162 self.maxCoverage = max(self.maxCoverage, sumValue)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
163 dataFile.write("%d\n" % (sumValue))
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
164 if found:
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
165 outputCoverage += 1
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
166 dataFile.close()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
167 self.log += "%s (%d nt):\n - both strands: %d (%.0f%%)\n - (+) strand: %d (%.0f%%)\n - (-) strand: %d (%.0f%%)\n" % (self.name, self.size, outputCoverage, float(outputCoverage) / self.size * 100, outputCoveragePerStrand[1], float(outputCoveragePerStrand[1]) / self.size * 100, outputCoveragePerStrand[-1], float(outputCoveragePerStrand[-1]) / self.size * 100)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
168 self.coverageScript += "data = scan(\"%s\", list(pos = -666, minus = -666, plus = -666, sumValue = -666), sep=\"\t\")\n" % (os.path.abspath(self.coverageDataFileName))
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
169 self.coverageScript += "lines(x = data$pos, y = data$minus, col = \"%s\")\n" % (colors[-1])
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
170 self.coverageScript += "lines(x = data$pos, y = data$plus, col = \"%s\")\n" % (colors[1])
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
171 self.coverageScript += "lines(x = data$pos, y = data$sumValue, col = \"%s\")\n" % (colors[0])
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
172
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
173 def setOverlapData(self, overlap):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
174 height = 1
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
175 self.maxOverlap = (len(overlap) + 1) * height
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
176 thisElement = SimpleTranscript(self.transcript, self.transcript, "black")
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
177 self.overlapScript += thisElement.getRScript(0, height)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
178 for cpt, transcript in enumerate(sorted(overlap, cmp=lambda c1, c2: c1.start - c2.start if c1.start != c2.start else c1.end - c2.end)):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
179 self.overlapScript += transcript.getRScript((cpt + 1) * height, height)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
180
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
181 def getFirstLine(self, suffix = None):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
182 return "png(file = \"%s_%s%s.png\", width = %d, height = %d, bg = \"white\")\n" % (self.outputFileName, self.name, "" if suffix == None or self.merge else "_%s" % (suffix), self.width, self.height)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
183
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
184 def getLastLine(self):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
185 return "dev.off()\n"
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
186
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
187 def startR(self, fileName, script):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
188 scriptFile = open(fileName, "w")
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
189 scriptFile.write(script)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
190 scriptFile.close()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
191 command = "R CMD BATCH %s" % (fileName)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
192 status = subprocess.call(command, shell=True)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
193 if status != 0:
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
194 raise Exception("Problem with the execution of script file %s, status is: %s" % (fileName, status))
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
195
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
196 def plot(self):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
197 if self.merge:
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
198 fileName = "%s_%d_%s.R" % (self.outputFileName, self.seed, self.index)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
199 plotLine = "plot(x = NA, y = NA, xlab=\"%s\", ylab=\"%s\", panel.first = grid(lwd = 1.0), xlim = c(0, %d), ylim = c(0, %d), cex.axis = 2, cex.lab = 2, cex.main=2, main = \"%s\")\n" % (self.xLabel, self.yLabel, self.size, max(self.maxCoverage, self.maxOverlap), self.title)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
200 script = self.getFirstLine() + plotLine + self.overlapScript + self.coverageScript + self.getLastLine()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
201 self.startR(fileName, script)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
202 else:
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
203 fileName = "%s_%d_%s_overlap.R" % (self.outputFileName, self.seed, self.index)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
204 plotLine = "plot(x = NA, y = NA, xlab=\"%s\", ylab=\"%s\", panel.first = grid(lwd = 1.0), xlim = c(0, %d), ylim = c(0, %d), cex.axis = 2, cex.lab = 2, cex.main=2, main = \"%s\")\n" % (self.xLabel, self.yLabel, self.size, self.maxOverlap, self.title)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
205 script = self.getFirstLine("overlap") + plotLine + self.overlapScript + self.getLastLine()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
206 self.startR(fileName, script)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
207 fileName = "%s_%d_%s_coverage.R" % (self.outputFileName, self.seed, self.index)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
208 plotLine = "plot(x = NA, y = NA, xlab=\"%s\", ylab=\"%s\", panel.first = grid(lwd = 1.0), xlim = c(0, %d), ylim = c(0, %d), cex.axis = 2, cex.lab = 2, cex.main=2, main = \"%s\")\n" % (self.xLabel, self.yLabel, self.size, self.maxCoverage, self.title)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
209 script = self.getFirstLine("coverage") + plotLine + self.coverageScript + self.getLastLine()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
210 self.startR(fileName, script)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
211
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
212
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
213 class PlotParser(object):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
214
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
215 def __init__(self, verbosity):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
216 self.verbosity = verbosity
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
217 self.parsers = [None, None]
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
218 self.sequenceParser = None
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
219 self.seed = random.randint(0, 10000)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
220 self.title = ""
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
221 self.merge = False
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
222
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
223 def __del__(self):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
224 for fileName in glob.glob("tmpFile_%d*.dat" % (self.seed)):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
225 os.remove(fileName)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
226 for fileName in glob.glob("%s*.R" % (os.path.abspath(self.outputFileName))):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
227 os.remove(fileName)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
228 for fileName in glob.glob("%s*.Rout" % (os.path.abspath(self.outputFileName))):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
229 os.remove(fileName)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
230
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
231 def addInput(self, inputNb, fileName, fileFormat):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
232 if fileName == None:
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
233 return
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
234 chooser = ParserChooser(self.verbosity)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
235 chooser.findFormat(fileFormat)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
236 self.parsers[inputNb] = chooser.getParser(fileName)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
237 if inputNb == 0:
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
238 self.parsers[1] = self.parsers[0]
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
239
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
240 def addSequence(self, fileName):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
241 if fileName == None:
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
242 return
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
243 self.sequenceParser = FastaParser(fileName, self.verbosity)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
244
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
245 def setOutput(self, fileName):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
246 self.outputFileName = fileName
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
247
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
248 def setPlotSize(self, width, height):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
249 self.width = width
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
250 self.height = height
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
251
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
252 def setLabels(self, xLabel, yLabel):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
253 self.xLabel = xLabel
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
254 self.yLabel = yLabel
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
255
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
256 def setTitle(self, title):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
257 self.title = title
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
258
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
259 def setMerge(self, merge):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
260 self.merge = merge
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
261
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
262 def initializeDataFromSequences(self):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
263 self.sizes = {}
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
264 self.coverage = {}
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
265 self.overlap = {}
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
266 for region in self.sequenceParser.getRegions():
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
267 self.sizes[region] = self.sequenceParser.getSizeOfRegion(region)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
268 self.coverage[region] = {}
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
269 self.overlap[region] = []
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
270 for strand in strands:
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
271 self.coverage[region][strand] = {}
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
272 self.coverage[region][strand][1] = 0
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
273 self.coverage[region][strand][self.sizes[region]] = 0
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
274
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
275 def initializeDataFromTranscripts(self):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
276 self.coverage = dict([i, None] for i in range(self.parsers[1].getNbTranscripts()))
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
277 self.overlap = dict([i, None] for i in range(self.parsers[1].getNbTranscripts()))
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
278 self.sizes = dict([i, 0] for i in range(self.parsers[1].getNbTranscripts()))
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
279 progress = Progress(self.parsers[1].getNbTranscripts(), "Reading regions", self.verbosity)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
280 for cpt, transcript in enumerate(self.parsers[1].getIterator()):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
281 self.coverage[cpt] = {}
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
282 self.overlap[cpt] = []
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
283 for strand in strands:
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
284 self.coverage[cpt][strand] = {}
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
285 self.coverage[cpt][strand][0] = 0
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
286 self.coverage[cpt][strand][transcript.getEnd() - transcript.getStart()] = 0
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
287 for exon in transcript.getExons():
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
288 self.sizes[cpt] += exon.getSize()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
289 progress.inc()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
290 progress.done()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
291
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
292 def initialize(self):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
293 if self.sequenceParser == None:
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
294 self.initializeDataFromTranscripts()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
295 else:
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
296 self.initializeDataFromSequences()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
297
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
298 def computeCoverage(self, transcript1, transcript2, id):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
299 strand = transcript1.getDirection() * transcript2.getDirection()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
300 for exon1 in transcript1.getExons():
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
301 for exon2 in transcript2.getExons():
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
302 if exon1.overlapWith(exon2):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
303 for position in range(max(exon1.getStart(), exon2.getStart()), min(exon1.getEnd(), exon2.getEnd()) + 1):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
304 relativePosition = position - transcript2.getStart() + 1
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
305 self.coverage[id][strand][relativePosition] = self.coverage[id][strand].get(relativePosition, 0) + 1
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
306
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
307 def computeOverlap(self, transcript1, transcript2, id):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
308 simpleTranscript = SimpleTranscript(transcript1, transcript2)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
309 self.overlap[id].append(simpleTranscript)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
310
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
311 def compute2TranscriptFiles(self):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
312 progress = Progress(self.parsers[1].getNbTranscripts(), "Comparing regions", self.verbosity)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
313 for cpt2, transcript2 in enumerate(self.parsers[1].getIterator()):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
314 for transcript1 in self.parsers[0].getIterator():
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
315 if transcript1.overlapWithExon(transcript2):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
316 self.computeCoverage(transcript1, transcript2, cpt2)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
317 self.computeOverlap(transcript1, transcript2, cpt2)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
318 progress.inc()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
319 progress.done()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
320
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
321 def extractReferenceQueryMapping(self, mapping):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
322 queryTranscript = mapping.getTranscript()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
323 referenceTranscript = Transcript()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
324 referenceTranscript.setChromosome(queryTranscript.getChromosome())
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
325 referenceTranscript.setName(queryTranscript.getChromosome())
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
326 referenceTranscript.setDirection("+")
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
327 referenceTranscript.setEnd(self.sizes[queryTranscript.getChromosome()])
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
328 referenceTranscript.setStart(1)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
329 return (referenceTranscript, queryTranscript)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
330
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
331 def extractReferenceQuery(self, inputTranscript):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
332 if "Target" not in inputTranscript.getTagNames():
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
333 raise Exception("Cannot extract Target field in line '%s'." % (inputTranscript))
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
334 id, start, end, strand = parseTargetField(inputTranscript.getTagValue("Target"))
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
335 if id not in self.sizes:
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
336 raise Exception("Target id '%s' of transcript '%s' does not correspond to anything in FASTA file." % (id, inputTranscript))
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
337 referenceTranscript = Transcript()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
338 referenceTranscript.setChromosome(id)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
339 referenceTranscript.setName(id)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
340 referenceTranscript.setDirection("+")
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
341 referenceTranscript.setEnd(self.sizes[id])
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
342 referenceTranscript.setStart(1)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
343 queryTranscript = Transcript()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
344 queryTranscript.setChromosome(id)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
345 queryTranscript.setName(id)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
346 queryTranscript.setStart(start)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
347 queryTranscript.setEnd(end)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
348 queryTranscript.setDirection(strand)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
349 if inputTranscript.getNbExons() > 1:
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
350 factor = float(end - start) / (inputTranscript.getEnd() - inputTranscript.getStart())
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
351 for exon in inputTranscript.getExons():
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
352 newExon = Interval()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
353 newExon.setChromosome(id)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
354 newExon.setDirection(strand)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
355 if "Target" in inputTranscript.getTagNames():
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
356 id, start, end, strand = parseTargetField(exon.getTagValue("Target"))
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
357 newExon.setStart(start)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
358 newExon.setEnd(end)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
359 else:
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
360 newExon.setStart(int(round((exon.getStart() - inputTranscript.getStart()) * factor)) + start)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
361 newExon.setEnd( int(round((exon.getEnd() - inputTranscript.getStart()) * factor)) + start)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
362 queryTranscript.addExon(newExon)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
363 return (referenceTranscript, queryTranscript)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
364
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
365 def compute1TranscriptFiles(self):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
366 progress = Progress(self.parsers[1].getNbItems(), "Comparing regions", self.verbosity)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
367 for transcript in self.parsers[1].getIterator():
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
368 if transcript.__class__.__name__ == "Mapping":
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
369 referenceTranscript, queryTranscript = self.extractReferenceQueryMapping(transcript)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
370 else:
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
371 referenceTranscript, queryTranscript = self.extractReferenceQuery(transcript)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
372 self.computeCoverage(queryTranscript, referenceTranscript, referenceTranscript.getName())
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
373 self.computeOverlap(queryTranscript, referenceTranscript, referenceTranscript.getName())
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
374 progress.inc()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
375 progress.done()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
376
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
377 def compute(self):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
378 if self.sequenceParser == None:
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
379 self.compute2TranscriptFiles()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
380 else:
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
381 self.compute1TranscriptFiles()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
382
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
383 def plotTranscript(self, index, transcript):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
384 plotter = Plotter(self.seed, index, self.verbosity)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
385 plotter.setOutputFileName(self.outputFileName)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
386 plotter.setTranscript(transcript)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
387 plotter.setTitle(self.title)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
388 plotter.setLabels(self.xLabel, self.yLabel)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
389 plotter.setPlotSize(self.width, self.height)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
390 plotter.setCoverageData(self.coverage[index])
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
391 plotter.setOverlapData(self.overlap[index])
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
392 plotter.setMerge(self.merge)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
393 plotter.plot()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
394 output = plotter.log
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
395 return output
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
396
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
397 def plot1TranscriptFile(self):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
398 self.outputCoverage = {}
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
399 self.outputCoveragePerStrand = {}
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
400 output = ""
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
401 progress = Progress(len(self.sequenceParser.getRegions()), "Plotting regions", self.verbosity)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
402 for cpt2, region in enumerate(self.sequenceParser.getRegions()):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
403 transcript = Transcript()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
404 transcript.setName(region)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
405 transcript.setDirection("+")
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
406 transcript.setEnd(self.sizes[region])
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
407 transcript.setStart(1)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
408 output += self.plotTranscript(region, transcript)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
409 progress.inc()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
410 progress.done()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
411 if self.verbosity > 0:
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
412 print output
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
413
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
414 def plot2TranscriptFiles(self):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
415 self.outputCoverage = [0] * self.parsers[1].getNbTranscripts()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
416 self.outputCoveragePerStrand = [None] * self.parsers[1].getNbTranscripts()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
417 for cpt in range(self.parsers[1].getNbTranscripts()):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
418 self.outputCoveragePerStrand[cpt] = dict([strand, 0] for strand in strands)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
419 progress = Progress(self.parsers[1].getNbTranscripts(), "Plotting regions", self.verbosity)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
420 output = ""
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
421 for cpt2, transcript2 in enumerate(self.parsers[1].getIterator()):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
422 output += self.plotTranscript(cpt2, transcript2)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
423 progress.inc()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
424 progress.done()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
425 if self.verbosity > 0:
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
426 print output
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
427
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
428 def plot(self):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
429 if self.sequenceParser == None:
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
430 self.plot2TranscriptFiles()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
431 else:
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
432 self.plot1TranscriptFile()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
433
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
434 def start(self):
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
435 self.initialize()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
436 self.compute()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
437 self.plot()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
438
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
439
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
440 if __name__ == "__main__":
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
441
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
442 # parse command line
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
443 description = "Plot Coverage v1.0.1: Plot the coverage of the first data with respect to the second one. [Category: Visualization]"
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
444
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
445 parser = OptionParser(description = description)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
446 parser.add_option("-i", "--input1", dest="inputFileName1", action="store", type="string", help="input file 1 [compulsory] [format: file in transcript or mapping format given by -f]")
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
447 parser.add_option("-f", "--inputFormat1", dest="inputFormat1", action="store", type="string", help="format of input file 1 [compulsory] [format: transcript or mapping file format]")
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
448 parser.add_option("-j", "--input2", dest="inputFileName2", action="store", type="string", help="input file 2 [compulsory] [format: file in transcript format given by -g]")
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
449 parser.add_option("-g", "--inputFormat2", dest="inputFormat2", action="store", type="string", help="format of input file 2 [compulsory] [format: transcript file format]")
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
450 parser.add_option("-q", "--sequence", dest="inputSequence", action="store", default=None, type="string", help="input sequence file [format: file in FASTA format] [default: None]")
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
451 parser.add_option("-o", "--output", dest="outputFileName", action="store", type="string", help="output file [compulsory] [format: output file in PNG format]")
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
452 parser.add_option("-w", "--width", dest="width", action="store", default=1500, type="int", help="width of the plots (in px) [format: int] [default: 1500]")
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
453 parser.add_option("-e", "--height", dest="height", action="store", default=1000, type="int", help="height of the plots (in px) [format: int] [default: 1000]")
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
454 parser.add_option("-t", "--title", dest="title", action="store", default="", type="string", help="title of the plots [format: string]")
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
455 parser.add_option("-x", "--xlab", dest="xLabel", action="store", default="", type="string", help="label on the x-axis [format: string]")
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
456 parser.add_option("-y", "--ylab", dest="yLabel", action="store", default="", type="string", help="label on the y-axis [format: string]")
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
457 parser.add_option("-p", "--plusColor", dest="plusColor", action="store", default="red", type="string", help="color for the elements on the plus strand [format: string] [default: red]")
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
458 parser.add_option("-m", "--minusColor", dest="minusColor", action="store", default="blue", type="string", help="color for the elements on the minus strand [format: string] [default: blue]")
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
459 parser.add_option("-s", "--sumColor", dest="sumColor", action="store", default="black", type="string", help="color for 2 strands coverage line [format: string] [default: black]")
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
460 parser.add_option("-l", "--lineColor", dest="lineColor", action="store", default="black", type="string", help="color for the lines [format: string] [default: black]")
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
461 parser.add_option("-1", "--merge", dest="merge", action="store_true", default=False, help="merge the 2 plots in 1 [format: boolean] [default: false]")
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
462 parser.add_option("-D", "--directory", dest="working_Dir", action="store", default=os.getcwd(), type="string", help="the directory to store the results [format: directory]")
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
463 parser.add_option("-v", "--verbosity", dest="verbosity", action="store", default=1, type="int", help="trace level [format: int]")
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
464 (options, args) = parser.parse_args()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
465
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
466 colors[1] = options.plusColor
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
467 colors[-1] = options.minusColor
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
468 colors[0] = options.sumColor
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
469 colorLine = options.lineColor
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
470
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
471 pp = PlotParser(options.verbosity)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
472 pp.addInput(0, options.inputFileName1, options.inputFormat1)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
473 pp.addInput(1, options.inputFileName2, options.inputFormat2)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
474 pp.addSequence(options.inputSequence)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
475 pp.setOutput(options.outputFileName if os.path.isabs(options.outputFileName) else os.path.join(options.working_Dir, options.outputFileName))
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
476 pp.setPlotSize(options.width, options.height)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
477 pp.setLabels(options.xLabel, options.yLabel)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
478 pp.setTitle(options.title)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
479 pp.setMerge(options.merge)
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
480 pp.start()
44d5973c188c Uploaded
m-zytnicki
parents:
diff changeset
481