#!/usr/bin/perl

$file=shift;

open in,$file;
while(<in>){
chop $_;
if($.<3){next;}
@vals=split(/\t/,$_);



$seq1=$vals[$#vals-1];
@seq1array=split(//,$seq1);
$seq2=$vals[$#vals];
@seq2array=split(//,$seq2);

$alnsubststring="";

for($i=0;$i<$#seq1array;$i++){
	if($seq1array[$i] ne $seq2array[$i] && $seq1array[$i] ne "-" && $seq2array[$i] ne "-"){
		$alnsubststring.="1";
	}else{
		$alnsubststring.="0";
	}
}


@triplets = ($alnsubststring =~ m/.../g); 


# reset Contatori mismatch:
$C1=0;
$C2=0;
$C3=0;


# Conteggio sostituzioni:
foreach $t (@triplets){
	@tripletarray=split(//,$t);
	$C1+=$tripletarray[0];
	$C2+=$tripletarray[1];
	$C3+=$tripletarray[2];
}

# Calcolo approssimazione potenza segnale sostituzioni a periodo 1/3

$SCP = (1/2) * ( ($C1-$C2)**2 + ($C2-$C3)**2 + ($C3-$C1)**2 );
$SCPnorm = $SCP/$vals[1];

# Stampa risultato:

print "$vals[0]\t$vals[1]\t$SCP\t$SCPnorm\t( C1: $C1 C2: $C2 C3: $C3 )\n";

}
close in;
