use strict;
use Bio::Seq;
use Bio::SeqIO;
use Bio::AlignIO;

open(OUT,">99.03hotmut.txt");

open(NAME168,"Bacillus_subtilis_subsp_subtilis_str_168.GCA_000009045.1.17.gtf");
my %gene1;
while(my $line=<NAME168>){
	chomp $line;
	my @data=split("\t",$line);
	my @data2=split(";",$data[8]);
	$data2[1]=~/transcript_id \"(.+)\"/;
	my $geneid1=$1;#CAB11777
	$data2[3]=~/gene_name \"(.+)\"/;
	my $geneid2=$1;	#dnaA
	$gene1{$geneid1}=$geneid2;
	$gene1{$geneid2}=$geneid1;
}

open(GENE,"99.00gene_nature.txt");
my %gene2;
while(my $line=<GENE>){
	chomp $line;
	my @data=split("\t",$line);
	if(exists $gene1{$data[0]}){
		$gene2{$gene1{$data[0]}}=$data[1];
	}
}

my %seq;
open(SPE,"99.00species.txt");
while(my $line=<SPE>){
	chomp $line;
	my $seqhu=new Bio::SeqIO(-file=>$line,-format=>"fasta");
	while(my $this_seq=$seqhu->next_seq){
		$seq{$this_seq->id}=$this_seq->seq;
	}
}

my %orth;
open(ORTH,"99.02orth.txt");
while(my $line=<ORTH>){
	chomp $line;
	my @data=split("\t",$line);
	if(exists $gene2{$data[0]}){
		my $i=0;
		my $temfile=$data[0].".txt";
		my $alnfile=$data[0].".aln";
		open(TEM,">".$temfile);
		foreach my $j (@data){
			print TEM ">".$j."\n".$seq{$j}."\n";
			if($seq{$j} eq ''){
				$i=1;
			}
		}
		close TEM;
		if($i==0){
			system<clustalw $temfile>;
			my $alnfilein=Bio::AlignIO->new('-file'=>$alnfile);
			my $aln=$alnfilein->next_aln();
			$aln=$aln->remove_gaps;
			my @match=split("",$aln->match_line);
			my $num_unmat=0;
			my $i=0;#frist and last aa cound
			while($i<scalar(@match)){
				if ($match[$i] ne "*"){
					$num_unmat+=1;
				}
				$i+=1;
			}
			print OUT $data[0]."\t".$gene1{$data[0]}."\t".$gene2{$data[0]}."\t".scalar(@match)."\t".$num_unmat."\n";
		}
	}
}

