Mercurial > repos > yufei-luo > s_mart

comparison smart_toolShed/commons/core/seq/FastaUtils.py @ 0:e0f8dcca02ed

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