#! usr/bin/perl

#2009/1/2

sub usage {
    print STDERR "usage: combines forward and reverse sequences  <standard fasta alignment containing forward and reverse sequences>\n\n";
    print STDERR "#program reads in full alignment including forward and reverse sequences and outputs only one sequence per sample\n";
    print STDERR "Creator: Laurie Stevison\n";
    exit;
} #end sub

if ($#ARGV  != 0) { # zero means one argument
    &usage;
}

$fasta = $ARGV[0];

if($fasta =~ /\.fas/) {
    $fasta_name = $`;
    $output = $fasta_name . "_consensus.fas";
} else {
    print "File not saved in fasta format. Please resave SeqMan output and try conversion again. Need \".fas\" extension.\n";
    exit;
} #end if

open(READFAS, $fasta);
open(WRITEOUT, ">$output");

print STDERR "Reading input file...\n";

while (<READFAS>) {
    chomp;
    if (/>/) {
        $seq_name = $'; #'
        push(@seq_names, $seq_name);
        if ($newseq==1) {push(@sequences, uc($sequence));}
        $sequence = "";
        $newseq=0;
        next;
    } else {
        $sequence .= $_;
        $newseq=1;
} #end else
} #end while

push(@sequences, $sequence);

print STDOUT "Please indicate labeling for forward sequences...\n";
chomp($forward_seqs = (<STDIN>));

print STDOUT "Please indicate labeling for reverse sequences...\n";
chomp($reverse_seqs = (<STDIN>));

print STDOUT "Please indicate maximum length of final sequence.\n";
chomp($seq_length = (<STDIN>));

for ($i = 0; $i <= $#seq_names; $i++) {

    $forward_sequence="";
    $seq_name = $seq_names[$i];
    $forward_sequence = $sequences[$i];
    @forward = ();
    @forward = split("", $forward_sequence);

    if($seq_name =~ /$forward_seqs/) {
	@sequence_title = split($forward_seqs, $seq_name);
    } else {
	next;
    } #end if                           

    $sample_name = $sequence_title[0];
    
    for  ($k=0; $k <= $#seq_names; $k++) {
	$rev_seq_name = $seq_names[$k];
        if($seq_names[$k] =~ /$reverse_seqs/) {
	    @reverse_sequence_title = split($reverse_seqs, $rev_seq_name);
	} else {
	    next;
	} #end if

	$reverse_sequence="";
	$reverse_sequence = $sequences[$k];
	@reverse = ();
	@reverse = split("", $reverse_sequence);

	if ($reverse_sequence_title[0] eq $sample_name) {
	    #print STDOUT "Reverse sequence name matches sample name\n";
	    @finalseqarray=();
	    for($j=0; $j <= $seq_length; $j++) {
		#print STDOUT "made into j for loop\n";
		if($forward[$j] eq $reverse[$j]) {
		    push (@finalseqarray, $forward[$j]);
		} elsif($forward[$j] eq undef) { 
		    push (@finalseqarray, lc($reverse[$j]));
		} elsif($reverse[$j] eq undef) {
		    push (@finalseqarray, lc($forward[$j]));
		} elsif($forward[$j]=~/-/) { 
		    push (@finalseqarray, lc($reverse[$j]));
		} elsif($reverse[$j]=~/-/) {
		    push (@finalseqarray, lc($forward[$j]));
		} elsif($forward[$j]=~/N/) { 
		    push (@finalseqarray, lc($reverse[$j]));
		} elsif($reverse[$j]=~/N/) {
		    push (@finalseqarray, lc($forward[$j]));
		} elsif  ($forward[$j]=~/C|A/ && $reverse[$j]=~/C|A/) {
		    push (@finalseqarray,M);
	        } elsif  ($forward[$j]=~/G|A/ && $reverse[$j]=~/G|A/) {
		    push (@finalseqarray,R);
        	} elsif  ($forward[$j]=~/T|A/ && $reverse[$j]=~/T|A/) {
		    push (@finalseqarray,W);
		} elsif  ($forward[$j]=~/C|G/ && $reverse[$j]=~/C|G/) {
		    push (@finalseqarray,S);
	        } elsif  ($forward[$j]=~/C|T/ && $reverse[$j]=~/C|T/) {
		    push (@finalseqarray,Y);
	        } elsif  ($forward[$j]=~/G|T/ && $reverse[$j]=~/G|T/) {
		    push (@finalseqarray,K);
		} #end if
	    } #end j for
	} else {
	    next;
        } #end if
    } #end k for
    
    $final_sequence = join("",@finalseqarray);
    $sample_name=~s/\'//;
    print WRITEOUT ">$sample_name" . "_consensus\n";
    print WRITEOUT "$final_sequence\n";

} #end i for

print STDERR "\n\nConversion done.\n\n";