Mercurial > repos > rdaveau > gfap

comparison gfapts/gfap_r1.0_samvcf_data_parser.pl @ 0:f753b30013e6 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author rdaveau
date Fri, 29 Jun 2012 10:20:55 -0400
parents
children 028f435b6cfb
comparison
equal deleted inserted replaced
-1:000000000000 0:f753b30013e6
1 #!/usr/bin/perl
2
3 use strict;
4 use lib 'inc/perlmod';
5 use ngsutil qw[ :DEFAULT &explode_varcall ];
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
16 my $annovar_dir = 'inc/annovar';
17 my $rdep = 'inc/R';
18
19 my($varfile, $outdir, $outfile, $i, @DP4, @buffer, @Temp, @previous, @fnames, %opts, %chr);
20
21 GetOptions(\%opts, "varfile=s", "outdir=s", "outfile=s");
22 $varfile = $opts{varfile};
23 $outdir = $opts{outdir};
24 $outfile = $opts{outfile};
25
26 my $fname = readlink($varfile) || $varfile;
27 $fname = basename($fname);
28
29 my %fh=(
30 'chr1' => *chr1, 'chr2' => *chr2, 'chr3' => *chr3, 'chr4' => *chr4, 'chr5' => *chr5,
31 'chr6' => *chr6, 'chr7' => *chr7, 'chr8' => *chr8, 'chr9' => *chr9, 'chr10' => *chr10,
32 'chr11' => *chr11, 'chr12' => *chr12, 'chr13' => *chr13, 'chr14' => *chr14, 'chr15' => *chr15,
33 'chr16' => *chr16, 'chr17' => *chr17, 'chr18' => *chr18, 'chr19' => *chr19, 'chr20' => *chr20,
34 'chr21' => *chr21, 'chr22' => *chr22, 'chrX' => *chrX, 'chrY' => *chrY, 'chrM' => *chrM
35 );
36
37 `${annovar_dir}/convert2annovar.pl -format vcf4 $varfile -includeinfo > ${outdir}/${fname}_Temp-00 2> /dev/null` and die $!;
38
39 open($fh{$_}, ">${outdir}/${fname}_${_}.Temp-00") or die $! foreach keys %fh;
40 open IN, "<${outdir}/${fname}_Temp-00" or die $!;
41 while(<IN>){
42 /^(\S+)\s+(?:\S+\s+){2}(\S+)\s+(\S+)/;
43 next if !exists $fh{$1};
44 if(min(length($2), length($3))!=1){
45 chomp;
46 @buffer=split /\s+/, $_;
47 @Temp=explode_varcall(@buffer[1,3..4]);
48 for($i=0; $i<$#{$Temp[0]}; $i+=2){
49 print{ $fh{$buffer[0]} } join("\t", $buffer[0], @{$Temp[0]}[$i..$i+1], @{$Temp[1]}[$i..$i+1], @buffer[6..$#buffer]), "\n";
50 }
51 next;
52 }
53 print{ $fh{$1} } $_;
54 $chr{$1}++;
55 }
56 close IN;
57 foreach (keys %fh){
58 close($fh{$_});
59 next if !exists $chr{$_};
60 `sort -k2,2n -k3,3n ${outdir}/${fname}_${_}.Temp-00 > ${outdir}/${fname}_${_}.Temp-01` and die $!;
61 open IN, "<${outdir}/${fname}_${_}.Temp-01" or die $!;
62 open OUT, ">${outdir}/${fname}_${_}.Temp-02" or die $!;
63 $_=readline(IN);
64 /^((?:\S+\s+){7})(?:\S+\s+){8}(\S+\s+\S+)/;
65 @buffer=split /\s+/, $1.$2;
66 ($_=pop(@buffer))=~s/.+DP4=([^;]+).+/$1/;
67 @DP4=split /,/, $_;
68 push @buffer, @DP4;
69 @previous=@buffer;
70 MAINLOOP: while(<IN>){
71 /^((?:\S+\s+){7})(?:\S+\s+){8}(\S+\s+\S+)/;
72 @buffer=split /\s+/, $1.$2;
73 ($_=pop(@buffer))=~s/.+DP4=([^;]+).+/$1/;
74 @DP4=split /,/, $_;
75 push @buffer, @DP4;
76 while(($previous[0] eq $buffer[0]) && ($buffer[2]==$previous[2]+1) && (join('', @previous[3..4]) !~ /-/) && (join('', @buffer[3..4]) !~ /-/)){
77 $previous[2]=$buffer[2];
78 $previous[$_].=$buffer[$_] for 3..4;
79 $previous[5]='unk' if $previous[5] ne $buffer[5];
80 $previous[7]='SKIP' if $previous[7] ne $buffer[7];
81 for (6,8..11){
82 $previous[$_]+=$buffer[$_];
83 $previous[$_]/=2;
84 }
85 next MAINLOOP;
86 }
87 $previous[7]='NONE' if $previous[7] eq '.';
88 $previous[$_]=sprintf("%.0f", $previous[$_]) for (6,8..11);
89 print OUT join("\t", @previous[0..6,8..11,7]), "\n";
90 @Temp=@previous if eof;
91 @previous=@buffer;
92 }
93 $previous[7]='NONE' if $previous[7] eq '.';
94 $previous[$_]=sprintf("%.0f", $previous[$_]) for (6,8..11);
95 print OUT join("\t", @previous[0..6,8..11,7]), "\n" if(join('_', @Temp[1..2]) ne join('_', @previous[1..2]));
96 close IN;
97 close OUT;
98 }
99 foreach (1..22, 'X', 'Y', 'M'){
100 push @fnames, "${outdir}/${fname}_chr${_}.Temp-02" if exists $chr{"chr$_"};
101 }
102 system join(' ', 'cat', @fnames, '>', "${outdir}/${fname}.Temp.2R") and die $!;
103 `${rbin} --vanilla --slave --args ${outdir}/${fname}.Temp.2R < ${rdep}/samvcf_data_parser.R` and die $!;
104 system "rm ${outdir}/${fname}*Temp* $outfile; ln -s ${outdir}/${fname}.var $outfile" and die $!;