|
0
|
1 $ofile=shift;
|
|
|
2 $plot1=shift;
|
|
|
3 $plot2=shift;
|
|
|
4 @genomes=();
|
|
|
5 foreach $var_file (@ARGV)
|
|
|
6 {
|
|
|
7 open(IN,$var_file);
|
|
|
8 $genome="genome";
|
|
|
9 while(<IN>)
|
|
|
10 {
|
|
|
11 ($pos,$ref,$alt,$gen)=(split(/\s+/))[1,2,3,-1];
|
|
|
12 next unless $ref=~/[ACTG]/ && $alt=~/[ACTG]/;
|
|
|
13 #print $pos $ref $alt\n;
|
|
|
14 if ($genome eq "genome")
|
|
|
15 {
|
|
|
16 push(@genomes,$gen);
|
|
|
17 $genome=$gen;
|
|
|
18 }
|
|
|
19 $is_present{"$pos\_$ref|$alt"}{$genome}=1;
|
|
|
20 }
|
|
|
21 }
|
|
|
22
|
|
|
23 @test=qw(241_C|T 514_T|C 1059_C|T 1397_G|A 1440_G|A 1605_ATG|... 2416_C|T 2480_A|G 2558_C|T 2891_G|A 3037_C|T 8782_C|T 9477_T|A 10097_G|A 11083_G|T 11916_C|T 14408_C|T 14805_C|T 15324_C|T 17247_T|C 17747_C|T 17858_A|G 18060_C|T 18877_C|T 18998_C|T 20268_A|G 23403_A|G 23731_C|T 24034_C|T 25429_G|T 25563_G|T 25979_G|T 26144_G|T 27046_C|T 27964_C|T 28144_T|C 28311_C|T 28657_C|T 28688_T|C 28851_G|T 28854_C|T 28863_C|T 28881_GGG|AAC 29540_G|A 29553_G|A 29742_G|T);
|
|
|
24
|
|
|
25 open(OUT,">$ofile");
|
|
|
26 print OUT " @test\n";
|
|
|
27
|
|
|
28 foreach $genome (@genomes)
|
|
|
29 {
|
|
|
30 $ostring="$genome ";
|
|
|
31 foreach $t (@test)
|
|
|
32 {
|
|
|
33 $val=$is_present{$t}{$genome} ? 1 : 0;
|
|
|
34 $ostring.="$val ";
|
|
|
35 }
|
|
|
36 chop($ostring);
|
|
|
37 print OUT "$ostring\n";
|
|
|
38 }
|
|
|
39 system("Rscript --vanilla bubbleplot.R $ofile $plot1 $plot2")==0||die("no plot");
|
