#!/usr/bin/perl -w
use strict;

########################################################################
# this script will run three tests to see whether components of a unigene
# are flipped:
# 	1. check for C/N in database for contigs, EST orientation of singletons
# 	2. BLAST against arabidopsis
# 	3. check for Poly T 'tails' at the start of sequences
# this data will be assembled into a tab delimited file as such:
# 	<est>	<0|1 - test 1>	<0|1 - test 2>	<0|1 - test 3>
########################################################################

use db_link;
use projects;
use runtime;

my $blast_program = "blastall";
my $ath1_loc = "/data/local/blast/ath1/ATH1_pep";

my @arg_pairs = split (/\-/, (join ' ', @ARGV));
my %args=();
my @unigene_build_ids;
foreach (@arg_pairs){
	$_ or next;
	my ($flag, @val)=split /\s+/;
	if ($flag eq 'b') {
		@unigene_build_ids=@val;
	}
	$args{$flag}=$val[0];
}
my $project=$args{'p'};
my $set=$args{'b'};

$project or die "You must set the project with the -p flag.\n";
$set or die "You must specify the unigene builds with the -b flag.\n";

my ($db,$usr) = @{projects::get_db_info($project)};
$db or die "No known database for project $project";

my $dbh = db_link::connect_db($db,$usr) or die "couldn't open database link!\n";

my $total_time = time;
my $likely_flipped;

foreach my $id (@unigene_build_ids) {
	my $outfile = "${project}_build_${id}_flipped_components.txt";

########################################################################
# gat data
########################################################################

	my %ests = ();
	my ($est, $assembly_orientation, $read_direction, $trimmed_sequence);
	my ($stm, $sth, $rv, $rc);

# 	1. check for C/N in database for contigs, EST orientation of singletons
########################################################################
	my $test1_time = time;
	my $test1_hits = 0;
	
	$stm = "select o.external_id, c.assembly_orientation from unigene_assembly_component as c left join unigenes as u on c.unigene_assembly_id=u.unigene_element_id left join other_identifier as o on c.seq_id=o.local_db_id where u.unigene_build_id=? and u.unigene_element_type_id='1' and o.external_id_type=1 and c.assembly_orientation='C';";
	$sth = $dbh->prepare($stm) or die "Can't prepare statement: $DBI::errstr";
	$rv = $sth->execute($id) or die "Can't execute statement: $DBI::errstr";
	$rc = $sth->bind_columns(\$est, \$assembly_orientation);
	while ($sth->fetch){
		$ests{$est}{1}++;
		$test1_hits++;
	}

# every read_direction in the database is 5... this query is a waste of time.
# 	$stm = "select o.external_id, e.read_direction from unigene_singleton as s left join unigenes as u on s.unigene_singleton_id=u.unigene_element_id left join est_info as e on e.seq_id=s.seq_id left join other_identifier as o on s.seq_id=o.local_db_id where u.unigene_build_id=? and u.unigene_element_type_id='2' and o.external_id_type=1 and e.read_direction='3';";
# 	$sth = $dbh->prepare($stm) or die "Can't prepare statement: $DBI::errstr";
# 	$rv = $sth->execute($id) or die "Can't execute statement: $DBI::errstr";
# 	$rc = $sth->bind_columns(\$est, \$read_direction);
# 	while ($sth->fetch){
# 		$ests{$est}{1}++;
#		$test1_hits++;
# 	}

	runtime::runtime_print($test1_time, "Test 1 for build $id");
	print "\t$test1_hits found\n";

# 	2. BLAST against arabidopsis
########################################################################
	my $test2_time = time;
	my $test2_hits = 0;
	
	# dump all relevant ests into a text file
	my $blast_temp = "${project}_build_${id}_est_sequences.fasta";
	open (SEQ_OUT, ">$blast_temp");
	
	# get contigs
	$stm = "select o.external_id, ts.sequence_data from unigene_assembly_component as c left join unigenes as u on c.unigene_assembly_id=u.unigene_element_id left join other_identifier as o on c.seq_id=o.local_db_id left join trimmed_sequence as ts on c.seq_id=ts.seq_id where u.unigene_build_id=? and u.unigene_element_type_id='1' and o.external_id_type=1;";
	$sth = $dbh->prepare($stm) or die "Can't prepare statement: $DBI::errstr";
	$rv = $sth->execute($id) or die "Can't execute statement: $DBI::errstr";
	$rc = $sth->bind_columns(\$est, \$trimmed_sequence);
	while ($sth->fetch){
		print SEQ_OUT ">$est\n$trimmed_sequence\n";
	}

	# get singletons
	$stm = "select o.external_id, ts.sequence_data from unigenes as u left join unigene_singleton as s on s.unigene_singleton_id=u.unigene_element_id left join trimmed_sequence as ts on s.seq_id=ts.seq_id left join other_identifier as o on s.seq_id=o.local_db_id where u.unigene_build_id=? and u.unigene_element_type_id='2' and o.external_id_type=1;";
	$sth = $dbh->prepare($stm) or die "Can't prepare statement: $DBI::errstr";
	$rv = $sth->execute($id) or die "Can't execute statement: $DBI::errstr";
	$rc = $sth->bind_columns(\$est, \$trimmed_sequence);
	while ($sth->fetch){
		print SEQ_OUT ">$est\n$trimmed_sequence\n";
	}

	close(SEQ_OUT);

	# blast the text file against arabidopsis
	my $blast_command = "$blast_program -p blastx -e 1e-10 -i $blast_temp -o $blast_temp.result -d $ath1_loc -m 8";
	system($blast_command);
	
	# parse the results of the blast
	open(RESULTS, "<$blast_temp.result") or die "Couldn't open $blast_temp.result: $!\n";
	while(<RESULTS>) {
		my ($query_idnt, $sbjt_idnt, $pct_idnt, $algn_lngth, $num_mis, $num_gaps, $start_query, $end_query, $start_sbjt, $end_subjt, $evalue, $bit_score) = split(/\t/, $_);
		if (($end_query<$start_query)&&($end_subjt>$start_sbjt)) {
			$ests{$query_idnt}{2}++;
			$ests{$query_idnt}{2} == 1 and $test2_hits++;
		}
	}
	
	# clean up
	unlink($blast_temp);
	unlink("$blast_temp.result");
	
	runtime::runtime_print($test2_time, "Test 2 for build $id");
	print "\t$test2_hits found\n";

# 	3. check for Poly T 'tails' at the start of sequences
########################################################################
	my $test3_time = time;
	my $test3_hits = 0;

	# number of T's in a row to consider odd
	my $num_t = 10;
	
	# get contigs
	$stm = "select o.external_id from unigene_assembly_component as c left join unigenes as u on c.unigene_assembly_id=u.unigene_element_id left join other_identifier as o on c.seq_id=o.local_db_id left join raw_sequence as rs on c.seq_id=rs.seq_id where unigene_build_id=? and u.unigene_element_type_id='1' and o.external_id_type=1 and rs.sequence_data regexp 'T{$num_t}\$';";
	$sth = $dbh->prepare($stm) or die "Can't prepare statement: $DBI::errstr";
	$rv = $sth->execute($id) or die "Can't execute statement: $DBI::errstr";
	$rc = $sth->bind_columns(\$est);
	while ($sth->fetch){
		$ests{$est}{3}++;
		$test3_hits++;
	}
	
	# get singletons
	$stm = "select o.external_id from unigene_singleton as s left join unigenes as u on s.unigene_singleton_id=u.unigene_element_id left join raw_sequence as rs on rs.seq_id=s.seq_id left join other_identifier as o on s.seq_id=o.local_db_id where u.unigene_build_id=? and u.unigene_element_type_id='2' and o.external_id_type=1 and rs.sequence_data regexp 'T{$num_t}\$';";
	$sth = $dbh->prepare($stm) or die "Can't prepare statement: $DBI::errstr";
	$rv = $sth->execute($id) or die "Can't execute statement: $DBI::errstr";
	$rc = $sth->bind_columns(\$est);
	while ($sth->fetch){
		$ests{$est}{3}++;
		$test3_hits++;
	}
	
	runtime::runtime_print($test3_time, "Test 3 for build $id");
	print "\t$test3_hits found\n";

########################################################################
# generate output
# 	<est>	<0|1 - test 1>	<0|1 - test 2>	<0|1 - test 3>
########################################################################
	
	print "\n\n";

	open(OUT, ">$outfile") or die "Couldn't open $outfile: $!\n";
	foreach (sort { uc($b) cmp uc($a) } keys %ests) {
		print OUT $_."\t";
		$ests{$_}{1}? print OUT "1\t" : print OUT "0\t";
		$ests{$_}{2}? print OUT "1\t" : print OUT "0\t";
		$ests{$_}{3}? print OUT "1\t" : print OUT "0\t";
		print OUT "\n";

		if ($ests{$_}{1} && $ests{$_}{2} && $ests{$_}{3}) {
			$likely_flipped .= "$_ positive for all three tests.\n";
		}
		elsif ($ests{$_}{1} && $ests{$_}{2}) {
			$likely_flipped .=  "$_ positive for tests 1 and 2.\n";
		}
		elsif ($ests{$_}{1} && $ests{$_}{3}) {
			$likely_flipped .= "$_ positive for tests 1 and 3.\n";
		}
		elsif ($ests{$_}{2} && $ests{$_}{3}) {
			$likely_flipped .= "$_ positive for tests 2 and 3.\n";
		}

	}
	close OUT;
	
}

runtime::runtime_print($total_time, "Orientation screens of $project builds $set");

db_link::disconnect_db($dbh);

print "\n\n$likely_flipped\n";