Mercurial > repos > yufei-luo > s_mart

comparison SMART/DiffExpAnal/compareOverlapping_parallel_unSQL.py @ 18:94ab73e8a190

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author m-zytnicki
date Mon, 29 Apr 2013 03:20:15 -0400
parents
children
comparison
equal deleted inserted replaced
17:b0e8584489e6 18:94ab73e8a190
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__()