perl/etc/WoodyMiaoLin/Tiger/print_svg_linkage_map.pl

   1 #!/usr/bin/perl
   2 use strict;
   3 use warnings;
   4
   5 die "Usage: $0 <input_bcf> <scaffold> <output_svg>\n" if @ARGV<3;
   6
   7 my @samples = qw/ JHH001 GZXJ03 GZXJ05 GZXJ26 GZXJ27 GZXJ28 GZXJ29 GZXJ30 GZXJ33 BHX011 BHX019 GZXJ04 GZXJ06 GZXJ31 GZXJ32 GZXJ36 GZXJ37 GZXJ38 /;
   8 my $ns = @samples;
   9 my $in = shift;
  10 my $sc = shift;
  11 my $out = shift;
  12
  13 open I, "-|", "bcftools view -I $in";
  14 open O, ">", "$out";
  15
  16 # Read genotype from bcf file, and push into @gt;
  17 my @gt;
  18 while (<I>) {
  19         next unless /^$sc\t/;
  20         chomp;
  21         my @line = split /\t/;
  22         next if $line[5] < 20;
  23         $line[7] =~ /;FQ=([0-9\-.]+)/;
  24         next unless $1 > 0;
  25         my $a = 0;
  26         my @a;
  27         foreach (9 .. ($ns+8)) {
  28                 my @sample = split /:/, $line[$_];
  29                 if ($sample[4] >= 20 and $sample[2] >= 1) {
  30                         ++$a;
  31                         push @a, $sample[0];
  32                 }
  33         }
  34         next if $a < $ns;
  35         next if (!grep {$_ ne $a[0]} @a);
  36         my @b;
  37         $b[0] = $line[0];
  38         $b[1] = $line[1];
  39         my $b0 = 0;
  40         my $b1 = 0;
  41         foreach (0 .. ($ns-1)) {
  42                 my ($a1, $a2) = split /\//, $a[$_];
  43                 if (($a1 eq "0") and ($a2 eq "0")) {
  44                         ++$b0;
  45                         push @b, $line[3];
  46                 } elsif (($a1 eq "1") and ($a2 eq "1")) {
  47                         ++$b1;
  48                         push @b, $line[4];
  49                 } else {
  50                         push @b, "$line[3]$line[4]";
  51                 }
  52         }
  53         if ($b0 > $b1) {
  54                 push @b, $line[3];
  55         } elsif ($b0 < $b1) {
  56                 push @b, $line[4];
  57         } else {
  58                 push @b, "equal";
  59         }
  60         push @gt, \@b;
  61 }
  62 close I;
  63
  64 # SVG attribute definitions;
  65 my $w = 30*$ns + 130;
  66 my $h = 10*@gt + 20;
  67 print O "<?xml version=\"1.0\"?>\n";
  68 print O "<svg xmlns=\"http://www.w3.org/2000/svg\" width=\"$w\" height=\"$h\">\n";
  69 print O "<g transform=\"translate(10,10)\">\n\n";
  70
  71 # Print title;
  72 my $x0 = 112;
  73 foreach (@samples) {
  74         print O "<text x=\"$x0\" y=\"0\" font-family=\"Courier\" font-size=\"7\" fill=\"black\">$_</text>\n";
  75         $x0 += 30;
  76 }
  77
  78 # Print genotype;
  79 my $yy;
  80 foreach my $a (@gt) {
  81         ++$yy;
  82         my $y = 10*$yy;
  83         print O "<text x=\"0\" y=\"$y\" font-family=\"Courier\" font-size=\"8\" fill=\"black\">${$a}[0]</text>\n";
  84         print O "<text x=\"70\" y=\"$y\" font-family=\"Courier\" font-size=\"8\" fill=\"black\">${$a}[1]</text>\n";
  85         foreach (2 .. ($ns+1)) {
  86                 my $xt = 30*$_ + 60;
  87                 my $yt = 10*$yy;
  88                 my $xr = 30*$_ + 50;
  89                 my $yr = 10*$yy - 8;
  90                 my $h = length ${$a}[$_];
  91                 my $color;
  92                 if ($h == 2) {
  93                         $color = "green";
  94                 } elsif ($h == 1) {
  95                         if (${$a}[$ns+2] eq "equal") {
  96                                 $color = "pink";
  97                         } elsif (${$a}[$ns+2] eq ${$a}[$_]) {
  98                                 $color = "yellow";
  99                         } else {
 100                                 $color = "blue";
 101                         }
 102                 } else {
 103                         $color = "white";
 104                 }
 105         print O "<rect x=\"$xr\" y=\"$yr\" width=\"30\" height=\"10\" fill=\"$color\" stroke-width=\"0.5\" stroke=\"black\"/>\n";
 106         print O "<text x=\"$xt\" y=\"$yt\" font-family=\"Courier\" font-size=\"8\" fill=\"black\">${$a}[$_]</text>\n";
 107         }
 108 }
 109 print O "\n</g>\n</svg>\n";
 110 close O;