perl/etc/WoodyMiaoLin/structure_histogram/structure_histogram.pl

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