Mercurial > repos > rdaveau > gfap

comparison gfapts/gfap_r1.0_samvcf_data_parser.pl @ 1:028f435b6cfb draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author rdaveau
date Fri, 03 Aug 2012 05:50:41 -0400
parents f753b30013e6
children
comparison
equal deleted inserted replaced
0:f753b30013e6 1:028f435b6cfb
1 #!/usr/bin/perl 1 #!/usr/bin/perl
2 2
3 use strict; 3 use strict;
4 use lib 'inc/perlmod'; 4 # use lib 'inc/perlmod';
5 use ngsutil qw[ :DEFAULT &explode_varcall ]; 5 # use ngsutil qw[ :DEFAULT &explode_varcall ];
6 use warnings FATAL => qw[ numeric uninitialized ]; 6 use warnings FATAL => qw[ numeric uninitialized ];
7 use List::Util qw[ sum min max ]; 7 use List::Util qw[ sum min max ];
8 use File::Basename; 8 use File::Basename;
9 use Getopt::Long; 9 use Getopt::Long;
10 10
11 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 11 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
12 # PATH TO YOUR R-bin DIRECTORY 12 # PATH TO YOUR R-bin DIRECTORY
13 my $rbin = '/usr/bin/R'; 13 my $rbin = '/usr/bin/R';
14 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
15
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 }
14 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 93 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
15 94
16 my $annovar_dir = 'inc/annovar'; 95 my $annovar_dir = 'inc/annovar';
17 my $rdep = 'inc/R'; 96 my $rdep = 'inc/R';
18 97