#! /usr/bin/perl

open(SNP,$ARGV[0]) || die "Unable to open file";
while(<SNP>)
 {
   $line = $_;
   chomp $line;
   $ind[$i] = $line;
   $i++;
 }

$numstrain = 22;
for($ind=1;$ind<$numstrain;$ind++)
 {
   $prea1 = 1 / $ind;
   $a1 = $a1 + $prea1;
   $prea2 = 1 / ($ind * $ind);
   $a2 = $a2 + $prea2;
   $a3 = $a3 + $ind;
 }

$b1 = ($numstrain + 1) / ( 3 * ($numstrain - 1));

$b2 = (2 * (($numstrain * $numstrain) + $numstrain + 3)) / (9 * $numstrain * ($numstrain-1));

$c1 = $b1 - (1 / $a1);

$c2 = $b2 - (($numstrain + 2) / ( $a1 * $numstrain)) + ($a2 / ($a1 * $a1));

$e1 = $c1 / $a1;

$e2 = $c2 / (($a1 * $a1) + $a2);

print "a1=$a1\ta2=$a2\tb1=$b1\tb2=$b2\tc1=$c1\tc2=$c2\te1=$e1\te2=$e2\ta=$a3\n";

open(FILE1,$ARGV[1]) || die "Unable to open file";

while(<FILE1>)
 {
   @b = split;
   $chr = $b[0];
   $start = $b[1];
   $end = $b[2];
   $base = $b[3];
   $pairwise=0;
   $polyloci=0;
   $tajimaD=0;
   $taji1=0;
   $taji2=0;
   $taji3=0;
   $pi=0;
   $p1=0;
    for($j=0;$j<$i;$j++)
     {
      @c = split(/\t/,$ind[$j]);
      $ind_chr = $c[0];
      $ind_pos = $c[1];
      next if $ind_chr ne $chr;
      if($ind_chr eq $chr && $ind_pos >= $start && $ind_pos <= $end)
        {
          for($k=5;$k<=$#c;$k++)
            {
              for($l=$k+1;$l<=$#c;$l++)
                {
                   $comb++;
                   if(($c[$k] == 0 && $c[$l] == 2)  || ($c[$k] == 0 && $c[$l] == 1) || ($c[$k] == 2 && $c[$l] == 0) || ($c[$k] == 1 && $c[$l] == 0))
                    {
                      $pairwise++;
                    }
                 }
              }
            #print "$comb\n";
            $comb=0;
            $polyloci++;
           }
       }
     $pi = $pairwise / $a3;

     if($polyloci > 0 && $pairwise > 0)
      {
        $p1= $polyloci / $a1;
        $taji1 = $pi - $p1;
        $taji2 = (($e1 * $polyloci) + (($e2 * $polyloci) * ($polyloci - 1)));
        $taji3 = sqrt(abs($taji2));
        #$taji3 = sqrt($taji2);
        $tajimaD= $taji1 / $taji3;
        print "$chr\t$start\t$end\t$polyloci\t$pairwise\t$p1\t$pi\t$taji1\t$taji3\t$tajimaD\n";
       }
     else 
      { 
       print "$chr\t$start\t$end\t$polyloci\t$pairwise\t$p1\t$pi\t$taji1\t$taji3\t$tajimaD\n";
      }
  }