|
0
|
1 #!/usr/bin/perl
|
|
|
2
|
|
|
3 use strict;
|
|
|
4 use lib 'inc/perlmod';
|
|
|
5 use ngsutil qw[ :DEFAULT &varscan ];
|
|
|
6 use warnings FATAL => qw[ numeric uninitialized ];
|
|
|
7 use File::Basename;
|
|
|
8 use Getopt::Long;
|
|
|
9
|
|
|
10 my($varfile, $buildver, $outdir, $dir_1000g, $dir_dbsnp, $dir_cosmic, $release_1000g, $release_dbsnp, $release_cosmic, $outfile, $k, @buffer, @varlist, %opts, %varlist);
|
|
|
11
|
|
|
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");
|
|
|
13 $varfile = $opts{varfile};
|
|
|
14 $buildver = $opts{buildver};
|
|
|
15 $outdir = $opts{outdir};
|
|
|
16 $dir_1000g = $opts{dir_1000g};
|
|
|
17 $dir_dbsnp = $opts{dir_dbsnp};
|
|
|
18 $dir_cosmic = $opts{dir_cosmic};
|
|
|
19 $release_1000g = $opts{release_1000g};
|
|
|
20 $release_dbsnp = $opts{release_dbsnp};
|
|
|
21 $release_cosmic = $opts{release_cosmic};
|
|
|
22 $outfile = $opts{outfile};
|
|
|
23
|
|
|
24 my $fname = readlink($varfile) || $varfile;
|
|
|
25 $fname = basename($fname);
|
|
|
26
|
|
|
27 my %k=(
|
|
|
28 '1000g' => {
|
|
|
29 'dir' => $dir_1000g, 'release' => $release_1000g, 'value' => join(':', ('0.00000')x5), 'header' => join(':', 'AF_ALL', 'AF_AFR', 'AF_AMR', 'AF_ASN', 'AF_EUR')
|
|
|
30 }, 'dbsnp' => {
|
|
|
31 'dir' => $dir_dbsnp, 'release' => $release_dbsnp, 'value' => join(':', ('na')x2), 'header' => join(':', 'rs', 'dbsnp')
|
|
|
32 }, 'cosmic_var' => {
|
|
|
33 'dir' => $dir_cosmic, 'release' => $release_cosmic, 'value' => join(':', '0.00000', 'na'), 'header' => join(':', 'AF_COS', 'cid')
|
|
|
34 }
|
|
|
35 );
|
|
|
36
|
|
|
37 my %legend=(
|
|
|
38 'chr' => 'chromosome identifier',
|
|
|
39 'start' => "${buildver} 1-based start position",
|
|
|
40 'end' => "${buildver} 1-based end position",
|
|
|
41 'ref' => 'reference allele',
|
|
|
42 'alt' => 'alternate allele',
|
|
|
43 'QC' => 'Phred-scaled call quality',
|
|
|
44 'NRF' => '#reads consistent w/ the reference allele on the F-strand',
|
|
|
45 'NRR' => '#reads consistent w/ the reference allele on the R-strand',
|
|
|
46 'NAF' => '#reads consistent w/ the alternate allele on the F-strand',
|
|
|
47 'NAR' => '#reads consistent w/ the alternate allele on the R-strand',
|
|
|
48 'DP' => 'total #reads in call ie. NRF+NRR+NAF+NAR',
|
|
|
49 'AD' => 'total #reads consistent w/ the alternate allele ie. NAF+NAR',
|
|
|
50 'AF' => 'alternate allele ratio ie. AD/DP',
|
|
|
51 'VCF.FILTER' => 'FILTER field from the input vcf file',
|
|
|
52 'DPT.FILTER' => 'check for heterogeneous depth in substituted blocks',
|
|
|
53 'VAR.FILTER' => 'GFAP default FILTER to discriminate between TP and FP variants',
|
|
|
54 'P.str' => 'NRF+NAF vs. NRR+NAR binomial test P-value ie. total strand bias',
|
|
|
55 'P.ref' => 'NRF vs. NRR binomial test P-value ie. reference allele strand bias',
|
|
|
56 'P.alt' => 'NAF vs. NAR binomial test P-value ie. alternate allele strand bias',
|
|
|
57 'AF_ALL' => "global AF in ${release_1000g} 1000g data",
|
|
|
58 'AF_AFR' => "AF in AFR ${release_1000g} 1000g data",
|
|
|
59 'AF_AMR' => "AF in AMR ${release_1000g} 1000g data",
|
|
|
60 'AF_ASN' => "AF in ASN ${release_1000g} 1000g data",
|
|
|
61 'AF_EUR' => "AF in EUR ${release_1000g} 1000g data",
|
|
|
62 'AF_COS' => "AF in ${release_cosmic} cosmic data",
|
|
|
63 'rs' => "dbsnp rs identifier(s) from ${release_dbsnp} release",
|
|
|
64 'dbsnp' => "dbsnp build version(s) from ${release_dbsnp} release",
|
|
|
65 'cid' => "cosmic mutation identifier from ${release_cosmic} release"
|
|
|
66 );
|
|
|
67 my @header=('chr', 'start', 'end', 'ref', 'alt', 'DPT.FILTER', 'QC', 'NRF', 'NRR', 'NAF', 'NAR', 'VCF.FILTER', 'P.str', 'P.ref', 'P.alt', 'DP', 'AD', 'AF', 'VAR.FILTER');
|
|
|
68 my @k=qw[ 1000g dbsnp cosmic_var ];
|
|
|
69
|
|
|
70 open IN, "<$varfile" or die $!;
|
|
|
71 while(<IN>){
|
|
|
72 chomp;
|
|
|
73 @buffer=split /\s+/, $_;
|
|
|
74 $buffer[0]=~s/^chr(.+)$/$1/;
|
|
|
75 push @varlist, ($k=join(':', @buffer[0..2]));
|
|
|
76 shift(@buffer) for 0..2;
|
|
|
77 $varlist{$k}->{$_}=shift(@buffer) foreach qw[ ref alt ];
|
|
|
78 $varlist{$k}->{cov}=join(':', (($buffer[0] eq 'unk')?'SKIP':'PASS'), @buffer[1..$#buffer]);
|
|
|
79 }
|
|
|
80 close IN;
|
|
|
81
|
|
|
82 foreach $k (@k){
|
|
|
83 push @header, split(/:/, $k{$k}->{header});
|
|
|
84 varscan($k, $k{$k}->{file}, \%varlist);
|
|
|
85 }
|
|
|
86
|
|
|
87 my @idx=(0..4,7..10,15..17,6,12..14,11,5,18..23,26..27,24..25);
|
|
|
88 open OUT, ">${outdir}/${fname}.dbi" or die $!;
|
|
|
89 print OUT '#', join(' = ', $_, $legend{$_}), "\n" foreach @header[@idx];
|
|
|
90 print OUT '#', join("\t", @header[@idx]), "\n";
|
|
|
91 foreach $k (@varlist){
|
|
|
92 @buffer=(split(/:/, 'chr'.$k), $varlist{$k}->{ref}, $varlist{$k}->{alt});
|
|
|
93 push @buffer, split(/:/, ($varlist{$k}->{$_} || $k{$_}->{value})) foreach ('cov', @k);
|
|
|
94 print OUT join("\t", @buffer[@idx]), "\n";
|
|
|
95 }
|
|
|
96 close OUT;
|
|
|
97
|
|
|
98 system "rm $outfile; ln -s ${outdir}/${fname}.dbi $outfile" and die $!; |
