#!/usr/bin/perl -w

###########################################################
# redup.pl
# John Garbe
# August 2012
#
# Remove exact duplicate reads from paired-end fastq files, printing out
# unique reads to .unique files. Prints out top N most duplicated sequences 
#
###########################################################

=head1 NAME

redup.pl - Remove exact duplicate reads from paired-end fastq files, printing out unique reads to .unique files. Prints out top N most duplicated sequences in fasta format

=head1 SYNOPSIS

redup.pl [-n N] sample1_R1.fastq sample1_R2.fastq > topdups.fasta

=head1 DESCRIPTION

This script removes duplicate paired-end reads from the input files sample1_R1.fastq and sample1_R2.fastq and prints out unique reads to the files sample1_R1.fastq.unique and sample2_R2.fastq.unique. Reads must have the exact same sequence to be called duplicates, quality scores are ignored. The top N (default 20) most duplicated sequences are printed out in fasta format, making it convenient for using BLAST to identify them.

=cut

use Getopt::Std;

# determine number of most duplicated sequences to return
$opt_n = 20;
$ok = getopts('n:');
$k = $opt_n - 1; # subtract one to work with Perl's zero-based arrays

die "USAGE: redup.pl [-n N] sample1_R1.fastq sample1_R2.fastq\n" unless ((($#ARGV == 1) || ($#ARGV == 3)) && $ok);

# open up input files
open F1, "<$ARGV[0]" or die "cannot open $ARGV[0]\n";
open F2, "<$ARGV[1]" or die "cannot open $ARGV[1]\n";

# open up output files
open O1, ">$ARGV[0].unique" or die "cannot open $ARGV[0].unique\n";
open O2, ">$ARGV[1].unique" or die "cannot open $ARGV[1].unique\n";

# open up output files
if ($#ARGV < 3) {
    open O1, ">$ARGV[0].unique" or die "cannot open $ARGV[0].unique\n";
    open O2, ">$ARGV[1].unique" or die "cannot open $ARGV[1].unique\n";
} else {
    open O1, ">$ARGV[2]" or die "cannot open $ARGV[2]\n";
    open O2, ">$ARGV[3]" or die "cannot open $ARGV[3]\n";
}


# run through the input files
$readcount = 0;
$dupcount = 0;
while ($f1line1 = <F1>) {
    $readcount++;
    if ($readcount % 1000000 == 0) { # print progress update
	print STDERR "$readcount reads processed, $dupcount duplicates removed\n";
    }
    $f1line2 = <F1>;
    $f1line3 = <F1>;
    $f1line4 = <F1>;

    $f2line1 = <F2>;
    $f2line2 = <F2>;
    $f2line3 = <F2>;
    $f2line4 = <F2>;

    $combo = $f1line2 . $f2line2; # concatenate seqs from reads 1 and 2
    if (defined($seen{$combo})) { # if seq is a dup skip it
	$seen{$combo}++;
	$dupcount++;
	next;
    } else { # seq is not a dup, print reads back out
	$seen{$combo} = 1;
	print O1 "$f1line1";
	print O1 "$f1line2";
	print O1 "$f1line3";
	print O1 "$f1line4";

	print O2 "$f2line1";
	print O2 "$f2line2";
	print O2 "$f2line3";
	print O2 "$f2line4";
    }
}

close O1;
close O2;
print STDERR "$readcount reads processed, $dupcount duplicates removed\n";

print STDERR "Identifying most frequent sequences (kill program to skip this step)\n";

### print out the top k most duplicated reads in fasta format
# It would be better to use a proper partial sort algorithm here

# initialize first k elements
for $i (0..$k) {
    $result[$i][0] = "seed";
    $result[$i][1] = -1;
}

# find reads with more duplicates than least duplicated sequence in results
for $key (keys %seen) {

    if ($seen{$key} > $result[$#result][1]) {
	push @result, [ ( $key, $seen{$key} ) ];
	@result = sort { $b->[1] <=> $a->[1] } @result;
	$#result = $k;
    }
}

# print out the top k most duplicated reads in fasta format
for $i (0..$k) {
    next if ($result[$i][1] < 0);
    $rank = $i + 1;
    print ">Rank:$rank:Count:$result[$i][1]\n";
    ($r1, $r2) = split /\n/, $result[$i][0];
    print $r1 . "N" . "$r2\n";
}

# print out total again in case it was missed
print STDERR "$readcount reads processed, $dupcount duplicates removed\n";

exit;