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

sub preprocess_alignments($$){
    my (
	$ref_align1,
	$ref_align2
	) = @_;

    my $align1 = uc($$ref_align1);
    my $align2 = uc($$ref_align2);

    my @align1_array = split(//,$align1);
    my @align2_array = split(//,$align2);

    my $align1_res = "";
    my $align2_res = "";

    for (my $i = 0;
	    $i < length($align1) && 
	    $i < length($align2) ;
	    $i++) {
	my $v1 = index("CFYWMLIVGPATSNHQEDRKBZX",$align1_array[$i]);
	my $v2 = index("CFYWMLIVGPATSNHQEDRKBZX",$align2_array[$i]);
	if ($v1 > 0 || $v2 > 0) {
	    $align1_res  .= $align1_array[$i];
	    $align2_res  .= $align2_array[$i];
	}
    }
    $$ref_align1 = $align1_res;
    $$ref_align2 = $align2_res; 
}


# Chech if the blast seed matches an alignment
sub blastp_seed_matches_alignment($$$$$){
    my (
	$align1,
	$align2,
	$refSM, #hash table ref of "x,x" where x is 1 letter aminoAcid symbol (-> gives score)
	$scorethreehold,
	$span,
	) = @_;

    my $j     = 0;
    my $score = 0;
    for (my $i = 0;
	    $i < (length $align1) - $span && 
	    $i < (length $align2) - $span ;
            $i++) {
      nextj:
	if ($j < $span) {
	    my $v1 = index("CFYWMLIVGPATSNHQEDRKBZX\*",substr($align1,$i+$j,1));
	    my $v2 = index("CFYWMLIVGPATSNHQEDRKBZX\*",substr($align2,$i+$j,1));
	    if (($v1 >= 0) && ($v2 >= 0)) {
		my $dscore = $$refSM{substr($align1,$i+$j,1).",".substr($align2,$i+$j,1)};
		$score += $dscore;
		$j++;
		goto nextj;
	    } else {
		$j     = 0;
		$score = 0;
	    }
	} else {
	    if ($score >= $scorethreehold){
		return 1;
	    }
	    $j     = 0;
	    $score = 0; 
	}
    } # ($i++)
    return 0;
}




# Chech if several seed tables match an alignment
sub subseed_matches_alignment($$$){
    my (
	$align1,
	$align2,
	$r_seedtable,
	) = @_;

    my $j = 0;
    for (
	my $i = 0 ; 
	   $i < (length $align1) - $#$r_seedtable && 
	   $i < (length $align2) - $#$r_seedtable ; 
	   $i++) {
      nextj:
	if ($j <= $#$r_seedtable) {
	    my $r_seed = $$r_seedtable[$j];	    
	    my $v1 = index("CFYWMLIVGPATSNHQEDRKBZX\*",substr($align1,$i+$j,1));
	    my $v2 = index("CFYWMLIVGPATSNHQEDRKBZX\*",substr($align2,$i+$j,1));
	    if (($v1 >= 0) && ($v2 >= 0)){
		for (my $k = 0; $k <= $#$r_seed ; $k++){
		    if ((index($$r_seed[$k],substr($align1,$i+$j,1)) >= 0) && (index($$r_seed[$k],substr($align2,$i+$j,1)) >= 0)) {
			$j++;
			goto nextj;
		    }
		}
		$j = 0;
	    } else {
		$j = 0;
	    }
	} else {
	    return 1;
	}
    } # ($i++)
    return 0;
}



# Chech if the seed table matches an alignment
sub subseeds_match_alignment($$$){
    my (
	$align1,
	$align2,
	$r_seeds,
	) = @_;

    for (my $i = 0 ; $i <= $#$r_seeds ; $i++) {
	#&print_table($$r_seeds[$i]);
	if (&subseed_matches_alignment($align1,$align2,$$r_seeds[$i])){
	    return 1;
	}
    }
    return 0;
}




# Read seed definition from a formatted file
sub read_subseed_alphabet($){
    my (
	$filename
	) = @_;
 
    my @seedalphabettable   = ();
    my $i = 0;
    open(IN, $filename) or die("Cant open \"$filename\"");
    while (<IN>){
	if ($_ =~ /^[CFYWMLIVGPATSNHQEDRK \t]+$/ ){ 
 	    chomp($_);
	    my @seedalphabetline = split(/[ \t]+/,$_);
	    for (my $j = 0 ; $j <= $#seedalphabetline ; $j++){
		$seedalphabettable[$i][$j] = $seedalphabetline[$j];
	    }
	    $i++;
	}
    }
    return \@seedalphabettable;
}

# read the nth Subset seeds of sensitivity at least ...
sub read_subseeds($$$$){
    my ($filename,$r_alphabet,$minselectivity,$n) = @_;
    my $c = 0;
    my @seedtable  = ();
    open(IN, $filename) or die("Cant open \"$filename\"");
    # search for a seed of selectivity at least ...
    while (<IN>){
	chomp($_);
	my @splitline = split (/[ \t]+/,$_); 
	if ($splitline[1] > $minselectivity && $c >= $n){
	    #read seeds (first column)
	    my @seeds = split (/,/,$splitline[0]);
	    for (my $i=0 ; $i <= $#seeds ; $i++) {
		my $seed = $seeds[$i];
		for (my $j=0 ; $j<length($seed); $j++) {
		    my $numberletter = substr($seed,$j,1);
		    #FIXME conversion des lettres 
		    if ( $numberletter =~ /^[A-Z]/) {
			$numberletter = ord($numberletter) - ord("A")  + 10;
		    }
		    #print($numberletter."<");
		    my $r_col = $$r_alphabet[$numberletter];
		    # for (my $k = 0; $k <= $#$r_col; $k++){
                    #   print $$r_col[$k].",";
                    # }
		    # print "\n";		   
		    $seedtable[$i][$j] = $r_col;
		}
	    }
	    return \@seedtable;
	} else {
	    $c++;
	}
    }
    return \@seedtable;
}



# Read blast formatted scoring matrix
sub read_blast_matrix($){
    my ($filename) = @_ ;
    
    my %dhashtable = ();
    my @letters = ();
    my $count = 0;
    open(IN, $filename) or die("Cant open \"$filename\"");
    while (<IN>){
	if ($_ =~ /^[ \t]*([CFYWMLIVGPATSNHQEDRKBZX\*][ \t]*)+$/ ){ 
	    chomp($_);
	    @letters = split(/[ \t]+/," ".$_);	    
	} 
	if ($_ =~ /^([ \t]*[+-]*[0-9]+[ \t]*)+$/ ){
	    chomp($_);
	    my @scores = split(/[ \t]+/," ".$_);
	    $count++;
	    for(my $i = 0 ;$i <= $#scores ; $i++){
		$dhashtable{$letters[$count].",".$letters[$i]} = $scores[$i];
		$dhashtable{$letters[$i].",".$letters[$count]} = $scores[$i];
	    }
	} 
    }
    return \%dhashtable;
}





# Print 2d table
sub print_table($){
    my ($r_T) = @_;
    for (my $i = 0 ; $i <= $#$r_T; $i++){
	my $r_L  = $$r_T[$i];
	for (my $j = 0 ; $j <= $#$r_L; $j++){
	    print $$r_L[$j].",";
	}
	print "\n";
    }
}

# Print 3d table
sub print_3table($){
    my ($r_T) = @_;
    for (my $i = 0 ; $i <= $#$r_T; $i++){
	my $r_X  = $$r_T[$i];
	for (my $j = 0 ; $j <= $#$r_X; $j++){
	    my $r_Y = $$r_X[$j];
	    for (my $k = 0 ; $k <= $#$r_Y; $k++){
		print $$r_Y[$k].",";	
	    }
	    print "\n";
	}
	print "--\n";
    }
}







#------#
# MAIN #
#------#

my @SELECTIVITIES = (0,0,0,0,0,
		     0.96234,0.975559,0.984249,0.99002,0.993846,
		     0.996342,0.997909,0.998827,0.999355,0.999644,
		     0.99981);

my $SCORINGMATRIX = &read_blast_matrix("blosum62.bla");

# 1) read alignment files names
foreach my $dir (
		 "BALIBASE3","HOMSTRAD2008","IRMBASE1","OXBENCH1.3","PFAM22","PREFAB4","SMART4"
		 ) {

    my $dirmfas = $dir."/mfas";

    print "% Base ".$dir."\n";

    opendir(DIR,$dirmfas);
    my @files = readdir(DIR);
    closedir(DIR);
    
    # 2) select a BLAST scoring threshold
    for my $SCORE_THRESHOLD (11,12,13) {

	my $SELECTIVITY = $SELECTIVITIES[$SCORE_THRESHOLD];

	# 3) select alphabet and scoring matrix  
	foreach my $ALPHABET_NAME  ("transitive-predefined","transitive-ab-initio","non-tree-transitive") {
	    # read alphabet
	    my $ALPHABET = &read_subseed_alphabet("Results-alphabets/Alphabet.".$ALPHABET_NAME.".txt");
	    # print alphabet
	    print "% Alphabet ".$ALPHABET_NAME."\n";
	    ##&print_table($ALPHABET);
	    
	    
	    # 4) select the n-th seed of selectivity >= BLASTP 
	    foreach my $N  (0) {
		
		my $Bonly  = 0;
		my $Sonly  = 0;
		my $togh   = 0;
		my $none   = 0;
		
		
		# 5) read subset seed file
		my $SUBSEEDS = &read_subseeds("Results-models/results.".$ALPHABET_NAME.".B62.L16.dat",$ALPHABET,$SELECTIVITY,$N);
		##&print_3table($SUBSEEDS);
		
		# 6) read and check alignment files
		foreach my $file (@files) {	
		    
		    # 6.1) memorize full alignment
		    my @align = ();
		    open(IN,"${dirmfas}/${file}") or die("Cant open ${dirmfas}/${file}");
		    my $alignment = "";
		    while(<IN>) {
			chomp($_);   
			if ($_ =~ /^>.*/ ) {			
			    if (!($alignment eq "")) {
				push(@align,$alignment);
				$alignment = "";
				
			    }
			} else {
			    $alignment .= $_;
			}
		    }
		    if (!($alignment eq "")) {
			push(@align,$alignment);
			$alignment = "";
		    }
		    close(IN);
		    
		    srand(0);
		    
		    # 6.2) test the alignment at most x-times
		    my $nb = $#align+1;
		    if ($nb >= 2) {
			
			# nb runs
			my $runs = $nb;
			if ($nb == 2) {
			    $runs = 1;			
			} else {
			    if ($nb < 50) {
				$runs = int($nb * $nb / 2);
			    }
			    if ($nb > 5000) {
				$runs = int(sqrt($nb));
			    }
			    
			}
			
			# tests
			for (my $r = 0 ; $r < $runs ; $r++) {			
			    my $i = int(rand($nb));
			    my $j = int(rand($nb));			
			    while ($i == $j) {
				$j = int(rand($nb));
			    }
			    
			    my $align1 = $align[$i];
			    my $align2 = $align[$j];
			    
			    &preprocess_alignments(\$align1,\$align2);
			    
			    # test
			    if (!($align1 eq "" || $align2 eq "")) {				
				
				my $S = &subseeds_match_alignment($align1,$align2,$SUBSEEDS);
				my $B = &blastp_seed_matches_alignment($align1,$align2,$SCORINGMATRIX,$SCORE_THRESHOLD,3);
				if ($S) {
				    if ($B) {
					$togh++;
				    } else {
					$Sonly++;
				    }
				} else {
				    if ($B) {
					$Bonly++;
				    } else {
					$none++;
				    }
				}
				$align1=""; 
				$align2="";
			    }
			} #$r
		    } #$nb
		} # files
		open OUT,">> results.${ALPHABET_NAME}.SC${SCORE_THRESHOLD}.dat" or die "Cant open results.${ALPHABET_NAME}.SC${SCORE_THRESHOLD}.dat";	
		print OUT "# ".$dir." : ".$Sonly." found only by subset seed (total=".(100*($Sonly+$togh)/($Sonly+$Bonly+$togh+$none))."%),".
		    " vs ".$Bonly." found only by blastp (total=".(100*($Bonly+$togh)/($Sonly+$Bonly+$togh+$none))."%)".
		    "(both ".$togh.", missed ".$none.")\n";
		print OUT "".$ALPHABET_NAME."\t"."x"."\t".$Sonly."\t".$Bonly."\t".$none."\t".$togh."\n";
		close OUT;
	    } # N
	} # Score-threshold
    } # Seed-Alphabet
} # Database
    





