use strict;
use warnings;

#Bioperl moduler
use Bio::Tools::Run::StandAloneBlast;
use Bio::SeqIO;
use Bio::SearchIO;
use Bio::SearchIO::Writer::HTMLResultWriter;


#Hovedprogram

my @params = (program => 'blastp', database => 'YAL001CDB', Filterquerysequence => 'F');	#parametre til StandAloneBlast
my $blast_obj =  Bio::Tools::Run::StandAloneBlast->new(@params);	#danner blast objekt

my $seqio_obj_input = Bio::SeqIO->new(-file => "YAL001C", -format => "fasta");	#input objekt (query) laves og parametre til den gives

my $seqio_obj_output = Bio::SeqIO->new(-file => '', -format => "fasta");	#output objekt laves og parametre til output filen gives (STDIN)

while (my $seq_obj = $seqio_obj_input -> next_seq) {	#én sekvens (query) trækkes ud af input filen af gangen vha. next_seq metoden

	my $report_obj = $blast_obj->blastall($seq_obj);	#denne sekvens (query) blastes med præferencerne sat ovenover, og et rapport objekt laves

	my $result_obj = $report_obj->next_result;	#ét resultat trækkes ud af rapporten af gangen vha. next_result metoden, og danner 									#resultat objektet. Objektet indeholder alle hits som igen indeholder alle hsp.

	my $parse_blast_results = parse_blast_result(\$result_obj);	#resultatet sendes til parse_blast_result subrutinen. Dvs. ét resultat 											#for hvert blastet input (query) sekvens parses af gangen.
	my ($unique_hash_key) = keys %$parse_blast_results;
	
	write_blast_report_html(\$result_obj, \$unique_hash_key);	#blast resultatet for hvert query protein visualiseres vha. denne 										#subrutine
}

exit;



################################################################

#Subrutiner

#################################################################

sub parse_blast_result {	#hovedsubrutine der parser blast resultatet primær ved at kalde andre subrutiner

	my ($result_obj) = @_;	#resultat input objekt fra hovedprogrammet. 

	#initialisering
	my $percent_identity = 100;
	my $count_hit = 0;
	my %result_hash = ();
	my $number_of_hits = $$result_obj->num_hits;
		
	my $total_query_protein_length = $$result_obj->query_length;	#længden af hele query protein input sekvensen

	#headeren for query sekvensen parses
	my ($query_protein_acc_ID, $query_chromosome_nr, $query_chromosome_position, $query_gene_orientation, 
	$number_of_query_introns) = parse_query_protein_header(\$$result_obj);	

	while (my $hit = $$result_obj->next_hit) {	#ét hit resultat trækkes ud af resultat objektet af gangen vha. next_hit metoden.
							#Hit resultatet indeholder den eller de alignet sekvens(er) der er mellem query sekvensen og 
							#denne ene hit (subject) sekvens.
		#initialisering
		my $count_hsp = 0;
		my $query_alignment_length_result = '';
		my $subject_alignment_length_result = '';
		my @aligned_query_region = ();
		my @aligned_subject_region = ();
		my %query_nomatch_hash = ();
		my %subject_nomatch_hash = ();
		my %query_conserved_hash = ();
		my %subject_conserved_hash = ();
		my $unique_hash_key = '';
		
		#angiver nummeret af hittet for query sekvensen
		$count_hit++;	
	
		#bestemmer antal HSP'ere for hittet
		my $number_of_hsps_in_hit = $hit->num_hsps;
			
		#headeren for hit (subjekt) sekvensen parses
		my ($subject_protein_acc_ID, $subject_chromosome_nr, $subject_chromosome_position, $subject_gene_orientation, 
		$number_of_subject_introns) = parse_subject_protein_header(\$hit);

		#Tjekker om query og hit sekvensen findes på samme kromosom
		my $check_on_same_cromosome = check_on_same_cromosome(\$query_chromosome_nr, \$subject_chromosome_nr);

		while (my $hsp = $hit->next_hsp) {	#ét hsp resultat trækkes ud af hit objektet af gangen vha. next_hsp metoden.

			#initialisering
			my $unique_hash_key = '';

			#angiver nummeret af hsp'en for query sekvensens givne hit
			$count_hsp++;

			#percent_identity metoden bestemmer den procentvise identitet mellem de alignet sekvenser 
			$percent_identity = $hsp->percent_identity;

			#her kan der sættes et lavere mål for identitetsprocenten af query-hsp. Dvs. man kan sortere dårlige alignments fra 
			if ($hsp->percent_identity < 90) {next};			
		
			#tjekker om hele query input protein sekvensen er med i alignmenten					
			$query_alignment_length_result = check_if_full_length_query_aligned(\$hsp, \$$result_obj);

			#tjekker om hele subjekt input protein sekvensen er med i alignmenten	
			$subject_alignment_length_result = check_if_full_length_subject_aligned(\$hsp, \$hit);

			#angiver intervallet af query regionen der er alignet
			@aligned_query_region = retrieve_aligned_query_region(\$hsp);

			#angiver intervallet af subjekt regionen der er alignet
			@aligned_subject_region = retrieve_aligned_subject_region(\$hsp);

			#linker query positionen til Aminosyren for de anminosyrer der ikke matcher en Aa i subjekt sekvensen. Dvs. en key i hashen er  			#en position og dens value er den Aa i query sekvensen der ikke er et match til i subjekt sekvensen.
			%query_nomatch_hash = query_seq_alignment_nomatch_residues(\$hsp);

			#linker subjekt positionen til Aminosyren for de anminosyrer der ikke matcher en Aa i query sekvensen. Dvs. en key i hashen er 				#en position og dens value er den Aa i subject sekvensen der ikke er et match til i query sekvensen.
			%subject_nomatch_hash = subject_seq_alignment_nomatch_residues(\$hsp);
	
			#linker query positionen til Aminosyren for de anminosyrer der er konserverede men ikke identiske til en Aa i subjekt sekvensen. 				#Dvs. en key i hashen er en position og dens value er den Aa i query sekvensen der er konserveret men ikke identisk til en Aa i 			#subjekt sekvensen.
			%query_conserved_hash = query_seq_alignment_conserved_residues(\$hsp);

			#linker subjekt positionen til Aminosyren for de anminosyrer der er konserverede men ikke identiske til en Aa i query sekvensen. 				#Dvs. en key i hashen er en position og dens value er den Aa i subjekt sekvensen der er konserveret men ikke identisk til en Aa i 				#query sekvensen.
			%subject_conserved_hash = subject_seq_alignment_conserved_residues(\$hsp);
	
			#laver et unikt ID for hver hsp. Medtager query ID, hit nummer og hsp nummer 
			$unique_hash_key = $query_protein_acc_ID."hit".$count_hit."hsp".$count_hsp;

			#hasher alle de ovenover fundne værdier til den unikke hash key i result hashen
			$result_hash{$unique_hash_key} = [$query_protein_acc_ID,  			
				       $query_chromosome_nr, $query_chromosome_position, $query_gene_orientation,
				       $number_of_query_introns, $subject_protein_acc_ID, $subject_chromosome_nr, 
	 			       $subject_chromosome_position, $subject_gene_orientation, $number_of_subject_introns, 				               $query_alignment_length_result, $subject_alignment_length_result, $percent_identity, 					       \%query_nomatch_hash, \%subject_nomatch_hash, \%subject_conserved_hash,      
				       \%query_conserved_hash, $number_of_hsps_in_hit, $count_hsp, \@aligned_query_region, 
				       \@aligned_subject_region, $check_on_same_cromosome,$total_query_protein_length];

			#printer de parsede resultater ud for hver query blastet mod databasen
			open(FH, ">$unique_hash_key") or die "cannot save $unique_hash_key to textfile\n";

			print FH ($unique_hash_key) = keys %result_hash, "\n";
			print FH $result_hash{$unique_hash_key}[0], "\n";	#query_protein_acc_ID,  			
			print FH $result_hash{$unique_hash_key}[1], "\n"; #query_chromosome_nr, 
			print FH $result_hash{$unique_hash_key}[2], "\n"; #query_chromosome_position, 
			print FH $result_hash{$unique_hash_key}[3], "\n"; #query_gene_orientation,
			print FH $result_hash{$unique_hash_key}[4], "\n"; #number_of_query_introns, 
			print FH $result_hash{$unique_hash_key}[5], "\n"; #subject_protein_acc_ID, 
			print FH $result_hash{$unique_hash_key}[6], "\n"; #subject_chromosome_nr,
			print FH $result_hash{$unique_hash_key}[7], "\n"; #subject_chromosome_position, 
			print FH $result_hash{$unique_hash_key}[8], "\n"; #subject_gene_orientation, 
			print FH $result_hash{$unique_hash_key}[9], "\n"; #number_of_subject_introns, 	               
			print FH $result_hash{$unique_hash_key}[10], "\n"; #query_alignment_length_result, 
			print FH $result_hash{$unique_hash_key}[11], "\n"; #subject_alignment_length_result, 
			print FH $result_hash{$unique_hash_key}[12], "%", "\n"; #percent_identity, 
			print FH "Query nomatches:\n";
			#query_nomatches 
			while (my ($key, $value) = each %{$result_hash{$unique_hash_key}[13]}) {print FH "$key => $value\n";} 	 
			print FH "Subject nomatches:\n";
			#subject_nomatches
			while (my ($key, $value) = each %{$result_hash{$unique_hash_key}[14]}) {print FH "$key => $value\n";} 
			print FH "Subject conserved:\n";
			#subject_conserved
			while (my ($key, $value) = each %{$result_hash{$unique_hash_key}[15]}) {print FH "$key => $value\n";} 
			print FH "Query conserved:\n";
			#query_conserved 
			while (my ($key, $value) = each %{$result_hash{$unique_hash_key}[16]}) {print FH "$key => $value\n";}  
			print FH $result_hash{$unique_hash_key}[17], "\n"; #number_of_hsps_in_hit, 
			print FH $result_hash{$unique_hash_key}[18], "\n"; #count_hsp, 
			print FH my $query_alignment_start_nr = ${$result_hash{$unique_hash_key}[19]}[0];	#start nr for query alignment 
			print FH "-";
			print FH my $query_alignment_end_nr = ${$result_hash{$unique_hash_key}[19]}[1], "\n"; #slut nr for query 																   #alignment 
			print FH my $subject_alignment_start_nr = ${$result_hash{$unique_hash_key}[20]}[0]; #start nr for subject 															          #alignment 
			print FH "-";
			print FH my $subject_alignment_end_nr = ${$result_hash{$unique_hash_key}[20]}[1], "\n"; #slut nr for subject 																      #alignment   
			print FH $result_hash{$unique_hash_key}[21], "\n"; #check_on_same_cromosom

			print FH $total_query_protein_length, "\n"; #længden af hele query protein input sekvensen		

			close FH;	

			my $database_name = $unique_hash_key.'_result_hash';	#generere navnet for HSP'ens dbm fil

			#gemmer hver resultat-hash som er genereret i parse_blast_result subrutinen for hvert $result_obj til biolinux_result_hash dbm'en
			dbmopen(%result_hash, $database_name, 0666) or die "cannot save $database_name to dbm\n";
			dbmclose(%result_hash); 

		}	

	}

	return \%result_hash;

}
				

#################################################################

sub check_if_full_length_query_aligned {	#tjekker om hele query input protein sekvensen er med i alignmenten

	my ($hsp, $result_obj) = @_;

	my $total_query_protein_length = $$result_obj->query_length;		#længden af hele query protein input sekvensen
	
	my $alignment_query_protein_length = $$hsp->length('query');	#længden af hele query protein sekvensen alignet, 											#eksklusiv gaps

	if ($total_query_protein_length == $alignment_query_protein_length) {

		return 'QUERY_FULL_LENGTH';

	} else {

	return 'QUERY_NOT_FULL_LENGTH';

	}

}

########################################################################

sub check_on_same_cromosome {	#tjekker om query og subject matchet blev fundet på det samme kromosom

	my ($query_chromosome_nr, $subject_chromosome_nr) = @_;	

	if ($$query_chromosome_nr eq $$subject_chromosome_nr) {

		return 'wright_cromosome'; 
	
	} else {return 'wrong_cromosome'}
	
}


##########################################################################

sub parse_query_protein_header {	#parser query proteinets accession ID, kromosom lokalisation, position på kromosom,  
					#gen's orientation på kromosom og antallet af gen introns, og returnere dem
				
	my ($result_obj) = @_;

	my $query_protein_header = $$result_obj->query_description;		#query_description metoden returnere query headeren

	my @query_protein_header_elements = split / /, $query_protein_header;	#splitter headeren

	my $query_protein_acc_ID = $$result_obj->query_accession;		#query_accession metoden returnere query accession ID'et
	my $query_chromosome_nr = $query_protein_header_elements[3];		#kromosome nr trækkes ud af query header
	my $query_chromosome_position = $query_protein_header_elements[5];	#kromosome position trækkes ud af query header	
	my $query_gene_orientation = $query_protein_header_elements[6];	#query gen orientation trækkes ud af query header
	
	$query_protein_acc_ID =~ s/\.//g;						#ryder op i query accession ID'et

	if ($query_gene_orientation =~ /reverse/ig) { $query_gene_orientation = 'revcom'};	#reverse oversættes til revcom

	my $number_of_introns = ($query_chromosome_position =~ tr/,//) - 1;			#antal introns bestemmes ud fra headerens 														#udtrukne kromosom position

	return ($query_protein_acc_ID, $query_chromosome_nr, $query_chromosome_position, $query_gene_orientation, $number_of_introns);

}


##########################################################################

sub parse_subject_protein_header {	#parser subject proteinets accession ID, kromosom lokalisation, position på kromosom,  
					#gen's orientation på kromosom og antallet af gen introns, og returnere dem til sub parse_blast_result 
				
	my ($hit) = @_;

	my $subject_protein_header = $$hit->description;					#description metoden returnere subjekt headeren

	my @subject_protein_header_elements = split / /, $subject_protein_header; 	#splitter headeren

	my $subject_protein_acc_ID = $$hit->accession;					#accession metoden returnere subjekt accession ID'et
	my $subject_chromosome_nr = $subject_protein_header_elements[3];		#kromosome nr trækkes ud af subjekt header
	my $subject_chromosome_position = $subject_protein_header_elements[5];	#kromosome position trækkes ud af subjekt header
	my $subject_gene_orientation = $subject_protein_header_elements[6];		#query gen orientation trækkes ud af subjekt header
	
	$subject_protein_acc_ID =~ s/\.//g;							#ryder op i subjekt accession ID'et

	if ($subject_gene_orientation =~ /reverse/ig) { $subject_gene_orientation = 'revcom'};	#reverse oversættes til revcom

	my $number_of_subject_introns = ($subject_chromosome_position =~ tr/,//) - 1;		#antal introns bestemmes ud fra headerens 														#udtrukne kromosom position

	return ($subject_protein_acc_ID, $subject_chromosome_nr, $subject_chromosome_position, $subject_gene_orientation, 
		$number_of_subject_introns);
}


#####################################################################################

sub write_blast_report_html {	#subrutine der udprinter blast rapporten i html format til filen "searchio2.html"

	my ($result_obj, $unique_hash_key) = @_;

	my $writerhtml = new Bio::SearchIO::Writer::HTMLResultWriter();

	my $outhtml = new Bio::SearchIO(-writer => $writerhtml, -file => ">$$unique_hash_key.html");	#ouput er kaldt den unikke hash key.html

	$outhtml->write_result($$result_obj);

}


#############################################################################

sub query_seq_alignment_nomatch_residues {		#linker query positionen til Aminosyren i en hash for de anminosyrer der ikke matcher en Aa i 								#subjekt sekvensen.  	
							
	my ($hsp_obj) = @_;
	my %query_nomatch_hash = ();
	my @new_query_string = ();

	my @query_string = split "", $$hsp_obj->query_string;	#query_string metoden trækker den alignet del af query sekvensen ud
 
	foreach (@query_string) {					#klipper gaps ud af query sekvensen. Ellers kan man få forkerte nomatch positioner
									#i visse tilfælde
		if ($_ ne '-') {push @new_query_string, $_};

	}

	my $start_query_number = $$hsp_obj->start('query');	#start('query') metoden trækker query alignment start positionen ud 
									#Dvs. den position hvor query'ets alignment med subjektet starter

	$start_query_number = $start_query_number - 1;		#nødvendig modifikation af query start positionen

	foreach ($$hsp_obj->seq_inds('query', 'nomatch')) {	#seq_inds('query', 'nomatch') metoden trækker positionerne ud på alle de Aa i query 									#sekvensen der ikke har et match på subjekt sekvensen 

		$query_nomatch_hash{$_} = $new_query_string[$_ -1 -$start_query_number];	#nødvendig modifikation af query nomatch 														#positionen, trækker derefter den tilsvarende 														#aminosyre ud af query sekvensen og hasher den 														#til positionen.
	}

	return %query_nomatch_hash;

}	


################################################################


sub subject_seq_alignment_nomatch_residues {	#linker subject positionen til Aminosyren i en hash for de anminosyrer der ikke matcher en Aa i 						#query sekvensen. 

	my ($hsp_obj) = @_;
	my %subject_nomatch_hash = ();
	my @new_subject_string = ();

	my @subject_string = split "", $$hsp_obj->hit_string;	#hit_string metoden trækker den alignet del af subjekt sekvensen ud

	foreach (@subject_string) {					#klipper gaps ud af subject sekvensen. Ellers kan man få forkerte nomatch 										#positioner i visse tilfælde
		if ($_ ne '-') {push @new_subject_string, $_};

	}

	my $start_subject_number = $$hsp_obj->start('hit');	#start('hit') metoden trækker subjekt alignment start positionen ud 
									#Dvs. den position hvor subjektets alignment med query'en starter
	$start_subject_number = $start_subject_number - 1;	#nødvendig modifikation af subject start positionen

	foreach ($$hsp_obj->seq_inds('hit', 'nomatch')) {	#seq_inds('hit', 'nomatch') metoden trækker positionerne ud på alle de Aa i subjekt 									#sekvensen der ikke har et match på query sekvensen 

		$subject_nomatch_hash{$_} = $new_subject_string[$_ -1 -$start_subject_number];#nødvendig modifikation af subjekt nomatch 														#positionen, trækker derefter den tilsvarende 														#aminosyre ud af subjekt sekvensen og hasher den 														#til positionen.
	}

	return %subject_nomatch_hash;

}	



###################################################################


sub check_if_full_length_subject_aligned {	#tjekker om hele subject input protein sekvensen er med i alignmenten

	my ($hsp, $hit) = @_;

	my $total_subject_protein_length = $$hit->length;		#længden af hele subject protein input sekvensen

	my $alignment_subject_protein_length = $$hsp->length('hit');	#længden af hele subject protein sekvensen 											#alignet, eksklusiv gaps

	if ($total_subject_protein_length == $alignment_subject_protein_length) {

		return 'SUBJECT_FULL_LENGTH';

	} else {

	return 'SUBJECT_NOT_FULL_LENGTH';

	}

}


#############################################################################

sub retrieve_aligned_query_region {	#udtrækker regionen af query input protein sekvensen der er med i alignmenten

	my ($hsp) = @_;
	my @query_aligned_region = $$hsp->range('query');		#range('query') metoden giver regionen af hele query protein sekvensen alignet
									#F.eks. (1, 1160)
	return @query_aligned_region;

}


#############################################################################

sub retrieve_aligned_subject_region {	#udtrækker regionen af subject input protein sekvensen der er med i alignmenten

	my ($hsp) = @_;

	my @subject_aligned_region = $$hsp->range('hit');		#range('query') metoden giver regionen af hele subject protein sekvensen alignet
									#F.eks. (1, 1160)
	return @subject_aligned_region;

}


#############################################################################


sub query_seq_alignment_conserved_residues {		#linker query positionen til Aminosyren i en hash for de anminosyrer der er konserverede men ikke 								#identiske til en Aa i subjekt sekvensen.

	my ($hsp_obj) = @_;
	my %query_conserved_hash = ();
	my @new_query_string = ();

	my @query_string = split "", $$hsp_obj->query_string;	#query_string metoden trækker den alignet del af query sekvensen ud

	foreach (@query_string) {					#klipper gaps ud af query sekvensen. Ellers kan man få forkerte  										#positioner for de konservede Aa. i visse tilfælde
		if ($_ ne '-') {push @new_query_string, $_};

	}

	my $start_query_number = $$hsp_obj->start('query');	#start('query') metoden trækker query alignment start positionen ud 
									#Dvs. den position hvor query'ets alignment med subjektet starter

	foreach ($$hsp_obj->seq_inds('query', 'conserved-not-identical')) {	#seq_inds('query', 'conserved-not-identical') metoden trækker 											#positionerne ud på alle de Aa i query sekvensen der er konserverede men 											#ikke identiske med et match på subjekt sekvensen 		

		$query_conserved_hash{$_} = $new_query_string[$_ -1];	#nødvendig modifikation af query conserved-not-identical positionen, 
										#trækker derefter den tilsvarende aminosyre ud af query 										#sekvensen og hasher den til positionen.	
	}

	return %query_conserved_hash;

}	


################################################################


sub subject_seq_alignment_conserved_residues {		#linker subject positionen til Aminosyren i en hash for de anminosyrer der er konserverede men 								#ikke identiske til en Aa i query sekvensen.
	my ($hsp_obj) = @_;
	my %subject_conserved_hash = ();
	my @new_subject_string = ();

	my @subject_string = split "", $$hsp_obj->hit_string;	#hit_string metoden trækker den alignet del af subjekt sekvensen ud

	foreach (@subject_string) {					#klipper gaps ud af subject sekvensen. Ellers kan man få forkerte  										#positioner for de konservede Aa. i visse tilfælde

		if ($_ ne '-') {push @new_subject_string, $_};

	}

	my $start_subject_number = $$hsp_obj->start('hit');	#start('hit') metoden trækker subjekt alignment start positionen ud 
									#Dvs. den position hvor subjektets alignment med query'en starter

	foreach ($$hsp_obj->seq_inds('hit', 'conserved-not-identical')) {	#seq_inds('hit', 'conserved-not-identical') metoden trækker 											#positionerne ud på alle de Aa i subjekt sekvensen der er konserverede  										#men ikke identiske med et match på subjekt sekvensen 		

		$subject_conserved_hash{$_} = $new_subject_string[$_ -1];	#nødvendig modifikation af subjekt conserved-not-identical 												#positionen, trækker derefter den tilsvarende 												#aminosyre ud af subjekt sekvensen og hasher den 												#til positionen.
	}

	return %subject_conserved_hash;

}	



###################################################################