|
0
|
1 #!/usr/bin/perl
|
|
|
2
|
|
|
3 use strict;
|
|
1
|
4 # use lib 'inc/perlmod';
|
|
|
5 # use ngsutil qw[ :DEFAULT &explode_varcall ];
|
|
0
|
6 use warnings FATAL => qw[ numeric uninitialized ];
|
|
|
7 use List::Util qw[ sum min max ];
|
|
|
8 use File::Basename;
|
|
|
9 use Getopt::Long;
|
|
|
10
|
|
|
11 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
12 # PATH TO YOUR R-bin DIRECTORY
|
|
|
13 my $rbin = '/usr/bin/R';
|
|
|
14 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
15
|
|
1
|
16 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
17 # TEMP include ngsutil.pm
|
|
|
18 sub explode_varcall{
|
|
|
19 my $N=0;
|
|
|
20 $_=shift @_ foreach my($POS, $REF, $ALT);
|
|
|
21 $_=$POS foreach my($START, $END);
|
|
|
22 my(@length, @range, @idx, @VAR, @POS);
|
|
|
23 @{$_}=() foreach (\@length, \@range, \@idx, \@VAR, \@POS);
|
|
|
24 push @length, length($_) foreach ($REF, $ALT);
|
|
|
25 @range=sort{ $a<=>$b } @length;
|
|
|
26 if($range[0]==1){
|
|
|
27 if($range[1]!=1){
|
|
|
28 foreach ($REF, $ALT){
|
|
|
29 $_=substr($_, 1);
|
|
|
30 $_=~s/^$/-/;
|
|
|
31 }
|
|
|
32 if($length[0]!=1){
|
|
|
33 $END+=$length[0]-1;
|
|
|
34 $START++;
|
|
|
35 }
|
|
|
36 }
|
|
|
37 push @POS, $START, $END;
|
|
|
38 push @VAR, $REF, $ALT;
|
|
|
39 }else{
|
|
|
40 my @N=();
|
|
|
41 undef $_ foreach my ($i, $VAR);
|
|
|
42 $_-=2 foreach (@length, @range);
|
|
|
43 $_++ foreach ($START, $END);
|
|
|
44 $_=substr($_, 1) foreach ($REF, $ALT);
|
|
|
45 my $indel='-' x ($range[1]-$range[0]);
|
|
|
46 $VAR.=($_>$range[0])?
|
|
|
47 ('-'):((substr($REF, $_, 1) ne substr($ALT, $_, 1))?
|
|
|
48 0:1) for 0 .. $range[1];
|
|
|
49 $N++ while $VAR =~ /0/g;
|
|
|
50 if($length[0]<$length[1]){
|
|
|
51 @VAR=($VAR);
|
|
|
52 @N=($N);
|
|
|
53 $N=0;
|
|
|
54 undef($VAR);
|
|
|
55 $VAR.=($_>$range[0])?
|
|
|
56 ('-'):((substr($REF, $length[0]-$_, 1) ne substr($ALT, $length[1]-$_, 1))?
|
|
|
57 0:1) for reverse 0 .. $range[1];
|
|
|
58 $N++ while $VAR =~ /0/g;
|
|
|
59 if($N>=$N[0]){ $N=shift(@N); $VAR=shift(@VAR); }
|
|
|
60 else{ $REF=$indel . $REF; }
|
|
|
61 }else{ $ALT.=$indel; }
|
|
|
62 foreach (qw[ 0 \- ]){
|
|
|
63 push @idx, [ $-[0], $+[0]-$-[0] ] while ($VAR =~ /$_+/g);
|
|
|
64 }
|
|
|
65 @{$_}=() foreach (\@VAR, \@POS);
|
|
|
66 foreach my $k (@idx){
|
|
|
67 push @VAR, substr($_, ${$k}[0], ${$k}[1]) || '-' foreach ($REF, $ALT);
|
|
|
68 push @POS, ${$k}[0], sum(@{$k})-1;
|
|
|
69 }
|
|
|
70 $_+=$START foreach @POS;
|
|
|
71 $_=~s/\-+/\-/ foreach @VAR;
|
|
|
72 for($i=0; $i<$#POS; $i+=2){ $POS[$i+1]=$POS[$i] if $VAR[$i] eq '-'; }
|
|
|
73 }
|
|
|
74 return(\@POS, \@VAR);
|
|
|
75 }
|
|
|
76
|
|
|
77 sub varscan{
|
|
|
78 $_=shift @_ foreach my($kname, $fpath, $href);
|
|
|
79 my($k, @buffer);
|
|
|
80 open IN, "<$fpath" or die $!;
|
|
|
81 while(<IN>){
|
|
|
82 next if /^#/;
|
|
|
83 chomp;
|
|
|
84 @buffer=split /\s+/, $_;
|
|
|
85 next if !exists $$href{($k=join(':', @buffer[0..2]))};
|
|
|
86 next if $$href{$k}->{ref} !~ $buffer[3];
|
|
|
87 next if $$href{$k}->{alt} !~ $buffer[4];
|
|
|
88 splice(@buffer, 0, 5);
|
|
|
89 $$href{$k}->{$kname}=join(':', @buffer);
|
|
|
90 }
|
|
|
91 close IN;
|
|
|
92 }
|
|
|
93 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
|
94
|
|
0
|
95 my $annovar_dir = 'inc/annovar';
|
|
|
96 my $rdep = 'inc/R';
|
|
|
97
|
|
|
98 my($varfile, $outdir, $outfile, $i, @DP4, @buffer, @Temp, @previous, @fnames, %opts, %chr);
|
|
|
99
|
|
|
100 GetOptions(\%opts, "varfile=s", "outdir=s", "outfile=s");
|
|
|
101 $varfile = $opts{varfile};
|
|
|
102 $outdir = $opts{outdir};
|
|
|
103 $outfile = $opts{outfile};
|
|
|
104
|
|
|
105 my $fname = readlink($varfile) || $varfile;
|
|
|
106 $fname = basename($fname);
|
|
|
107
|
|
|
108 my %fh=(
|
|
|
109 'chr1' => *chr1, 'chr2' => *chr2, 'chr3' => *chr3, 'chr4' => *chr4, 'chr5' => *chr5,
|
|
|
110 'chr6' => *chr6, 'chr7' => *chr7, 'chr8' => *chr8, 'chr9' => *chr9, 'chr10' => *chr10,
|
|
|
111 'chr11' => *chr11, 'chr12' => *chr12, 'chr13' => *chr13, 'chr14' => *chr14, 'chr15' => *chr15,
|
|
|
112 'chr16' => *chr16, 'chr17' => *chr17, 'chr18' => *chr18, 'chr19' => *chr19, 'chr20' => *chr20,
|
|
|
113 'chr21' => *chr21, 'chr22' => *chr22, 'chrX' => *chrX, 'chrY' => *chrY, 'chrM' => *chrM
|
|
|
114 );
|
|
|
115
|
|
|
116 `${annovar_dir}/convert2annovar.pl -format vcf4 $varfile -includeinfo > ${outdir}/${fname}_Temp-00 2> /dev/null` and die $!;
|
|
|
117
|
|
|
118 open($fh{$_}, ">${outdir}/${fname}_${_}.Temp-00") or die $! foreach keys %fh;
|
|
|
119 open IN, "<${outdir}/${fname}_Temp-00" or die $!;
|
|
|
120 while(<IN>){
|
|
|
121 /^(\S+)\s+(?:\S+\s+){2}(\S+)\s+(\S+)/;
|
|
|
122 next if !exists $fh{$1};
|
|
|
123 if(min(length($2), length($3))!=1){
|
|
|
124 chomp;
|
|
|
125 @buffer=split /\s+/, $_;
|
|
|
126 @Temp=explode_varcall(@buffer[1,3..4]);
|
|
|
127 for($i=0; $i<$#{$Temp[0]}; $i+=2){
|
|
|
128 print{ $fh{$buffer[0]} } join("\t", $buffer[0], @{$Temp[0]}[$i..$i+1], @{$Temp[1]}[$i..$i+1], @buffer[6..$#buffer]), "\n";
|
|
|
129 }
|
|
|
130 next;
|
|
|
131 }
|
|
|
132 print{ $fh{$1} } $_;
|
|
|
133 $chr{$1}++;
|
|
|
134 }
|
|
|
135 close IN;
|
|
|
136 foreach (keys %fh){
|
|
|
137 close($fh{$_});
|
|
|
138 next if !exists $chr{$_};
|
|
|
139 `sort -k2,2n -k3,3n ${outdir}/${fname}_${_}.Temp-00 > ${outdir}/${fname}_${_}.Temp-01` and die $!;
|
|
|
140 open IN, "<${outdir}/${fname}_${_}.Temp-01" or die $!;
|
|
|
141 open OUT, ">${outdir}/${fname}_${_}.Temp-02" or die $!;
|
|
|
142 $_=readline(IN);
|
|
|
143 /^((?:\S+\s+){7})(?:\S+\s+){8}(\S+\s+\S+)/;
|
|
|
144 @buffer=split /\s+/, $1.$2;
|
|
|
145 ($_=pop(@buffer))=~s/.+DP4=([^;]+).+/$1/;
|
|
|
146 @DP4=split /,/, $_;
|
|
|
147 push @buffer, @DP4;
|
|
|
148 @previous=@buffer;
|
|
|
149 MAINLOOP: while(<IN>){
|
|
|
150 /^((?:\S+\s+){7})(?:\S+\s+){8}(\S+\s+\S+)/;
|
|
|
151 @buffer=split /\s+/, $1.$2;
|
|
|
152 ($_=pop(@buffer))=~s/.+DP4=([^;]+).+/$1/;
|
|
|
153 @DP4=split /,/, $_;
|
|
|
154 push @buffer, @DP4;
|
|
|
155 while(($previous[0] eq $buffer[0]) && ($buffer[2]==$previous[2]+1) && (join('', @previous[3..4]) !~ /-/) && (join('', @buffer[3..4]) !~ /-/)){
|
|
|
156 $previous[2]=$buffer[2];
|
|
|
157 $previous[$_].=$buffer[$_] for 3..4;
|
|
|
158 $previous[5]='unk' if $previous[5] ne $buffer[5];
|
|
|
159 $previous[7]='SKIP' if $previous[7] ne $buffer[7];
|
|
|
160 for (6,8..11){
|
|
|
161 $previous[$_]+=$buffer[$_];
|
|
|
162 $previous[$_]/=2;
|
|
|
163 }
|
|
|
164 next MAINLOOP;
|
|
|
165 }
|
|
|
166 $previous[7]='NONE' if $previous[7] eq '.';
|
|
|
167 $previous[$_]=sprintf("%.0f", $previous[$_]) for (6,8..11);
|
|
|
168 print OUT join("\t", @previous[0..6,8..11,7]), "\n";
|
|
|
169 @Temp=@previous if eof;
|
|
|
170 @previous=@buffer;
|
|
|
171 }
|
|
|
172 $previous[7]='NONE' if $previous[7] eq '.';
|
|
|
173 $previous[$_]=sprintf("%.0f", $previous[$_]) for (6,8..11);
|
|
|
174 print OUT join("\t", @previous[0..6,8..11,7]), "\n" if(join('_', @Temp[1..2]) ne join('_', @previous[1..2]));
|
|
|
175 close IN;
|
|
|
176 close OUT;
|
|
|
177 }
|
|
|
178 foreach (1..22, 'X', 'Y', 'M'){
|
|
|
179 push @fnames, "${outdir}/${fname}_chr${_}.Temp-02" if exists $chr{"chr$_"};
|
|
|
180 }
|
|
|
181 system join(' ', 'cat', @fnames, '>', "${outdir}/${fname}.Temp.2R") and die $!;
|
|
|
182 `${rbin} --vanilla --slave --args ${outdir}/${fname}.Temp.2R < ${rdep}/samvcf_data_parser.R` and die $!;
|
|
|
183 system "rm ${outdir}/${fname}*Temp* $outfile; ln -s ${outdir}/${fname}.var $outfile" and die $!; |
