###################################################################################
#	Written by Justin Gardin November 2014
#	Bruce Futcher Lab, Stony Brook University
#	
#	Code accepts paired end fastq files and splits the files by 5' barcodes
#
###################################################################################
#!/usr/bin/perl

use strict; use warnings;
use Getopt::Long;
use Symbol;

my $BC_filename='';
my $read_1_filename='';
my $read_2_filename='';
my $mismatches_allowed=0;
my $barcode_length=7;
my @barcodes;
my $suffix='out.fastq';
my @prefixes;

GetOptions( 	"b=s" => \$BC_filename,
		"l=i" => \$barcode_length,
		"m=i" => \$mismatches_allowed,
		"o=s" => \$suffix,
		"1=s" => \$read_1_filename,
		"2=s" => \$read_2_filename
);


print "\nStarting barcode seperation: Loading files\n\n";
if($BC_filename eq ''|| !(open(BCFILE,'<',$BC_filename))){die "Please provide a valid file name for the barcode file with the -b option\n";
}
my $count=0;
while(my $line=<BCFILE>){
	chomp $line;
	my @a=split(/\t/,$line);
	my @bc_string=split(//,$a[1]);
	$prefixes[$count]=$a[0];
	$barcodes[$count]= \@bc_string;
	++$count;
}
close BCFILE;

if($read_1_filename eq '' || !(open(READ1,'<',$read_1_filename))) {die "Please provide a valid file name for the read 2 file with the -1 option\n"; }
if($read_2_filename eq '' || !(open(READ2,'<',$read_2_filename))) {die "Please provide a valid file name for the read 2 file with the -2 option\n"; }

my @file_handles_r1;
my @file_handles_r2;
for(my $i=0;$i<@barcodes;++$i){
	my $fh1=gensym;
	my $filename1=$prefixes[$i]."_read1_".$suffix;
	open($fh1,'>',$filename1)||die "Cannot open $filename1\n";
	push(@file_handles_r1,$fh1);
	my $fh2=gensym;
	my $filename2=$prefixes[$i]."_read2_".$suffix;
	open($fh2,'>',$filename2)||die "Cannot open $filename2\n";
	push(@file_handles_r2,$fh2);
}

print "Files loaded successfully\n\nParsing Read files\n\n";

my $read_discarded=0;
while(my $a1=<READ1>){
	my $a2=<READ2>;
	my $b1=<READ1>; my $b2=<READ2>;
	my $c1=<READ1>; my $c2=<READ2>;
	my $d1=<READ1>; my $d2=<READ2>;
	$a1=~/(.*?)\s.*?$/;
	my $loc1=$1;	
	$a2=~/(.*?)\s.*?$/;
	my $loc2=$1;
	unless($loc1 eq $loc2){
		die "Fastq read 1 or Fastq read 2 files appear to be unsorted. Please sort the Fastq files first\n$a1 does not match $a2\n"; 
	}
	my @bc_c1=split(//,substr($b1,0,$barcode_length));
	my @bc_c2=split(//,substr($b2,0,$barcode_length));
	my @mismatches_r1=((0)x$count);
	my @mismatches_r2=((0)x$count);
	for(my $i=0;$i<$count;++$i){
		for(my $j=0;$j<$barcode_length;++$j){
			++$mismatches_r1[$i] unless ($bc_c1[$j] eq $barcodes[$i][$j]);
			++$mismatches_r2[$i] unless ($bc_c2[$j] eq $barcodes[$i][$j]);
		}
	}
	my $num_accepted_r1=0;
	my $match_r1;
	for(my $i=0;$i<$count;++$i){
		if($mismatches_r1[$i] <= $mismatches_allowed){
			++$num_accepted_r1;
			$match_r1=$i;
		}
	}
	if($num_accepted_r1 == 1){
		my $handle1=$file_handles_r1[$match_r1];
		my $handle2=$file_handles_r2[$match_r1];
		print $handle1 $a1,$b1,$c1,$d1;
		print $handle2 $a2,$b2,$c2,$d2;
	}
	else{
		++$read_discarded;
	}
}
		
print "Finished: $read_discarded reads could not be matched to a barcode\n";
close READ1;
close READ2;