perl/etc/WoodyMiaoLin/PbeBefore2015/Structure_histogram.pl

   1 #!/usr/bin/perl
   2 use strict;
   3 use warnings;
   4
   5 my $in = shift;
   6 open I, "<", $in;
   7 while (<I>) {
   8         last if /\(%Miss\) :  Inferred clusters/;
   9 }
  10
  11 my @color = ("#FF0000", "#00FF00", "#0000FF", "#00FFFF", "#FF00FF", "#FFFF00", "#800000", "#008000", "#000080", "#008080", "#800080", "#808000");
  12 my @PO; # Probability of inferred clusters with original order
  13 my @pg; # Probability of inferred clusters grouped by numerical order of the max probability
  14 my $num_clust; # Number of clusters
  15 my $num_indiv; # Number of individuals
  16 while (<I>) {
  17         chomp;
  18         last unless $_;
  19         my @a = split /:  /;
  20         $a[1] =~ s/ $//;
  21         my @b = split / /, $a[1];
  22         $num_clust = @b unless $num_clust;
  23         ++$num_indiv;
  24         unshift @b, $num_indiv;
  25         my $max_order = 1; # The numerical order of the max probability
  26         foreach (1 .. $num_clust) {
  27                 if ($b[$_] > $b[$max_order]) {
  28                         $max_order = $_;
  29                 }
  30         }
  31         push @{$pg[$max_order]}, \@b;
  32         push @PO, \@b;
  33 }
  34 close I;
  35
  36 my @PQ; # Probability of inferred clusters sort by Q
  37 foreach (1 .. $num_clust) {
  38         next unless $pg[$_];
  39         push @PQ, sort {${$b}[$_] <=> ${$a}[$_]} @{$pg[$_]};
  40 }
  41
  42 open O, ">", "$in.svg";
  43 # SVG attribute definitions
  44 my $w = 10 * $num_indiv + 20;
  45 my $h = 240;
  46 print O "<?xml version=\"1.0\"?>\n";
  47 print O "<svg xmlns=\"http://www.w3.org/2000/svg\" width=\"$w\" height=\"$h\">\n";
  48 print O "<g transform=\"translate(10,10)\">\n\n";
  49
  50 # Print histogram
  51 foreach my $indiv (0 .. $num_indiv-1) {
  52         my $xr = 10*($indiv);
  53         my $yrPO = 0;
  54         my $yrPQ = 115;
  55         my $xt = $xr + 1;
  56         foreach my $clust (1 .. $num_clust) {
  57                 my $hPO = 100*$PO[$indiv][$clust];
  58                 my $hPQ = 100*$PQ[$indiv][$clust];
  59                 print O "<rect x=\"$xr\" y=\"$yrPO\" width=\"10\" height=\"$hPO\" fill=\"$color[$clust-1]\"/>\n";
  60                 print O "<rect x=\"$xr\" y=\"$yrPQ\" width=\"10\" height=\"$hPQ\" fill=\"$color[$clust-1]\"/>\n";
  61                 $yrPO += $hPO;
  62                 $yrPQ += $hPQ
  63         }
  64         my $numPO = sprintf("%03d", $PO[$indiv][0]);
  65         my $numPQ = sprintf("%03d", $PQ[$indiv][0]);
  66         print O "<text x=\"$xt\" y=\"105\" font-family=\"Courier\" font-size=\"4\" fill=\"black\">$numPO</text>\n";
  67         print O "<text x=\"$xt\" y=\"220\" font-family=\"Courier\" font-size=\"4\" fill=\"black\">$numPQ</text>\n";
  68 }
  69
  70 print O "\n</g>\n</svg>\n";
  71 close O;