perl/etc/WoodyMiaoLin/Tiger/drawSynteny2.pl

   1 #!/usr/bin/perl
   2 use strict;
   3 use warnings;
   4
   5 die "Usage: $0 <gene> <SNP> <output>\n" if @ARGV<3;
   6 my $in1=shift;
   7 my $in2=shift;
   8 my $out=shift;
   9
  10 open I1, "<", $in1 or die "Error opening $in1 : $!\n";
  11 open I2, "<", $in2 or die "Error opening $in2 : $!\n";
  12 open O, ">", $out or die "Error opening $out : $!\n";
  13
  14
  15 my $nSpeci = 4; # 物种数
  16 my $nGene = 22; # 基因数
  17 my $yScale = 0.0002; # y轴方向的缩放系数，等于基因的高度除以DNA碱基数
  18 my $genWid = 10; # 基因矩形宽度的一半
  19 my $spaChr = 200; # 染色体间距
  20
  21 # 读取染色体与基因信息
  22 my @yCoo;
  23 for my $a (0..$nSpeci) {
  24   for my $b (0..$nGene) {
  25     my $c = <I1>;
  26     chomp $c;
  27     @{$yCoo[$a][$b]} = split /\t/, $c;
  28   }
  29 }
  30 close I1;
  31
  32 # 读取SNP信息
  33 my %SNPcoo;
  34 while (<I2>) {
  35   chomp;
  36   my @a = split /\t/;
  37   $SNPcoo{$a[0]} = $a[1];
  38 }
  39 close I2;
  40
  41 print O "<?xml version=\"1.0\"?>\n\n<svg xmlns=\"http://www.w3.org/2000/svg\">\n\n";
  42
  43 # 画染色体 Scaffold75长5658748
  44 for (2..$nSpeci) {
  45   my $x = $spaChr*$_;
  46   my $y1 = 0;
  47   my $y2 = 9302904*$yScale;
  48   print O "<line x1=\"$x\" y1=\"$y1\" x2=\"$x\" y2=\"$y2\" stroke=\"black\" stroke-width=\"5\"/>\n";
  49 }
  50 print O "<line x1=\"$spaChr\" y1=\"0\" x2=\"$spaChr\" y2=\"", 5658748*$yScale, "\" stroke=\"black\" stroke-width=\"5\"/>\n";
  51 print O "<line x1=\"$spaChr\" y1=\"", 5758748*$yScale, "\" x2=\"$spaChr\" y2=\"", 9302904*$yScale, "\" stroke=\"black\" stroke-width=\"5\"/>\n";
  52
  53 # 画基因
  54 for my $a (1..$nSpeci) {
  55   for my $b (1..$nGene) {
  56     my $x = $spaChr*$a - $genWid;
  57     my $y = $yScale*$yCoo[$a][$b][0];
  58     my $w = 2*$genWid;
  59     my $h = $yScale*($yCoo[$a][$b][1]-$yCoo[$a][$b][0]);
  60     print O "<rect x=\"$x\" y=\"$y\" width=\"$w\" height=\"$h\" fill=\"black\"/>\n";
  61   }
  62 }
  63
  64 # 画基因连线
  65 for my $b (1..$nGene) {
  66   for my $a (1..$nSpeci-1) {
  67     for my $c (0..1) {
  68       my $x1 = $spaChr*$a + $genWid;
  69       my $y1 = $yScale*$yCoo[$a][$b][$c];
  70       my $x2 = $spaChr*($a+1) - $genWid;
  71       my $y2 = $yScale*$yCoo[$a+1][$b][$c];
  72       print O "<line x1=\"$x1\" y1=\"$y1\" x2=\"$x2\" y2=\"$y2\" stroke=\"black\" stroke-width=\"0.5\"/>\n";
  73     }
  74   }
  75 }
  76
  77 # 打印基因名
  78 for (1..$nGene) {
  79   my $x = $nSpeci*$spaChr + 2*$genWid;
  80   my $y = 0.5*$yScale*($yCoo[$nSpeci][$_][0] + $yCoo[$nSpeci][$_][1]) + 2;
  81   my $gene = $yCoo[0][$_][0];
  82   print O "<text x=\"$x\" y=\"$y\" font-family=\"Courier\" font-size=\"8\" fill=\"black\">$gene</text>\n";
  83 }
  84
  85 # 打印染色体名
  86 for (1..$nSpeci) {
  87   my $x = ($_-0.4)*$spaChr;
  88   my $scaffold = $yCoo[$_][0][0];
  89   print O "<text x=\"$x\" y=\"-10\" font-family=\"Courier\" font-size=\"12\" fill=\"black\">$scaffold</text>\n";
  90 }
  91
  92 # 标记SNP
  93 for (sort {$a <=> $b} keys %SNPcoo) {
  94   my $x1 = $spaChr - 3*$genWid;
  95   my $x2 = $spaChr;
  96   my $y = $yScale*$_;
  97   print O "<line x1=\"$x1\" y1=\"$y\" x2=\"$x2\" y2=\"$y\" stroke=\"black\" stroke-width=\"1\"/>\n";
  98   my $SNP = "$_($SNPcoo{$_})";
  99   my $ty = $yScale*$_ + 2;
 100   print O "<text x=\"120\" y=\"$ty\" font-family=\"Courier\" font-size=\"8\" fill=\"black\">$SNP</text>\n";
 101 }
 102
 103 print O "\n</svg>\n";
 104 close O;