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