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

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 &varscan ]; 5 #use ngsutil qw[ :DEFAULT &varscan ];
6 use warnings FATAL => qw[ numeric uninitialized ]; 6 use warnings FATAL => qw[ numeric uninitialized ];
7 use File::Basename; 7 use File::Basename;
8 use Getopt::Long; 8 use Getopt::Long;
9
10 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
11 # TEMP include ngsutil.pm
12 sub explode_varcall{
13 my $N=0;
14 $_=shift @_ foreach my($POS, $REF, $ALT);
15 $_=$POS foreach my($START, $END);
16 my(@length, @range, @idx, @VAR, @POS);
17 @{$_}=() foreach (\@length, \@range, \@idx, \@VAR, \@POS);
18 push @length, length($_) foreach ($REF, $ALT);
19 @range=sort{ $a<=>$b } @length;
20 if($range[0]==1){
21 if($range[1]!=1){
22 foreach ($REF, $ALT){
23 $_=substr($_, 1);
24 $_=~s/^$/-/;
25 }
26 if($length[0]!=1){
27 $END+=$length[0]-1;
28 $START++;
29 }
30 }
31 push @POS, $START, $END;
32 push @VAR, $REF, $ALT;
33 }else{
34 my @N=();
35 undef $_ foreach my ($i, $VAR);
36 $_-=2 foreach (@length, @range);
37 $_++ foreach ($START, $END);
38 $_=substr($_, 1) foreach ($REF, $ALT);
39 my $indel='-' x ($range[1]-$range[0]);
40 $VAR.=($_>$range[0])?
41 ('-'):((substr($REF, $_, 1) ne substr($ALT, $_, 1))?
42 0:1) for 0 .. $range[1];
43 $N++ while $VAR =~ /0/g;
44 if($length[0]<$length[1]){
45 @VAR=($VAR);
46 @N=($N);
47 $N=0;
48 undef($VAR);
49 $VAR.=($_>$range[0])?
50 ('-'):((substr($REF, $length[0]-$_, 1) ne substr($ALT, $length[1]-$_, 1))?
51 0:1) for reverse 0 .. $range[1];
52 $N++ while $VAR =~ /0/g;
53 if($N>=$N[0]){ $N=shift(@N); $VAR=shift(@VAR); }
54 else{ $REF=$indel . $REF; }
55 }else{ $ALT.=$indel; }
56 foreach (qw[ 0 \- ]){
57 push @idx, [ $-[0], $+[0]-$-[0] ] while ($VAR =~ /$_+/g);
58 }
59 @{$_}=() foreach (\@VAR, \@POS);
60 foreach my $k (@idx){
61 push @VAR, substr($_, ${$k}[0], ${$k}[1]) || '-' foreach ($REF, $ALT);
62 push @POS, ${$k}[0], sum(@{$k})-1;
63 }
64 $_+=$START foreach @POS;
65 $_=~s/\-+/\-/ foreach @VAR;
66 for($i=0; $i<$#POS; $i+=2){ $POS[$i+1]=$POS[$i] if $VAR[$i] eq '-'; }
67 }
68 return(\@POS, \@VAR);
69 }
70
71 sub varscan{
72 $_=shift @_ foreach my($kname, $fpath, $href);
73 my($k, @buffer);
74 open IN, "<$fpath" or die $!;
75 while(<IN>){
76 next if /^#/;
77 chomp;
78 @buffer=split /\s+/, $_;
79 next if !exists $$href{($k=join(':', @buffer[0..2]))};
80 next if $$href{$k}->{ref} !~ $buffer[3];
81 next if $$href{$k}->{alt} !~ $buffer[4];
82 splice(@buffer, 0, 5);
83 $$href{$k}->{$kname}=join(':', @buffer);
84 }
85 close IN;
86 }
87 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
9 88
10 my($varfile, $buildver, $outdir, $dir_1000g, $dir_dbsnp, $dir_cosmic, $release_1000g, $release_dbsnp, $release_cosmic, $outfile, $k, @buffer, @varlist, %opts, %varlist); 89 my($varfile, $buildver, $outdir, $dir_1000g, $dir_dbsnp, $dir_cosmic, $release_1000g, $release_dbsnp, $release_cosmic, $outfile, $k, @buffer, @varlist, %opts, %varlist);
11 90
12 GetOptions(\%opts, "varfile=s", "buildver=s", "outdir=s", "dir_1000g=s", "dir_dbsnp=s", "dir_cosmic=s", "release_1000g=s", "release_dbsnp=s", "release_cosmic=s", "outfile=s"); 91 GetOptions(\%opts, "varfile=s", "buildver=s", "outdir=s", "dir_1000g=s", "dir_dbsnp=s", "dir_cosmic=s", "release_1000g=s", "release_dbsnp=s", "release_cosmic=s", "outfile=s");
13 $varfile = $opts{varfile}; 92 $varfile = $opts{varfile};