I save this bit of code in a file called removesmalls, made it executable chmod +x
#!/bin/awk -f
!/^>/ {next}
{getline s}
length(s) >= i { print $0 “\n” s }
and then i used it with the following command
awk removesmalls i=200 contigs.fa > contigs-l200.fa
Result :Now its generating the out put file also but without any data
Please help me out.
Your original awk script assumes that the FASTA sequence following the header line spans only one line, which is not normally the case. As suggested by @SES above, bioawk is definitely a neat alternative.
A traditional Perl based approach would be like so:
## removesmalls.pl
#!/usr/bin/perl
use strict;
use warnings;
my $minlen = shift or die "Error: `minlen` parameter not provided\n";
{
local $/=">";
while(<>) {
chomp;
next unless /\w/;
s/>$//gs;
my @chunk = split /\n/;
my $header = shift @chunk;
my $seqlen = length join "", @chunk;
print ">$_" if($seqlen >= $minlen);
}
local $/="\n";
}
You would then invoke the script like so:
perl removesmalls.pl 200 contigs.fasta > contigs-l200.fasta
This solution is making use of the Perl record separator (which is a "\n" newline by default) and switching it to ">" so that when iterating through the FASTA file, we encounter whole FASTA records. This addresses the issue I mentioned above.
Hi this script is really helpful. I have run it with my fasta files (Assembled transcriptomes) removing contigs less than 300 and the output files do not have these contigs.
Is there a way to print an output file for contigs with a minimum and maximum lenght for example for 300 to 600 bp. I need to do this task to separate (filter) my fasta files by contig length: 1) to see how many contigs of a specific lenght are (>300, 300 to 600, 600 to 1200, 1200 to 1800,,, etc) 2)These filtered files i will blast against a Uniprot dataset to see how many of these hit off?
Thanks for your help
You can extend the script above to include a maxlen parameter to extract sequences within a range of lengths, like so:
# filterfastarange.pl
#!/usr/bin/perl
use strict;
use warnings;
my $minlen = shift or die "Error: `minlen` parameter not provided\n";
my $maxlen = shift or die "Error: `maxlen` parameter not provided\n";
{
local $/=">";
while(<>) {
chomp;
next unless /\w/;
s/>$//gs;
my @chunk = split /\n/;
my $header = shift @chunk;
my $seqlen = length join "", @chunk;
print ">$_" if($seqlen >= $minlen and $seqlen <= maxlen);
}
local $/="\n";
}
You would the invoke this script like so:
perl filterfastarange.pl 301 600 contigs.fasta > contigs-gt300-lte600.fasta
Please note, the if() condition above assumes a closed interval [minlen, maxlen], i.e. sequence length is greater than or equal to minlen and less than or equal to maxlen. So your filter intervals should be defined appropriately: [301, 600], [601, 1200], [1201, 1800] and so on, so as to not cause any overlaps between the consecutive sequence sets.
You could easily extend this script to loop through your desired ranges.
Hope this helps!
The above script works for me but there is a problem, i have 109 different fasta files with different file names. i can run the script for individual file but i want to run the script at once for all files and the result file should be individually different for each.
file names are like SRR8224532.fasta, SRR8224533.fasta, SRR8224534.fasta, and so on i want the result files after removing the contigs (i.e., for me less than 1000) something like SRR8224532-out.fasta, SRR8224533-out.fasta, and so on.
Any help or suggestion would be helpfull.
@subrahmanya, you should be able to achieve this via the following loop:
for fasta_file in $(find . -name "SRR*.fasta" -type f); do
output_fasta_file="$(basename $fasta_file .fasta)-out.fasta"
perl removesmalls.pl 1000 $fasta_file > $output_fasta_file
done
In the above code, the find command assumes that all your FASTA files are located in the current directory, as signified by the dot (.)
I recommend you try bioawk, which is a modified version of awk that will parse some common sequence formats. In my opinion, there is no reason to mess with trying to parse files if someone has already figured it out and created some easy to use methods. Also, I would argue against writing a script for simple tasks when a single line of code will do the job. Here is a one-liner to get sequences over 50 bp (will take fasta or fastq):
bioawk -c fastx '{ if(length($seq) > 50) { print ">"$name; print $seq }}' some_seq.fastq
bioawk i tried earlier : I found some problem to install it. Once i am doing like this
bioawk -c fastx '{ if(length($seq) > 1000) { print ">"$name; print $seq }}' list4a.fasta
bioawk: command not found
Then i tried to install bioawk like this
make
make: `awk' is up to date.
I could not solve it. Could you please check it once?
I see, you will need to compile bioawk first, then create a link to awk and name it bioawk. This is not strictly necessary, but I do this so bioawk does not conflict with the system awk (both are named 'awk'). After you type make to compile it, just create a link ln -s awk bioawk and try again. Your shell will not know it's there so you'll have to add it to your PATH, or specify the full path to bioawk.
Here is a PHP solution for your problem using Regular Expression. I have used a limit of length 30 here, you change it as per your wish. The sequences will be in fasta format in seqs.txt file and the resulting screened sequences will be written to screened_seqs.txt:
<?php
$readfile = "seqs.txt";
$writefile = "screened_seqs.txt";
$screenlen = 30;
$fr = fopen($readfile, "r");
$text = fread($fr, filesize($readfile));
$fw = fopen($writefile, "w");
$ntext = preg_replace("/\>[^\n]{1,}\n[ATGCatgc\n]{1,".($screenlen-1)."}\>/", ">", $text);
fwrite($fw, $ntext);
fclose($fr);
fclose($fw);
?>
Hope this helps you.
In the text you give you use typographers quotes (the curly double quotes), which may be stuffing things up also
could it be that your missing the -f between awk and removesmalls. Copying your text and using the following command worked for me:
awk -f removesmalls i=200 test.fa >someoutput.fa
also you should be able to do:
./removesmalls i=200 test.fa >someoutput.fa
awk -f removesmalls i=200 list4a.fasta > list4aout.fasta
i tried with this command but same result it is generating file without any data
Here is the link from where i got the code http://bioinfofly.wordpress.com/2012/07/09/binawk-f-nex/ colud you please ckeck once in your system
Simply copying the code from your question or from the wordpress article produces an error for me. As I say in my answer above, the text in your question contains typographers quotes (also called smart quotes or curly quotes) around the \n in your print statement. You need to remove them and replace them with normal quotes. This can be simply done by deleting those two characters and replacing them with a standard ". With this modification and with the -f change the code you give works for me
Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
version 2.3.6
Hi, you can change the record separator
RSinto>and turn your script into a short one-liner:If you want the minimum size to be an external variable, use the
-voption:Thank you so much