|
36
|
1 # Copyright INRA (Institut National de la Recherche Agronomique)
|
|
|
2 # http://www.inra.fr
|
|
|
3 # http://urgi.versailles.inra.fr
|
|
|
4 #
|
|
|
5 # This software is governed by the CeCILL license under French law and
|
|
|
6 # abiding by the rules of distribution of free software. You can use,
|
|
|
7 # modify and/ or redistribute the software under the terms of the CeCILL
|
|
|
8 # license as circulated by CEA, CNRS and INRIA at the following URL
|
|
|
9 # "http://www.cecill.info".
|
|
|
10 #
|
|
|
11 # As a counterpart to the access to the source code and rights to copy,
|
|
|
12 # modify and redistribute granted by the license, users are provided only
|
|
|
13 # with a limited warranty and the software's author, the holder of the
|
|
|
14 # economic rights, and the successive licensors have only limited
|
|
|
15 # liability.
|
|
|
16 #
|
|
|
17 # In this respect, the user's attention is drawn to the risks associated
|
|
|
18 # with loading, using, modifying and/or developing or reproducing the
|
|
|
19 # software by the user in light of its specific status of free software,
|
|
|
20 # that may mean that it is complicated to manipulate, and that also
|
|
|
21 # therefore means that it is reserved for developers and experienced
|
|
|
22 # professionals having in-depth computer knowledge. Users are therefore
|
|
|
23 # encouraged to load and test the software's suitability as regards their
|
|
|
24 # requirements in conditions enabling the security of their systems and/or
|
|
|
25 # data to be ensured and, more generally, to use and operate it in the
|
|
|
26 # same conditions as regards security.
|
|
|
27 #
|
|
|
28 # The fact that you are presently reading this means that you have had
|
|
|
29 # knowledge of the CeCILL license and that you accept its terms.
|
|
|
30
|
|
|
31
|
|
|
32 import os
|
|
|
33 import sys
|
|
|
34 import string
|
|
|
35 import math
|
|
|
36 import shutil
|
|
|
37 import re
|
|
|
38 import glob
|
|
|
39 from operator import itemgetter
|
|
|
40 from commons.core.seq.BioseqDB import BioseqDB
|
|
|
41 from commons.core.seq.Bioseq import Bioseq
|
|
|
42 from commons.core.coord.MapUtils import MapUtils
|
|
|
43 from commons.core.coord.Range import Range
|
|
|
44 from commons.core.checker.CheckerUtils import CheckerUtils
|
|
|
45 from commons.core.launcher.LauncherUtils import LauncherUtils
|
|
|
46 from commons.core.coord.ConvCoord import ConvCoord
|
|
|
47 from commons.core.parsing.FastaParser import FastaParser
|
|
|
48
|
|
|
49
|
|
|
50 ## Static methods for fasta file manipulation
|
|
|
51 #
|
|
|
52 class FastaUtils( object ):
|
|
|
53
|
|
|
54 ## Count the number of sequences in the input fasta file
|
|
|
55 #
|
|
|
56 # @param inFile name of the input fasta file
|
|
|
57 #
|
|
|
58 # @return integer number of sequences in the input fasta file
|
|
|
59 #
|
|
|
60 @staticmethod
|
|
|
61 def dbSize( inFile ):
|
|
|
62 nbSeq = 0
|
|
|
63 inFileHandler = open( inFile, "r" )
|
|
|
64 line = inFileHandler.readline()
|
|
|
65 while line:
|
|
|
66 if line[0] == ">":
|
|
|
67 nbSeq = nbSeq + 1
|
|
|
68 line = inFileHandler.readline()
|
|
|
69 inFileHandler.close()
|
|
|
70
|
|
|
71 return nbSeq
|
|
|
72
|
|
|
73
|
|
|
74 ## Compute the cumulative sequence length in the input fasta file
|
|
|
75 #
|
|
|
76 # @param inFile handler of the input fasta file
|
|
|
77 #
|
|
|
78 @staticmethod
|
|
|
79 def dbCumLength( inFile ):
|
|
|
80 cumLength = 0
|
|
|
81 line = inFile.readline()
|
|
|
82 while line:
|
|
|
83 if line[0] != ">":
|
|
|
84 cumLength += len(string.rstrip(line))
|
|
|
85 line = inFile.readline()
|
|
|
86
|
|
|
87 return cumLength
|
|
|
88
|
|
|
89
|
|
|
90 ## Return a list with the length of each sequence in the input fasta file
|
|
|
91 #
|
|
|
92 # @param inFile string name of the input fasta file
|
|
|
93 #
|
|
|
94 @staticmethod
|
|
|
95 def dbLengths( inFile ):
|
|
|
96 lLengths = []
|
|
|
97 inFileHandler = open( inFile, "r" )
|
|
|
98 currentLength = 0
|
|
|
99 line = inFileHandler.readline()
|
|
|
100 while line:
|
|
|
101 if line[0] == ">":
|
|
|
102 if currentLength != 0:
|
|
|
103 lLengths.append( currentLength )
|
|
|
104 currentLength = 0
|
|
|
105 else:
|
|
|
106 currentLength += len(line[:-1])
|
|
|
107 line = inFileHandler.readline()
|
|
|
108 lLengths.append( currentLength )
|
|
|
109 inFileHandler.close()
|
|
|
110 return lLengths
|
|
|
111
|
|
|
112
|
|
|
113 ## Retrieve the sequence headers present in the input fasta file
|
|
|
114 #
|
|
|
115 # @param inFile string name of the input fasta file
|
|
|
116 # @param verbose integer level of verbosity
|
|
|
117 #
|
|
|
118 # @return list of sequence headers
|
|
|
119 #
|
|
|
120 @staticmethod
|
|
|
121 def dbHeaders( inFile, verbose=0 ):
|
|
|
122 lHeaders = []
|
|
|
123
|
|
|
124 inFileHandler = open( inFile, "r" )
|
|
|
125 line = inFileHandler.readline()
|
|
|
126 while line:
|
|
|
127 if line[0] == ">":
|
|
|
128 lHeaders.append( string.rstrip(line[1:]) )
|
|
|
129 if verbose > 0:
|
|
|
130 print string.rstrip(line[1:])
|
|
|
131 line = inFileHandler.readline()
|
|
|
132 inFileHandler.close()
|
|
|
133
|
|
|
134 return lHeaders
|
|
|
135
|
|
|
136
|
|
|
137 ## Cut a data bank into chunks according to the input parameters
|
|
|
138 # If a sequence is shorter than the threshold, it is only renamed (not cut)
|
|
|
139 #
|
|
|
140 # @param inFileName string name of the input fasta file
|
|
|
141 # @param chkLgth string chunk length (in bp, default=200000)
|
|
|
142 # @param chkOver string chunk overlap (in bp, default=10000)
|
|
|
143 # @param wordN string N stretch word length (default=11, 0 for no detection)
|
|
|
144 # @param outFilePrefix string prefix of the output files (default=inFileName + '_chunks.fa' and '_chunks.map')
|
|
|
145 # @param clean boolean remove 'cut' and 'Nstretch' files
|
|
|
146 # @param verbose integer (default = 0)
|
|
|
147 #
|
|
|
148 @staticmethod
|
|
|
149 def dbChunks( inFileName, chkLgth="200000", chkOver="10000", wordN="11", outFilePrefix="", clean=False, verbose=0 ):
|
|
|
150 nbSeq = FastaUtils.dbSize( inFileName )
|
|
|
151 if verbose > 0:
|
|
|
152 print "cut the %i input sequences with cutterDB..." % ( nbSeq )
|
|
|
153 sys.stdout.flush()
|
|
|
154
|
|
|
155 prg = "cutterDB"
|
|
|
156 cmd = prg
|
|
|
157 cmd += " -l %s" % ( chkLgth )
|
|
|
158 cmd += " -o %s" %( chkOver )
|
|
|
159 cmd += " -w %s" % ( wordN )
|
|
|
160 cmd += " %s" % ( inFileName )
|
|
|
161 returnStatus = os.system( cmd )
|
|
|
162 if returnStatus != 0:
|
|
|
163 msg = "ERROR: '%s' returned '%i'" % ( prg, returnStatus )
|
|
|
164 sys.stderr.write( "%s\n" % ( msg ) )
|
|
|
165 sys.exit(1)
|
|
|
166
|
|
|
167 nbChunks = FastaUtils.dbSize( "%s_cut" % ( inFileName ) )
|
|
|
168 if verbose > 0:
|
|
|
169 print "done (%i chunks)" % ( nbChunks )
|
|
|
170 sys.stdout.flush()
|
|
|
171
|
|
|
172 if verbose > 0:
|
|
|
173 print "rename the headers..."
|
|
|
174 sys.stdout.flush()
|
|
|
175
|
|
|
176 if outFilePrefix == "":
|
|
|
177 outFastaName = inFileName + "_chunks.fa"
|
|
|
178 outMapName = inFileName + "_chunks.map"
|
|
|
179 else:
|
|
|
180 outFastaName = outFilePrefix + ".fa"
|
|
|
181 outMapName = outFilePrefix + ".map"
|
|
|
182
|
|
|
183 inFile = open( "%s_cut" % ( inFileName ), "r" )
|
|
|
184 line = inFile.readline()
|
|
|
185
|
|
|
186 outFasta = open( outFastaName, "w" )
|
|
|
187 outMap = open( outMapName, "w" )
|
|
|
188
|
|
|
189 # read line after line (no need for big RAM) and change the sequence headers
|
|
|
190 while line:
|
|
|
191
|
|
|
192 if line[0] == ">":
|
|
|
193 if verbose > 1:
|
|
|
194 print "rename '%s'" % ( line[:-1] ); sys.stdout.flush()
|
|
|
195 data = line[:-1].split(" ")
|
|
|
196 seqID = data[0].split(">")[1]
|
|
|
197 newHeader = "chunk%s" % ( str(seqID).zfill( len(str(nbChunks)) ) )
|
|
|
198 oldHeader = data[2]
|
|
|
199 seqStart = data[4].split("..")[0]
|
|
|
200 seqEnd = data[4].split("..")[1]
|
|
|
201 outMap.write( "%s\t%s\t%s\t%s\n" % ( newHeader, oldHeader, seqStart, seqEnd ) )
|
|
|
202 outFasta.write( ">%s\n" % ( newHeader ) )
|
|
|
203
|
|
|
204 else:
|
|
|
205 outFasta.write( line.upper() )
|
|
|
206
|
|
|
207 line = inFile.readline()
|
|
|
208
|
|
|
209 inFile.close()
|
|
|
210 outFasta.close()
|
|
|
211 outMap.close()
|
|
|
212
|
|
|
213 if clean == True:
|
|
|
214 os.remove(inFileName + "_cut")
|
|
|
215 os.remove(inFileName + ".Nstretch.map")
|
|
|
216
|
|
|
217
|
|
|
218 ## Split the input fasta file in several output files
|
|
|
219 #
|
|
|
220 # @param inFile string name of the input fasta file
|
|
|
221 # @param nbSeqPerBatch integer number of sequences per output file
|
|
|
222 # @param newDir boolean put the sequences in a new directory called 'batches'
|
|
|
223 # @param useSeqHeader boolean use sequence header (only if 'nbSeqPerBatch=1')
|
|
|
224 # @param prefix prefix in output file name
|
|
|
225 # @param verbose integer verbosity level (default = 0)
|
|
|
226 #
|
|
|
227 @staticmethod
|
|
|
228 def dbSplit( inFile, nbSeqPerBatch, newDir, useSeqHeader=False, prefix="batch", verbose=0 ):
|
|
|
229 if not os.path.exists( inFile ):
|
|
|
230 msg = "ERROR: file '%s' doesn't exist" % ( inFile )
|
|
|
231 sys.stderr.write( "%s\n" % ( msg ) )
|
|
|
232 sys.exit(1)
|
|
|
233
|
|
|
234 nbSeq = FastaUtils.dbSize( inFile )
|
|
|
235
|
|
|
236 nbBatches = int( math.ceil( nbSeq / float(nbSeqPerBatch) ) )
|
|
|
237 if verbose > 0:
|
|
|
238 print "save the %i input sequences into %i batches" % ( nbSeq, nbBatches )
|
|
|
239 sys.stdout.flush()
|
|
|
240
|
|
|
241 if nbSeqPerBatch > 1 and useSeqHeader:
|
|
|
242 useSeqHeader = False
|
|
|
243
|
|
|
244 if newDir == True:
|
|
|
245 if os.path.exists( "batches" ):
|
|
|
246 shutil.rmtree( "batches" )
|
|
|
247 os.mkdir( "batches" )
|
|
|
248 os.chdir( "batches" )
|
|
|
249 os.system( "ln -s ../%s ." % ( inFile ) )
|
|
|
250
|
|
|
251 inFileHandler = open( inFile, "r" )
|
|
|
252 inFileHandler.seek( 0, 0 )
|
|
|
253 countBatch = 0
|
|
|
254 countSeq = 0
|
|
|
255 line = inFileHandler.readline()
|
|
|
256 while line:
|
|
|
257 if line == "":
|
|
|
258 break
|
|
|
259 if line[0] == ">":
|
|
|
260 countSeq += 1
|
|
|
261 if nbSeqPerBatch == 1 or countSeq % nbSeqPerBatch == 1:
|
|
|
262 if "outFile" in locals():
|
|
|
263 outFile.close()
|
|
|
264 countBatch += 1
|
|
|
265 if nbSeqPerBatch == 1 and useSeqHeader:
|
|
|
266 outFileName = "%s.fa" % ( line[1:-1].replace(" ","_") )
|
|
|
267 else:
|
|
|
268 outFileName = "%s_%s.fa" % ( prefix, str(countBatch).zfill( len(str(nbBatches)) ) )
|
|
|
269 outFile = open( outFileName, "w" )
|
|
|
270 if verbose > 1:
|
|
|
271 print "saving seq '%s' in file '%s'..." % ( line[1:40][:-1], outFileName )
|
|
|
272 sys.stdout.flush()
|
|
|
273 outFile.write( line )
|
|
|
274 line = inFileHandler.readline()
|
|
|
275 inFileHandler.close()
|
|
|
276
|
|
|
277 if newDir == True:
|
|
|
278 os.remove( os.path.basename( inFile ) )
|
|
|
279 os.chdir( ".." )
|
|
|
280
|
|
|
281
|
|
|
282 ## Split the input fasta file in several output files
|
|
|
283 #
|
|
|
284 # @param inFileName string name of the input fasta file
|
|
|
285 # @param maxSize integer max cumulative length for each output file
|
|
|
286 #
|
|
|
287 @staticmethod
|
|
|
288 def splitFastaFileInBatches(inFileName, maxSize = 200000):
|
|
|
289 iBioseqDB = BioseqDB(inFileName)
|
|
|
290 lHeadersSizeTuples = []
|
|
|
291 for iBioseq in iBioseqDB.db:
|
|
|
292 lHeadersSizeTuples.append((iBioseq.getHeader(), iBioseq.getLength()))
|
|
|
293
|
|
|
294 lHeadersList = LauncherUtils.createHomogeneousSizeList(lHeadersSizeTuples, maxSize)
|
|
|
295 os.mkdir("batches")
|
|
|
296 os.chdir("batches")
|
|
|
297
|
|
|
298 iterator = 0
|
|
|
299 for lHeader in lHeadersList :
|
|
|
300 iterator += 1
|
|
|
301 with open("batch_%s.fa" % iterator, 'w') as f :
|
|
|
302 for header in lHeader :
|
|
|
303 iBioseq = iBioseqDB.fetch(header)
|
|
|
304 iBioseq.write(f)
|
|
|
305 os.chdir("..")
|
|
|
306
|
|
|
307
|
|
|
308 ## Split the input fasta file in several output files according to their cluster identifier
|
|
|
309 #
|
|
|
310 # @param inFileName string name of the input fasta file
|
|
|
311 # @param clusteringMethod string name of the clustering method (Grouper, Recon, Piler, Blastclust)
|
|
|
312 # @param simplifyHeader boolean simplify the headers
|
|
|
313 # @param createDir boolean put the sequences in different directories
|
|
|
314 # @param outPrefix string prefix of the output files (default='seqCluster')
|
|
|
315 # @param verbose integer (default = 0)
|
|
|
316 #
|
|
|
317 @staticmethod
|
|
|
318 def splitSeqPerCluster( inFileName, clusteringMethod, simplifyHeader, createDir, outPrefix="seqCluster", verbose=0 ):
|
|
|
319 if not os.path.exists( inFileName ):
|
|
|
320 print "ERROR: %s doesn't exist" % ( inFileName )
|
|
|
321 sys.exit(1)
|
|
|
322
|
|
|
323 inFile = open( inFileName, "r" )
|
|
|
324
|
|
|
325 line = inFile.readline()
|
|
|
326 if line:
|
|
|
327 name = line.split(" ")[0]
|
|
|
328 if "Cluster" in name:
|
|
|
329 clusterID = name.split("Cluster")[1].split("Mb")[0]
|
|
|
330 seqID = name.split("Mb")[1]
|
|
|
331 else:
|
|
|
332 clusterID = name.split("Cl")[0].split("Gr")[1] # the notion of 'group' in Grouper corresponds to 'cluster' in Piler, Recon and Blastclust
|
|
|
333 if "Q" in name.split("Gr")[0]:
|
|
|
334 seqID = name.split("Gr")[0].split("MbQ")[1]
|
|
|
335 elif "S" in name:
|
|
|
336 seqID = name.split("Gr")[0].split("MbS")[1]
|
|
|
337 sClusterIDs = set( [ clusterID ] )
|
|
|
338 if simplifyHeader == True:
|
|
|
339 header = "%s_Cluster%s_Seq%s" % ( clusteringMethod, clusterID, seqID )
|
|
|
340 else:
|
|
|
341 header = line[1:-1]
|
|
|
342 if createDir == True:
|
|
|
343 if not os.path.exists( "%s_cluster_%s" % ( inFileName, clusterID ) ):
|
|
|
344 os.mkdir( "%s_cluster_%s" % ( inFileName, clusterID ) )
|
|
|
345 os.chdir( "%s_cluster_%s" % ( inFileName, clusterID ) )
|
|
|
346 outFileName = "%s%s.fa" % ( outPrefix, clusterID )
|
|
|
347 outFile = open( outFileName, "w" )
|
|
|
348 outFile.write( ">%s\n" % ( header ) )
|
|
|
349 prevClusterID = clusterID
|
|
|
350
|
|
|
351 line = inFile.readline()
|
|
|
352 while line:
|
|
|
353 if line[0] == ">":
|
|
|
354 name = line.split(" ")[0]
|
|
|
355 if "Cluster" in name:
|
|
|
356 clusterID = name.split("Cluster")[1].split("Mb")[0]
|
|
|
357 seqID = name.split("Mb")[1]
|
|
|
358 else:
|
|
|
359 clusterID = name.split("Cl")[0].split("Gr")[1]
|
|
|
360 if "Q" in name.split("Gr")[0]:
|
|
|
361 seqID = name.split("Gr")[0].split("MbQ")[1]
|
|
|
362 elif "S" in name:
|
|
|
363 seqID = name.split("Gr")[0].split("MbS")[1]
|
|
|
364
|
|
|
365 if clusterID != prevClusterID:
|
|
|
366 outFile.close()
|
|
|
367
|
|
|
368 if simplifyHeader == True:
|
|
|
369 header = "%s_Cluster%s_Seq%s" % ( clusteringMethod, clusterID, seqID )
|
|
|
370 else:
|
|
|
371 header = line[1:-1]
|
|
|
372
|
|
|
373 if createDir == True:
|
|
|
374 os.chdir( ".." )
|
|
|
375 if not os.path.exists( "%s_cluster_%s" % ( inFileName, clusterID ) ):
|
|
|
376 os.mkdir( "%s_cluster_%s" % ( inFileName, clusterID ) )
|
|
|
377 os.chdir( "%s_cluster_%s" % ( inFileName, clusterID ) )
|
|
|
378
|
|
|
379 outFileName = "%s%s.fa" % ( outPrefix, clusterID )
|
|
|
380 if not os.path.exists( outFileName ):
|
|
|
381 outFile = open( outFileName, "w" )
|
|
|
382 else:
|
|
|
383 if clusterID != prevClusterID:
|
|
|
384 outFile.close()
|
|
|
385 outFile = open( outFileName, "a" )
|
|
|
386 outFile.write( ">%s\n" % ( header ) )
|
|
|
387 prevClusterID = clusterID
|
|
|
388 sClusterIDs.add( clusterID )
|
|
|
389
|
|
|
390 else:
|
|
|
391 outFile.write( line )
|
|
|
392
|
|
|
393 line = inFile.readline()
|
|
|
394
|
|
|
395 outFile.close()
|
|
|
396 if verbose > 0:
|
|
|
397 print "number of clusters: %i" % ( len(sClusterIDs) ); sys.stdout.flush()
|
|
|
398
|
|
|
399 if createDir == True:
|
|
|
400 os.chdir("..")
|
|
|
401 else:
|
|
|
402 print "WARNING: empty input file - no cluster found"; sys.stdout.flush()
|
|
|
403
|
|
|
404
|
|
|
405 ## Filter a fasta file in two fasta files using the length of each sequence as a criteron
|
|
|
406 #
|
|
|
407 # @param len_min integer length sequence criterion to filter
|
|
|
408 # @param inFileName string name of the input fasta file
|
|
|
409 # @param verbose integer (default = 0)
|
|
|
410 #
|
|
|
411 @staticmethod
|
|
|
412 def dbLengthFilter( len_min, inFileName, verbose=0 ):
|
|
|
413 file_db = open( inFileName, "r" )
|
|
|
414 file_dbInf = open( inFileName+".Inf"+str(len_min), "w" )
|
|
|
415 file_dbSup = open( inFileName+".Sup"+str(len_min), "w" )
|
|
|
416 seq = Bioseq()
|
|
|
417 numseq = 0
|
|
|
418 nbsave = 0
|
|
|
419
|
|
|
420 seq.read( file_db )
|
|
|
421 while seq.sequence:
|
|
|
422 l = seq.getLength()
|
|
|
423 numseq = numseq + 1
|
|
|
424 if l >= len_min:
|
|
|
425 seq.write( file_dbSup )
|
|
|
426 if verbose > 0:
|
|
|
427 print 'sequence #',numseq,'=',l,'[',seq.header[0:40],'...] Sup !!'
|
|
|
428 nbsave=nbsave+1
|
|
|
429 else:
|
|
|
430 seq.write( file_dbInf )
|
|
|
431 if verbose > 0:
|
|
|
432 print 'sequence #',numseq,'=',l,'[',seq.header[0:40],'...] Inf !!'
|
|
|
433 nbsave=nbsave+1
|
|
|
434 seq.read( file_db )
|
|
|
435
|
|
|
436 file_db.close()
|
|
|
437 file_dbInf.close()
|
|
|
438 file_dbSup.close()
|
|
|
439 if verbose > 0:
|
|
|
440 print nbsave,'saved sequences in ',inFileName+".Inf"+str(len_min)," and ", inFileName+".Sup"+str(len_min)
|
|
|
441
|
|
|
442
|
|
|
443 ## Extract the longest sequences from a fasta file
|
|
|
444 #
|
|
|
445 # @param num integer maximum number of sequences in the output file
|
|
|
446 # @param inFileName string name of the input fasta file
|
|
|
447 # @param outFileName string name of the output fasta file
|
|
|
448 # @param minThresh integer minimum length threshold (default=0)
|
|
|
449 # @param verbose integer (default = 0)
|
|
|
450 #
|
|
|
451 @staticmethod
|
|
|
452 def dbLongestSequences( num, inFileName, outFileName="", verbose=0, minThresh=0 ):
|
|
|
453 bsDB = BioseqDB( inFileName )
|
|
|
454 if verbose > 0:
|
|
|
455 print "nb of input sequences: %i" % ( bsDB.getSize() )
|
|
|
456
|
|
|
457 if outFileName == "":
|
|
|
458 outFileName = inFileName + ".best" + str(num)
|
|
|
459 outFile = open( outFileName, "w" )
|
|
|
460
|
|
|
461 if bsDB.getSize()==0:
|
|
|
462 return 0
|
|
|
463
|
|
|
464 num = int(num)
|
|
|
465 if verbose > 0:
|
|
|
466 print "keep the %i longest sequences" % ( num )
|
|
|
467 if minThresh > 0:
|
|
|
468 print "with length > %i bp" % ( minThresh )
|
|
|
469 sys.stdout.flush()
|
|
|
470
|
|
|
471 # retrieve the length of each input sequence
|
|
|
472 tmpLSeqLgth = []
|
|
|
473 seqNum = 0
|
|
|
474 for bs in bsDB.db:
|
|
|
475 seqNum += 1
|
|
|
476 tmpLSeqLgth.append( bs.getLength() )
|
|
|
477 if verbose > 1:
|
|
|
478 print "%d seq %s : %d bp" % ( seqNum, bs.header[0:40], bs.getLength() )
|
|
|
479 sys.stdout.flush()
|
|
|
480
|
|
|
481 # sort the lengths
|
|
|
482 tmpLSeqLgth.sort()
|
|
|
483 tmpLSeqLgth.reverse()
|
|
|
484
|
|
|
485 # select the longest
|
|
|
486 lSeqLgth = []
|
|
|
487 for i in xrange( 0, min(num,len(tmpLSeqLgth)) ):
|
|
|
488 if tmpLSeqLgth[i] >= minThresh:
|
|
|
489 lSeqLgth.append( tmpLSeqLgth[i] )
|
|
|
490 if verbose > 0:
|
|
|
491 print "selected max length: %i" % ( max(lSeqLgth) )
|
|
|
492 print "selected min length: %i" % ( min(lSeqLgth) )
|
|
|
493 sys.stdout.flush()
|
|
|
494
|
|
|
495 # save the longest
|
|
|
496 inFile = open( inFileName )
|
|
|
497 seqNum = 0
|
|
|
498 nbSave = 0
|
|
|
499 for bs in bsDB.db:
|
|
|
500 seqNum += 1
|
|
|
501 if bs.getLength() >= min(lSeqLgth) and bs.getLength() >= minThresh:
|
|
|
502 bs.write( outFile )
|
|
|
503 if verbose > 1:
|
|
|
504 print "%d seq %s : saved !" % ( seqNum, bs.header[0:40] )
|
|
|
505 sys.stdout.flush()
|
|
|
506 nbSave += 1
|
|
|
507 if nbSave == num:
|
|
|
508 break
|
|
|
509 inFile.close()
|
|
|
510 outFile.close()
|
|
|
511 if verbose > 0:
|
|
|
512 print nbSave, "saved sequences in ", outFileName
|
|
|
513 sys.stdout.flush()
|
|
|
514
|
|
|
515 return 0
|
|
|
516
|
|
|
517
|
|
|
518 ## Extract all the sequence headers from a fasta file and write them in a new fasta file
|
|
|
519 #
|
|
|
520 # @param inFileName string name of the input fasta file
|
|
|
521 # @param outFileName string name of the output file recording the headers (default = inFileName + '.headers')
|
|
|
522 #
|
|
|
523 @staticmethod
|
|
|
524 def dbExtractSeqHeaders( inFileName, outFileName="" ):
|
|
|
525 lHeaders = FastaUtils.dbHeaders( inFileName )
|
|
|
526
|
|
|
527 if outFileName == "":
|
|
|
528 outFileName = inFileName + ".headers"
|
|
|
529
|
|
|
530 outFile = open( outFileName, "w" )
|
|
|
531 for i in lHeaders:
|
|
|
532 outFile.write( i + "\n" )
|
|
|
533 outFile.close()
|
|
|
534
|
|
|
535 return 0
|
|
|
536
|
|
|
537
|
|
|
538 ## Extract sequences and their headers selected by a given pattern from a fasta file and write them in a new fasta file
|
|
|
539 #
|
|
|
540 # @param pattern regular expression to search in headers
|
|
|
541 # @param inFileName string name of the input fasta file
|
|
|
542 # @param outFileName string name of the output file recording the selected bioseq (default = inFileName + '.extracted')
|
|
|
543 # @param verbose integer verbosity level (default = 0)
|
|
|
544 #
|
|
|
545 @staticmethod
|
|
|
546 def dbExtractByPattern( pattern, inFileName, outFileName="", verbose=0 ):
|
|
|
547 if pattern == "":
|
|
|
548 return
|
|
|
549
|
|
|
550 if outFileName == "":
|
|
|
551 outFileName = inFileName + '.extracted'
|
|
|
552 outFile = open( outFileName, 'w' )
|
|
|
553
|
|
|
554 patternTosearch = re.compile( pattern )
|
|
|
555 bioseq = Bioseq()
|
|
|
556 bioseqNb = 0
|
|
|
557 savedBioseqNb = 0
|
|
|
558 inFile = open( inFileName, "r" )
|
|
|
559 bioseq.read( inFile )
|
|
|
560 while bioseq.sequence:
|
|
|
561 bioseqNb = bioseqNb + 1
|
|
|
562 m = patternTosearch.search( bioseq.header )
|
|
|
563 if m:
|
|
|
564 bioseq.write( outFile )
|
|
|
565 if verbose > 1:
|
|
|
566 print 'sequence num',bioseqNb,'matched on',m.group(),'[',bioseq.header[0:40],'...] saved !!'
|
|
|
567 savedBioseqNb = savedBioseqNb + 1
|
|
|
568 bioseq.read( inFile )
|
|
|
569 inFile.close()
|
|
|
570
|
|
|
571 outFile.close()
|
|
|
572
|
|
|
573 if verbose > 0:
|
|
|
574 print "%i sequences saved in file '%s'" % ( savedBioseqNb, outFileName )
|
|
|
575
|
|
|
576
|
|
|
577 ## Extract sequences and their headers selected by patterns contained in a file, from a fasta file and write them in a new fasta file
|
|
|
578 #
|
|
|
579 # @param patternFileName string file containing regular expression to search in headers
|
|
|
580 # @param inFileName string name of the input fasta file
|
|
|
581 # @param outFileName string name of the output file recording the selected bioseq (default = inFileName + '.extracted')
|
|
|
582 # @param verbose integer verbosity level (default = 0)
|
|
|
583 #
|
|
|
584 @staticmethod
|
|
|
585 def dbExtractByFilePattern( patternFileName, inFileName, outFileName="", verbose=0 ):
|
|
|
586
|
|
|
587 if patternFileName == "":
|
|
|
588 print "ERROR: no file of pattern"
|
|
|
589 sys.exit(1)
|
|
|
590
|
|
|
591 bioseq = Bioseq()
|
|
|
592 bioseqNb = 0
|
|
|
593 savedBioseqNb = 0
|
|
|
594 lHeaders = []
|
|
|
595
|
|
|
596 inFile = open( inFileName, "r" )
|
|
|
597 bioseq.read( inFile )
|
|
|
598 while bioseq.sequence != None:
|
|
|
599 lHeaders.append( bioseq.header )
|
|
|
600 bioseq.read( inFile )
|
|
|
601 inFile.close()
|
|
|
602
|
|
|
603 lHeadersToKeep = []
|
|
|
604 patternFile = open( patternFileName, "r" )
|
|
|
605 for pattern in patternFile:
|
|
|
606 if verbose > 0:
|
|
|
607 print "pattern: ",pattern[:-1]; sys.stdout.flush()
|
|
|
608
|
|
|
609 patternToSearch = re.compile(pattern[:-1])
|
|
|
610 for h in lHeaders:
|
|
|
611 if patternToSearch.search(h):
|
|
|
612 lHeadersToKeep.append(h)
|
|
|
613 patternFile.close()
|
|
|
614
|
|
|
615 if outFileName == "":
|
|
|
616 outFileName = inFileName + ".extracted"
|
|
|
617 outFile=open( outFileName, "w" )
|
|
|
618
|
|
|
619 inFile = open( inFileName, "r" )
|
|
|
620 bioseq.read(inFile)
|
|
|
621 while bioseq.sequence:
|
|
|
622 bioseqNb += 1
|
|
|
623 if bioseq.header in lHeadersToKeep:
|
|
|
624 bioseq.write(outFile)
|
|
|
625 if verbose > 1:
|
|
|
626 print 'sequence num',bioseqNb,'[',bioseq.header[0:40],'...] saved !!'; sys.stdout.flush()
|
|
|
627 savedBioseqNb += 1
|
|
|
628 bioseq.read(inFile)
|
|
|
629 inFile.close()
|
|
|
630
|
|
|
631 outFile.close()
|
|
|
632
|
|
|
633 if verbose > 0:
|
|
|
634 print "%i sequences saved in file '%s'" % ( savedBioseqNb, outFileName )
|
|
|
635
|
|
|
636
|
|
|
637 ## Extract sequences and their headers not selected by a given pattern from a fasta file and write them in a new fasta file
|
|
|
638 #
|
|
|
639 # @param pattern regular expression to search in headers
|
|
|
640 # @param inFileName string name of the input fasta file
|
|
|
641 # @param outFileName string name of the output file recording the selected bioseq (default = inFileName + '.extracted')
|
|
|
642 # @param verbose integer verbosity level (default = 0)
|
|
|
643 #
|
|
|
644 @staticmethod
|
|
|
645 def dbCleanByPattern( pattern, inFileName, outFileName="", verbose=0 ):
|
|
|
646 if pattern == "":
|
|
|
647 return
|
|
|
648
|
|
|
649 patternToSearch = re.compile(pattern)
|
|
|
650
|
|
|
651 if outFileName == "":
|
|
|
652 outFileName = inFileName + '.cleaned'
|
|
|
653 outFile = open(outFileName,'w')
|
|
|
654
|
|
|
655 bioseq = Bioseq()
|
|
|
656 bioseqNb = 0
|
|
|
657 savedBioseqNb = 0
|
|
|
658 inFile = open(inFileName)
|
|
|
659 bioseq.read(inFile)
|
|
|
660 while bioseq.sequence != None:
|
|
|
661 bioseqNb += 1
|
|
|
662 if not patternToSearch.search(bioseq.header):
|
|
|
663 bioseq.write(outFile)
|
|
|
664 if verbose > 1:
|
|
|
665 print 'sequence num',bioseqNb,'[',bioseq.header[0:40],'...] saved !!'
|
|
|
666 savedBioseqNb += 1
|
|
|
667 bioseq.read(inFile)
|
|
|
668 inFile.close()
|
|
|
669
|
|
|
670 outFile.close()
|
|
|
671
|
|
|
672 if verbose > 0:
|
|
|
673 print "%i sequences saved in file '%s'" % ( savedBioseqNb, outFileName )
|
|
|
674
|
|
|
675
|
|
|
676 ## Extract sequences and their headers not selected by patterns contained in a file, from a fasta file and write them in a new fasta file
|
|
|
677 #
|
|
|
678 # @param patternFileName string file containing regular expression to search in headers
|
|
|
679 # @param inFileName string name of the input fasta file
|
|
|
680 # @param outFileName string name of the output file recording the selected bioseq (default = inFileName + '.extracted')
|
|
|
681 # @param verbose integer verbosity level (default = 0)
|
|
|
682 #
|
|
|
683 @staticmethod
|
|
|
684 def dbCleanByFilePattern( patternFileName, inFileName, outFileName="", verbose=0 ):
|
|
|
685 if patternFileName == "":
|
|
|
686 print "ERROR: no file of pattern"
|
|
|
687 sys.exit(1)
|
|
|
688
|
|
|
689 bioseq = Bioseq()
|
|
|
690 bioseqNb = 0
|
|
|
691 savedBioseqNb = 0
|
|
|
692 lHeaders = []
|
|
|
693 inFile = open( inFileName, "r" )
|
|
|
694 bioseq.read( inFile )
|
|
|
695 while bioseq.sequence != None:
|
|
|
696 bioseqNb += 1
|
|
|
697 lHeaders.append( bioseq.header )
|
|
|
698 bioseq.read( inFile )
|
|
|
699 inFile.close()
|
|
|
700
|
|
|
701 patternFile = open( patternFileName, "r")
|
|
|
702 lHeadersToRemove = []
|
|
|
703 for pattern in patternFile:
|
|
|
704 if verbose > 0:
|
|
|
705 print "pattern: ",pattern[:-1]; sys.stdout.flush()
|
|
|
706
|
|
|
707 patternToSearch = re.compile( pattern[:-1] )
|
|
|
708 for h in lHeaders:
|
|
|
709 if patternToSearch.search(h):
|
|
|
710 lHeadersToRemove.append(h)
|
|
|
711 patternFile.close()
|
|
|
712
|
|
|
713 if outFileName == "":
|
|
|
714 outFileName = inFileName + '.cleaned'
|
|
|
715 outFile = open( outFileName, 'w' )
|
|
|
716
|
|
|
717 bioseqNum = 0
|
|
|
718 inFile=open( inFileName )
|
|
|
719 bioseq.read( inFile )
|
|
|
720 while bioseq.sequence != None:
|
|
|
721 bioseqNum += 1
|
|
|
722 if bioseq.header not in lHeadersToRemove:
|
|
|
723 bioseq.write( outFile )
|
|
|
724 if verbose > 1:
|
|
|
725 print 'sequence num',bioseqNum,'/',bioseqNb,'[',bioseq.header[0:40],'...] saved !!'; sys.stdout.flush()
|
|
|
726 savedBioseqNb += 1
|
|
|
727 bioseq.read( inFile )
|
|
|
728 inFile.close()
|
|
|
729
|
|
|
730 outFile.close()
|
|
|
731
|
|
|
732 if verbose > 0:
|
|
|
733 print "%i sequences saved in file '%s'" % ( savedBioseqNb, outFileName )
|
|
|
734
|
|
|
735
|
|
|
736 ## Find sequence's ORFs from a fasta file and write them in a tab file
|
|
|
737 #
|
|
|
738 # @param inFileName string name of the input fasta file
|
|
|
739 # @param orfMaxNb integer Select orfMaxNb best ORFs
|
|
|
740 # @param orfMinLength integer Keep ORFs with length > orfMinLength
|
|
|
741 # @param outFileName string name of the output fasta file (default = inFileName + '.orf.map')
|
|
|
742 # @param verbose integer verbosity level (default = 0)
|
|
|
743 #
|
|
|
744 @staticmethod
|
|
|
745 def dbORF( inFileName, orfMaxNb = 0, orfMinLength = 0, outFileName = "", verbose=0 ):
|
|
|
746 if outFileName == "":
|
|
|
747 outFileName = inFileName + ".ORF.map"
|
|
|
748 outFile = open( outFileName, "w" )
|
|
|
749
|
|
|
750 bioseq = Bioseq()
|
|
|
751 bioseqNb = 0
|
|
|
752
|
|
|
753 inFile = open( inFileName )
|
|
|
754 bioseq.read( inFile )
|
|
|
755 while bioseq.sequence != None:
|
|
|
756 bioseq.upCase()
|
|
|
757 bioseqNb += 1
|
|
|
758 if verbose > 0:
|
|
|
759 print 'sequence num',bioseqNb,'=',bioseq.getLength(),'[',bioseq.header[0:40],'...]'
|
|
|
760
|
|
|
761 orf = bioseq.findORF()
|
|
|
762 bestOrf = []
|
|
|
763 for i in orf.keys():
|
|
|
764 orfLen = len(orf[i])
|
|
|
765 for j in xrange(1, orfLen):
|
|
|
766 start = orf[i][j-1] + 4
|
|
|
767 end = orf[i][j] + 3
|
|
|
768 if end - start >= orfMinLength:
|
|
|
769 bestOrf.append( ( end-start, i+1, start, end ) )
|
|
|
770
|
|
|
771 bioseq.reverseComplement()
|
|
|
772
|
|
|
773 orf = bioseq.findORF()
|
|
|
774 seqLen = bioseq.getLength()
|
|
|
775 for i in orf.keys():
|
|
|
776 orfLen = len(orf[i])
|
|
|
777 for j in xrange(1, orfLen):
|
|
|
778 start = seqLen - orf[i][j-1] - 3
|
|
|
779 end = seqLen - orf[i][j] - 2
|
|
|
780 if start - end >= orfMinLength:
|
|
|
781 bestOrf.append( ( start-end, (i+1)*-1, start, end ) )
|
|
|
782
|
|
|
783 bestOrf.sort()
|
|
|
784 bestOrf.reverse()
|
|
|
785 bestOrfNb = len(bestOrf)
|
|
|
786 if orfMaxNb != 0 and orfMaxNb < bestOrfNb:
|
|
|
787 bestOrfNb = orfMaxNb
|
|
|
788 for i in xrange(0, bestOrfNb):
|
|
|
789 if verbose > 0:
|
|
|
790 print bestOrf[i]
|
|
|
791 outFile.write("%s\t%s\t%d\t%d\n"%("ORF|"+str(bestOrf[i][1])+\
|
|
|
792 "|"+str(bestOrf[i][0]),bioseq.header,
|
|
|
793 bestOrf[i][2],bestOrf[i][3]))
|
|
|
794 bioseq.read( inFile )
|
|
|
795
|
|
|
796 inFile.close()
|
|
|
797 outFile.close()
|
|
|
798
|
|
|
799 return 0
|
|
|
800
|
|
|
801
|
|
|
802 ## Sort sequences by increasing length (write a new file)
|
|
|
803 #
|
|
|
804 # @param inFileName string name of the input fasta file
|
|
|
805 # @param outFileName string name of the output fasta file
|
|
|
806 # @param verbose integer verbosity level
|
|
|
807 #
|
|
|
808 @staticmethod
|
|
|
809 def sortSequencesByIncreasingLength(inFileName, outFileName, verbose=0):
|
|
|
810 if verbose > 0:
|
|
|
811 print "sort sequences by increasing length"
|
|
|
812 sys.stdout.flush()
|
|
|
813 if not os.path.exists( inFileName ):
|
|
|
814 print "ERROR: file '%s' doesn't exist" % ( inFileName )
|
|
|
815 sys.exit(1)
|
|
|
816
|
|
|
817 # read each seq one by one
|
|
|
818 # save them in distinct temporary files
|
|
|
819 # with their length in the name
|
|
|
820 inFileHandler = open( inFileName, "r" )
|
|
|
821 countSeq = 0
|
|
|
822 bs = Bioseq()
|
|
|
823 bs.read( inFileHandler )
|
|
|
824 while bs.header:
|
|
|
825 countSeq += 1
|
|
|
826 tmpFile = "%ibp_%inb" % ( bs.getLength(), countSeq )
|
|
|
827 bs.appendBioseqInFile( tmpFile )
|
|
|
828 if verbose > 1:
|
|
|
829 print "%s (%i bp) saved in '%s'" % ( bs.header, bs.getLength(), tmpFile )
|
|
|
830 bs.header = ""
|
|
|
831 bs.sequence = ""
|
|
|
832 bs.read( inFileHandler )
|
|
|
833 inFileHandler.close()
|
|
|
834
|
|
|
835 # sort temporary file names
|
|
|
836 # concatenate them into the output file
|
|
|
837 if os.path.exists( outFileName ):
|
|
|
838 os.remove( outFileName )
|
|
|
839 lFiles = glob.glob( "*bp_*nb" )
|
|
|
840 lFiles.sort( key=lambda s:int(s.split("bp_")[0]) )
|
|
|
841 for fileName in lFiles:
|
|
|
842 cmd = "cat %s >> %s" % ( fileName, outFileName )
|
|
|
843 returnValue = os.system( cmd )
|
|
|
844 if returnValue != 0:
|
|
|
845 print "ERROR while concatenating '%s' with '%s'" % ( fileName, outFileName )
|
|
|
846 sys.exit(1)
|
|
|
847 os.remove( fileName )
|
|
|
848
|
|
|
849 return 0
|
|
|
850
|
|
|
851
|
|
|
852 ## Sort sequences by header
|
|
|
853 #
|
|
|
854 # @param inFileName string name of the input fasta file
|
|
|
855 # @param outFileName string name of the output fasta file
|
|
|
856 # @param verbose integer verbosity level
|
|
|
857 #
|
|
|
858 @staticmethod
|
|
|
859 def sortSequencesByHeader(inFileName, outFileName = "", verbose = 0):
|
|
|
860 if outFileName == "":
|
|
|
861 outFileName = "%s_sortByHeaders.fa" % os.path.splitext(inFileName)[0]
|
|
|
862 iBioseqDB = BioseqDB(inFileName)
|
|
|
863 f = open(outFileName, "w")
|
|
|
864 lHeaders = sorted(iBioseqDB.getHeaderList())
|
|
|
865 for header in lHeaders:
|
|
|
866 iBioseq = iBioseqDB.fetch(header)
|
|
|
867 iBioseq.write(f)
|
|
|
868 f.close()
|
|
|
869
|
|
|
870
|
|
|
871 ## Return a dictionary which keys are the headers and values the length of the sequences
|
|
|
872 #
|
|
|
873 # @param inFile string name of the input fasta file
|
|
|
874 # @param verbose integer verbosity level
|
|
|
875 #
|
|
|
876 @staticmethod
|
|
|
877 def getLengthPerHeader( inFile, verbose=0 ):
|
|
|
878 dHeader2Length = {}
|
|
|
879
|
|
|
880 inFileHandler = open( inFile, "r" )
|
|
|
881 currentSeqHeader = ""
|
|
|
882 currentSeqLength = 0
|
|
|
883 line = inFileHandler.readline()
|
|
|
884 while line:
|
|
|
885 if line[0] == ">":
|
|
|
886 if currentSeqHeader != "":
|
|
|
887 dHeader2Length[ currentSeqHeader ] = currentSeqLength
|
|
|
888 currentSeqLength = 0
|
|
|
889 currentSeqHeader = line[1:-1]
|
|
|
890 if verbose > 0:
|
|
|
891 print "current header: %s" % ( currentSeqHeader )
|
|
|
892 sys.stdout.flush()
|
|
|
893 else:
|
|
|
894 currentSeqLength += len( line.replace("\n","") )
|
|
|
895 line = inFileHandler.readline()
|
|
|
896 dHeader2Length[ currentSeqHeader ] = currentSeqLength
|
|
|
897 inFileHandler.close()
|
|
|
898
|
|
|
899 return dHeader2Length
|
|
|
900
|
|
|
901
|
|
|
902 ## Convert headers from a fasta file having chunk coordinates
|
|
|
903 #
|
|
|
904 # @param inFile string name of the input fasta file
|
|
|
905 # @param mapFile string name of the map file with the coordinates of the chunks on the chromosomes
|
|
|
906 # @param outFile string name of the output file
|
|
|
907 #
|
|
|
908 @staticmethod
|
|
|
909 def convertFastaHeadersFromChkToChr(inFile, mapFile, outFile):
|
|
|
910 inFileHandler = open(inFile, "r")
|
|
|
911 outFileHandler = open(outFile, "w")
|
|
|
912 dChunk2Map = MapUtils.getDictPerNameFromMapFile(mapFile)
|
|
|
913 iConvCoord = ConvCoord()
|
|
|
914 line = inFileHandler.readline()
|
|
|
915 while line:
|
|
|
916 if line[0] == ">":
|
|
|
917 if "{Fragment}" in line:
|
|
|
918 chkName = line.split(" ")[1]
|
|
|
919 chrName = dChunk2Map[chkName].seqname
|
|
|
920 lCoordPairs = line.split(" ")[3].split(",")
|
|
|
921 lRangesOnChk = []
|
|
|
922 for i in lCoordPairs:
|
|
|
923 iRange = Range(chkName, int(i.split("..")[0]), int(i.split("..")[1]))
|
|
|
924 lRangesOnChk.append(iRange)
|
|
|
925 lRangesOnChr = []
|
|
|
926 for iRange in lRangesOnChk:
|
|
|
927 lRangesOnChr.append(iConvCoord.getRangeOnChromosome(iRange, dChunk2Map))
|
|
|
928 newHeader = line[1:-1].split(" ")[0]
|
|
|
929 newHeader += " %s" % chrName
|
|
|
930 newHeader += " {Fragment}"
|
|
|
931 newHeader += " %i..%i" % (lRangesOnChr[0].start, lRangesOnChr[0].end)
|
|
|
932 for iRange in lRangesOnChr[1:]:
|
|
|
933 newHeader += ",%i..%i" % (iRange.start, iRange.end)
|
|
|
934 outFileHandler.write(">%s\n" % newHeader)
|
|
|
935 else:
|
|
|
936 chkName = line.split("_")[1].split(" ")[0]
|
|
|
937 chrName = dChunk2Map[chkName].seqname
|
|
|
938 coords = line[line.find("[")+1 : line.find("]")]
|
|
|
939 start = int(coords.split(",")[0])
|
|
|
940 end = int(coords.split(",")[1])
|
|
|
941 iRangeOnChk = Range(chkName, start, end)
|
|
|
942 iRangeOnChr = iConvCoord.getRangeOnChromosome(iRangeOnChk, dChunk2Map)
|
|
|
943 newHeader = line[1:-1].split("_")[0]
|
|
|
944 newHeader += " %s" % chrName
|
|
|
945 newHeader += " %s" % line[line.find("(") : line.find(")")+1]
|
|
|
946 newHeader += " %i..%i" % (iRangeOnChr.getStart(), iRangeOnChr.getEnd())
|
|
|
947 outFileHandler.write(">%s\n" % newHeader)
|
|
|
948 else:
|
|
|
949 outFileHandler.write(line)
|
|
|
950 line = inFileHandler.readline()
|
|
|
951 inFileHandler.close()
|
|
|
952 outFileHandler.close()
|
|
|
953
|
|
|
954
|
|
|
955 ## Convert a fasta file to a length file
|
|
|
956 #
|
|
|
957 # @param inFile string name of the input fasta file
|
|
|
958 # @param outFile string name of the output file
|
|
|
959 #
|
|
|
960 @staticmethod
|
|
|
961 def convertFastaToLength(inFile, outFile = ""):
|
|
|
962 if outFile == "":
|
|
|
963 outFile = "%s.length" % inFile
|
|
|
964
|
|
|
965 if inFile != "":
|
|
|
966 with open(inFile, "r") as inFH:
|
|
|
967 with open(outFile, "w") as outFH:
|
|
|
968 bioseq = Bioseq()
|
|
|
969 bioseq.read(inFH)
|
|
|
970 while bioseq.sequence != None:
|
|
|
971 seqLen = bioseq.getLength()
|
|
|
972 outFH.write("%s\t%d\n" % (bioseq.header.split()[0], seqLen))
|
|
|
973 bioseq.read(inFH)
|
|
|
974
|
|
|
975
|
|
|
976 ## Convert a fasta file to a seq file
|
|
|
977 #
|
|
|
978 # @param inFile string name of the input fasta file
|
|
|
979 # @param outFile string name of the output file
|
|
|
980 #
|
|
|
981 @staticmethod
|
|
|
982 def convertFastaToSeq(inFile, outFile = ""):
|
|
|
983 if outFile == "":
|
|
|
984 outFile = "%s.seq" % inFile
|
|
|
985
|
|
|
986 if inFile != "":
|
|
|
987 with open(inFile, "r") as inFH:
|
|
|
988 with open(outFile, "w") as outFH:
|
|
|
989 bioseq = Bioseq()
|
|
|
990 bioseq.read(inFH)
|
|
|
991 while bioseq.sequence != None:
|
|
|
992 seqLen = bioseq.getLength()
|
|
|
993 outFH.write("%s\t%s\t%s\t%d\n" % (bioseq.header.split()[0], \
|
|
|
994 bioseq.sequence, bioseq.header, seqLen))
|
|
|
995 bioseq.read(inFH)
|
|
|
996
|
|
|
997
|
|
|
998 ## Splice an input fasta file using coordinates in a Map file
|
|
|
999 #
|
|
|
1000 # @note the coordinates should be merged beforehand!
|
|
|
1001 #
|
|
|
1002 @staticmethod
|
|
|
1003 def spliceFromCoords( genomeFile, coordFile, obsFile ):
|
|
|
1004 genomeFileHandler = open( genomeFile, "r" )
|
|
|
1005 obsFileHandler = open( obsFile, "w" )
|
|
|
1006 dChr2Maps = MapUtils.getDictPerSeqNameFromMapFile( coordFile )
|
|
|
1007
|
|
|
1008 bs = Bioseq()
|
|
|
1009 bs.read( genomeFileHandler )
|
|
|
1010 while bs.sequence:
|
|
|
1011 if dChr2Maps.has_key( bs.header ):
|
|
|
1012 lCoords = MapUtils.getMapListSortedByIncreasingMinThenMax( dChr2Maps[ bs.header ] )
|
|
|
1013 splicedSeq = ""
|
|
|
1014 currentSite = 0
|
|
|
1015 for iMap in lCoords:
|
|
|
1016 minSplice = iMap.getMin() - 1
|
|
|
1017 if minSplice > currentSite:
|
|
|
1018 splicedSeq += bs.sequence[ currentSite : minSplice ]
|
|
|
1019 currentSite = iMap.getMax()
|
|
|
1020 splicedSeq += bs.sequence[ currentSite : ]
|
|
|
1021 bs.sequence = splicedSeq
|
|
|
1022 bs.write( obsFileHandler )
|
|
|
1023 bs.read( genomeFileHandler )
|
|
|
1024
|
|
|
1025 genomeFileHandler.close()
|
|
|
1026 obsFileHandler.close()
|
|
|
1027
|
|
|
1028
|
|
|
1029 ## Shuffle input sequences (single file or files in a directory)
|
|
|
1030 #
|
|
|
1031 @staticmethod
|
|
|
1032 def dbShuffle( inData, outData, verbose=0 ):
|
|
|
1033 if CheckerUtils.isExecutableInUserPath("esl-shuffle"):
|
|
|
1034 prg = "esl-shuffle"
|
|
|
1035 else : prg = "shuffle"
|
|
|
1036 genericCmd = prg + " -d INPUT > OUTPUT"
|
|
|
1037 if os.path.isfile( inData ):
|
|
|
1038 if verbose > 0:
|
|
|
1039 print "shuffle input file '%s'" % inData
|
|
|
1040 cmd = genericCmd.replace("INPUT",inData).replace("OUTPUT",outData)
|
|
|
1041 print cmd
|
|
|
1042 returnStatus = os.system( cmd )
|
|
|
1043 if returnStatus != 0:
|
|
|
1044 sys.stderr.write( "ERROR: 'shuffle' returned '%i'\n" % returnStatus )
|
|
|
1045 sys.exit(1)
|
|
|
1046
|
|
|
1047 elif os.path.isdir( inData ):
|
|
|
1048 if verbose > 0:
|
|
|
1049 print "shuffle files in input directory '%s'" % inData
|
|
|
1050 if os.path.exists( outData ):
|
|
|
1051 shutil.rmtree( outData )
|
|
|
1052 os.mkdir( outData )
|
|
|
1053 lInputFiles = glob.glob( "%s/*.fa" %( inData ) )
|
|
|
1054 nbFastaFiles = 0
|
|
|
1055 for inputFile in lInputFiles:
|
|
|
1056 nbFastaFiles += 1
|
|
|
1057 if verbose > 1:
|
|
|
1058 print "%3i / %3i" % ( nbFastaFiles, len(lInputFiles) )
|
|
|
1059 fastaBaseName = os.path.basename( inputFile )
|
|
|
1060 prefix, extension = os.path.splitext( fastaBaseName )
|
|
|
1061 cmd = genericCmd.replace("INPUT",inputFile).replace("OUTPUT","%s/%s_shuffle.fa"%(outData,prefix))
|
|
|
1062 returnStatus = os.system( cmd )
|
|
|
1063 if returnStatus != 0:
|
|
|
1064 sys.stderr.write( "ERROR: 'shuffle' returned '%i'\n" % returnStatus )
|
|
|
1065 sys.exit(1)
|
|
|
1066
|
|
|
1067
|
|
|
1068 ## Convert a cluster file (one line = one cluster = one headers list) into a fasta file with cluster info in headers
|
|
|
1069 #
|
|
|
1070 # @param inClusterFileName string input cluster file name
|
|
|
1071 # @param inFastaFileName string input fasta file name
|
|
|
1072 # @param outFileName string output file name
|
|
|
1073 # @param verbosity integer verbosity
|
|
|
1074 #
|
|
|
1075 @staticmethod
|
|
|
1076 def convertClusterFileToFastaFile(inClusterFileName, inFastaFileName, outFileName, clusteringTool = "", verbosity = 0):
|
|
|
1077 dHeader2ClusterClusterMember, clusterIdForSingletonCluster = FastaUtils._createHeader2ClusterMemberDict(inClusterFileName, verbosity)
|
|
|
1078 iFastaParser = FastaParser(inFastaFileName)
|
|
|
1079 with open(outFileName, "w") as f:
|
|
|
1080 for iSequence in iFastaParser.getIterator():
|
|
|
1081
|
|
|
1082 header = iSequence.getName()
|
|
|
1083 if dHeader2ClusterClusterMember.get(header):
|
|
|
1084 cluster = dHeader2ClusterClusterMember[header][0]
|
|
|
1085 member = dHeader2ClusterClusterMember[header][1]
|
|
|
1086 else:
|
|
|
1087 clusterIdForSingletonCluster += 1
|
|
|
1088 cluster = clusterIdForSingletonCluster
|
|
|
1089 member = 1
|
|
|
1090
|
|
|
1091 newHeader = "%sCluster%sMb%s_%s" % (clusteringTool, cluster, member, header)
|
|
|
1092 iSequence.setName(newHeader)
|
|
|
1093 f.write(iSequence.printFasta())
|
|
|
1094
|
|
|
1095 @staticmethod
|
|
|
1096 def _createHeader2ClusterMemberDict(inClusterFileName, verbosity = 0):
|
|
|
1097 dHeader2ClusterClusterMember = {}
|
|
|
1098 clusterId = 0
|
|
|
1099 with open(inClusterFileName) as f:
|
|
|
1100 line = f.readline()
|
|
|
1101 while line:
|
|
|
1102 lineWithoutLastChar = line.rstrip()
|
|
|
1103 lHeaders = lineWithoutLastChar.split("\t")
|
|
|
1104 clusterId += 1
|
|
|
1105 if verbosity > 0:
|
|
|
1106 print "%i sequences in cluster %i" % (len(lHeaders), clusterId)
|
|
|
1107 memberId = 0
|
|
|
1108 for header in lHeaders:
|
|
|
1109 memberId += 1
|
|
|
1110 dHeader2ClusterClusterMember[header] = (clusterId, memberId)
|
|
|
1111 line = f.readline()
|
|
|
1112 if verbosity > 0:
|
|
|
1113 print "%i clusters" % clusterId
|
|
|
1114 return dHeader2ClusterClusterMember, clusterId
|
|
|
1115
|
|
|
1116 @staticmethod
|
|
|
1117 def convertClusteredFastaFileToMapFile(fastaFileNameFromClustering, outMapFileName = ""):
|
|
|
1118 """
|
|
|
1119 Write a map file from fasta output of clustering tool.
|
|
|
1120 Warning: only works if input fasta headers are formated like LTRharvest fasta output.
|
|
|
1121 """
|
|
|
1122 if not outMapFileName:
|
|
|
1123 outMapFileName = "%s.map" % (os.path.splitext(fastaFileNameFromClustering)[0])
|
|
|
1124
|
|
|
1125 fileDb = open(fastaFileNameFromClustering , "r")
|
|
|
1126 fileMap = open(outMapFileName, "w")
|
|
|
1127 seq = Bioseq()
|
|
|
1128 numseq = 0
|
|
|
1129 while 1:
|
|
|
1130 seq.read(fileDb)
|
|
|
1131 if seq.sequence == None:
|
|
|
1132 break
|
|
|
1133 numseq = numseq + 1
|
|
|
1134 ID = seq.header.split(' ')[0].split('_')[0]
|
|
|
1135 chunk = seq.header.split(' ')[0].split('_')[1]
|
|
|
1136 start = seq.header.split(' ')[-1].split(',')[0][1:]
|
|
|
1137 end = seq.header.split(' ')[-1].split(',')[1][:-1]
|
|
|
1138 line = '%s\t%s\t%s\t%s' % (ID, chunk, start, end)
|
|
|
1139 fileMap.write(line + "\n")
|
|
|
1140
|
|
|
1141 fileDb.close()
|
|
|
1142 fileMap.close()
|
|
|
1143 print "saved in %s" % outMapFileName
|
|
|
1144
|
|
|
1145 @staticmethod
|
|
|
1146 def writeNstreches(fastaFileName, nbN = 2, outFileName = "", outFormat = "map"):
|
|
|
1147 outFormat = outFormat.lower()
|
|
|
1148 if outFormat in ["gff", "gff3"]:
|
|
|
1149 outFormat = "gff3"
|
|
|
1150 else:
|
|
|
1151 outFormat = "map"
|
|
|
1152
|
|
|
1153 lTNstretches = []
|
|
|
1154 if nbN != 0:
|
|
|
1155 iBSDB = BioseqDB(fastaFileName)
|
|
|
1156 for iBS in iBSDB.db:
|
|
|
1157 nbNFound = 0
|
|
|
1158 start = 1
|
|
|
1159 pos = 1
|
|
|
1160 lastPos = 0
|
|
|
1161
|
|
|
1162 while pos <= iBS.getLength():
|
|
|
1163 if nbNFound == 0:
|
|
|
1164 start = pos
|
|
|
1165
|
|
|
1166 while pos <= iBS.getLength() and iBS.getNtFromPosition(pos).lower() in ['n', 'x']:
|
|
|
1167 nbNFound += 1
|
|
|
1168 pos += 1
|
|
|
1169 lastPos = pos
|
|
|
1170
|
|
|
1171 if pos - lastPos >= nbN:
|
|
|
1172 if nbNFound >= nbN:
|
|
|
1173 lTNstretches.append((iBS.getHeader(), start, lastPos - 1))
|
|
|
1174 nbNFound = 0
|
|
|
1175 pos += 1
|
|
|
1176
|
|
|
1177 if nbNFound >= nbN:
|
|
|
1178 lTNstretches.append((iBS.getHeader(), start, lastPos - 1))
|
|
|
1179
|
|
|
1180 lTNstretches.sort(key = itemgetter(0, 1, 2))
|
|
|
1181
|
|
|
1182 if outFileName == "":
|
|
|
1183 outFileName = "%s_Nstretches.%s" % (os.path.splitext(os.path.split(fastaFileName)[1])[0], outFormat)
|
|
|
1184
|
|
|
1185 with open(outFileName, "w") as fH:
|
|
|
1186 if outFormat == "gff3":
|
|
|
1187 fH.write("##gff-version 3\n")
|
|
|
1188 for item in lTNstretches:
|
|
|
1189 seq = item[0]
|
|
|
1190 start = item[1]
|
|
|
1191 end = item[2]
|
|
|
1192 if outFormat == "gff3":
|
|
|
1193 fH.write("%s\tFastaUtils\tN_stretch\t%s\t%s\t.\t.\t.\tName=N_stretch_%s-%s\n" % (seq, start, end, start, end))
|
|
|
1194 else:
|
|
|
1195 fH.write("N_stretch\t%s\t%s\t%s\n" % (seq, start, end))
|
|
|
1196
|
|
|
1197 |
