diff SMART/Java/Python/plotCoverage.py @ 46:169d364ddd91

Uploaded
author m-zytnicki
date Mon, 30 Sep 2013 03:19:26 -0400
parents 44d5973c188c
children
line wrap: on
line diff
--- a/SMART/Java/Python/plotCoverage.py	Wed Sep 18 08:51:22 2013 -0400
+++ b/SMART/Java/Python/plotCoverage.py	Mon Sep 30 03:19:26 2013 -0400
@@ -32,7 +32,7 @@
 from optparse import OptionParser
 from SMART.Java.Python.structure.Interval import Interval
 from SMART.Java.Python.structure.Transcript import Transcript
-from commons.core.parsing.ParserChooser import ParserChooser
+from SMART.Java.Python.structure.TranscriptContainer import TranscriptContainer
 from SMART.Java.Python.misc.RPlotter import RPlotter
 from SMART.Java.Python.misc.Progress import Progress
 from commons.core.parsing.FastaParser import FastaParser
@@ -42,440 +42,430 @@
 colorLine = "black"
 
 def parseTargetField(field):
-	strand             = "+"
-	splittedFieldSpace = field.split()
-	splittedFieldPlus  = field.split("+", 4)
-	if len(splittedFieldSpace) == 3:
-		id, start, end = splittedFieldSpace
-	elif len(splittedFieldSpace) == 4:
-		id, start, end, strand = splittedFieldSpace
-	elif len(splittedFieldPlus) == 3:
-		id, start, end = splittedFieldPlus
-	elif len(splittedFieldPlus) == 4:
-		id, start, end, strand = splittedFieldPlus
-	else:
-		raise Exception("Cannot parse Target field '%s'." % (field))
-	return (id, int(start), int(end), strand)
+    strand             = "+"
+    splittedFieldSpace = field.split()
+    splittedFieldPlus  = field.split("+", 4)
+    if len(splittedFieldSpace) == 3:
+        id, start, end = splittedFieldSpace
+    elif len(splittedFieldSpace) == 4:
+        id, start, end, strand = splittedFieldSpace
+    elif len(splittedFieldPlus) == 3:
+        id, start, end = splittedFieldPlus
+    elif len(splittedFieldPlus) == 4:
+        id, start, end, strand = splittedFieldPlus
+    else:
+        raise Exception("Cannot parse Target field '%s'." % (field))
+    return (id, int(start), int(end), strand)
 
 
 class SimpleTranscript(object):
-	def __init__(self, transcript1, transcript2, color = None):
-		self.start  = max(0, transcript1.getStart() - transcript2.getStart())
-		self.end    = min(transcript2.getEnd() - transcript2.getStart(), transcript1.getEnd() - transcript2.getStart())
-		self.strand = transcript1.getDirection() * transcript2.getDirection()
-		self.exons  = []
-		for exon in transcript1.getExons():
-			if exon.getEnd() >= transcript2.getStart() and exon.getStart() <= transcript2.getEnd():
-				start = max(0, exon.getStart() - transcript2.getStart())
-				end   = min(transcript2.getEnd() - transcript2.getStart(), exon.getEnd() - transcript2.getStart())
-				self.addExon(start, end, self.strand, color)
+    def __init__(self, transcript1, transcript2, color = None):
+        self.start  = max(0, transcript1.getStart() - transcript2.getStart())
+        self.end    = min(transcript2.getEnd() - transcript2.getStart(), transcript1.getEnd() - transcript2.getStart())
+        self.strand = transcript1.getDirection() * transcript2.getDirection()
+        self.exons  = []
+        for exon in transcript1.getExons():
+            if exon.getEnd() >= transcript2.getStart() and exon.getStart() <= transcript2.getEnd():
+                start = max(0, exon.getStart() - transcript2.getStart())
+                end   = min(transcript2.getEnd() - transcript2.getStart(), exon.getEnd() - transcript2.getStart())
+                self.addExon(start, end, self.strand, color)
 
-	def addExon(self, start, end, strand, color):
-		exon = SimpleExon(start, end, strand, color)
-		self.exons.append(exon)
+    def addExon(self, start, end, strand, color):
+        exon = SimpleExon(start, end, strand, color)
+        self.exons.append(exon)
 
-	def getRScript(self, yOffset, height):
-		rString     = ""
-		previousEnd = None
-		for exon in sorted(self.exons, key=lambda exon: exon.start):
-			if previousEnd != None:
-				rString += "segments(%.1f, %.1f, %.1f, %.1f, col = \"%s\")\n" % (previousEnd, yOffset + height / 4.0, exon.start, yOffset + height / 4.0, colorLine)
-			rString    += exon.getRScript(yOffset, height)
-			previousEnd = exon.end
-		return rString
+    def getRScript(self, yOffset, height):
+        rString     = ""
+        previousEnd = None
+        for exon in sorted(self.exons, key=lambda exon: exon.start):
+            if previousEnd != None:
+                rString += "segments(%.1f, %.1f, %.1f, %.1f, col = \"%s\")\n" % (previousEnd, yOffset + height / 4.0, exon.start, yOffset + height / 4.0, colorLine)
+            rString    += exon.getRScript(yOffset, height)
+            previousEnd = exon.end
+        return rString
 
 
 class SimpleExon(object):
-	def __init__(self, start, end, strand, color = None):
-		self.start  = start
-		self.end    = end
-		self.strand = strand
-		self.color  = color
-		
-	def getRScript(self, yOffset, height):
-		color = self.color if self.color != None else colors[self.strand]
-		return "rect(%.1f, %.1f, %.1f, %.1f, col=\"%s\", border = \"%s\")\n" % (self.start, yOffset, self.end, yOffset + height / 2.0, color, colorLine)
+    def __init__(self, start, end, strand, color = None):
+        self.start  = start
+        self.end    = end
+        self.strand = strand
+        self.color  = color
+        
+    def getRScript(self, yOffset, height):
+        color = self.color if self.color != None else colors[self.strand]
+        return "rect(%.1f, %.1f, %.1f, %.1f, col=\"%s\", border = \"%s\")\n" % (self.start, yOffset, self.end, yOffset + height / 2.0, color, colorLine)
 
 
 class Plotter(object):
-	
-	def __init__(self, seed, index, verbosity):
-		self.seed        = seed
-		self.index       = index
-		self.verbosity   = verbosity
-		self.maxCoverage = 0
-		self.maxOverlap  = 0
-		self.log         = ""
-		self.merge       = False
-		self.width       = 1500
-		self.heigth      = 1000
-		self.xLabel      = ""
-		self.yLabel      = ""
-		self.title       = None
-		self.absPath     = os.getcwd()
-		self.coverageDataFileName    = "tmpFile_%d_%s.dat" % (seed, index)
-		self.coverageScript          = ""
-		self.overlapScript           = ""
-		self.outputFileName          = None
+    
+    def __init__(self, seed, index, verbosity):
+        self.seed        = seed
+        self.index       = index
+        self.verbosity   = verbosity
+        self.maxCoverage = 0
+        self.maxOverlap  = 0
+        self.log         = ""
+        self.merge       = False
+        self.width       = 1500
+        self.heigth      = 1000
+        self.xLabel      = ""
+        self.yLabel      = ""
+        self.title       = None
+        self.absPath     = os.getcwd()
+        self.coverageDataFileName    = "tmpFile_%d_%s.dat" % (seed, index)
+        self.coverageScript          = ""
+        self.overlapScript           = ""
+        self.outputFileName          = None
 
-	def setOutputFileName(self, fileName):
-		self.outputFileName = fileName
+    def setOutputFileName(self, fileName):
+        self.outputFileName = fileName
 
-	def setTranscript(self, transcript):
-		self.transcript = transcript
-		self.name       = transcript.getName()
-		self.size       = transcript.getEnd() - transcript.getStart() + 1
-		if self.title == None:
-			self.title = self.name
-		else:
-			self.title += " " + self.name
+    def setTranscript(self, transcript):
+        self.transcript = transcript
+        self.name       = transcript.getName()
+        self.size       = transcript.getEnd() - transcript.getStart() + 1
+        if self.title == None:
+            self.title = self.name
+        else:
+            self.title += " " + self.name
 
-	def setTitle(self, title):
-		self.title = title + " " + self.name
+    def setTitle(self, title):
+        self.title = title + " " + self.name
 
-	def setPlotSize(self, width, height):
-		self.width  = width
-		self.height = height
+    def setPlotSize(self, width, height):
+        self.width  = width
+        self.height = height
 
-	def setLabels(self, xLabel, yLabel):
-		self.xLabel = xLabel
-		self.yLabel = yLabel
+    def setLabels(self, xLabel, yLabel):
+        self.xLabel = xLabel
+        self.yLabel = yLabel
 
-	def setMerge(self, merge):
-		self.merge = merge
+    def setMerge(self, merge):
+        self.merge = merge
 
-	def setCoverageData(self, coverage):
-		outputCoveragePerStrand = dict([strand, 0] for strand in strands)
-		outputCoverage          = 0
-		dataFile = open(os.path.abspath(self.coverageDataFileName), "w")
-		for position in range(self.size+1):
-			sumValue = 0
-			found    = False
-			dataFile.write("%d\t" % (position))
-			for strand in strands:
-				value     = coverage[strand].get(position, 0)
-				sumValue += value
-				dataFile.write("%d\t" % (value))
-				if value > 0:
-					found = True
-					outputCoveragePerStrand[strand] += 1
-			self.maxCoverage = max(self.maxCoverage, sumValue)
-			dataFile.write("%d\n" % (sumValue))
-			if found:
-				outputCoverage += 1
-		dataFile.close()
-		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) 
-		self.coverageScript += "data = scan(\"%s\", list(pos = -666, minus = -666, plus = -666, sumValue = -666), sep=\"\t\")\n" % (os.path.abspath(self.coverageDataFileName))
-		self.coverageScript += "lines(x = data$pos, y = data$minus,    col = \"%s\")\n" % (colors[-1])
-		self.coverageScript += "lines(x = data$pos, y = data$plus,     col = \"%s\")\n" % (colors[1])
-		self.coverageScript += "lines(x = data$pos, y = data$sumValue, col = \"%s\")\n" % (colors[0])
+    def setCoverageData(self, coverage):
+        outputCoveragePerStrand = dict([strand, 0] for strand in strands)
+        outputCoverage          = 0
+        dataFile = open(os.path.abspath(self.coverageDataFileName), "w")
+        for position in range(self.size+1):
+            sumValue = 0
+            found    = False
+            dataFile.write("%d\t" % (position))
+            for strand in strands:
+                value     = coverage[strand].get(position, 0)
+                sumValue += value
+                dataFile.write("%d\t" % (value))
+                if value > 0:
+                    found = True
+                    outputCoveragePerStrand[strand] += 1
+            self.maxCoverage = max(self.maxCoverage, sumValue)
+            dataFile.write("%d\n" % (sumValue))
+            if found:
+                outputCoverage += 1
+        dataFile.close()
+        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) 
+        self.coverageScript += "data = scan(\"%s\", list(pos = -666, minus = -666, plus = -666, sumValue = -666), sep=\"\t\")\n" % (os.path.abspath(self.coverageDataFileName))
+        self.coverageScript += "lines(x = data$pos, y = data$minus,    col = \"%s\")\n" % (colors[-1])
+        self.coverageScript += "lines(x = data$pos, y = data$plus,     col = \"%s\")\n" % (colors[1])
+        self.coverageScript += "lines(x = data$pos, y = data$sumValue, col = \"%s\")\n" % (colors[0])
 
-	def setOverlapData(self, overlap):
-		height              = 1
-		self.maxOverlap     = (len(overlap) + 1) * height
-		thisElement         = SimpleTranscript(self.transcript, self.transcript, "black")
-		self.overlapScript += thisElement.getRScript(0, height)
-		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)):
-			self.overlapScript += transcript.getRScript((cpt + 1) * height, height)
+    def setOverlapData(self, overlap):
+        height              = 1
+        self.maxOverlap     = (len(overlap) + 1) * height
+        thisElement         = SimpleTranscript(self.transcript, self.transcript, "black")
+        self.overlapScript += thisElement.getRScript(0, height)
+        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)):
+            self.overlapScript += transcript.getRScript((cpt + 1) * height, height)
 
-	def getFirstLine(self, suffix = None):
-		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)
+    def getFirstLine(self, suffix = None):
+        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)
 
-	def getLastLine(self):
-		return "dev.off()\n"
+    def getLastLine(self):
+        return "dev.off()\n"
 
-	def startR(self, fileName, script):
-		scriptFile = open(fileName, "w")
-		scriptFile.write(script)
-		scriptFile.close()
-		command = "R CMD BATCH %s" % (fileName)
-		status  = subprocess.call(command, shell=True)
-		if status != 0:
-			raise Exception("Problem with the execution of script file %s, status is: %s" % (fileName, status))
+    def startR(self, fileName, script):
+        scriptFile = open(fileName, "w")
+        scriptFile.write(script)
+        scriptFile.close()
+        command = "R CMD BATCH %s" % (fileName)
+        status  = subprocess.call(command, shell=True)
+        if status != 0:
+            raise Exception("Problem with the execution of script file %s, status is: %s" % (fileName, status))
 
-	def plot(self):
-		if self.merge:
-			fileName = "%s_%d_%s.R" % (self.outputFileName, self.seed, self.index)
-			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)
-			script   = self.getFirstLine() + plotLine + self.overlapScript + self.coverageScript + self.getLastLine()
-			self.startR(fileName, script)
-		else:
-			fileName = "%s_%d_%s_overlap.R" % (self.outputFileName, self.seed, self.index)
-			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)
-			script   = self.getFirstLine("overlap") + plotLine + self.overlapScript + self.getLastLine()
-			self.startR(fileName, script)
-			fileName = "%s_%d_%s_coverage.R" % (self.outputFileName, self.seed, self.index)
-			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)
-			script   = self.getFirstLine("coverage") + plotLine + self.coverageScript + self.getLastLine()
-			self.startR(fileName, script)
+    def plot(self):
+        print "outputfileName is written in :", self.outputFileName
+        if self.merge:
+            fileName = "%s_%d_%s.R" % (self.outputFileName, self.seed, self.index)
+            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)
+            script   = self.getFirstLine() + plotLine + self.overlapScript + self.coverageScript + self.getLastLine()
+            self.startR(fileName, script)
+        else:
+            fileName = "%s_%d_%s_overlap.R" % (self.outputFileName, self.seed, self.index)
+            print "overlap file is written in :", fileName
+            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)
+            script   = self.getFirstLine("overlap") + plotLine + self.overlapScript + self.getLastLine()
+            self.startR(fileName, script)
+            fileName = "%s_%d_%s_coverage.R" % (self.outputFileName, self.seed, self.index)
+            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)
+            script   = self.getFirstLine("coverage") + plotLine + self.coverageScript + self.getLastLine()
+            self.startR(fileName, script)
 
 
 class PlotParser(object):
 
-	def __init__(self, verbosity):
-		self.verbosity      = verbosity
-		self.parsers        = [None, None]
-		self.sequenceParser = None
-		self.seed           = random.randint(0, 10000)
-		self.title          = ""
-		self.merge          = False
+    def __init__(self, verbosity):
+        self.verbosity      = verbosity
+        self.parsers        = [None, None]
+        self.sequenceParser = None
+        self.seed           = random.randint(0, 10000)
+        self.title          = ""
+        self.merge          = False
 
-	def __del__(self):
-		for fileName in glob.glob("tmpFile_%d*.dat" % (self.seed)):
-			os.remove(fileName)
-		for fileName in glob.glob("%s*.R" % (os.path.abspath(self.outputFileName))):
-			os.remove(fileName)
-		for fileName in glob.glob("%s*.Rout" % (os.path.abspath(self.outputFileName))):
-			os.remove(fileName)
+    def __del__(self):
+        for fileName in glob.glob("tmpFile_%d*.dat" % (self.seed)):
+            os.remove(fileName)
+        for fileName in glob.glob("%s*.R" % (os.path.abspath(self.outputFileName))):
+            os.remove(fileName)
+        for fileName in glob.glob("%s*.Rout" % (os.path.abspath(self.outputFileName))):
+            os.remove(fileName)
 
-	def addInput(self, inputNb, fileName, fileFormat):
-		if fileName == None:
-			return
-		chooser = ParserChooser(self.verbosity)
-		chooser.findFormat(fileFormat)
-		self.parsers[inputNb] = chooser.getParser(fileName)
-		if inputNb == 0:
-			self.parsers[1] = self.parsers[0]
+    def addInput(self, inputNb, fileName, fileFormat):
+        if fileName == None:
+            return
+        self.parsers[inputNb] = TranscriptContainer(fileName, fileFormat, self.verbosity)
+        if inputNb == 0:
+            self.parsers[1] = self.parsers[0]
 
-	def addSequence(self, fileName):
-		if fileName == None:
-			return
-		self.sequenceParser = FastaParser(fileName, self.verbosity)
+    def addSequence(self, fileName):
+        if fileName == None:
+            return
+        self.sequenceParser = FastaParser(fileName, self.verbosity)
 
-	def setOutput(self, fileName):
-		self.outputFileName = fileName
+    def setOutput(self, fileName):
+        self.outputFileName = fileName
 
-	def setPlotSize(self, width, height):
-		self.width  = width
-		self.height = height
+    def setPlotSize(self, width, height):
+        self.width  = width
+        self.height = height
 
-	def setLabels(self, xLabel, yLabel):
-		self.xLabel = xLabel
-		self.yLabel = yLabel
+    def setLabels(self, xLabel, yLabel):
+        self.xLabel = xLabel
+        self.yLabel = yLabel
 
-	def setTitle(self, title):
-		self.title = title
+    def setTitle(self, title):
+        self.title = title
 
-	def setMerge(self, merge):
-		self.merge = merge
+    def setMerge(self, merge):
+        self.merge = merge
 
-	def initializeDataFromSequences(self):
-		self.sizes    = {}
-		self.coverage = {}
-		self.overlap  = {}
-		for region in self.sequenceParser.getRegions():
-			self.sizes[region]    = self.sequenceParser.getSizeOfRegion(region)
-			self.coverage[region] = {}
-			self.overlap[region]  = []
-			for strand in strands:
-				self.coverage[region][strand] = {}
-				self.coverage[region][strand][1] = 0
-				self.coverage[region][strand][self.sizes[region]] = 0
+    def initializeDataFromSequences(self):
+        self.sizes    = {}
+        self.coverage = {}
+        self.overlap  = {}
+        for region in self.sequenceParser.getRegions():
+            self.sizes[region]    = self.sequenceParser.getSizeOfRegion(region)
+            self.coverage[region] = {}
+            self.overlap[region]  = []
+            for strand in strands:
+                self.coverage[region][strand] = {}
+                self.coverage[region][strand][1] = 0
+                self.coverage[region][strand][self.sizes[region]] = 0
+
 
-	def initializeDataFromTranscripts(self):
-		self.coverage = dict([i, None] for i in range(self.parsers[1].getNbTranscripts()))
-		self.overlap  = dict([i, None] for i in range(self.parsers[1].getNbTranscripts()))
-		self.sizes    = dict([i, 0]    for i in range(self.parsers[1].getNbTranscripts()))
-		progress = Progress(self.parsers[1].getNbTranscripts(), "Reading regions", self.verbosity)
-		for cpt, transcript in enumerate(self.parsers[1].getIterator()):
-			self.coverage[cpt] = {}
-			self.overlap[cpt]  = []
-			for strand in strands:
-				self.coverage[cpt][strand] = {}
-				self.coverage[cpt][strand][0] = 0
-				self.coverage[cpt][strand][transcript.getEnd() - transcript.getStart()] = 0
-			for exon in transcript.getExons():
-				self.sizes[cpt] += exon.getSize()
-			progress.inc()
-		progress.done()
-
-	def initialize(self):
-		if self.sequenceParser == None:
-			self.initializeDataFromTranscripts()
-		else:
-			self.initializeDataFromSequences()
+    def initializeDataFromTranscripts(self):
+        self.coverage = dict([i, None] for i in range(self.parsers[1].getNbTranscripts()))
+        self.overlap  = dict([i, None] for i in range(self.parsers[1].getNbTranscripts()))
+        self.sizes    = dict([i, 0]    for i in range(self.parsers[1].getNbTranscripts()))
+        self.parsers[0].findData()
+        progress = Progress(self.parsers[1].getNbTranscripts(), "Reading regions", self.verbosity)
+        for cpt, transcript in enumerate(self.parsers[1].getIterator()):
+            self.coverage[cpt] = {}
+            self.overlap[cpt]  = []
+            for strand in strands:
+                self.coverage[cpt][strand] = {}
+                self.coverage[cpt][strand][0] = 0
+                self.coverage[cpt][strand][transcript.getEnd() - transcript.getStart()] = 0
+            for exon in transcript.getExons():
+                self.sizes[cpt] += exon.getSize()
+            progress.inc()
+        progress.done()
 
-	def computeCoverage(self, transcript1, transcript2, id):
-		strand = transcript1.getDirection() * transcript2.getDirection()
-		for exon1 in transcript1.getExons():
-			for exon2 in transcript2.getExons():
-				if exon1.overlapWith(exon2):
-					for position in range(max(exon1.getStart(), exon2.getStart()), min(exon1.getEnd(), exon2.getEnd()) + 1):
-						relativePosition = position - transcript2.getStart() + 1
-						self.coverage[id][strand][relativePosition] = self.coverage[id][strand].get(relativePosition, 0) + 1
+    def initialize(self):
+        if self.sequenceParser == None:
+            self.initializeDataFromTranscripts()
+        else:
+            self.initializeDataFromSequences()
 
-	def computeOverlap(self, transcript1, transcript2, id):
-		simpleTranscript = SimpleTranscript(transcript1, transcript2)
-		self.overlap[id].append(simpleTranscript)
-		
-	def compute2TranscriptFiles(self):
-		progress = Progress(self.parsers[1].getNbTranscripts(), "Comparing regions", self.verbosity)
-		for cpt2, transcript2 in enumerate(self.parsers[1].getIterator()):
-			for transcript1 in self.parsers[0].getIterator():
-				if transcript1.overlapWithExon(transcript2):
-					self.computeCoverage(transcript1, transcript2, cpt2)
-					self.computeOverlap(transcript1, transcript2, cpt2)
-			progress.inc()
-		progress.done()
+    def computeCoverage(self, transcript1, transcript2, id):
+        strand = transcript1.getDirection() * transcript2.getDirection()
+        for exon1 in transcript1.getExons():
+            for exon2 in transcript2.getExons():
+                if exon1.overlapWith(exon2):
+                    for position in range(max(exon1.getStart(), exon2.getStart()), min(exon1.getEnd(), exon2.getEnd()) + 1):
+                        relativePosition = position - transcript2.getStart() + 1
+                        self.coverage[id][strand][relativePosition] = self.coverage[id][strand].get(relativePosition, 0) + 1
 
-	def extractReferenceQueryMapping(self, mapping):
-		queryTranscript = mapping.getTranscript()
-		referenceTranscript = Transcript()
-		referenceTranscript.setChromosome(queryTranscript.getChromosome())
-		referenceTranscript.setName(queryTranscript.getChromosome())
-		referenceTranscript.setDirection("+")
-		referenceTranscript.setEnd(self.sizes[queryTranscript.getChromosome()])
-		referenceTranscript.setStart(1)
-		return (referenceTranscript, queryTranscript)
+    def computeOverlap(self, transcript1, transcript2, id):
+        simpleTranscript = SimpleTranscript(transcript1, transcript2)
+        self.overlap[id].append(simpleTranscript)
+        
+    def compute2TranscriptFiles(self):
+        progress = Progress(self.parsers[1].getNbTranscripts(), "Comparing regions", self.verbosity)
+        for cpt2, transcript2 in enumerate(self.parsers[1].getIterator()):
+            for transcript1 in self.parsers[0].getIterator():
+                if transcript1.overlapWithExon(transcript2):
+                    self.computeCoverage(transcript1, transcript2, cpt2)
+                    self.computeOverlap(transcript1, transcript2, cpt2)
+            progress.inc()
+        progress.done()
 
-	def extractReferenceQuery(self, inputTranscript):
-		if "Target" not in inputTranscript.getTagNames():
-			raise Exception("Cannot extract Target field in line '%s'." % (inputTranscript))
-		id, start, end, strand = parseTargetField(inputTranscript.getTagValue("Target"))
-		if id not in self.sizes:
-			raise Exception("Target id '%s' of transcript '%s' does not correspond to anything in FASTA file." % (id, inputTranscript))
-		referenceTranscript = Transcript()
-		referenceTranscript.setChromosome(id)
-		referenceTranscript.setName(id)
-		referenceTranscript.setDirection("+")
-		referenceTranscript.setEnd(self.sizes[id])
-		referenceTranscript.setStart(1)
-		queryTranscript = Transcript()
-		queryTranscript.setChromosome(id)
-		queryTranscript.setName(id)
-		queryTranscript.setStart(start)
-		queryTranscript.setEnd(end)
-		queryTranscript.setDirection(strand)
-		if inputTranscript.getNbExons() > 1:
-			factor = float(end - start) / (inputTranscript.getEnd() - inputTranscript.getStart())
-			for exon in inputTranscript.getExons():
-				newExon = Interval()
-				newExon.setChromosome(id)
-				newExon.setDirection(strand)
-				if "Target" in inputTranscript.getTagNames():
-					id, start, end, strand = parseTargetField(exon.getTagValue("Target"))
-					newExon.setStart(start)
-					newExon.setEnd(end)
-				else:
-					newExon.setStart(int(round((exon.getStart() - inputTranscript.getStart()) * factor)) + start)
-					newExon.setEnd(  int(round((exon.getEnd() -   inputTranscript.getStart()) * factor)) + start)
-				queryTranscript.addExon(newExon)
-		return (referenceTranscript, queryTranscript)
+    def extractReferenceQuery(self, inputTranscript):
+        if "Target" not in inputTranscript.getTagNames():
+            raise Exception("Cannot extract Target field in line '%s'." % (inputTranscript))
+        id, start, end, strand = parseTargetField(inputTranscript.getTagValue("Target"))
+        if id not in self.sizes:
+            raise Exception("Target id '%s' of transcript '%s' does not correspond to anything in FASTA file." % (id, inputTranscript))
+        referenceTranscript = Transcript()
+        referenceTranscript.setChromosome(id)
+        referenceTranscript.setName(id)
+        referenceTranscript.setDirection("+")
+        referenceTranscript.setEnd(self.sizes[id])
+        referenceTranscript.setStart(1)
+        queryTranscript = Transcript()
+        queryTranscript.setChromosome(id)
+        queryTranscript.setName(id)
+        queryTranscript.setStart(start)
+        queryTranscript.setEnd(end)
+        queryTranscript.setDirection(strand)
+        if inputTranscript.getNbExons() > 1:
+            factor = float(end - start) / (inputTranscript.getEnd() - inputTranscript.getStart())
+            for exon in inputTranscript.getExons():
+                newExon = Interval()
+                newExon.setChromosome(id)
+                newExon.setDirection(strand)
+                if "Target" in inputTranscript.getTagNames():
+                    id, start, end, strand = parseTargetField(exon.getTagValue("Target"))
+                    newExon.setStart(start)
+                    newExon.setEnd(end)
+                else:
+                    newExon.setStart(int(round((exon.getStart() - inputTranscript.getStart()) * factor)) + start)
+                    newExon.setEnd(  int(round((exon.getEnd() -   inputTranscript.getStart()) * factor)) + start)
+                queryTranscript.addExon(newExon)
+        return (referenceTranscript, queryTranscript)
 
-	def compute1TranscriptFiles(self):
-		progress = Progress(self.parsers[1].getNbItems(), "Comparing regions", self.verbosity)
-		for transcript in self.parsers[1].getIterator():
-			if transcript.__class__.__name__ == "Mapping":
-				referenceTranscript, queryTranscript = self.extractReferenceQueryMapping(transcript)
-			else:
-				referenceTranscript, queryTranscript = self.extractReferenceQuery(transcript)
-			self.computeCoverage(queryTranscript, referenceTranscript, referenceTranscript.getName())
-			self.computeOverlap(queryTranscript, referenceTranscript, referenceTranscript.getName())
-			progress.inc()
-		progress.done()
+    def compute1TranscriptFiles(self):
+        progress = Progress(self.parsers[1].getNbTranscripts(), "Comparing regions", self.verbosity)
+        for transcript in self.parsers[1].getIterator():
+            referenceTranscript, queryTranscript = self.extractReferenceQuery(transcript)
+            self.computeCoverage(queryTranscript, referenceTranscript, referenceTranscript.getName())
+            self.computeOverlap(queryTranscript, referenceTranscript, referenceTranscript.getName())
+            progress.inc()
+        progress.done()
 
-	def compute(self):
-		if self.sequenceParser == None:
-			self.compute2TranscriptFiles()
-		else:
-			self.compute1TranscriptFiles()
+    def compute(self):
+        if self.sequenceParser == None:
+            self.compute2TranscriptFiles()
+        else:
+            self.compute1TranscriptFiles()
 
-	def plotTranscript(self, index, transcript):
-		plotter = Plotter(self.seed, index, self.verbosity)
-		plotter.setOutputFileName(self.outputFileName)
-		plotter.setTranscript(transcript)
-		plotter.setTitle(self.title)
-		plotter.setLabels(self.xLabel, self.yLabel)
-		plotter.setPlotSize(self.width, self.height)
-		plotter.setCoverageData(self.coverage[index])
-		plotter.setOverlapData(self.overlap[index])
-		plotter.setMerge(self.merge)
-		plotter.plot()
-		output = plotter.log
-		return output
-		
-	def plot1TranscriptFile(self):
-		self.outputCoverage          = {}
-		self.outputCoveragePerStrand = {}
-		output   = ""
-		progress = Progress(len(self.sequenceParser.getRegions()), "Plotting regions", self.verbosity)
-		for cpt2, region in enumerate(self.sequenceParser.getRegions()):
-			transcript = Transcript()
-			transcript.setName(region)
-			transcript.setDirection("+")
-			transcript.setEnd(self.sizes[region])
-			transcript.setStart(1)
-			output += self.plotTranscript(region, transcript)
-			progress.inc()
-		progress.done()
-		if self.verbosity > 0:
-			print output
+    def plotTranscript(self, index, transcript):
+        plotter = Plotter(self.seed, index, self.verbosity)
+        plotter.setOutputFileName(self.outputFileName)
+        plotter.setTranscript(transcript)
+        plotter.setTitle(self.title)
+        plotter.setLabels(self.xLabel, self.yLabel)
+        plotter.setPlotSize(self.width, self.height)
+        plotter.setCoverageData(self.coverage[index])
+        plotter.setOverlapData(self.overlap[index])
+        plotter.setMerge(self.merge)
+        plotter.plot()
+        output = plotter.log
+        return output
+        
+    def plot1TranscriptFile(self):
+        self.outputCoverage          = {}
+        self.outputCoveragePerStrand = {}
+        output   = ""
+        progress = Progress(len(self.sequenceParser.getRegions()), "Plotting regions", self.verbosity)
+        for cpt2, region in enumerate(self.sequenceParser.getRegions()):
+            transcript = Transcript()
+            transcript.setName(region)
+            transcript.setDirection("+")
+            transcript.setEnd(self.sizes[region])
+            transcript.setStart(1)
+            output += self.plotTranscript(region, transcript)
+            progress.inc()
+        progress.done()
+        if self.verbosity > 0:
+            print output
 
-	def plot2TranscriptFiles(self):
-		self.outputCoverage          = [0] * self.parsers[1].getNbTranscripts()
-		self.outputCoveragePerStrand = [None] * self.parsers[1].getNbTranscripts()
-		for cpt in range(self.parsers[1].getNbTranscripts()):
-			self.outputCoveragePerStrand[cpt] = dict([strand, 0] for strand in strands)
-		progress = Progress(self.parsers[1].getNbTranscripts(), "Plotting regions", self.verbosity)
-		output = ""
-		for cpt2, transcript2 in enumerate(self.parsers[1].getIterator()):
-			output += self.plotTranscript(cpt2, transcript2)
-			progress.inc()
-		progress.done()
-		if self.verbosity > 0:
-			print output
+    def plot2TranscriptFiles(self):
+        self.outputCoverage          = [0] * self.parsers[1].getNbTranscripts()
+        self.outputCoveragePerStrand = [None] * self.parsers[1].getNbTranscripts()
+        for cpt in range(self.parsers[1].getNbTranscripts()):
+            self.outputCoveragePerStrand[cpt] = dict([strand, 0] for strand in strands)
+        progress = Progress(self.parsers[1].getNbTranscripts(), "Plotting regions", self.verbosity)
+        output = ""
+        for cpt2, transcript2 in enumerate(self.parsers[1].getIterator()):
+            output += self.plotTranscript(cpt2, transcript2)
+            progress.inc()
+        progress.done()
+        if self.verbosity > 0:
+            print output
 
-	def plot(self):
-		if self.sequenceParser == None:
-			self.plot2TranscriptFiles()
-		else:
-			self.plot1TranscriptFile()
+    def plot(self):
+        if self.sequenceParser == None:
+            self.plot2TranscriptFiles()
+        else:
+            self.plot1TranscriptFile()
 
-	def start(self):
-		self.initialize()
-		self.compute()
-		self.plot()
+    def start(self):
+        self.initialize()
+        self.compute()
+        self.plot()
 
 
 if __name__ == "__main__":
-	
-	# parse command line
-	description = "Plot Coverage v1.0.1: Plot the coverage of the first data with respect to the second one. [Category: Visualization]"
+    
+    # parse command line
+    description = "Plot Coverage v1.0.1: Plot the coverage of the first data with respect to the second one. [Category: Visualization]"
 
-	parser = OptionParser(description = description)
-	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]")
-	parser.add_option("-f", "--inputFormat1", dest="inputFormat1",   action="store",                       type="string", help="format of input file 1 [compulsory] [format: transcript or mapping file format]")
-	parser.add_option("-j", "--input2",       dest="inputFileName2", action="store",                       type="string", help="input file 2 [compulsory] [format: file in transcript format given by -g]")
-	parser.add_option("-g", "--inputFormat2", dest="inputFormat2",   action="store",                       type="string", help="format of input file 2 [compulsory] [format: transcript file format]")
-	parser.add_option("-q", "--sequence",     dest="inputSequence",  action="store",      default=None,    type="string", help="input sequence file [format: file in FASTA format] [default: None]")
-	parser.add_option("-o", "--output",       dest="outputFileName", action="store",                       type="string", help="output file [compulsory] [format: output file in PNG format]")
-	parser.add_option("-w", "--width",        dest="width",          action="store",      default=1500,    type="int",    help="width of the plots (in px) [format: int] [default: 1500]")
-	parser.add_option("-e", "--height",       dest="height",         action="store",      default=1000,    type="int",    help="height of the plots (in px) [format: int] [default: 1000]")
-	parser.add_option("-t", "--title",        dest="title",          action="store",      default="",      type="string", help="title of the plots [format: string]")
-	parser.add_option("-x", "--xlab",         dest="xLabel",         action="store",      default="",      type="string", help="label on the x-axis [format: string]")
-	parser.add_option("-y", "--ylab",         dest="yLabel",         action="store",      default="",      type="string", help="label on the y-axis [format: string]")
-	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]")
-	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]")
-	parser.add_option("-s", "--sumColor",     dest="sumColor",       action="store",      default="black", type="string", help="color for 2 strands coverage line [format: string] [default: black]")
-	parser.add_option("-l", "--lineColor",    dest="lineColor",      action="store",      default="black", type="string", help="color for the lines [format: string] [default: black]")
-	parser.add_option("-1", "--merge",        dest="merge",          action="store_true", default=False,                  help="merge the 2 plots in 1 [format: boolean] [default: false]")
-	parser.add_option("-D", "--directory",    dest="working_Dir",    action="store",      default=os.getcwd(), type="string", help="the directory to store the results [format: directory]")
-	parser.add_option("-v", "--verbosity",    dest="verbosity",      action="store",      default=1,       type="int",    help="trace level [format: int]")
-	(options, args) = parser.parse_args()
+    parser = OptionParser(description = description)
+    parser.add_option("-i", "--input1",       dest="inputFileName1", action="store",                       type="string", help="input file 1 [compulsory] [format: file in transcript format given by -f]")
+    parser.add_option("-f", "--inputFormat1", dest="inputFormat1",   action="store",                       type="string", help="format of input file 1 [compulsory] [format: transcript file format]")
+    parser.add_option("-j", "--input2",       dest="inputFileName2", action="store",                       type="string", help="input file 2 [compulsory] [format: file in transcript format given by -g]")
+    parser.add_option("-g", "--inputFormat2", dest="inputFormat2",   action="store",                       type="string", help="format of input file 2 [compulsory] [format: transcript file format]")
+    parser.add_option("-q", "--sequence",     dest="inputSequence",  action="store",      default=None,    type="string", help="input sequence file [format: file in FASTA format] [default: None]")
+    parser.add_option("-o", "--output",       dest="outputFileName", action="store",                       type="string", help="output file [compulsory] [format: output file in PNG format]")
+    parser.add_option("-w", "--width",        dest="width",          action="store",      default=1500,    type="int",    help="width of the plots (in px) [format: int] [default: 1500]")
+    parser.add_option("-e", "--height",       dest="height",         action="store",      default=1000,    type="int",    help="height of the plots (in px) [format: int] [default: 1000]")
+    parser.add_option("-t", "--title",        dest="title",          action="store",      default="",      type="string", help="title of the plots [format: string]")
+    parser.add_option("-x", "--xlab",         dest="xLabel",         action="store",      default="",      type="string", help="label on the x-axis [format: string]")
+    parser.add_option("-y", "--ylab",         dest="yLabel",         action="store",      default="",      type="string", help="label on the y-axis [format: string]")
+    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]")
+    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]")
+    parser.add_option("-s", "--sumColor",     dest="sumColor",       action="store",      default="black", type="string", help="color for 2 strands coverage line [format: string] [default: black]")
+    parser.add_option("-l", "--lineColor",    dest="lineColor",      action="store",      default="black", type="string", help="color for the lines [format: string] [default: black]")
+    parser.add_option("-1", "--merge",        dest="merge",          action="store_true", default=False,                  help="merge the 2 plots in 1 [format: boolean] [default: false]")
+    parser.add_option("-D", "--directory",    dest="working_Dir",    action="store",      default=os.getcwd(), type="string", help="the directory to store the results [format: directory]")
+    parser.add_option("-v", "--verbosity",    dest="verbosity",      action="store",      default=1,       type="int",    help="trace level [format: int]")
+    (options, args) = parser.parse_args()
 
-	colors[1]  = options.plusColor
-	colors[-1] = options.minusColor
-	colors[0]  = options.sumColor
-	colorLine  = options.lineColor
+    colors[1]  = options.plusColor
+    colors[-1] = options.minusColor
+    colors[0]  = options.sumColor
+    colorLine  = options.lineColor
 
-	pp = PlotParser(options.verbosity)
-	pp.addInput(0, options.inputFileName1, options.inputFormat1)
-	pp.addInput(1, options.inputFileName2, options.inputFormat2)
-	pp.addSequence(options.inputSequence)
-	pp.setOutput(options.outputFileName if os.path.isabs(options.outputFileName) else os.path.join(options.working_Dir, options.outputFileName))
-	pp.setPlotSize(options.width, options.height)
-	pp.setLabels(options.xLabel, options.yLabel)
-	pp.setTitle(options.title)
-	pp.setMerge(options.merge)
-	pp.start()
+    pp = PlotParser(options.verbosity)
+    pp.addInput(0, options.inputFileName1, options.inputFormat1)
+    pp.addInput(1, options.inputFileName2, options.inputFormat2)
+    pp.addSequence(options.inputSequence)
+    pp.setOutput(options.outputFileName if os.path.isabs(options.outputFileName) else os.path.join(options.working_Dirpath, options.outputFileName))
+    pp.setPlotSize(options.width, options.height)
+    pp.setLabels(options.xLabel, options.yLabel)
+    pp.setTitle(options.title)
+    pp.setMerge(options.merge)
+    pp.start()
 
+