This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Filtering Fastq Sequences Based On Lengths

As the question says,

I have a fastq file from small RNA sequencing with sequence lengths between 15 - 30. I wanted to filter sequence lengths between 21-25 and write to another file. how can i do that?

awk perl unix
• 35,853 views • link
updated by komjinhubuio • komjinhubuio • written by empyrean999 empyrean999

5 answers 

cat your.fastq | paste - - - - | awk 'length($2)  >= 21 && length($2) <= 25' | sed 's/\t/\n/g' > filtered.fastq
• 323 views • link
written by Wen.Huang Wen.Huang

If I may, a pure Awk command is twice faster:

awk 'BEGIN {FS = "\t" ; OFS = "\n"} {header = $0 ; getline seq ; getline qheader ; getline qseq ; if (length(seq) >= 21 && length(seq) <= 25) {print header, seq, qheader, qseq}}' < your.fastq > filtered.fastq

Like your solution with paste, it assumes that a fastq record takes exactly 4 lines.

Edit: deal with spaces in sequence names, as suggested by brianpenghe.

• 0 views • link
written by Frédéric Mahé Frédéric Mahé

How would this be modified for gzip's fastqs?

• 1 views • link
written by shanasabri • shanasabri

zcat decompresses the data of all the input files, and writes the result on the standard output. zcat concatenates the data in the same way cat do

zcat your.fastq.gz | ...
• 1 views • link
written by Medhat Medhat

....................removed....................

• 0 views • link
written by shanasabri • shanasabri

Like, a normal for loop?

for fastq in *.fastq
do
awk ... < $fastq > filtered_$fastq
done
• 0 views • link
written by WouterDeCoster WouterDeCoster

One note: This command doesn't work when the read names contain spaces. better use awk -F"\t" instead of awk

• 0 views • link
written by brianpenghe brianpenghe

This is just printing zero lines! Note I changed constraints to either >=16 only, or >=16 && <= 500.

• 0 views • link
written by caverill • caverill

I would like to introduce you to a powerful software: seqkit (https://bioinf.shenwei.me/seqkit/usage/), with which you can easily manipulate fastq/fasta format seq. seqkit seq youseq.fastq -m 21 -M 25 > result.fq

• 0 views • link
written by komjinhubuio • komjinhubuio

Using Biopieces www.biopieces.org)

read_fastq -i in.fq | grab -e 'SEQ_LEN>=21' | grab -e 'SEQ_LEN<=25' | write_fastq -o out.fq -x

And when you realize that you want to do a lot of extra things besides filtering on sequence length you will find lots of useful tools in Biopieces.

• 809 views • link
written by Martin A Hansen Martin A Hansen

Edit: If you're coming via Google, my answer is very, very old. Consider komjinhubuio's answer, seqkit is the 'modern' way to go

Using the awesome readfq-library in perl, and their modified example:

  my @aux = undef; # this is for keeping intermediate data
  while (my ($name, $seq, $qual) = readfq(\*STDIN, \@aux)) { 
     if( (length($seq) >= 21) && (length($seq) <= 25) ) { 
         print "@$name\n";
         print "$seq\n"; 
         print "+\n";
         print "$qual\n";
     }
  }

(Beware: Haven't tested this yet)

• 0 views • link
written by Philipp Bayer Philipp Bayer

You are missing a closing bracket, a "n" in the last line (there's also an unnecessary comma), and use strict; use warnings; which will tell you such things :)

• 303 views • link
written by SES SES

:) Thank you! I usually never use Perl.

• 0 views • link
written by Philipp Bayer Philipp Bayer

You can easily do this with prinseq-lite:

FILTER OPTIONS 

-min_len <integer>
        Filter sequence shorter than min_len.

-max_len <integer>
        Filter sequence longer than max_len.

prinseq-lite.pl -fastq yourfile.fastq -out_format 4 -out_good seqs_good -min_len 21 -trim_to_len 25

http://prinseq.sourceforge.net/manual.html

• 1 views • link
written by etienne.rifa • etienne.rifa

Log in to answer this question.

More posts like this