fastq-filter

#!/usr/bin/perl -w

# Copyright © 2015-16, Owen Marshall

# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or (at
# your option) any later version. 
# 
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# General Public License for more details. 
# 
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 
# USA

# Only works on Illumina v1.8 or greater fastq data ...

use strict;
$|++;

my %vars = (
	'trim_ads' => 1,
	'trim_dimers' => 1,
	'Q' => 30,
);

my %vars_details = (
	'trim_ads' => 'Trim DamID adaptors',
	'trim_dimers' => 'Trim DamID adaptor dimers',
	'Q' => 'Trim bases below this minimum quality score',
);

unless (@ARGV) {
	help();
}

my @in_files;
process_cli();

my $min_len = 10;
my $max_len = 40;

my $gatc_re = qr/^GATC/;
my $ad_re = qr/^[G|N]GTCGCGGCCGAGGATC/;
my $dimer_re = qr/^[C|N]TAATACGACTCACTATAGGGCAGCGT/;

foreach my $fn (@ARGV) {
	print STDERR "Reading file $fn ...\n";
	print STDERR "Trimming to Q$vars{'Q'}\n";
	open(IN, "gunzip -c $fn |") || die "ERROR: Cannot open input for reading: $!";
	
	my ($fhead) = $fn =~ m/(.*).fastq/; 
	open(OUT, "| /bin/gzip -c > $fhead.trimmed.fastq.gz") || die "Cannot start gzip process: $!\n";
	
	print STDERR "Writing data to $fhead.trimmed.fastq.gz ...\n";
	
	my $lines=0;
	my ($id, $seq, $desc, $qual, $total, $inc, $len20, $len30);
	$inc = $len20 = $len30 = 0;
	while (<IN>) {
		chomp;
		my $lnum = $lines%4 +1;
		
		# Switch is too slow ...
		$lines++;		
		if ($lnum == 1) {$id = $_; $seq =""; $desc=""; $qual=""; next};
		if ($lnum == 2) {$seq = $_; next};
		if ($lnum == 3) {$desc = $_; next} ;
		$qual = $_;
		
		if ($qual) {
			# full set; process
			$total ++;
			if ($total%20000 == 0) {
				my ($reads_in, $reads_out, $pc, $pc20, $pc30) = format_output($total, $inc, $len20, $len30);
				print STDERR "[ $reads_in"."M in | $reads_out"."M out ($pc%) | $pc20% <20bp | $pc30% <30bp ]\r";
			}
			
			# ignore DamID adapter dimers
			next if ($vars{'trim_dimers'} && $seq =~ $dimer_re);
			
			# process quality
			my $len = length($qual);
			
			# it's *much* faster to process like this than slurp into an array ...
			my ($count, $trim_start, $trim_end);
			$count = $trim_start = $trim_end = 0;
			foreach (split(//,$qual)) {
				$count++;
				my $val = ord("$_")-33; # Illumina v1.8 or greater
				unless ($trim_start) {
					next if $val < $vars{'Q'};
					$trim_start = $count;
				} elsif ($count > 4) {
					$trim_end = $count if $val > $vars{'Q'};
				} else {
					$trim_start = $count+1 if $val < $vars{'Q'};
				}
			}
			
			# Remove DamID adapter seqs if present
			if ($vars{'trim_ads'} && $seq =~ $ad_re) {
				$trim_start = $trim_start > 13 ? $trim_start : 13;
			}
			
			$trim_end = min($trim_end,$max_len);
			
			my $trimmed_len = $trim_end - $trim_start;
			
			next if $trimmed_len < $min_len;
			$len20++ if $trimmed_len >= 20;
			$len30++ if $trimmed_len >= 30;
			
			my $seq_ed = substr($seq, $trim_start-1, $trim_end-$trim_start+1);
			my $qual_ed = substr($qual, $trim_start-1,$trim_end-$trim_start+1);
			
			# output trimmed
			print OUT "$id\n$seq_ed\n$desc\n$qual_ed\n";
			$inc++;
		}
	}
	
	close OUT;
	close IN;
	
	my ($reads_in, $reads_out, $pc, $pc20, $pc30) = format_output($total, $inc, $len20, $len30);
	print STDOUT "[ $reads_in"."M in | $reads_out"."M out ($pc%) | $pc20% <20bp | $pc30% <30bp ]   \n\n"; 
}

sub format_output {
	my ($total, $inc, $len20, $len30) = @_;
	my $reads_in = sprintf("%0.2f",$total/1000000);
	my $reads_out = sprintf("%0.2f",$inc/1000000);
	my $pc = sprintf("%0.2f",$inc/$total*100);
	my $pc20 = sprintf("%0.2f",($inc-$len20)/$total*100);
	my $pc30 = sprintf("%0.2f",($inc-$len30)/$total*100);
	return ($reads_in, $reads_out, $pc, $pc20, $pc30);
}


sub max {
    my ($max, @vars) = @_;
	my $index=0;
	$max||=0;
    for my $i (0..$#vars) {
        ($max, $index) = ($vars[$i], $i+1) if $vars[$i] > $max;
    }
    return ($index, $max);
}

sub min {
    my ($min, @vars) = @_;
	my $index=0;
	$min||=0;
    for my $i (0..$#vars) {
        ($min, $index) = ($vars[$i],$i+1) if $vars[$i] < $min;
    }
    return ($index, $min);
}


sub process_cli {
	foreach (@ARGV) {
		if (/--(.*)=(.*)/) {
			unless (defined($vars{$1})) {
				print STDERR "Did not understand $_ ...\n";
				help();
			}
			my ($v, $opt) = ($1,$2);
			$vars{$v} = $opt;
			next;
		} elsif (/--h[elp]*/) {
			help();
		} elsif (/--(.*)/) {
			if (defined($vars{$1})) {
				$vars{$1} = 1;
			} else {
				print STDERR "Did not understand $_ ...\n";
				help();
			}
			next;
		}
		push @in_files, $_;
	}
}


sub help {
	print STDOUT <<EOT;

fastq-filter: removes damid adaptor sequences and trims to a minimum read quality
	
Options:
EOT
	
	my $opt_len = 0;
	foreach (keys %vars) {
		my $l = length($_);
		$opt_len = $l if $l > $opt_len;
	}
	
	$opt_len+=2;
	
	my $cols= `tput cols` || 80;
	
	my ($v, $val, $def, $def_format);
	my $help_format = "format STDOUT =\n"
		.' '.'^'.'<'x$opt_len . ' '. '^' . '<'x($cols-$opt_len-4) . "\n"
		.'$v, $def_format'."\n"
		.' '.'^'.'<'x$opt_len . '   '. '^' . '<'x($cols-$opt_len-6) . "~~\n"
		.'$v, $def_format'."\n"
		.".\n";
		
	eval $help_format;
	die $@ if $@;
	
	foreach my $k (sort (keys %vars)) {
		($v, $val, $def) = ($k, $vars{$k}, $vars_details{$k});
		$def||="";
		$def_format = $val ? "$def\n\r[Current value: $val]" : $def;
		$v = "--$v";
#		format =
# ^<<<<<<<<<<<<<<<<<<<< ^<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
#$v, $def_format
# ^<<<<<<<<<<<<<<<<<<<<   ^<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ~~
#$v, $def_format
#.

		write();
		
	}
	print STDOUT "\n";
	exit 1;
}