#!/usr/bin/perl -w
##############################################
### Huang He                               ###
### take_taxo_from_bwa.sam                 ###
##############################################

($bwa,$silva)=@ARGV;
open(File,"$bwa")||die "can't open the $bwa\n";
while(<File>){
chomp;
s/\c@*//g;
s/\cM*$//;
$mle=0;
	if(/^\@SQ/){
		next;
	}
	if(/^\@PG/){
		next;
	}
	@temp=split;
	if($temp[2] eq "*"){
		next;
	}
	push(@{$haxi{$temp[2]}},$temp[0]);
	@temp1=split(/M/,$temp[5]);
	$tlen=$temp1[0];
	if($temp[5] ne "101M"){
		$tlen=101;
	}
#		print"$tlen\t";
	for($a=0;$a<@temp;$a++){
		if($temp[$a]=~/MD\:Z\:/){
			@temp2=split(/Z\:/,$temp[$a]);
			$mlen=$temp2[1];
			@temp3=split(/[A|T|G|C|\^]/,$mlen);
			for($b=0;$b<@temp3;$b++){
				if($temp3[$b]=~/[\d]/){
					$mle+=$temp3[$b];
#					print"$mlen\t$temp3[$b]\t";
				}	
			}
		}
	}
	$per=($mle/$tlen)*100;
	push(@{$perc{$temp[2]}},$per);
#	print"$per\n";
#print"\n";
}
close File;

#foreach(keys %haxi){
#	$reads=$_;
#	$num=@{$haxi{$reads}};
#	print"$reads\t$num\n";
#}
open(File,"$silva")||die "can't open the $silva\n";
while(<File>){
chomp;
s/\c@*//g;
s/\cM*$//;
	if(/^\>/){
		@temp=();
		$line=$_;
		@temp=split(/\s/,$line);
		$id=$temp[0];
		$id=~s/\>//g;
##		print"$id\n";
		if(! exists $haxi{$id}){
			next;
		}
##		print"$id\n";
		$idl{$id}=$line;
	}
}
close File;

print"ID\tCount\tAvg.% Identity\tGenome\tTaxonomy\n";
$t_cou=0;
foreach(keys %idl){
	$kkk=0;
	@temp=();
	$refn=$_;
	$cou=@{$haxi{$refn}};
	$t_cou+=$cou;
	@num=@{$perc{$refn}};
	$nu=@num;
	for($c=0;$c<@num;$c++){
		$kkk+=$num[$c];
	}
	$per=$kkk/$nu;
	@temp=split(/\;/,$idl{$refn});
	@temp1=split(/\s/,$temp[0]);
	$fir=$temp1[1];
	$gen=$temp[-1];
	print"$refn\t$cou\t$per\t$gen\t$fir\t";
	for($a=1;$a<@temp;$a++){
		print"$temp[$a]\t";
	}
	print"\n";
}
#print"\t$t_cou\n";