filter sequence by length

[NicoBxl]

Hi,

To filter sequencesby length from a fasta file, I'm using now a bioperl script. But is there an another method (faster) ?

here's my little script

#!/usr/bin/perl

use warnings;
use strict;
use Bio::SeqIO;


my $file = $ARGV[0]; 
my $min = $ARGV[1];
my $max = $ARGV[2];
my $out = $ARGV[3];

open (FILE, ">>$out") or die ("Error : Cannot open file $out for writing..!\n");

my $seq_in  = Bio::SeqIO->new( -format => 'fasta',-file => $file);

while( my $seq1 = $seq_in->next_seq() ) {	
	
	my $id  = $seq1->primary_id;
	chomp $id;
	my $seq = $seq1->seq;
	chomp $seq;
	my $lseq = length($seq);
	if($lseq>=$min && $lseq <=$max){
		print FILE ">",$id,"\n",$seq,"\n";	
	}
}

[Giorgio C]

You could try Galaxy:

http://main.g2.bx.psu.edu/root?tool_...lter_by_length

It's very fast and you can do different type of analysis.

[NicoBxl]

Thanks,

is it possible to have this tool in command line ?

[Giorgio C]

From Galaxy i don't think is possible. But probably if you search you may find a similar script.

[maasha]

With Biopieces (www.biopieces.org) you can do:

read_fasta -i test.fna | grab -e 'SEQ_LEN > 10' | grab -e 'SEQ_LEN <= 100' | write_fasta -x

Cheers,


Martin

[tnabtaf]

Originally Posted by Giorgio C

From Galaxy i don't think is possible. But probably if you search you may find a similar script.

This is a command line tool that is included in the Galaxy distribution. I've included the code below.

To get any tool that is included with Galaxy, follow the instructions at http://getgalaxy.org/. Tools are defined in the tools directory. Many will have dependencies (although the one below does not).

#!/usr/bin/env python
"""
Input: fasta, minimal length, maximal length
Output: fasta
Return sequences whose lengths are within the range.
"""

import sys, os

assert sys.version_info[:2] >= ( 2, 4 )

def stop_err( msg ):
    sys.stderr.write( msg )
    sys.exit()

def __main__():
    input_filename = sys.argv[1]
    try:
        min_length = int( sys.argv[2] )
    except:
        stop_err( "Minimal length of the return sequence requires a numerical value." )
    try:
        max_length = int( sys.argv[3] )
    except:
        stop_err( "Maximum length of the return sequence requires a numerical value." )
    output_filename = sys.argv[4]
    output_handle = open( output_filename, 'w' )
    tmp_size = 0 #-1
    tmp_buf = ''
    at_least_one = 0
    for line in file(input_filename):
        if not line or line.startswith('#'):
            continue
        if line[0] == '>':
            if min_length <= tmp_size <= max_length or (min_length <= tmp_size and max_length == 0):
                output_handle.write(tmp_buf)
                at_least_one = 1
            tmp_buf = line
            tmp_size = 0                                                       
        else:
            if max_length == 0 or tmp_size < max_length:
                tmp_size += len(line.rstrip('\r\n'))
                tmp_buf += line
    # final flush of buffer
    if min_length <= tmp_size <= max_length or (min_length <= tmp_size and max_length == 0):
        output_handle.write(tmp_buf.rstrip('\r\n'))
        at_least_one = 1
    output_handle.close()
    if at_least_one == 0:
        print "There is no sequence that falls within your range."

if __name__ == "__main__" : __main__()

[Giorgio C]

Very useful. Thanks