18
+ − 1 #! /usr/bin/env python
+ − 2 #This program is a wrapp for CompareOverlapping.py.
+ − 3 import os, sys, tarfile, optparse
+ − 4 from commons.core.launcher.Launcher import Launcher
+ − 5 from commons.core.sql.TableJobAdaptatorFactory import TableJobAdaptatorFactory
+ − 6 from optparse import OptionParser
+ − 7 from commons.core.utils.FileUtils import FileUtils
+ − 8 from commons.core.parsing.ParserChooser import ParserChooser
+ − 9 from SMART.Java.Python.structure.TranscriptList import TranscriptList
+ − 10 from commons.core.writer.WriterChooser import WriterChooser
+ − 11
+ − 12 def stop_err( msg ):
+ − 13 sys.stderr.write( "%s\n" % msg )
+ − 14 sys.exit()
+ − 15
+ − 16 def toTar(tarFileName, overlapOutputNames):
+ − 17 dir = os.path.dirname(tarFileName)
+ − 18 tfile = tarfile.open(tarFileName + ".tmp.tar", "w")
+ − 19 currentPath = os.getcwd()
+ − 20 os.chdir(dir)
+ − 21 for file in overlapOutputNames:
+ − 22 relativeFileName = os.path.basename(file)
+ − 23 tfile.add(relativeFileName)
+ − 24 os.system("mv %s %s" % (tarFileName + ".tmp.tar", tarFileName))
+ − 25 tfile.close()
+ − 26 os.chdir(currentPath)
+ − 27
+ − 28 def _createCompareOverlappingCmd(iLauncher, options, inputFileName, annotationFile, overlapOutputName):
+ − 29 lArgs = []
+ − 30 lArgs.append("-i %s" % annotationFile)
+ − 31 lArgs.append("-f %s" % options.format1)
+ − 32 lArgs.append("-j %s" % inputFileName)
+ − 33 lArgs.append("-g %s" % options.format2)
+ − 34 lArgs.append("-o %s" % overlapOutputName)
+ − 35 if options.notOverlapping:
+ − 36 lArgs.append("-O")
+ − 37 if options.exclude:
+ − 38 lArgs.append("-x")
+ − 39 if options.distance != None:
+ − 40 lArgs.append("-d %s" % options.distance)
+ − 41 return(iLauncher.getSystemCommand("python %s/SMART/Java/Python/CompareOverlappingSmallQuery.py" % os.environ["REPET_PATH"], lArgs))
+ − 42
+ − 43 def _map(iLauncher, cmd, cmdStart, cmdFinish ):
+ − 44 lCmds = []
+ − 45 lCmds.append(cmd)
+ − 46 lCmdStart = []
+ − 47 lCmdStart.append(cmdStart)
+ − 48 lCmdFinish = []
+ − 49 lCmdFinish.append(cmdFinish)
+ − 50 return(iLauncher.prepareCommands_withoutIndentation(lCmds, lCmdStart, lCmdFinish))
+ − 51
+ − 52 def split(fileName, nbOfSeqPerBatch):
+ − 53 filePrefix, fileExt = os.path.splitext(os.path.basename(fileName))
+ − 54 resDir = os.path.dirname(fileName)
+ − 55 lInputName = []
+ − 56 fileNb = 1
+ − 57 SeqNb = 0
+ − 58 outFileName = "%s/%s-%s%s" %(resDir, filePrefix, fileNb, fileExt)
+ − 59 lInputName.append(outFileName)
+ − 60 outFile = open(outFileName, "w")
+ − 61 f = open(fileName, "r")
+ − 62 line = f.readline()
+ − 63 previousRefName = ""
+ − 64 while line != "":
+ − 65 if not line.startswith('@SQ'):
+ − 66 if SeqNb == nbOfSeqPerBatch:
+ − 67 SeqNb = 0
+ − 68 fileNb += 1
+ − 69 outFile.close()
+ − 70 outFileName = "%s/%s-%s%s" %(resDir, filePrefix, fileNb, fileExt)
+ − 71 lInputName.append(outFileName)
+ − 72 outFile = open(outFileName, "w")
+ − 73 refName = line.split("\t")[2]
+ − 74 if previousRefName != refName:
+ − 75 SeqNb += 1
+ − 76 outFile.write(line)
+ − 77 else:
+ − 78 previousRefName = refName
+ − 79 outFile.write(line)
+ − 80 line = f.readline()
+ − 81 return lInputName
+ − 82
+ − 83 def join(dCutOut2Out, options):
+ − 84 chooser = ParserChooser()
+ − 85 chooser.findFormat("gtf")
+ − 86 gtfParser = chooser.getParser(options.inputFileName1)
+ − 87 ref = {}
+ − 88 for transcript in gtfParser.getIterator():
+ − 89 ref[transcript.getTagValue("ID")] = transcript
+ − 90 for key in dCutOut2Out.keys():
+ − 91 writerChooser = WriterChooser()
+ − 92 writerChooser.findFormat("gff3")
+ − 93 for inputFile in dCutOut2Out[key]:
+ − 94 chooser = ParserChooser()
+ − 95 chooser.findFormat("gff")
+ − 96 gffParser = chooser.getParser(inputFile)
+ − 97 for transcript in gffParser.getIterator():
+ − 98 finalTranscript = ref[transcript.getTagValue("ID")]
+ − 99 if finalTranscript.getTagValue("nbOverlaps"):
+ − 100 nbOverlap = int(finalTranscript.getTagValue("nbOverlaps")) + int(transcript.getTagValue("nbOverlaps"))
+ − 101 finalTranscript.setTagValue("nbOverlaps", nbOverlap)
+ − 102 else:
+ − 103 finalTranscript.setTagValue("nbOverlaps", transcript.getTagValue("nbOverlaps"))
+ − 104
+ − 105 if finalTranscript.getTagValue("overlapsWith") and transcript.getTagValue("overlapsWith") != None:
+ − 106 overlapName = "--".join([finalTranscript.getTagValue("overlapsWith"), transcript.getTagValue("overlapsWith")])
+ − 107 finalTranscript.setTagValue("overlapsWith", overlapName)
+ − 108 else:
+ − 109 if transcript.getTagValue("overlapsWith") != None:
+ − 110 finalTranscript.setTagValue("overlapsWith", transcript.getTagValue("overlapsWith"))
+ − 111
+ − 112 gffWriter = writerChooser.getWriter(key)
+ − 113 gffWriter.setTitle("S-MART")
+ − 114 for transcript in ref.values():
+ − 115 gffWriter.addTranscript(transcript)
+ − 116 gffWriter.write()
+ − 117 transcript.deleteTag("nbOverlaps")
+ − 118 transcript.deleteTag("overlapsWith")
+ − 119 gffWriter.close()
+ − 120
+ − 121 def __main__():
+ − 122 description = "Compare Overlapping wrapp script: Get the a list of data which overlap with a reference set. [Category: Data Comparison]"
+ − 123 parser = OptionParser(description = description)
+ − 124 parser.add_option("-i", "--input1", dest="inputFileName1", action="store", type="string", help="input file 1 (for annotation) [compulsory] [format: file in transcript format given by -f]")
+ − 125 parser.add_option("-f", "--format1", dest="format1", action="store", type="string", help="format of file 1 [compulsory] [format: transcript file format]")
+ − 126 parser.add_option("", "--inputTxt", dest="inputTxt", action="store", type="string", help="input, a txt file for a list of input reads files. Should identify all reads files format, given by -g [compulsory]")
+ − 127 #parser.add_option("-j", "--input2", dest="inputFileName2", action="store", default="inputRead", type="string", help="input file 2 [compulsory] [format: file in transcript format given by -g]")
+ − 128 parser.add_option("-g", "--format2", dest="format2", action="store", type="string", help="format of file 2 [compulsory] [format: transcript file format]")
+ − 129 #parser.add_option("-o", "--output", dest="output", action="store", default=None, type="string", help="output file [compulsory] [format: output file in GFF3 format]")
+ − 130 parser.add_option("-S", "--start1", dest="start1", action="store", default=None, type="int", help="only consider the n first nucleotides of the transcripts in file 1 (do not use it with -U) [format: int]")
+ − 131 parser.add_option("-s", "--start2", dest="start2", action="store", default=None, type="int", help="only consider the n first nucleotides of the transcripts in file 2 (do not use it with -u) [format: int]")
+ − 132 parser.add_option("-U", "--end1", dest="end1", action="store", default=None, type="int", help="only consider the n last nucleotides of the transcripts in file 1 (do not use it with -S) [format: int]")
+ − 133 parser.add_option("-u", "--end2", dest="end2", action="store", default=None, type="int", help="only consider the n last nucleotides of the transcripts in file 2 (do not use it with -s) [format: int]")
+ − 134 parser.add_option("-t", "--intron", dest="introns", action="store_true", default=False, help="also report introns [format: bool] [default: false]")
+ − 135 parser.add_option("-E", "--5primeExtension1", dest="fivePrime1", action="store", default=None, type="int", help="extension towards 5' in file 1 [format: int]")
+ − 136 parser.add_option("-e", "--5primeExtension2", dest="fivePrime2", action="store", default=None, type="int", help="extension towards 5' in file 2 [format: int]")
+ − 137 parser.add_option("-N", "--3primeExtension1", dest="threePrime1", action="store", default=None, type="int", help="extension towards 3' in file 1 [format: int]")
+ − 138 parser.add_option("-n", "--3primeExtension2", dest="threePrime2", action="store", default=None, type="int", help="extension towards 3' in file 2 [format: int]")
+ − 139 parser.add_option("-c", "--colinear", dest="colinear", action="store_true", default=False, help="colinear only [format: bool] [default: false]")
+ − 140 parser.add_option("-a", "--antisense", dest="antisense", action="store_true", default=False, help="antisense only [format: bool] [default: false]")
+ − 141 parser.add_option("-d", "--distance", dest="distance", action="store", default=None, type="int", help="accept some distance between query and reference [format: int]")
+ − 142 parser.add_option("-k", "--included", dest="included", action="store_true", default=False, help="keep only elements from file 1 which are included in an element of file 2 [format: bool] [default: false]")
+ − 143 parser.add_option("-K", "--including", dest="including", action="store_true", default=False, help="keep only elements from file 2 which are included in an element of file 1 [format: bool] [default: false]")
+ − 144 parser.add_option("-m", "--minOverlap", dest="minOverlap", action="store", default=None, type="int", help="minimum number of nucleotides overlapping to declare an overlap [format: int] [default: 1]")
+ − 145 parser.add_option("-p", "--pcOverlap", dest="pcOverlap", action="store", default=None, type="int", help="minimum percentage of nucleotides to overlap to declare an overlap [format: int]")
+ − 146 parser.add_option("-O", "--notOverlapping", dest="notOverlapping", action="store_true", default=False, help="also output not overlapping data [format: bool] [default: false]")
+ − 147 parser.add_option("-x", "--exclude", dest="exclude", action="store_true", default=False, help="invert the match [format: bool] [default: false]")
+ − 148 parser.add_option("-v", "--verbosity", dest="verbosity", action="store", default=1, type="int", help="trace level [format: int]")
+ − 149 parser.add_option('', '--tar', dest='outputTar', default=None, help='output all SAM results in a tar file.' )
+ − 150 parser.add_option( '', '--outTxt', dest='outTxtFile', help='The output list of results files on txt format.[compulsory]' )
+ − 151 (options, args) = parser.parse_args()
+ − 152
+ − 153
+ − 154 #Parse the input txt file and read a list of BAM files.
+ − 155 file = open(options.inputTxt, "r")
+ − 156 lines = file.readlines()
+ − 157 inputFileNames = []
+ − 158 overlapOutputNames = []
+ − 159 outputName = options.outTxtFile
+ − 160 resDirName = os.path.dirname(outputName) + "/"
+ − 161 #Write output txt file and define all output sam file names.
+ − 162 out = open(outputName, "w")
+ − 163 for line in lines:
+ − 164 tab = line.split()
+ − 165 inputFileNames.append(tab[1])
+ − 166 overlapOutName = resDirName + tab[0] + '_overlapOut.gff3'
+ − 167 overlapOutputNames.append(overlapOutName)
+ − 168 out.write(tab[0] + '\t' + overlapOutName + '\n')
+ − 169 file.close()
+ − 170 out.close()
+ − 171
+ − 172 #Launch on nodes
+ − 173 acronym = "compareOverlapping"
+ − 174 jobdb = TableJobAdaptatorFactory.createJobInstance()
+ − 175 iLauncher = Launcher(jobdb, os.getcwd(), "", "", os.getcwd(), os.getcwd(), "jobs", "test", acronym, acronym, False, True)
+ − 176
+ − 177
+ − 178
+ − 179
+ − 180 #construction the commandes for each input file
+ − 181 lCmdsTuples = []
+ − 182 dCutOut2Out = {}
+ − 183 lAllFile2remove = []
+ − 184 for i in range(len(inputFileNames)):
+ − 185 lCutInputFile = split(inputFileNames[i], 20000)
+ − 186 lAllFile2remove.extend(lCutInputFile)
+ − 187 lCutOutput = []
+ − 188 for cutInput in lCutInputFile:
+ − 189 cutOutput = "%s_out" % cutInput
+ − 190 lCutOutput.append(cutOutput)
+ − 191 lAllFile2remove.extend(lCutOutput)
+ − 192 cmd2Launch = _createCompareOverlappingCmd(iLauncher, options, cutInput, options.inputFileName1, cutOutput)
+ − 193 lCmdsTuples.append(_map(iLauncher, cmd2Launch, "", ""))
+ − 194 chooser = ParserChooser()
+ − 195 chooser.findFormat(options.format2)
+ − 196 dCutOut2Out[overlapOutputNames[i]] = lCutOutput
+ − 197 iLauncher.runLauncherForMultipleJobs(acronym, lCmdsTuples, True)
+ − 198
+ − 199 join(dCutOut2Out, options)
+ − 200 FileUtils.removeFilesFromListIfExist(lAllFile2remove)
+ − 201
+ − 202 if options.outputTar != None:
+ − 203 toTar(options.outputTar, overlapOutputNames)
+ − 204
+ − 205 if __name__=="__main__": __main__()
