#!usr/bin/perl
open (IN, "HMOX2_1000genomes"); open (IN1, "positions_red"); open (IN2, "positions_green"); open(IN3, "positions_blue");
open (OUT, ">H_ind");
@a = ();
while (<IN>){ chomp $_; push (@a, $_);}
$x = 0;
for $n (0..$#a){
        $point = 0;
        if ($a[$n] =~m/^#CHROM/){ $start = $n; $point = 1;}
        last if $point == 1;
}
for $n ($start..$#a){
        @{'line'.$x} = split (/\t/,$a[$n]);
        $x++ unless $n == $#a;
}
$c = 1; @coord = @reverse_SNPs = @{'h'.$c} = @{'H'.$c} = @{'None'.$c} = ();
#while (<IN1>){  chomp $_ ; unless ($_ =~m/^#\d/){ push (@coord, $_); push (@reverse_SNPs, $_);}}
#while (<IN2>){  chomp $_ ; unless ($_ =~m/^#\d/){ push (@coord, $_); }}
while (<IN3>){  chomp $_ ; unless ($_ =~m/^#\d/){ push (@coord, $_); }}
for $n (0..$#coord){
        $flag = 0; $flag2 = 0;
        for $z (1..$x){
                if (${'line'.$z}[1] =~m/^($coord[$n])/){
                        for $k (0..$#reverse_SNPs){ if ($coord[$n] == $reverse_SNPs[$k]){$flag2 = 1;}}
                        if ($flag2 == 1){
                                for $y (0..$#{'line'.$z}){
                                        if (${'line'.$z}[$y] =~m/^(\d)\|(\d).*/){
                                                if ($1 == 0 && $2 == 0){ push ( @{'H'.$c}, $line0[$y]);}
                                                elsif ($1 == 0 && $2 == 1){ push ( @{'h'.$c}, $line0[$y]);}
                                                elsif ($1 == 1 && $2 == 0){ push ( @{'h'.$c}, $line0[$y]);}
                                                elsif ($1 == 1 && $2 == 1){ push ( @{'None'.$c}, $line0[$y]);}
                                                $flag = 1;
                                        }
                                }
                        }
                        else{
                                for $y (0..$#{'line'.$z}){
                                        if (${'line'.$z}[$y] =~m/^(\d)\|(\d).*/){
                                                if ($1 == 0 && $2 == 0){ push ( @{'None'.$c}, $line0[$y]);}
                                                elsif ($1 == 0 && $2 == 1){ push ( @{'h'.$c}, $line0[$y]);}
                                                elsif ($1 == 1 && $2 == 0){ push ( @{'h'.$c}, $line0[$y]);}
                                                elsif ($1 == 1 && $2 == 1){ push ( @{'H'.$c}, $line0[$y]);}
                                                $flag = 1;
                                        }
                                }
                        }
                }
        }
        if ($flag == 1){$c++ unless $c == $#coord+1 ;}
}
print "\tNone\th\tH\n";
for $v (1..$c){
        $p = $#{'None'.$v}+1; $q = $#{'h'.$v}+1; $r = $#{'H'.$v}+1;
        print $v,".","\t",$p,"\t",$q,"\t",$r,"\n";
}
%N = %h = %H = @homozygous = @heterozygous = @common_homo = @common_hetero = @none_haplotypes = @common_none = ();
for $n (1..$c){
        push (@none_haplotypes, @{'None'.$n});
        push (@heterozygous , @{'h'.$n});
        push (@homozygous, @{'H'.$n});
}
foreach $x (@none_haplotypes){ $N{$x}++;}
foreach $x (@heterozygous){ $h{$x}++;}
foreach $x (@homozygous){ $H{$x}++;}
foreach $k (keys %N){
        if ($N{$k} == $c){ push (@common_none , $k);}
}
foreach $k (keys %h){
        if ($h{$k} == $c){ push (@common_hetero , $k);}
}
foreach $k (keys %H){
        if ($H{$k} == $c){ push (@common_homo , $k);}
}
$N = $#common_none +1; $h = $#common_hetero +1; $H = $#common_homo +1;
print "Common ancestral homozygous = $N , Common heterozygous = $h , common mutated homozygous = $H\n";
for $f (0..$#common_homo){ print OUT $common_homo[$f],"\n";}