Mercurial > repos > elixir-it > vinyl_pca

comparison makeToM.pl @ 0:0011da72f65a draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author elixir-it
date Tue, 09 Jun 2020 16:11:33 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:0011da72f65a
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 }