|
0
|
1 #!/usr/bin/perl -w
|
|
|
2 #
|
|
|
3
|
|
|
4 use strict;
|
|
|
5
|
|
|
6
|
|
|
7 my $f1=shift; #varscan
|
|
|
8 my $f2=shift; #freebayes
|
|
|
9 my $f3=shift; #GATK
|
|
|
10 my $fout=shift;
|
|
|
11 my %ALL=();
|
|
|
12 my %L=();
|
|
|
13 my %H=();
|
|
|
14
|
|
|
15 die("non trovo il file $f1") unless -e $f1;
|
|
|
16 die("non trovo il file $f2") unless -e $f2;
|
|
|
17 die("non trovo il file $f2") unless -e $f3;
|
|
|
18
|
|
|
19 die ("non specificato il file di output") unless $fout;
|
|
|
20
|
|
|
21
|
|
|
22 open(OUT,">$fout.common.snps.vcf");
|
|
|
23 open(UN,">$fout.unique.snps.vcf");
|
|
|
24
|
|
|
25 open(IN,$f1);
|
|
|
26 while(<IN>)
|
|
|
27 {
|
|
|
28 next if $_=~/^\#/;
|
|
|
29 my ($chr,$st,$b1,$b2)=(split())[0,1,3,4];
|
|
|
30 my @all_vl=(split(/\t/));
|
|
|
31 next unless (length($b1)==length($b2));
|
|
|
32 next if ($b1=~/\./ || $b1=~/\-/ || $b2=~/\./ || $b2=~/\-/);
|
|
|
33 my @b1=(split('',$b1));
|
|
|
34 my @b2=(split('',$b2));
|
|
|
35 for (my $i=0;$i<=$#b1;$i++)
|
|
|
36 {
|
|
|
37 my $lb1=$b1[$i];
|
|
|
38 my $lb2=$b2[$i];
|
|
|
39 next if $lb1 eq $lb2;
|
|
|
40 next if ($lb1=~/\./ || $lb1=~/\-/ || $lb2=~/\./ || $lb2=~/\-/ || $lb2=~/\,/);
|
|
|
41 my $lst=$st+$i;
|
|
|
42 $all_vl[1]=$lst;
|
|
|
43 $all_vl[3]=$lb1;
|
|
|
44 $all_vl[4]=$lb2;
|
|
|
45 my $J=join("\t",@all_vl);
|
|
|
46 $H{$chr}{$lst}{$lb1}{$lb2}++;
|
|
|
47 $ALL{$chr}{$lst}{$lb1}{$lb2}=$J;
|
|
|
48 push(@{$L{$chr}{$lst}{$lb1}{$lb2}},"VS");
|
|
|
49 }
|
|
|
50 }
|
|
|
51 close(IN);
|
|
|
52
|
|
|
53 open(IN,$f2);
|
|
|
54 while(<IN>)
|
|
|
55 {
|
|
|
56 next if $_=~/^\#/;
|
|
|
57 my ($chr,$st,$b1,$b2)=(split())[0,1,3,4];
|
|
|
58 my @all_vl=(split(/\t/));
|
|
|
59 next unless (length($b1)==length($b2));
|
|
|
60 next if ($b1=~/\./ || $b1=~/\-/ || $b2=~/\./ || $b2=~/\-/);
|
|
|
61 my @b1=(split('',$b1));
|
|
|
62 my @b2=(split('',$b2));
|
|
|
63
|
|
|
64 for (my $i=0;$i<=$#b1;$i++)
|
|
|
65 {
|
|
|
66 my $lb1=$b1[$i];
|
|
|
67 my $lb2=$b2[$i];
|
|
|
68 next if $lb1 eq $lb2;
|
|
|
69 next if ($lb1=~/\./ || $lb1=~/\-/ || $lb2=~/\./ || $lb2=~/\-/ || $lb2=~/\,/);
|
|
|
70 my $lst=$st+$i;
|
|
|
71 $all_vl[1]=$lst;
|
|
|
72 $all_vl[3]=$lb1;
|
|
|
73 $all_vl[4]=$lb2;
|
|
|
74 my $J=join("\t",@all_vl);
|
|
|
75 $H{$chr}{$lst}{$lb1}{$lb2}++;
|
|
|
76 $ALL{$chr}{$lst}{$lb1}{$lb2}=$J;
|
|
|
77 push(@{$L{$chr}{$lst}{$lb1}{$lb2}},"FB");
|
|
|
78 }
|
|
|
79
|
|
|
80 }
|
|
|
81 close(IN);
|
|
|
82
|
|
|
83 open(IN,$f3);
|
|
|
84 while(<IN>)
|
|
|
85 {
|
|
|
86 if ($_=~/^\#/)
|
|
|
87 {
|
|
|
88 print OUT if $_=~/fileformat/;
|
|
|
89 print OUT if $_=~/contig/;
|
|
|
90 print OUT if $_=~/CHROM/;
|
|
|
91 print UN if $_=~/fileformat/;
|
|
|
92 print UN if $_=~/contig/;
|
|
|
93 print UN if $_=~/CHROM/;
|
|
|
94 next;
|
|
|
95 }
|
|
|
96 my ($chr,$st,$b1,$b2)=(split())[0,1,3,4];
|
|
|
97 my @all_vl=(split(/\t/));
|
|
|
98 next unless (length($b1)==length($b2));
|
|
|
99 next if ($b1=~/\./ || $b1=~/\-/ || $b2=~/\./ || $b2=~/\-/);
|
|
|
100 my @b1=(split('',$b1));
|
|
|
101 my @b2=(split('',$b2));
|
|
|
102
|
|
|
103 for (my $i=0;$i<=$#b1;$i++)
|
|
|
104 {
|
|
|
105 my $lb1=$b1[$i];
|
|
|
106 my $lb2=$b2[$i];
|
|
|
107 next if $lb1 eq $lb2;
|
|
|
108 next if ($lb1=~/\./ || $lb1=~/\-/ || $lb2=~/\./ || $lb2=~/\-/ || $lb2=~/\,/);
|
|
|
109 my $lst=$st+$i;
|
|
|
110 $all_vl[1]=$lst;
|
|
|
111 $all_vl[3]=$lb1;
|
|
|
112 $all_vl[4]=$lb2;
|
|
|
113 my $J=join("\t",@all_vl);
|
|
|
114 $H{$chr}{$lst}{$lb1}{$lb2}++;
|
|
|
115 $ALL{$chr}{$lst}{$lb1}{$lb2}=$J;
|
|
|
116 push(@{$L{$chr}{$lst}{$lb1}{$lb2}},"GA");
|
|
|
117 }
|
|
|
118 }
|
|
|
119 close(IN);
|
|
|
120
|
|
|
121 foreach my $chr (sort keys %H)
|
|
|
122 {
|
|
|
123 foreach my $s (sort{$a<=>$b} keys %{$H{$chr}})
|
|
|
124 {
|
|
|
125 foreach my $b1 (keys %{$H{$chr}{$s}})
|
|
|
126 {
|
|
|
127 foreach my $b2 (keys %{$H{$chr}{$s}{$b1}})
|
|
|
128 {
|
|
|
129
|
|
|
130 my $vl=$H{$chr}{$s}{$b1}{$b2};
|
|
|
131 unless ($H{$chr}{$s}{$b1}{$b2})
|
|
|
132 {
|
|
|
133 warn("non trovo il numero di occorrenze di $chr $s $b1 $b2\n");
|
|
|
134 next;
|
|
|
135 }
|
|
|
136 my $out_vl=$ALL{$chr}{$s}{$b1}{$b2};
|
|
|
137
|
|
|
138 unless ($ALL{$chr}{$s}{$b1}{$b2})
|
|
|
139 {
|
|
|
140 warn("non trovo le info di $chr $s $b1 $b2\n");
|
|
|
141 next;
|
|
|
142 }
|
|
|
143 my @vls=(split(/\s+/,$out_vl));
|
|
|
144 my $comm=$vls[7];
|
|
|
145
|
|
|
146 if ($vl>=2)
|
|
|
147 {
|
|
|
148 my $mm=join("_",@{$L{$chr}{$s}{$b1}{$b2}});
|
|
|
149 my $comm2=$comm;
|
|
|
150 $comm2.=";NM=$vl;LM=$mm";
|
|
|
151 $vls[7]=$comm2;
|
|
|
152 my $out_vl=join("\t",@vls);
|
|
|
153 print OUT "$out_vl\n";
|
|
|
154
|
|
|
155 }else{
|
|
|
156 my $mm=$L{$chr}{$s}{$b1}{$b2}[0];
|
|
|
157 my $comm2=$comm;
|
|
|
158 $comm2.=";NM=$vl;LM=$mm";
|
|
|
159 $vls[7]=$comm2;
|
|
|
160 my $out_vl=join("\t",@vls);
|
|
|
161 print UN "$out_vl\n";
|
|
|
162
|
|
|
163 }
|
|
|
164 }
|
|
|
165 }
|
|
|
166 }
|
|
|
167 }
|
|
|
168
|
|
|
169 #foreach $V (keys %V)
|
|
|
170 #{
|
|
|
171 # print "$V $V{$V}\n";
|
|
|
172 #}
|
