Mercurial > repos > rdaveau > gfap
diff 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 |
line wrap: on
line diff
--- a/gfapts/gfap_r1.0_known_var_finder.pl Fri Jun 29 10:20:55 2012 -0400 +++ b/gfapts/gfap_r1.0_known_var_finder.pl Fri Aug 03 05:50:41 2012 -0400 @@ -1,12 +1,91 @@ #!/usr/bin/perl use strict; -use lib 'inc/perlmod'; -use ngsutil qw[ :DEFAULT &varscan ]; +#use lib 'inc/perlmod'; +#use ngsutil qw[ :DEFAULT &varscan ]; use warnings FATAL => qw[ numeric uninitialized ]; use File::Basename; use Getopt::Long; +#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +# TEMP include ngsutil.pm +sub explode_varcall{ + my $N=0; + $_=shift @_ foreach my($POS, $REF, $ALT); + $_=$POS foreach my($START, $END); + my(@length, @range, @idx, @VAR, @POS); + @{$_}=() foreach (\@length, \@range, \@idx, \@VAR, \@POS); + push @length, length($_) foreach ($REF, $ALT); + @range=sort{ $a<=>$b } @length; + if($range[0]==1){ + if($range[1]!=1){ + foreach ($REF, $ALT){ + $_=substr($_, 1); + $_=~s/^$/-/; + } + if($length[0]!=1){ + $END+=$length[0]-1; + $START++; + } + } + push @POS, $START, $END; + push @VAR, $REF, $ALT; + }else{ + my @N=(); + undef $_ foreach my ($i, $VAR); + $_-=2 foreach (@length, @range); + $_++ foreach ($START, $END); + $_=substr($_, 1) foreach ($REF, $ALT); + my $indel='-' x ($range[1]-$range[0]); + $VAR.=($_>$range[0])? + ('-'):((substr($REF, $_, 1) ne substr($ALT, $_, 1))? + 0:1) for 0 .. $range[1]; + $N++ while $VAR =~ /0/g; + if($length[0]<$length[1]){ + @VAR=($VAR); + @N=($N); + $N=0; + undef($VAR); + $VAR.=($_>$range[0])? + ('-'):((substr($REF, $length[0]-$_, 1) ne substr($ALT, $length[1]-$_, 1))? + 0:1) for reverse 0 .. $range[1]; + $N++ while $VAR =~ /0/g; + if($N>=$N[0]){ $N=shift(@N); $VAR=shift(@VAR); } + else{ $REF=$indel . $REF; } + }else{ $ALT.=$indel; } + foreach (qw[ 0 \- ]){ + push @idx, [ $-[0], $+[0]-$-[0] ] while ($VAR =~ /$_+/g); + } + @{$_}=() foreach (\@VAR, \@POS); + foreach my $k (@idx){ + push @VAR, substr($_, ${$k}[0], ${$k}[1]) || '-' foreach ($REF, $ALT); + push @POS, ${$k}[0], sum(@{$k})-1; + } + $_+=$START foreach @POS; + $_=~s/\-+/\-/ foreach @VAR; + for($i=0; $i<$#POS; $i+=2){ $POS[$i+1]=$POS[$i] if $VAR[$i] eq '-'; } + } + return(\@POS, \@VAR); + } + +sub varscan{ + $_=shift @_ foreach my($kname, $fpath, $href); + my($k, @buffer); + open IN, "<$fpath" or die $!; + while(<IN>){ + next if /^#/; + chomp; + @buffer=split /\s+/, $_; + next if !exists $$href{($k=join(':', @buffer[0..2]))}; + next if $$href{$k}->{ref} !~ $buffer[3]; + next if $$href{$k}->{alt} !~ $buffer[4]; + splice(@buffer, 0, 5); + $$href{$k}->{$kname}=join(':', @buffer); + } + close IN; + } +#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + my($varfile, $buildver, $outdir, $dir_1000g, $dir_dbsnp, $dir_cosmic, $release_1000g, $release_dbsnp, $release_cosmic, $outfile, $k, @buffer, @varlist, %opts, %varlist); 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");