|
0
|
1 use strict;
|
|
|
2 use Cwd qw(cwd);
|
|
|
3
|
|
|
4
|
|
|
5 my $file_aff=shift;
|
|
|
6 my $file_cont=shift;
|
|
|
7 my $ofile=shift;
|
|
|
8
|
|
|
9 my $dir=cwd;
|
|
|
10
|
|
|
11 open(O,">$ofile.tmp");
|
|
|
12 my %Final_data=();
|
|
|
13 my @ids=();
|
|
|
14
|
|
|
15
|
|
|
16
|
|
|
17 my $ND=populate($file_aff);
|
|
|
18 my $NH=populate($file_cont);
|
|
|
19
|
|
|
20 my $head=join("\t",@ids);
|
|
|
21 print O "\t$head\n";
|
|
|
22 foreach my $gene (sort keys %Final_data)
|
|
|
23 {
|
|
|
24 print O "$gene\t";
|
|
|
25 foreach (my $j=0;$j<=$#ids;$j++)
|
|
|
26 {
|
|
|
27 my $individual=$ids[$j];
|
|
|
28 my $SCORE=$Final_data{$gene}{$individual} ? $Final_data{$gene}{$individual} : 0;
|
|
|
29 if ($j==$#ids)
|
|
|
30 {
|
|
|
31 print O "$SCORE\n";
|
|
|
32 }else{
|
|
|
33 print O "$SCORE\t";
|
|
|
34 }
|
|
|
35 }
|
|
|
36 }
|
|
|
37
|
|
|
38 #print "Rscript --vanilla $dir/PCA.R $ofile $ND $NH $ofile.png\n";
|
|
|
39 system ("Rscript --vanilla $dir/PCA.R $ofile.tmp $ND $NH $ofile")==0||die($!);
|
|
|
40
|
|
|
41 sub populate
|
|
|
42 {
|
|
|
43 my $file=$_[0];
|
|
|
44 open(IN,$file);
|
|
|
45 my @P=();
|
|
|
46 my $N=0;
|
|
|
47 while(<IN>)
|
|
|
48 {
|
|
|
49 if ($_=~/^#CHROM/)
|
|
|
50 {
|
|
|
51 my ($chr,$pos,$id,$ref,$alt,$qual,$filter,$info,$format,@Pids)=(split());
|
|
|
52 foreach my $P (@Pids)
|
|
|
53 {
|
|
|
54 push(@ids,$P);
|
|
|
55 push(@P,$P);
|
|
|
56 $N++;
|
|
|
57 }
|
|
|
58
|
|
|
59 }else{
|
|
|
60 my ($chr,$pos,$id,$ref,$alt,$qual,$filter,$info,$format,@data)=(split());
|
|
|
61 my @infos=split(/\;/,$info);
|
|
|
62 my $gene=".";
|
|
|
63 my $score=0;
|
|
|
64 foreach my $i (@infos)
|
|
|
65 {
|
|
|
66 if ($i=~/Gene.refGene/)
|
|
|
67 {
|
|
|
68 $gene=(split(/\=/,$i))[1];
|
|
|
69 }elsif ($i=~/VINYL_score/){
|
|
|
70 $score=(split(/\=/,$i))[1];
|
|
|
71 }
|
|
|
72 }
|
|
|
73 next if $gene eq ".";
|
|
|
74 foreach (my $j=0;$j<=$#data;$j++)
|
|
|
75 {
|
|
|
76 my $individual=$P[$j];
|
|
|
77 my $call=$data[$j];
|
|
|
78 next if $call eq "." || $call eq "0|0";
|
|
|
79 if ($Final_data{$gene}{$individual})
|
|
|
80 {
|
|
|
81 $Final_data{$gene}{$individual}=$score if $score>$Final_data{$gene}{$individual}
|
|
|
82 }else{
|
|
|
83 $Final_data{$gene}{$individual}=$score;
|
|
|
84 }
|
|
|
85
|
|
|
86 }
|
|
|
87 }
|
|
|
88 }
|
|
|
89 return($N);
|
|
|
90 }
|
