|
0
|
1 #!/usr/bin/env perl
|
|
|
2
|
|
|
3 use strict;
|
|
|
4 use warnings;
|
|
|
5
|
|
|
6 # This script adds sample names to the last column of an HGVS file if the coverage for that position was insufficient to make the call the other samples had.
|
|
|
7 @ARGV == 8 or die "Usage: $0 <hgvs_in.txt> <hgvs_out.txt> <sample name1> <sample1.bam> <homo FN coverage cutoff> <het coverage FN cutoff> <min map quality> <min base quality>\n";
|
|
|
8 my $hgvs_in_file = shift @ARGV;
|
|
|
9 my $hgvs_out_file = shift @ARGV;
|
|
|
10 my $sample_name = shift @ARGV;
|
|
|
11 my $bam_file = shift @ARGV;
|
|
|
12 my $homo_fn = shift @ARGV;
|
|
|
13 my $het_fn = shift @ARGV;
|
|
|
14 my $min_maq = shift @ARGV; # 30 (a.k.a. p < 0.001)
|
|
|
15 my $min_baq = shift @ARGV; # 20 (a.k.a. p < 0.01)
|
|
|
16
|
|
|
17 open(HGVS, $hgvs_in_file)
|
|
|
18 or die "Cannot open $hgvs_in_file for reading: $!\n";
|
|
|
19 my $header_line = <HGVS>;
|
|
|
20 chomp $header_line;
|
|
|
21 my @headers = split /\t/, $header_line;
|
|
|
22 my ($chr_column, $pos_column, $ref_column, $alt_column, $zygosity_column, $alt_cnt_column, $tot_cnt_column, $srcs_column);
|
|
|
23 for(my $i = 0; $i <= $#headers; $i++){
|
|
|
24 if($headers[$i] eq "Chr"){
|
|
|
25 $chr_column = $i;
|
|
|
26 }
|
|
|
27 elsif($headers[$i] eq "DNA From"){
|
|
|
28 $pos_column = $i;
|
|
|
29 }
|
|
|
30 elsif($headers[$i] eq "Ref base"){
|
|
|
31 $ref_column = $i;
|
|
|
32 }
|
|
|
33 elsif($headers[$i] eq "Obs base"){
|
|
|
34 $alt_column = $i;
|
|
|
35 }
|
|
|
36 elsif($headers[$i] eq "Zygosity"){
|
|
|
37 $zygosity_column = $i;
|
|
|
38 }
|
|
|
39 elsif($headers[$i] eq "Variant Reads"){
|
|
|
40 $alt_cnt_column = $i;
|
|
|
41 }
|
|
|
42 elsif($headers[$i] eq "Total Reads"){
|
|
|
43 $tot_cnt_column = $i;
|
|
|
44 }
|
|
|
45 elsif($headers[$i] eq "Sources"){
|
|
|
46 $srcs_column = $i;
|
|
|
47 }
|
|
|
48 }
|
|
|
49 die "Could not find header Chr in $hgvs_in_file, aborting." if not defined $chr_column;
|
|
|
50 die "Could not find header 'DNA From' in $hgvs_in_file, aborting." if not defined $pos_column;
|
|
|
51 die "Could not find header 'Ref base' in $hgvs_in_file, aborting." if not defined $ref_column;
|
|
|
52 die "Could not find header 'Obs base' in $hgvs_in_file, aborting." if not defined $alt_column;
|
|
|
53 die "Could not find header 'Variant Reads' in $hgvs_in_file, aborting." if not defined $alt_cnt_column;
|
|
|
54 die "Could not find header 'Total Reads' in $hgvs_in_file, aborting." if not defined $tot_cnt_column;
|
|
|
55 die "Could not find header Zygosity in $hgvs_in_file, aborting." if not defined $zygosity_column;
|
|
|
56 die "Could not find header Sources in $hgvs_in_file, aborting." if not defined $srcs_column;
|
|
|
57 open(OUT, ">$hgvs_out_file")
|
|
|
58 or die "Cannot open $hgvs_out_file for writing: $!\n";
|
|
|
59 print OUT $header_line, "\n"; # header line
|
|
|
60 my ($last_mpileup, $last_key);
|
|
|
61 while(<HGVS>){
|
|
|
62 chomp;
|
|
|
63 my @F = split /\t/, $_;
|
|
|
64 if(index($F[$#F], $sample_name) != -1){
|
|
|
65 print OUT "$_\n"; # keep line as-is
|
|
|
66 next;
|
|
|
67 }
|
|
|
68 else{
|
|
|
69 # Not in a call right now
|
|
|
70 # If there are reads supporting the call, add the sample name to the list (last column)
|
|
|
71 # If there are no reads supporting the call, but coverage is below threshold, put a ~ in front of the sample name
|
|
|
72 my $chr = $F[$chr_column];
|
|
|
73 my $pos = $F[$pos_column];
|
|
|
74 my $ref = $F[$ref_column];
|
|
|
75 my $alt = $F[$alt_column];
|
|
|
76
|
|
|
77 # some tools report SNPs with leading or trailing matched ref seqs, so trim as necessary
|
|
|
78 if(length($ref) > 1 and length($ref) == length($alt)){
|
|
|
79 while(length($ref) and substr($ref, 0, 1) eq substr($alt, 0, 1)){
|
|
|
80 substr($ref, 0, 1) = "";
|
|
|
81 substr($alt, 0, 1) = "";
|
|
|
82 $pos++;
|
|
|
83 }
|
|
|
84 while(length($ref) and substr($ref, -1) eq substr($alt, -1)){
|
|
|
85 substr($ref, -1) = "";
|
|
|
86 substr($alt, -1) = "";
|
|
|
87 }
|
|
|
88 }
|
|
|
89
|
|
|
90 # (todo: what about compound het?)
|
|
|
91 my $zygosity = $F[$zygosity_column];
|
|
|
92 # see if the coverage is below a called zygosity threshold
|
|
|
93 my $end = $pos;
|
|
|
94 # comment out next bit if we assume that reads have all indels left aligned
|
|
|
95 #if($alt ne "NA" and $ref ne "NA" and length($alt) > length($ref)){
|
|
|
96 # $end += length($alt) - length($ref);
|
|
|
97 #}
|
|
|
98 my $mpileup;
|
|
|
99 if(defined $last_key and $last_key eq "$chr:$pos:$ref:$alt"){ # used cached response if this is a repeated input line
|
|
|
100 $mpileup = $last_mpileup;
|
|
|
101 }
|
|
|
102 else{
|
|
|
103 $mpileup = `samtools mpileup -q $min_maq -Q $min_baq -r "$chr:$pos-$end" $bam_file 2>&1`;
|
|
|
104 $last_mpileup = $mpileup;
|
|
|
105 $last_key = "$chr:$pos:$ref:$alt";
|
|
|
106 }
|
|
|
107
|
|
|
108 $ref = "" if $ref eq "NA";
|
|
|
109 $alt = "" if $alt eq "NA";
|
|
|
110 # lines look something like chr10 74100774 N 164 G$ggGgg$GgGgGgggggG$ggggggggGggGGGgg...
|
|
|
111 # deletions are -1N, while insertions are +1g
|
|
|
112 for my $line (split /\n/, $mpileup){
|
|
|
113 my @M = split /\t/, $line;
|
|
|
114 next if $M[0] ne $chr; # header or message stuff
|
|
|
115 my $depth = $M[3];
|
|
|
116 if($depth == 0){ # No reads at all
|
|
|
117 $F[$srcs_column] .= "; ~$sample_name"; # add this sample to the last column
|
|
|
118 $F[$zygosity_column] .= "; low"; # zygosity unknown
|
|
|
119 $F[$alt_cnt_column] .= "; 0"; # alt count
|
|
|
120 $F[$tot_cnt_column] .= "; 0"; # tot count
|
|
|
121 #print OUT join("\t", @F), "\n";
|
|
|
122 last;
|
|
|
123 }
|
|
|
124 my %base_count;
|
|
|
125 while($M[4] =~ /[ACGT]|\^.|(-\d+|\+\d+)/gi){
|
|
|
126 my $call = uc($&); # upper case version of letter
|
|
|
127 if($1){ # call is -1n, or +3gcc, etc.
|
|
|
128 if(substr($1, 0, 1) eq "+"){ # insertion
|
|
|
129 my $offset = substr($1, 1);
|
|
|
130 $call .= uc(substr($M[4], pos($M[4]), $offset)); # append the letters (not captures in the regex above) to the number
|
|
|
131 pos($M[4]) += $offset; # skip the global search ahead of the letters just consumed
|
|
|
132 }
|
|
|
133 }
|
|
|
134 elsif($call eq '$' or length($call) > 1){ # the weird ^Y kind of call
|
|
|
135 next; # ignore
|
|
|
136 }
|
|
|
137 $base_count{$call}++;
|
|
|
138 }
|
|
|
139 my $base_count = 0;
|
|
|
140 if(length($ref) == 1 and length($alt) == 1){ #SNP
|
|
|
141 if(exists $base_count{$alt}){
|
|
|
142 $base_count = $base_count{$alt};
|
|
|
143 }
|
|
|
144 }
|
|
|
145 elsif(length($ref) > length($alt)){ #del
|
|
|
146 my $del_length = length($alt)-length($ref); # e.g. -2
|
|
|
147 if(exists $base_count{$del_length}){
|
|
|
148 $base_count = $base_count{$del_length};
|
|
|
149 }
|
|
|
150 }
|
|
|
151 elsif(length($ref) < length($alt)){ #ins
|
|
|
152 my $ins_call = "+".(length($alt)-length($ref)).substr($alt, length($ref)); # e.g. +3AGT
|
|
|
153 if(exists $base_count{$ins_call}){
|
|
|
154 $base_count = $base_count{$ins_call};
|
|
|
155 }
|
|
|
156 }
|
|
|
157 else{
|
|
|
158 # MNP variant?
|
|
|
159 #warn "Cannot check $ref -> $alt: not a SNP insertion or deletion\n";
|
|
|
160 next;
|
|
|
161 }
|
|
|
162 if(defined $base_count and $base_count and $depth/($base_count+1) < $het_fn){ # at least the min het prop (rounded up) of bases agreeing with the (missed in this sample) alt call
|
|
|
163 $F[$srcs_column] .= "; +$sample_name"; # add this sample to the last column
|
|
|
164 $F[$zygosity_column] .= $base_count/$depth < 0.78 ? "; heterozygote" : ($depth < $het_fn ? "; low" : "; homozygote"); # zygosity
|
|
|
165 $F[$alt_cnt_column] .= "; ".$base_count; # alt count
|
|
|
166 $F[$tot_cnt_column] .= "; $depth"; # tot count
|
|
|
167 }
|
|
|
168 elsif($zygosity =~ /het/ and $depth < $het_fn){
|
|
|
169 $F[$srcs_column].= "; ~$sample_name";
|
|
|
170 $F[$zygosity_column] .= $depth < $homo_fn ? "; low" : "; heterozygote"; # zygosity
|
|
|
171 $F[$alt_cnt_column] .= "; $base_count";
|
|
|
172 $F[$tot_cnt_column] .= "; $depth"; # tot count
|
|
|
173 }
|
|
|
174 elsif($depth < $homo_fn){
|
|
|
175 $F[$srcs_column].= "; ~$sample_name";
|
|
|
176 $F[$zygosity_column] .= "; low"; # zygosity
|
|
|
177 $F[$alt_cnt_column] .= "; $base_count";
|
|
|
178 $F[$tot_cnt_column] .= "; $depth"; # tot count
|
|
|
179 }
|
|
|
180 }
|
|
|
181 print OUT join("\t", @F), "\n";
|
|
|
182 }
|
|
|
183 }
|
|
|
184 close(HGVS);
|
|
|
185 close(OUT);
|
