perl/etc/WoodyMiaoLin/Tiger/compare.pl

   1 #!/usr/bin/perl
   2 use strict;
   3 use warnings;
   4
   5 open IP, "<", "phase.ht";
   6 open IM, "<", "manual.ht";
   7
   8 my (@id, @phase, @manual, @compare, @conflict, @out);
   9
  10 my $n = 0;
  11 while (<IP>) {
  12         chomp;
  13         s/^0 //;
  14         push @id, $_;
  15         my $h1 = <IP>;
  16         my $h2 = <IP>;
  17         chomp ($h1, $h2);
  18         $phase[$n][0] = [split / /, $h1];
  19         $phase[$n][1] = [split / /, $h2];
  20         ++$n;
  21 }
  22 close IP;
  23
  24 while (<IM>) {
  25         chomp;
  26         my @a = split /\t/;
  27         my $n = 0;
  28         foreach (@a) {
  29                 if ($n == int($n)) {
  30                         push @{$manual[int($n)][0]}, $_;
  31                 } else {
  32                         push @{$manual[int($n)][1]}, $_;
  33                 }
  34                 $n += 0.5;
  35         }
  36 }
  37 close IM;
  38
  39 foreach my $n (0 .. 11) {
  40         $conflict[$n][0] = 0;
  41         $conflict[$n][1] = 0;
  42         foreach my $p (0 .. 130) {
  43                 foreach my $a1 (0 .. 1) {
  44                         foreach my $a2 (0 .. 1) {
  45                                 if ($phase[$n][$a1][$p] eq $manual[$n][$a2][$p]) {
  46                                         $compare[$n]{"c$a1$a2"}[$p] = 0;
  47                                 } else {
  48                                         $compare[$n]{"c$a1$a2"}[$p] = 1;
  49                                 }
  50                         }
  51                 }
  52                 $conflict[$n][0] += $compare[$n]{c00}[$p] + $compare[$n]{c11}[$p];
  53                 $conflict[$n][1] += $compare[$n]{c01}[$p] + $compare[$n]{c10}[$p];
  54         }
  55         if ($conflict[$n][1] > $conflict[$n][0]) {
  56                 $out[$n][0][0] = $phase[$n][0];
  57                 $out[$n][0][1] = $manual[$n][0];
  58                 $out[$n][0][2] = $compare[$n]{c00};
  59                 $out[$n][1][0] = $phase[$n][1];
  60                 $out[$n][1][1] = $manual[$n][1];
  61                 $out[$n][1][2] = $compare[$n]{c11};
  62         } else {
  63                 $out[$n][0][0] = $phase[$n][0];
  64                 $out[$n][0][1] = $manual[$n][1];
  65                 $out[$n][0][2] = $compare[$n]{c01};
  66                 $out[$n][1][0] = $phase[$n][1];
  67                 $out[$n][1][1] = $manual[$n][0];
  68                 $out[$n][1][2] = $compare[$n]{c10};
  69         }
  70 }
  71
  72 open O, ">", "compare.ht";
  73
  74 foreach (0 .. 11) {
  75         print O "$id[$_]\n";
  76         print O "hyplotype0\tphase\t", (join " ", @{$out[$_][0][0]}), "\n";
  77         print O "hyplotype0\tmanual\t", (join " ", @{$out[$_][0][1]}), "\n";
  78         print O "hyplotype0\tconfl.\t", (join " ", @{$out[$_][0][2]}), "\n\n";
  79         print O "hyplotype1\tphase\t", (join " ", @{$out[$_][1][0]}), "\n";
  80         print O "hyplotype1\tmanual\t", (join " ", @{$out[$_][1][1]}), "\n";
  81         print O "hyplotype1\tconfl.\t", (join " ", @{$out[$_][1][2]}), "\n\n";
  82 }
  83 close O;