annotate SMART/Java/Python/plotCoverage.py @ 62:8c42a6d7ffd4

Added simple test BED file.

#! /usr/bin/env python
#
# Copyright INRA-URGI 2009-2010
# 
# This software is governed by the CeCILL license under French law and
# abiding by the rules of distribution of free software. You can use,
# modify and/ or redistribute the software under the terms of the CeCILL
# license as circulated by CEA, CNRS and INRIA at the following URL
# "http://www.cecill.info".
# 
# As a counterpart to the access to the source code and rights to copy,
# modify and redistribute granted by the license, users are provided only
# with a limited warranty and the software's author, the holder of the
# economic rights, and the successive licensors have only limited
# liability.
# 
# In this respect, the user's attention is drawn to the risks associated
# with loading, using, modifying and/or developing or reproducing the
# software by the user in light of its specific status of free software,
# that may mean that it is complicated to manipulate, and that also
# therefore means that it is reserved for developers and experienced
# professionals having in-depth computer knowledge. Users are therefore
# encouraged to load and test the software's suitability as regards their
# requirements in conditions enabling the security of their systems and/or
# data to be ensured and, more generally, to use and operate it in the
# same conditions as regards security.
# 
# The fact that you are presently reading this means that you have had
# knowledge of the CeCILL license and that you accept its terms.
#
import os, os.path, subprocess, glob, random
from optparse import OptionParser
from SMART.Java.Python.structure.Interval import Interval
from SMART.Java.Python.structure.Transcript import Transcript
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

strands = [-1, 1]
colors  = {-1: "blue", 1: "red", 0: "black"}
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)


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 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


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)


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 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 setTitle(self, title):
        self.title = title + " " + self.name

    def setPlotSize(self, width, height):
        self.width  = width
        self.height = height

    def setLabels(self, xLabel, yLabel):
        self.xLabel = xLabel
        self.yLabel = yLabel

    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 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 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 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 __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
        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 setOutput(self, fileName):
        self.outputFileName = fileName

    def setPlotSize(self, width, height):
        self.width  = width
        self.height = height

    def setLabels(self, xLabel, yLabel):
        self.xLabel = xLabel
        self.yLabel = yLabel

    def setTitle(self, title):
        self.title = title

    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 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 initialize(self):
        if self.sequenceParser == None:
            self.initializeDataFromTranscripts()
        else:
            self.initializeDataFromSequences()

    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 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 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 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 plot(self):
        if self.sequenceParser == None:
            self.plot2TranscriptFiles()
        else:
            self.plot1TranscriptFile()

    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]"

    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

    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()