#!/usr/bin/env perl

#gsprimer ver. 0.33
#
#This program is free software; you can redistribute it and/or
#modify it under the terms of the GNU General Public License.
#
#Naoki Shinoda
#

use warnings;
use strict;
use Data::Dumper;
use File::Copy;
use Getopt::Long qw(:config gnu_compat);

sub usage() {
	print "Usage: $0 [options]\n",
		"       $0 [options] <-a file> <-o outputfile>\n",
		"       $0 [options] <-s file> <-n file> <-o outputfile>\n",
		"Options:\n",
		"  -a <file>, --align=<file>          Alignment file (made by this program)\n",
		"  -t <file>, --target=<file>         target sequences file\n",
		"  -n <file>, --nontarget=<file>      non-target sequences file\n",
		"  -o <file>, --outfile=<file>        Output file (required)\n",
		"  -tm <int>,  --tm=<int>             Tm value (in integer)\n",
		"  -z <int>,  --zeronum=<int>         Allowed # of score zero species\n",
		"  -th <int>,  --threshold=<int>      Threshold of score\n",
		"  -dc <int>,  --discord=<int>        Allowed discordance(loose:1---3:strict)\n",
		"  -dm,  --dimer                      Modify homo-dimer\n"
}

our %opt;
my $dm = 0;
my $optcheck = GetOptions(\%opt,
		'align|a=s',
		'target|t=s',
		'nontarget|n=s',
		'outfile|o=s',
		'tm=i',
		'zeronum|z=i',
		'threshold|th=i',
		'discord|dc=i',
		'dimer|dm'=> \$dm
	);

if (!$optcheck) {
	&usage();
	exit 1;
}


print "Searches species(group)-specific primers.\n";

#
#If melting temp is not specify, 55.
#

if (!exists($opt{tm})) {
	$opt{tm} = 55;
}
if (!exists($opt{threshold})) {
	$opt{threshold} = 0;
}
if (!exists($opt{zeronum})) {
	$opt{zeronum} = 1;
}
if (!exists($opt{discord})) {
	$opt{discord} = 2;
}


#
#discordance between target sequences
#(1:loose --- strict: 3)
#

my $discord1;
my $discord2;

if ($opt{discord} == 1) {
	$discord1 = 3;
	$discord2 = 7;
}
elsif ($opt{discord} == 2) {
	$discord1 = 5;
	$discord2 = 4;
}
elsif ($opt{discord} == 3) {
	$discord1 = 7;
	$discord2 = 1;
}
else {
	print "-d # should be 1-3\n";
}



if (exists($opt{align}) && (exists($opt{target}) || exists($opt{nontarget}))) {
	print STDERR "Eather align file and seq file ware specified. Aborted.\n";
	&usage();
	exit 1;
}

	my $input_spfilename = $opt{target};
	my $input_othfilename = $opt{nontarget};
	my $output_filename = $opt{outfile};

if (!exists($opt{align})) {
	unless ($input_spfilename) {
		print "Input target species sequences: ";
		$input_spfilename = <>;
		chomp($input_spfilename);
		unless ($input_spfilename) {
			print "Aborted.\n";
			exit 1;
		}
	}
	unless (-f $input_spfilename) {
		print STDERR "target sequence file '$input_spfilename' not found\n";
		exit 1;
	}

	unless ($input_othfilename) {
		print "Input non-target species sequences: ";
		$input_othfilename = <>;
		chomp($input_othfilename);
		unless ($input_othfilename) {
			print "Aborted.\n";
			exit 1;
		}
	}
	unless (-f $input_othfilename) {
		print STDERR "non-target sequence file '$input_othfilename' not found\n";
		exit 1;
	}

	(my $df_outname = $input_spfilename) =~ s/(^.+)\.[^\.]+$/result\_$1\.txt/;
	unless ($output_filename) {
		print "output filename [$df_outname] : ";
		$output_filename = <>;
		chomp($output_filename);
		unless ($output_filename) {
			$output_filename = $df_outname;
		}
	}

	#count ">" : sequence number
	open(INPUT_FH, $input_spfilename);
	my $pnum = 0;
	while (my $line = <INPUT_FH>) {
		chomp($line);
		if (substr($line, 0, 1) eq ">") {
			$pnum++;
			next;
		}
	}
	close(INPUT_FH);

	if ( $pnum == 1 ) {
		open(FH, "$input_spfilename");
		open(CONS_FH, ">cons_est.fasta");
		open(PR_FH, ">cons_ref.fasta");

		print CONS_FH ">cons_est\n";
		print PR_FH ">cons_ref\n";

		#overwrite from ">" to ">cons_est"
		while (my $line = <FH>) {
			if (substr($line, 0, 1) ne ">") {
				print CONS_FH $line;
				print PR_FH $line;
				next;
			}
		}
		close(FH);
		close(CONS_FH);
		close(PR_FH);

	}
	else {

		#clustalw cons_est
		system("emma $input_spfilename -dendoutfile _dndfile0 -outseq consensus.aln");

		$pnum = $pnum - 1;

		#cons
		system("cons consensus.aln -outseq cons_est.fasta -plurality $pnum -name cons_est");
		system("cons consensus.aln -outseq cons_ref.fasta -plurality 0 -name cons_ref")

	}

	#sqseq + otherseq
	system("cat cons_est.fasta cons_ref.fasta $input_othfilename >emboss_tmp0");

	(my $aln = $input_spfilename) =~ s/(^.+)\.[^\.]+$/gsprimer\_$1\.aln/;
	#clustalw sqseq + otherseq
	system("emma emboss_tmp0 -dendoutfile _dndfile1 -outseq $aln");

	open(INPUT_FH, "$aln");

	print "Alighnment file '$aln' are created."

}
else {
	my $input_alnfilename = $opt{align};

	unless (-f $input_alnfilename) {
		print STDERR "Alignment file '$input_alnfilename' not found\n";
		exit 1;
	}

	(my $df_outname = $input_alnfilename) =~ s/(^.+)\.[^\.]+$/result\_$1\.txt/;
	unless ($output_filename) {
		print "output filename [$df_outname] : ";
		$output_filename = <>;
		chomp($output_filename);
		unless ($output_filename) {
			$output_filename = $df_outname;
		}
	}


	open(INPUT_FH, $input_alnfilename);
}

#
# read data to arrays
#

my %seq_data;
my $current_seq_name;
my $current_seq_idx = -1;
while (my $line = <INPUT_FH>) {
	chomp($line);
	if (substr($line, 0, 1) eq ">") { 
		# when a character of "pos 0, length 1"(i.e. lead-off character) is ">"
		$current_seq_name = substr($line, 1);
		unless ($current_seq_name) { next; }
		$seq_data{$current_seq_name} = "";
		next;
	}
	unless (defined($current_seq_name)) {
		print STDERR "Skipping Invalid Input line: $line\n";
		next;
	}

	$seq_data{$current_seq_name} .= $line;
}
close(INPUT_FH);

# check data
#print Dumper(\%seq_data);

#
#reference of tm value
#
our %check_tm = (
	'AA' => { 'h'=>-9.1, 's'=>-24.0 }, 
	'TT' => { 'h'=>-9.1, 's'=>-24.0 },
	'AT' => { 'h'=>-8.6, 's'=>-23.9 }, 
	'TA' => { 'h'=>-6.0, 's'=>-16.9 }, 
	'CA' => { 'h'=>-5.8, 's'=>-12.9 },
	'TG' => { 'h'=>-5.8, 's'=>-12.9 }, 
	'GT' => { 'h'=>-6.5, 's'=>-17.3 },
	'AC' => { 'h'=>-6.5, 's'=>-17.3 },
	'CT' => { 'h'=>-7.8, 's'=>-20.8 },
	'AG' => { 'h'=>-7.8, 's'=>-20.8 }, 
	'GA' => { 'h'=>-5.6, 's'=>-13.5 },
	'TC' => { 'h'=>-5.6, 's'=>-13.5 }, 
	'CG' => { 'h'=>-11.9, 's'=>-27.8 },
	'GC' => { 'h'=>-11.1, 's'=>-26.7 }, 
	'GG' => { 'h'=>-11.0, 's'=>-26.6 },
	'CC' => { 'h'=>-11.0, 's'=>-26.6 }
);


#
# scan
#
my $last_pos = length($seq_data{cons_est});

open(OUT_FH, ">$output_filename");
print OUT_FH "Forward primer\n";

my $current_pos = 0;
while($current_pos < $last_pos - 25) {
	&check_primer_cand(\%seq_data, $current_pos, 50);
	$current_pos++;
}

print OUT_FH "\nReverse primer\n\n";

$current_pos = 0;
while($current_pos < $last_pos - 25) {
	&check_primer_cand(\%seq_data, $current_pos, 50, 1);
	$current_pos++;
}

sub check_primer_cand($$$$) {
	my $seq_data_ref = $_[0]; # 뵉썽ㅞ과퓖E	my $pos = $_[1];
	my $len = $_[2];
	my $reverse_flag = $_[3];

	#print "pos ===> $pos\n";

	my @score_def = (8, 5, 3, 2, 1);
	#print Dumper(\%check_dat_ref);

	my %check_seqs;
	my $cur_seq_key;
	foreach $cur_seq_key (keys %$seq_data_ref) {
		$check_seqs{$cur_seq_key} = substr($seq_data_ref->{$cur_seq_key}, $pos, $len);

		if (!$reverse_flag) {
			$check_seqs{$cur_seq_key} = reverse($check_seqs{$cur_seq_key});
		}

		$check_seqs{$cur_seq_key} = {"str" => $check_seqs{$cur_seq_key},
			"score" => 0,};
		}
	#print Dumper(\%check_seqs);

	my $primer_cand = $check_seqs{cons_est}->{str};
	my $primer = $check_seqs{cons_ref}->{str};
	delete $check_seqs{cons_est};
	delete $check_seqs{cons_ref};

	if (substr($primer_cand, 0, 1) eq '-') {
		# When 3' terminal is "-", return.
		return "";
	}

	if (substr($primer_cand, 0, $discord1) =~ tr/Nn// ) {
		# When "N" is a lot in "0--$discord1", return
		return "";
	}

	#
	#When "N" is a lot in "$discord1--17", return
	#
	my $countn = substr($primer_cand, $discord1, 17 - $discord1) =~ s/n/N/gi;
	if ($countn >= $discord2) {
		return "";
	}

#
#calculate scores
#

	for (my $p = 0; $p < @score_def; $p++) {
		my $cur_base = substr($primer_cand, $p, 1);
		if ($cur_base eq 'N') {
			next;
		}

		foreach my $cur_seq_key (keys %check_seqs) {
			my $cur_base_oth = substr($check_seqs{$cur_seq_key}->{str}, $p, 1);
			if ($cur_base ne $cur_base_oth) {
				$check_seqs{$cur_seq_key}->{score} += $score_def[$p];
				if ($cur_base eq 'G' || $cur_base eq 'C') {
					$check_seqs{$cur_seq_key}->{score} += 0.5 * $score_def[$p];
				}
			}

		}
	}

	my $max_def = 0;
	foreach (@score_def) { $max_def += $_ };
	my $min_score = $max_def;

	my $zeronum = 0;
	my $score_2nd = $min_score;

	foreach $cur_seq_key (keys %check_seqs) {
		if ($check_seqs{$cur_seq_key}->{score} == 0) {
			$zeronum += 1;
		}
	       	if ($check_seqs{$cur_seq_key}->{score} > 0
			&& $check_seqs{$cur_seq_key}->{score} < $score_2nd) {
			$score_2nd = $check_seqs{$cur_seq_key}->{score};
		}
		if ($check_seqs{$cur_seq_key}->{score} < $min_score){
			$min_score = $check_seqs{$cur_seq_key}->{score};
		}
	}

	if ($min_score < $opt{threshold}) {
		return "";
	}

	#print "min: $min_score\n";
	if ($zeronum > $opt{zeronum}) {
		return "";
	}

	#
	#erase "-"
	#
	$primer =~ tr/-//d;
	$primer = substr($primer, 0, 25);
	#
	#modify dimer
	#
	my $primpos = $pos;
	if ($dm) {

		my $checkdm1 = substr($primer, 0, 3);
		$checkdm1 =~ tr/ATGC/TACG/;
		my $checkdm0 = $checkdm1;

		my $check_mer = substr($checkdm1, 0, 1);
		if ($check_mer eq "A") {
		substr($checkdm1, 0, 1) = "T";
		}
		elsif ($check_mer eq "T") {
			substr($checkdm1, 0, 1) = "A";
		}
		elsif ($check_mer eq "G") {
			substr($checkdm1, 0, 1) = "T";
		}
		elsif ($check_mer eq "C") {
			substr($checkdm1, 0, 1) = "A";
		}

		if (!$reverse_flag) {
			$primer = reverse($primer);
		}
		else {
			$checkdm0 = reverse($checkdm0);
			$checkdm1 = reverse($checkdm1);
		}
		while ($primer =~ s/$checkdm0/$checkdm1/g != 0){};
		if (!$reverse_flag) {
			$primer = reverse($primer);
		}
	}

	printf OUT_FH "score:%2d", $min_score;
	if ($min_score == 0) {
		if ($zeronum > 1) {
			print OUT_FH "($score_2nd:$zeronum)";
		}
		else {
			print OUT_FH "($score_2nd)";
		}
	}
	if (!$reverse_flag) {
		$primpos = $primpos + 50;
	}


	#
	#calculate Tm value
	#
	my $tm_old = 0;
	my $tm_cur = 0;
	my $tm_pos;
	my $primer_length = length($primer) - 2;
	for ($tm_pos = 15; $tm_pos <= $primer_length; $tm_pos++) {

		my $tm_h = 0;
		my $tm_s = 0;
		my $key;
		for (my $p = 0; $p <= $tm_pos; $p++) {
			$key = substr($primer, $p, 2);
			if (!$reverse_flag) {
				$key = reverse($key);
			}
			if ( !exists($check_tm{$key}) ) {
				next;
			}
			$tm_h = $tm_h + $check_tm{$key}->{'h'};
			$tm_s = $tm_s + $check_tm{$key}->{'s'};
		}

		$tm_old = $tm_cur;
		$tm_cur = 1000 * $tm_h / ($tm_s - 42.3833) - 294.747;
		if ($tm_cur >= $opt{tm}) {
			last;
		}
	}


	my $tm;
	if ($tm_cur < $opt{tm}) {
		$tm = $tm_cur;
		$tm_pos = length($primer);
	}
	elsif ($tm_cur - $opt{tm} > $opt{tm} - $tm_old) {
		$tm = $tm_old;
		$tm_pos += 1;
	}
	else {
		$tm = $tm_cur;
		$tm_pos += 2;
	}
	$primer = substr($primer, 0, $tm_pos);

	printf OUT_FH ("\tTm:%2.2f", $tm);

	print OUT_FH "\tpos:$primpos ";
	if ($reverse_flag) {
		print OUT_FH $primer;
		my $rprimer = reverse($primer);
		$rprimer =~ tr/ATGC/TACG/;
		print OUT_FH " ($rprimer)";
	}
	else {
		$primer = reverse($primer);
		print OUT_FH $primer;
	}
	print OUT_FH "\n";
}
close(OUT_FH);
print "\n";
