Removing :Y: flagged reads

[lukas1848]

Hi,
what is the best (i.e. quickest) tool to remove :Y: flagged reads from illumina data?

So far I used a perl script which is not very memory effective and I'd like to speed up this step.

[Rocketknight]

You could do it in Python (code at the bottom). It takes one argument, the path to the fastq file to process.

Since python only stores one line in memory at once with this implementation, it should be pretty efficient. You can save time by piping the output of this file directly to the downstream program, using a line like this:

readmapper <(/path/to/removebadreads.py /path/to/reads.fastq)

This will call program 'readmapper' and pass it the output from the python script below as it runs, so you don't need to waste time and disc space making a huge output file. Modify as required for whatever program you want to run yourself.

Also, this program won't work with gzipped fastq files, but you can use the same piping concept to call gunzip and pass unzipped output straight to the script:

./removebadreads.py <(gunzip -c /path/to/reads.fastq.gz)

You can even pipe the gunzip output to removebadreads.py and then pipe the output from that to your read mapper, all at once with no temporary files or delays! Linux is fun!

readmapper <(/path/to/removebadreads.py <(gunzip -c /path/to/reads.fastq.gz))

Python code:

#!/usr/bin/python
import sys

try:
	f = open(sys.argv[1])
except:
	exit('Supply an input file!')

linecounter = 0
goodread = False
for line in f:
	linecounter = 1 if linecounter == 4 else linecounter+1
	if linecounter == 1:
		if ':N:' in line:
			goodread = True
			print line,
		else:
			goodread = False
	elif goodread:
		print line,

[swbarnes2]

I always used zgrep. It did put '--' between the entries, but you can get rid of that, and it didn't seem to interfere with most programs anyway.

[HESmith]

Here's a Perl script that's quick and efficient (courtesy of James Platt):

#!/usr/bin perl

use warnings;

# usage: perl Illumina_fastq_filter-JP.pl filename > outname
my $n = 0; #line counter

open (FILENAME, "$ARGV[0]") or die "Can't find $ARGV[0]\n";

while (<FILENAME>) {
	chomp;

	if ($_ =~ /^@\S+\s\d\:N/) {
		$n = 4;
	}

	if ($n > 0) {
		print "$_\n";
		$n--;
	}
}

close FILENAME;