0
+ − 1 #!/usr/bin/env perl
+ − 2
+ − 3 BEGIN{
+ − 4 my $prog_dir = `dirname $0`;
+ − 5 chomp $prog_dir;
+ − 6 push @INC, $prog_dir; # so DisjointSets.pm can be found no matter the working directory
+ − 7 }
+ − 8
+ − 9 use DisjointSets; # homebrew module
+ − 10 use Bio::DB::Sam; # for FastA reference pulls
+ − 11 use Bio::SeqUtils;
+ − 12 use Bio::Tools::CodonTable;
+ − 13 use Statistics::Zed;
+ − 14 use Getopt::Long;
+ − 15 use Set::IntervalTree;
+ − 16 use strict;
+ − 17 use warnings;
+ − 18 use vars qw($min_prop $zed $codonTable $default_transl_table %transl_except %internal_prop %dbsnp_info %chr2variant_locs %chr2dbsnp_vcf_lines %chr2internal_vcf_lines %chr2caveats %chr2phase @snvs $fasta_index $max_args $quiet);
+ − 19
+ − 20 if(@ARGV == 1 and $ARGV[0] eq "-v"){
+ − 21 print "Version 1.0\n";
+ − 22 exit;
+ − 23 }
+ − 24
+ − 25 #$max_args = `getconf ARG_MAX`; # largest number of args you can send to a system command (enviroment included, see limits.h)
+ − 26 #chomp $max_args;
+ − 27 $max_args = 4096; # if not defined $max_args or $max_args < 1; # the minimum since System V
+ − 28 $max_args -= 50;
+ − 29
+ − 30 # find out if a variant appears in the user provided data
+ − 31 sub internal_prop($$$$){
+ − 32 my ($chr,$pos,$ref,$variant) = @_;
+ − 33
+ − 34 my $key = "$chr:$pos:$ref:$variant";
+ − 35 if(exists $internal_prop{$key}){
+ − 36 return $internal_prop{$key};
+ − 37 }
+ − 38
+ − 39 #print STDERR "Checking if internal_prop for $key exists: ";
+ − 40 if(exists $chr2internal_vcf_lines{$chr}->{$pos}){
+ − 41 for(@{$chr2internal_vcf_lines{$chr}->{$pos}}){
+ − 42 my @fields = split /\t/, $_;
+ − 43 if($pos == $fields[1] and length($fields[3]) == length($ref) and $fields[4] eq $variant){
+ − 44 #print STDERR "yes\n";
+ − 45 if(/MAF=(\d\.\d+)/){
+ − 46 $internal_prop{$key} = $1; # change from percent to proportion
+ − 47 return $1;
+ − 48 }
+ − 49 }
+ − 50 }
+ − 51 }
+ − 52 else{
+ − 53 #print STDERR "no\n";
+ − 54 }
+ − 55
+ − 56 $internal_prop{$key} = "NA";
+ − 57 return "NA";
+ − 58 }
+ − 59
+ − 60 # find out if a variant appears in the NCBI's dbSNP
+ − 61 sub dbsnp_info($$$$){
+ − 62 my ($chr,$pos,$ref,$variant) = @_;
+ − 63
+ − 64 my $key = "$chr:$pos:$ref:$variant";
+ − 65 if(exists $dbsnp_info{$key}){
+ − 66 return @{$dbsnp_info{$key}};
+ − 67 }
+ − 68
+ − 69 if(exists $chr2dbsnp_vcf_lines{$chr}->{$pos}){
+ − 70 #print STDERR "Checking existing SNP data for $chr:$pos -> ", join("\n", @{$chr2dbsnp_vcf_lines{$chr}->{$pos}}), "\n";
+ − 71 for(@{$chr2dbsnp_vcf_lines{$chr}->{$pos}}){
+ − 72 my @fields = split /\t/, $_;
+ − 73 for my $var (split /,/, $fields[4]){
+ − 74 # Allows for different reference seqs between dbSNP and input, assuming patches only
+ − 75 if(length($fields[3]) == length($ref) and ($var eq $variant or $ref eq $var and $variant eq $fields[3])){
+ − 76 my ($freq, $subpop) = ("","");
+ − 77 $freq = $1 if $fields[7] =~ /(?:\A|;)MMAF=(0\.\d+)(?:;|\Z)/;
+ − 78 $subpop = $1 if $fields[7] =~ /(?:\A|;)MMAF_SRC=(\S+?)(?:;|\Z)/;
+ − 79 $dbsnp_info{$key} = [$subpop, $freq || "NA", $fields[2]];
+ − 80 return @{$dbsnp_info{$key}};
+ − 81 }
+ − 82 }
+ − 83 }
+ − 84 }
+ − 85 $dbsnp_info{$key} = ["novel", "NA", "NA"];
+ − 86 return @{$dbsnp_info{$key}};
+ − 87 }
+ − 88
+ − 89 sub record_snv{
+ − 90 my $line = join("", @_);
+ − 91 push @snvs, $line;
+ − 92
+ − 93 my @fields = split /\t/, $line;
+ − 94 my $prop_info_key = $fields[9];
+ − 95 my ($chr,$pos,$ref,$variant) = split /:/, $prop_info_key;
+ − 96 $chr2variant_locs{$chr} = {} unless exists $chr2variant_locs{$chr};
+ − 97 return unless $ref; # ref not defined for CNVs
+ − 98 # Need to grab whole range for MNPs
+ − 99 for(my $i = 0; $i < length($ref); $i++){
+ − 100 $chr2variant_locs{$chr}->{$pos+$i} = 1;
+ − 101 }
+ − 102 }
+ − 103
+ − 104 sub retrieve_vcf_lines($$$){
+ − 105 my ($dbsnp_file, $internal_snp_file, $chr) = @_;
+ − 106
+ − 107 my (%dbsnp_lines, %internal_snp_lines);
+ − 108
+ − 109 if(not defined $dbsnp_file or not exists $chr2variant_locs{$chr}){
+ − 110 return ({}, {}, {}, {}); # no data requested for this chromosome
+ − 111 }
+ − 112
+ − 113 # build up the request
+ − 114 my @tabix_regions;
+ − 115 my @var_locs = keys %{$chr2variant_locs{$chr}};
+ − 116 # sort by variant start location
+ − 117 for my $var_loc (sort {$a <=> $b} @var_locs){
+ − 118 push @tabix_regions, $chr.":".$var_loc."-".$var_loc;
+ − 119 }
+ − 120 for(my $i = 0; $i <= $#tabix_regions; $i += $max_args){ # chunkify tabix request if too many for the system to handle
+ − 121 my $end = $i + $max_args > $#tabix_regions ? $#tabix_regions : $i + $max_args;
+ − 122 my $regions = "'".join("' '", @tabix_regions[$i..$end])."'";
+ − 123 # From file is very slow for some reason
+ − 124 #my $regions_file = "/tmp/vcf2hgvs_$$.bed";
+ − 125 #open(REQ_BED, ">$regions_file")
+ − 126 # or die "Cannot open $regions_file for writing: $!\n";
+ − 127 #print REQ_BED join("\n", @tabix_regions), "\n";
+ − 128 #close(REQ_BED);
+ − 129
+ − 130 # retrieve the data
+ − 131 die "Cannot find dbSNP VCF file $dbsnp_file\n" if not -e $dbsnp_file;
+ − 132
+ − 133 open(VCF, "tabix $dbsnp_file $regions |")
+ − 134 or die "Cannot run tabix on $dbsnp_file (args ".substr($regions, 0, length($regions)>100? 100 : length($regions))."): $!\n";
+ − 135 while(<VCF>){
+ − 136 #if(/^(\S+\t(\d+)(?:\t\S+){6})/ and grep {$_ eq $2} @var_locs){ # take only main columns to save room, if possible
+ − 137 if(/^(\S+\t(\d+)(?:\t\S+){6})/ and exists $chr2variant_locs{$chr}->{$2}){ # take only main columns to save room, if possible
+ − 138 $dbsnp_lines{$2} = [] unless exists $dbsnp_lines{$2};
+ − 139 push @{$dbsnp_lines{$2}}, $1;
+ − 140 }
+ − 141 }
+ − 142 close(VCF);
+ − 143
+ − 144 if($internal_snp_file){
+ − 145 die "Cannot find internal VCF file $internal_snp_file\n" if not -e $internal_snp_file;
+ − 146 open(VCF, "tabix $internal_snp_file $regions |")
+ − 147 or die "Cannot run tabix on $internal_snp_file: $!\n";
+ − 148 while(<VCF>){
+ − 149 #if(/^(\S+\t(\d+)(?:\t\S+){6})/ and grep {$_ eq $2} @var_locs){ # take only main columns to save room, if possible
+ − 150 if(/^(\S+\t(\d+)(?:\t\S+){5})/ and exists $chr2variant_locs{$chr}->{$2}){ # take only main columns to save room, if possible
+ − 151 $internal_snp_lines{$2} = [] unless exists $internal_snp_lines{$2};
+ − 152 push @{$internal_snp_lines{$2}}, $1;
+ − 153 }
+ − 154 }
+ − 155 close(VCF);
+ − 156 }
+ − 157 }
+ − 158
+ − 159 #unlink $regions_file;
+ − 160
+ − 161 return (\%dbsnp_lines, \%internal_snp_lines);
+ − 162 }
+ − 163
+ − 164 sub prop_info_key{
+ − 165 my($chr,$pos,$ref,$variant,$exon_edge_dist) = @_;
+ − 166
+ − 167 $chr =~ s/^chr//;
+ − 168 if($chr eq "M"){
+ − 169 $chr = "MT"; # NCBI uses different name for mitochondrial chromosome
+ − 170 $pos-- if $pos >= 3107; # also, doesn't keep the old positioning (historical)
+ − 171 }
+ − 172 return join(":", $chr,$pos,$ref,$variant, ($exon_edge_dist ? $exon_edge_dist : ""));
+ − 173 }
+ − 174
+ − 175 sub prop_info($$$){
+ − 176 my($snpfile,$internal_snps_file,$prop_info_key) = @_;
+ − 177
+ − 178 my($chr,$pos,$ref,$variant) = split /:/, $prop_info_key;
+ − 179
+ − 180 # is this the first call for this chromosome? If so, retrieve the VCF lines for it en masse
+ − 181 if(not exists $chr2dbsnp_vcf_lines{$chr}){
+ − 182 ($chr2dbsnp_vcf_lines{$chr}, $chr2internal_vcf_lines{$chr}) = retrieve_vcf_lines($snpfile,$internal_snps_file,$chr);
+ − 183 }
+ − 184 my $internal_maf = 0;
+ − 185 if($internal_snps_file){
+ − 186 $internal_maf = internal_prop($chr,$pos,$ref,$variant);
+ − 187 $internal_maf = 0 if $internal_maf eq "NA";
+ − 188 }
+ − 189
+ − 190 my @results = dbsnp_info($chr,$pos,$ref,$variant);
+ − 191
+ − 192 # Not all entries have a proportion in dbSNP
+ − 193 return $internal_snps_file ? ($ref, $variant, @results, $internal_maf) : ($ref, $variant, @results);
+ − 194 }
+ − 195
+ − 196 #offset a given HGVS nomenclature position (single position only) by a given number of bases
+ − 197 sub hgvs_plus($$){
+ − 198 my ($hgvs, $offset) = @_;
+ − 199 if($hgvs =~ /^(\S+)(-\d+)(.*)/){
+ − 200 # all negative
+ − 201 if($2+$offset<0){
+ − 202 return $1.($2+$offset).$3;
+ − 203 }
+ − 204 # switches to positive, need to mod
+ − 205 else{
+ − 206 return $1+($2+$offset);
+ − 207 }
+ − 208 }
+ − 209 elsif($hgvs =~ /^(\S+)\+(\d+)(.*)/){
+ − 210 # all positive
+ − 211 if($2+$offset>0){
+ − 212 return $1."+".($2+$offset).$3;
+ − 213 }
+ − 214 # switches to negative, need to mod
+ − 215 else{
+ − 216 return $1+($2+$offset);
+ − 217 }
+ − 218 }
+ − 219 elsif($hgvs =~ /^(-?\d+)(.*)/){
+ − 220 # special case if offset spans -/+ since there is no position 0
+ − 221 if($1 < 0 and $1+$offset >= 0){
+ − 222 $offset++;
+ − 223 }
+ − 224 elsif($1 > 0 and $1+$offset <= 0){
+ − 225 $offset--;
+ − 226 }
+ − 227 return ($1+$offset).$2;
+ − 228 }
+ − 229 else{
+ − 230 die "Cannot convert $hgvs to a new offset ($offset), only single base position nomenclature is currently supported\n";
+ − 231 }
+ − 232 }
+ − 233
+ − 234 # offset a given position by a given number of bases,
+ − 235 # taking into account that if the new offset crosses the threshold in the last argument,
+ − 236 # HGVS boundary nomenclature has to be introduced
+ − 237 sub hgvs_plus_exon($$$){
+ − 238 my ($pos, $offset, $boundary) = @_;
+ − 239
+ − 240 # special case if offset spans -/+ since there is no position 0
+ − 241 if($pos =~ /^(-?\d+)(.*)/){
+ − 242 if($1 < 0 and $1+$offset >= 0){
+ − 243 $offset++;
+ − 244 }
+ − 245 elsif($1 > 0 and $1+$offset <= 0){
+ − 246 $offset--;
+ − 247 }
+ − 248 }
+ − 249 my $new_pos = $pos + $offset;
+ − 250 if($new_pos > $boundary and $pos <= $boundary){
+ − 251 # just moved into an intron 3'
+ − 252 $new_pos = $boundary."+".($new_pos-$boundary);
+ − 253 }
+ − 254 elsif($new_pos < $boundary and $pos >= $boundary){
+ − 255 # just moved into an intron 5'
+ − 256 $new_pos = $boundary.($new_pos-$boundary);
+ − 257 }
+ − 258 return $new_pos;
+ − 259 }
+ − 260
+ − 261 # given a nucleotide position, calculates the AA there (assumes coding region)
+ − 262 sub getCodonFromSeq($$$$){
+ − 263 my ($chr_ref, $location, $frame_offset, $strand) = @_;
+ − 264
+ − 265 my $codon;
+ − 266 if($strand eq "+"){
+ − 267 $codon = substr($$chr_ref, $location-1-$frame_offset, 3);
+ − 268 }
+ − 269 else{
+ − 270 $codon = substr($$chr_ref, $location-3+$frame_offset, 3);
+ − 271 $codon = reverse($codon);
+ − 272 $codon =~ tr/ACGTacgt/TGCAtgca/;
+ − 273 }
+ − 274 return $codon;
+ − 275 }
+ − 276
+ − 277 sub getCodonFromSeqIndex($$$$){
+ − 278 my ($chr, $location, $frame_offset, $strand) = @_;
+ − 279
+ − 280 my $codon;
+ − 281 if($strand eq "+"){
+ − 282 $codon = $fasta_index->fetch($chr.":".($location-$frame_offset)."-".($location-$frame_offset+2));
+ − 283 }
+ − 284 else{
+ − 285 $codon = $fasta_index->fetch($chr.":".($location-2+$frame_offset)."-".($location+$frame_offset));
+ − 286 $codon = reverse($codon);
+ − 287 $codon =~ tr/ACGTacgt/TGCAtgca/;
+ − 288 }
+ − 289 return $codon;
+ − 290 }
+ − 291
+ − 292 sub getAAFromSeq($$$$$){
+ − 293 return $_[4]->translate(getCodonFromSeq($_[0], $_[1], $_[2], $_[3]));
+ − 294 }
+ − 295
+ − 296 sub getAAFromSeqIndex($$$$$){
+ − 297 # convert codon to AA
+ − 298 if(exists $transl_except{"$_[0]:$_[1]"}){
+ − 299 return $transl_except{"$_[0]:$_[1]"};
+ − 300 }
+ − 301 else{
+ − 302 return $_[4]->translate(getCodonFromSeqIndex($_[0], $_[1], $_[2], $_[3]));
+ − 303 }
+ − 304 }
+ − 305
+ − 306 sub hgvs_protein{
+ − 307 my ($chr, $location, $ref, $variant, $cdna_pos, $strand, $transl_table) = @_;
+ − 308
+ − 309 if(substr($ref,0,1) eq substr($variant,0,1)){
+ − 310 substr($ref,0,1) = "";
+ − 311 substr($variant,0,1) = "";
+ − 312 $location++;
+ − 313 if($strand eq "-"){
+ − 314 $cdna_pos--;
+ − 315 }
+ − 316 else{
+ − 317 $cdna_pos++;
+ − 318 }
+ − 319 }
+ − 320
+ − 321 if($cdna_pos !~ /^\d+/){
+ − 322 die "Aborting: got illegal cDNA position ($cdna_pos) for protein HGVS conversion of position ",
+ − 323 "$location, ref $ref, variant $variant. Please correct the program code.\n";
+ − 324 }
+ − 325 # Get the correct frame for the protein translation, to know what codons are affected
+ − 326 my $aapos = int(($cdna_pos-1)/3)+1;
+ − 327
+ − 328 # does it destroy the start codon?
+ − 329 if($cdna_pos < 4){ # assumes animal codon usage
+ − 330 return "p.0?"; # indicates start codon missing, unsure of effect
+ − 331 }
+ − 332
+ − 333 my $table = $transl_table ne $default_transl_table ? # non standard translation table requested
+ − 334 Bio::Tools::CodonTable->new(-id=>$transl_table) : $codonTable;
+ − 335
+ − 336 my $frame_offset = ($cdna_pos-1)%3;
+ − 337 my $origAA = getAAFromSeqIndex($chr, $location, $frame_offset, $strand, $table);
+ − 338 # take 100000 bp on either side for translation context of variant seq
+ − 339 my $five_prime_buffer = $location < 10000 ? $location-1 : 10000;
+ − 340 my $mutSeq = $fasta_index->fetch($chr.":".($location-$five_prime_buffer)."-".($location+10000));
+ − 341
+ − 342 # substitute all of the immediately adjacent variants in phase with this one to get the correct local effect
+ − 343 substr($mutSeq, $five_prime_buffer, length($ref)) = $variant;
+ − 344
+ − 345 # does it cause a frameshift?
+ − 346 my $length_diff = length($variant)-length($ref);
+ − 347 if($length_diff%3){ # insertion or deletion not a multiple of three
+ − 348 my $fs_codon = getCodonFromSeq(\$mutSeq, $five_prime_buffer+1, $frame_offset, $strand);
+ − 349 my $ext = 0;
+ − 350 my $newAA;
+ − 351 do{
+ − 352 $ext++;
+ − 353 # The "NA"s below make it so that we don't pick up any translation exceptions from the original reference annotation
+ − 354 if($strand eq "+"){
+ − 355 $newAA = getAAFromSeq(\$mutSeq, $five_prime_buffer+1+$ext*3, $frame_offset, $strand, $table);
+ − 356 }
+ − 357 else{
+ − 358 $newAA = getAAFromSeq(\$mutSeq, $five_prime_buffer+1-$ext*3, $frame_offset, $strand, $table);
+ − 359 }
+ − 360 } while($newAA ne "*");
+ − 361
+ − 362 return "p.".$origAA.$aapos.$table->translate($fs_codon)."fs*$ext";
+ − 363 }
+ − 364
+ − 365 # does it cause a stop codon to be lost?
+ − 366 if($origAA eq "*"){
+ − 367 my $stopChangeCodon = getCodonFromSeq(\$mutSeq, $five_prime_buffer+1, $frame_offset, $strand);
+ − 368 # still a stop after the mutation (ignore translation exceptions)
+ − 369 if($table->is_ter_codon($stopChangeCodon)){
+ − 370 return "p.*$aapos=";
+ − 371 }
+ − 372 # calculate the new stop, assuming there aren't mutations downstream in candidate stop codons
+ − 373 my $ext = 0;
+ − 374 my $newCodon;
+ − 375 do{
+ − 376 if($strand eq "+"){
+ − 377 $newCodon = getCodonFromSeq(\$mutSeq, $five_prime_buffer+1+(++$ext*3), $frame_offset, $strand);
+ − 378 }
+ − 379 else{
+ − 380 $newCodon = getCodonFromSeq(\$mutSeq, $five_prime_buffer+1-(++$ext*3), $frame_offset, $strand);
+ − 381 }
+ − 382 } while(not $table->is_ter_codon($newCodon));
+ − 383
+ − 384 return "p.*".$aapos.$table->translate($stopChangeCodon)."ext*".$ext;
+ − 385 }
+ − 386
+ − 387 # if we get this far, it's a "regular" AA level change
+ − 388 my $origAAs = "";
+ − 389 for(my $i = 0; $i < length($ref)+$frame_offset; $i+=3){
+ − 390 my $oldAA = getAAFromSeqIndex($chr, $location+$i, $frame_offset, $strand, $table);
+ − 391 if($strand eq "+"){
+ − 392 $origAAs .= $oldAA;
+ − 393 }
+ − 394 else{
+ − 395 $origAAs = $oldAA . $origAAs;
+ − 396 }
+ − 397 }
+ − 398 my $newAAs = "";
+ − 399 for(my $i = 0; $i < length($variant)+$frame_offset; $i+=3){
+ − 400 # NA means we don't take translation exceptions from the original
+ − 401 my $newAA = getAAFromSeq(\$mutSeq, $five_prime_buffer+1+$i, $frame_offset, $strand, $table);
+ − 402 if($strand eq "+"){
+ − 403 $newAAs .= $newAA;
+ − 404 }
+ − 405 else{
+ − 406 $newAAs = $newAA . $newAAs;
+ − 407 }
+ − 408 }
+ − 409
+ − 410 # silent
+ − 411 if($origAAs eq $newAAs){
+ − 412 return "p.".$origAAs.$aapos."=";
+ − 413 }
+ − 414
+ − 415 # minimize the difference if there are leading or trailing AAs the same
+ − 416 my $delLength = length($ref);
+ − 417 while(substr($newAAs, 0, 1) eq substr($origAAs, 0, 1)){
+ − 418 $newAAs = substr($newAAs, 1);
+ − 419 $origAAs = substr($origAAs, 1);
+ − 420 $location+=3;
+ − 421 $delLength-=3;
+ − 422 $aapos++;
+ − 423 }
+ − 424 while(substr($newAAs, -1) eq substr($origAAs, -1)){
+ − 425 $newAAs = substr($newAAs, 0, length($newAAs)-1);
+ − 426 $origAAs = substr($origAAs, 0, length($origAAs)-1);
+ − 427 }
+ − 428
+ − 429 # insertion
+ − 430 if(length($origAAs) == 0){
+ − 431 my $insAAs = getAAFromSeqIndex($chr,$location-3,$frame_offset,$strand,$table).($aapos-1)."_".
+ − 432 getAAFromSeqIndex($chr,$location,$frame_offset,$strand,$table);
+ − 433 return "p.".$insAAs.$aapos."ins".$newAAs;
+ − 434 }
+ − 435 # deletion
+ − 436 elsif(length($newAAs) == 0){
+ − 437 my $delAAs;
+ − 438 if(length($origAAs) == 1){
+ − 439 $delAAs = getAAFromSeqIndex($chr,$location,$frame_offset,$strand,$table).$aapos; # single AA deletion
+ − 440 }
+ − 441 else{ # deleting a stretch
+ − 442 if($strand eq "+"){
+ − 443 my $endPoint = $location+$delLength-1;
+ − 444 $delAAs = getAAFromSeqIndex($chr,$location,$frame_offset,$strand,$table).$aapos."_".
+ − 445 getAAFromSeqIndex($chr,$endPoint,$frame_offset,$strand,$table).($aapos+int(($delLength-1)/3));
+ − 446 }
+ − 447 else{
+ − 448 my $endPoint = $location-$delLength+1;
+ − 449 $delAAs = getAAFromSeqIndex($chr,$endPoint,$frame_offset,$strand,$table).($aapos-int(($delLength-1)/3))."_".
+ − 450 getAAFromSeqIndex($chr,$location,$frame_offset,$strand,$table).$aapos;
+ − 451 }
+ − 452 }
+ − 453 return "p.".$delAAs."del";
+ − 454 }
+ − 455 else{
+ − 456 # substitution
+ − 457 if(length($origAAs) == 1 and length($newAAs) == 1){
+ − 458 return "p.".$origAAs.$aapos.$newAAs;
+ − 459 }
+ − 460 # indel
+ − 461 elsif(length($origAAs) != 1){
+ − 462 # convert ref stretch into range syntax
+ − 463 if($strand eq "+"){
+ − 464 $origAAs = substr($origAAs, 0, 1).$aapos."_".substr($origAAs, -1).($aapos+length($origAAs)-1);
+ − 465 }
+ − 466 else{
+ − 467 $origAAs = substr($origAAs, 0, 1).($aapos-length($origAAs)+1)."_".substr($origAAs, -1).$aapos;
+ − 468 }
+ − 469 }
+ − 470 return "p.".$origAAs."delins".$newAAs;
+ − 471 }
+ − 472 return ("NA", "");
+ − 473 }
+ − 474
+ − 475 sub z2p{
+ − 476 if(not defined $zed){
+ − 477 $zed = new Statistics::Zed;
+ − 478 }
+ − 479 my $p = $zed->z2p(value => $_[0]);
+ − 480 return $p < 0.0000000001 ? 0 : $p;
+ − 481 }
+ − 482 sub gq2p{
+ − 483 return $_[0] > 200 ? 0 : 10**($_[0]/-10);
+ − 484 }
+ − 485
+ − 486 my ($multi_phased, $min_depth, $flanking_bases, $dbsnp, $internal_snp, $genename_bed_file, $dir_1000G, $dir_esp6500, $min_pvalue, $mappability_file, $reference_file, $samtools_phasing_file, $exons_file, $input_file, $output_file, $cnv_file, $dgv_file, $which_chr, $enrichment_regions_file, $rare_variant_prop);
+ − 487 $multi_phased = 0;
+ − 488 $min_depth = 2;
+ − 489 $flanking_bases = 30;
+ − 490 $min_pvalue = 0.01;
+ − 491 $min_prop = 0.14;
+ − 492 $rare_variant_prop = 0.05;
+ − 493 $input_file = "-"; # STDIN by default
+ − 494 $output_file = "-"; # STDOUT by default
+ − 495 $default_transl_table = "1"; # assumes NCBI 'Standard' table, unless it is an argument to the script...
+ − 496 &GetOptions("d=i" => \$min_depth,
+ − 497 "f=i" => \$flanking_bases,
+ − 498 "s=s" => \$dbsnp,
+ − 499 "t=s" => \$dir_1000G,
+ − 500 "n=s" => \$dir_esp6500,
+ − 501 "u=s" => \$internal_snp,
+ − 502 "q" => \$quiet,
+ − 503 "p=f" => \$min_pvalue,
+ − 504 "h=f" => \$min_prop,
+ − 505 "m=s" => \$mappability_file,
+ − 506 "r=s" => \$reference_file,
+ − 507 "z=s" => \$samtools_phasing_file,
+ − 508 "e=s" => \$exons_file,
+ − 509 "i=s" => \$input_file,
+ − 510 "c=s" => \$cnv_file,
+ − 511 "g=s" => \$dgv_file,
+ − 512 "b=s" => \$genename_bed_file,
+ − 513 "w=s" => \$which_chr,
+ − 514 "o=s" => \$output_file,
+ − 515 "a=i" => \$default_transl_table,
+ − 516 "v=f" => \$rare_variant_prop,
+ − 517 "x=s" => \$enrichment_regions_file); # if enrichment regions are specified, variants without a transcript model but in these ranges will be reported
+ − 518
+ − 519 if(($input_file ne "/dev/null" and not defined $reference_file) or
+ − 520 not defined $exons_file or
+ − 521 (defined $cnv_file and not defined $dgv_file)){
+ − 522 die "Usage: $0 [-v(ersion)] [-q(uiet)] [-w(hich) contig_to_report (default is all)] [-d(epth of variant reads req'd) #] [-v(ariant max freq to count as rare)] [-f(lanking exon bases to report) #] [-p(robability of random genotype, maximum to report) 0.#]\n",
+ − 523 " [-h(et proportion of variant reads, minimum to report) 0.#] [-c(opy number) variants_file.bed -g(enomic structural) variants_control_db.txt.gz] [-z file_containing_samtools_phase_output.txt]\n",
+ − 524 " [-t(housand) genomes_integrated_vcfs_gz_dir] [-n ESP6500_dir] [-u(ser) specified_population.vcf.gz] [-m(appability) crg_file.bed]\n",
+ − 525 " [-x enrichment_regions_file.bed] [-a(mino) acid translation table number from NCBI]\n",
+ − 526 " [-i(nput) genotypes.vcf <-r(eference) sequence_file.fasta>] [-o(utput) hgvs_file.tsv] [-s(np) database_from_ncbi.vcf.gz]\n",
+ − 527 " <-b(ed) file of named gene regions.bed> <-e(xons) file.gtf>\n\n",
+ − 528 "Input gz files must be indexed with Tabix.\nDefault input is STDIN, default output is STDOUT. Note: if -c is specified, polyploidies are are assume to be proximal. Other defaults: -d 2, -v 0.05, -f 30, -p 0.01, -h 0.14 -a 1\nReference sequence is not strictly necessary if only CNV are being annotated.\n";
+ − 529 }
+ − 530
+ − 531 print STDERR "Considering $flanking_bases flanking bases for variants as well\n" unless $quiet;
+ − 532
+ − 533 $codonTable = new Bio::Tools::CodonTable(id => $default_transl_table);
+ − 534
+ − 535 my %enrichment_regions;
+ − 536 # Note, we assume the regions are non-overlapping
+ − 537 if(defined $enrichment_regions_file){
+ − 538 print STDERR "Loading enrichment regions...\n" unless $quiet;
+ − 539 open(BED, $enrichment_regions_file)
+ − 540 or die "Cannot open $enrichment_regions_file for reading: $!\n";
+ − 541 while(<BED>){
+ − 542 chomp;
+ − 543 my @F = split /\t/, $_;
+ − 544 $enrichment_regions{$F[0]} = [] if not exists $enrichment_regions{$F[0]};
+ − 545 push @{$enrichment_regions{$F[0]}}, [$F[1], $F[2]];
+ − 546 }
+ − 547 close(BED);
+ − 548 }
+ − 549 for my $chr (keys %enrichment_regions){ # sort by start
+ − 550 $enrichment_regions{$chr} = [sort {$a->[0] <=> $b->[0]} @{$enrichment_regions{$chr}}];
+ − 551 }
+ − 552
+ − 553 if(defined $reference_file){
+ − 554 print STDERR "Scanning reference FastA info\n" unless $quiet;
+ − 555 if(not -e $reference_file){
+ − 556 die "Reference FastA file ($reference_file) does not exist.\n";
+ − 557 }
+ − 558 if(not -e $reference_file.".fai" and not -w dirname($reference_file)){
+ − 559 die "Reference FastA file ($reference_file) is not indexed, and the directory is not writable.\n";
+ − 560 }
+ − 561 $fasta_index = Bio::DB::Sam::Fai->load($reference_file);
+ − 562 }
+ − 563
+ − 564 my %chr2mappability;
+ − 565 if(defined $mappability_file){
+ − 566 print STDERR "Reading in mappability data\n" unless $quiet;
+ − 567 my ($nmer) = $mappability_file =~ /(\d+).*?$/;
+ − 568 die "Cannot determine nmer from nmer file name $mappability_file, aborting\n" unless $nmer;
+ − 569 open(MAP, $mappability_file)
+ − 570 or die "Cannot open mappability data file $mappability_file for reading: $!\n";
+ − 571 <MAP>; # header
+ − 572 while(<MAP>){
+ − 573 next if /^#/;
+ − 574 chomp;
+ − 575 my @F = split /\t/, $_;
+ − 576 my $x = int(1/$F[3]+0.5);
+ − 577 $chr2mappability{$F[0]} = Set::IntervalTree->new() if not exists $chr2mappability{$F[0]};
+ − 578 $chr2mappability{$F[0]}->insert("non-unique mapping region (x$x)", $F[1], $F[2]+$nmer-1);
+ − 579 }
+ − 580 close(MAP);
+ − 581 }
+ − 582
+ − 583 # Is phasing data provided?
+ − 584 if(defined $samtools_phasing_file){
+ − 585 print STDERR "Reading in phasing data\n" unless $quiet;
+ − 586 open(PHASE, $samtools_phasing_file)
+ − 587 or die "Cannot open phasing data file $samtools_phasing_file for reading: $!\n";
+ − 588 my $phase_range;
+ − 589 while(<PHASE>){
+ − 590 if(/^PS/){
+ − 591 chomp;
+ − 592 my @F = split /\t/, $_;
+ − 593 $phase_range = "$F[2]-$F[3]";
+ − 594 }
+ − 595 if(/^M[12]/){
+ − 596 chomp;
+ − 597 my @F = split /\t/, $_;
+ − 598 #ignore strange cases where haplotype reference has no cases (weird samtools call)
+ − 599 next if $F[9] == 0 or $F[7] == 0;
+ − 600 my $chr = $F[1];
+ − 601 next if defined $which_chr and not $chr eq $which_chr;
+ − 602 my $pos = $F[3];
+ − 603 #print STDERR "Recording phase for $chr:$pos:$F[4] , $chr:$pos:$F[5] as A-$chr:$phase_range and B-$chr:$phase_range\n" if $pos == 12907379;
+ − 604 if(($F[10]+$F[8])/($F[9]+$F[7]) >= $min_prop){ # error meets reporting threshold
+ − 605 $chr2caveats{"$chr:$pos"} .= "; " if exists $chr2caveats{"$chr:$pos"};
+ − 606 $chr2caveats{"$chr:$pos"} .= "inconsistent haplotype phasing";
+ − 607 }
+ − 608 else{ # appears to be a genuine phasing
+ − 609 $chr2phase{"$chr:$pos:$F[4]"} = "A-$chr:$phase_range"; # grouping for haplotype
+ − 610 $chr2phase{"$chr:$pos:$F[5]"} = "B-$chr:$phase_range"; # grouping for haplotype
+ − 611 }
+ − 612 }
+ − 613 }
+ − 614 close(PHASE);
+ − 615 }
+ − 616
+ − 617 # Check the VCF file to see if contains phase data
+ − 618 open(VCFIN, $input_file)
+ − 619 or die "Cannot open $input_file for reading: $!\n";
+ − 620 my $phase_chr = "";
+ − 621 my @phase_dataA;
+ − 622 my @phase_dataB;
+ − 623 while(<VCFIN>){
+ − 624 if(/^\s*(?:#|$)/){ # blank or hash comment
+ − 625 next;
+ − 626 }
+ − 627 my @F = split /\t/, $_;
+ − 628 next if exists $chr2caveats{"$F[0]:$F[1]"} and $chr2caveats{"$F[0]:$F[1]"} =~ /inconsistent haplotype phasing/;
+ − 629 # | indicates phased
+ − 630 if($F[8] =~ m(^(\d+)\|(\d+):)){
+ − 631 next if $1 eq $2; # not useful to us (actually would mess up phase combining later on), but is provided sometimes
+ − 632 # start of a phasing block
+ − 633 if($phase_chr eq ""){
+ − 634 $phase_chr = $F[0];
+ − 635 }
+ − 636 my @vars = split /,/, $F[4];
+ − 637 if($1 > @vars){
+ − 638 die "Invalid VCF file (line #$.): First haplotype listed as $1, but only ", scalar(@vars), " variants were provided (", join(",", @vars), "\n";
+ − 639 }
+ − 640 if($2 > @vars){
+ − 641 die "Invalid VCF file (line #$.): Second haplotype listed as $1, but only ", scalar(@vars), " variants were provided (", join(",", @vars), "\n";
+ − 642 }
+ − 643 unshift @vars, $F[3];
+ − 644 push @phase_dataA, [$F[1], $vars[$1]];
+ − 645 push @phase_dataB, [$F[1], $vars[$2]];
+ − 646 }
+ − 647 # non phased het call, ends any phasing block there might be
+ − 648 elsif($F[8] =~ m(^0/1)){
+ − 649 # Did we just finish a phased block? If so, output it.
+ − 650 if(@phase_dataA > 1){
+ − 651 my $phase_def = "G-$phase_chr:".$phase_dataA[0]->[0]."-".$phase_dataA[$#phase_dataA]->[0];
+ − 652 for my $d (@phase_dataA){
+ − 653 my ($p, $v) = @$d;
+ − 654 if(exists $chr2phase{"$phase_chr:$p:$v"}){
+ − 655 $chr2phase{"$phase_chr:$p:$v"} .= ",$phase_def";
+ − 656 $multi_phased ||= 1;
+ − 657 }
+ − 658 else{
+ − 659 $chr2phase{"$phase_chr:$p:$v"} = $phase_def;
+ − 660 }
+ − 661 }
+ − 662 $phase_def = "H-$phase_chr:".$phase_dataB[0]->[0]."-".$phase_dataB[$#phase_dataB]->[0];
+ − 663 for my $d (@phase_dataB){
+ − 664 my ($p, $v) = @$d;
+ − 665 if(exists $chr2phase{"$phase_chr:$p:$v"}){
+ − 666 $chr2phase{"$phase_chr:$p:$v"} = ",$phase_def";
+ − 667 $multi_phased ||= 1;
+ − 668 }
+ − 669 else{
+ − 670 $chr2phase{"$phase_chr:$p:$v"} = $phase_def;
+ − 671 }
+ − 672 }
+ − 673 }
+ − 674 if($phase_chr ne ""){
+ − 675 $phase_chr = "";
+ − 676 @phase_dataA = ();
+ − 677 @phase_dataB = ();
+ − 678 }
+ − 679 }
+ − 680 }
+ − 681
+ − 682 print STDERR "Reading in feature GTF data..." unless $quiet;
+ − 683 my %feature_range; # chr => transcript_id => [[genomic_exon_start,genomic_exon_end,cdna_start_pos],...]
+ − 684 my %feature_intervaltree; # chr => transcript_id => [[genomic_exon_start,genomic_exon_end,cdna_start_pos],...]
+ − 685 my %feature_strand; # transcript_id => +|-
+ − 686 my $feature_count = 0;
+ − 687 my %feature_min;
+ − 688 my %feature_max;
+ − 689 my %feature_cds_min;
+ − 690 my %feature_cds_max;
+ − 691 my %feature_contig;
+ − 692 my %feature_length;
+ − 693 my %feature_type;
+ − 694 my %feature_transl_table; # note alternate translation table usage
+ − 695 my %chr_read;
+ − 696 open(GTF, $exons_file)
+ − 697 or die "Cannot open $exons_file for reading: $!\n";
+ − 698 while(<GTF>){
+ − 699 next if /^\s*#/;
+ − 700 my @fields = split /\t/, $_;
+ − 701 next if defined $which_chr and $fields[0] ne $which_chr and "chr$fields[0]" ne $which_chr and $fields[0] ne "chr$which_chr";
+ − 702
+ − 703 if($fields[2] eq "exon" or $fields[2] eq "CDS"){
+ − 704 next unless $fields[$#fields] =~ /transcript_id \"(.*?)\"/o;
+ − 705 my $parent = $1;
+ − 706 if(not $quiet and not exists $chr_read{$fields[0]}){
+ − 707 print STDERR " $fields[0]";
+ − 708 $chr_read{$fields[0]} = 1;
+ − 709 }
+ − 710 if(not exists $feature_strand{$parent}){
+ − 711 $feature_strand{$parent} = $fields[6];
+ − 712 $feature_contig{$parent} = $fields[0];
+ − 713 if($fields[$#fields] =~ /transcript_type \"(.*?)\"/){
+ − 714 $feature_type{$parent} = $1;
+ − 715 }
+ − 716 else{
+ − 717 $feature_type{$parent} = "NA";
+ − 718 }
+ − 719 }
+ − 720 if($fields[2] eq "CDS"){
+ − 721 #print STDERR "CDS value for $parent is $fields[2]..$fields[3]\n";
+ − 722 if(not exists $feature_cds_min{$parent} or $fields[3] < $feature_cds_min{$parent}){
+ − 723 $feature_cds_min{$parent} = $fields[3];
+ − 724 }
+ − 725 if(not exists $feature_cds_max{$parent} or $fields[4] > $feature_cds_max{$parent}){
+ − 726 $feature_cds_max{$parent} = $fields[4];
+ − 727 }
+ − 728 if($fields[$#fields] =~ /transl_table \"(\d+)\"/){
+ − 729 $feature_transl_table{$parent} = $1; #assume one translation table per CDS, which should be reasonable
+ − 730 }
+ − 731 while($fields[$#fields] =~ /transl_except \"pos:(\S+?),aa:(\S+?)\"/g){
+ − 732 my $pos = $1;
+ − 733 my $new_aa = $2; # needs to change from three letter code to 1
+ − 734 if($new_aa =~ /^ter/i){ # can be funny so have special case (allows TERM, etc.)
+ − 735 $new_aa = "*";
+ − 736 }
+ − 737 elsif(length($new_aa) == 3){
+ − 738 $new_aa = Bio::SeqUtils->new()->seq3in($new_aa);
+ − 739 }
+ − 740 if($pos =~ /^(\d+)\.\.(\d+)/){
+ − 741 for my $p ($1..$2){
+ − 742 $transl_except{"$fields[0]:$p"} = $new_aa;
+ − 743 }
+ − 744 }
+ − 745 else{
+ − 746 $transl_except{"$fields[0]:$pos"} = $new_aa;
+ − 747 }
+ − 748 }
+ − 749 next;
+ − 750 }
+ − 751 if(not exists $feature_min{$parent} or $fields[3] < $feature_min{$parent}){
+ − 752 $feature_min{$parent} = $fields[3];
+ − 753 }
+ − 754 if(not exists $feature_max{$parent} or $fields[4] > $feature_max{$parent}){
+ − 755 $feature_max{$parent} = $fields[4];
+ − 756 }
+ − 757
+ − 758 $feature_count++;
+ − 759 if(not exists $feature_range{$fields[0]}){
+ − 760 $feature_range{$fields[0]} = {}; # Chr => {parentID => [start,stop]}
+ − 761 $feature_intervaltree{$fields[0]} = Set::IntervalTree->new();
+ − 762 }
+ − 763 if(not exists $feature_range{$fields[0]}->{$parent}){
+ − 764 $feature_range{$fields[0]}->{$parent} = [];
+ − 765 }
+ − 766 push @{$feature_range{$fields[0]}->{$parent}}, [$fields[3],$fields[4]];
+ − 767 $feature_intervaltree{$fields[0]}->insert($parent, $fields[3], $fields[4]+1); # ranges need to have positive length for module to work properly
+ − 768 $feature_length{$parent} += $fields[4]-$fields[3]+1;
+ − 769 }
+ − 770 }
+ − 771 close(GTF);
+ − 772 print STDERR "\nFound $feature_count exons on ", scalar(keys %feature_range), " contigs in the GTF file\n" unless $quiet;
+ − 773
+ − 774 for my $contig (keys %feature_range){
+ − 775 for my $parent (keys %{$feature_range{$contig}}){
+ − 776 # sort by subrange start
+ − 777 my @feature_ranges = sort {$a->[0] <=> $b->[0]} @{$feature_range{$contig}->{$parent}};
+ − 778 $feature_range{$contig}->{$parent} = \@feature_ranges;
+ − 779 $feature_range{"chr".$contig}->{$parent} = \@feature_ranges if not $contig =~ /^chr/;
+ − 780 $feature_range{$1}->{$parent} = \@feature_ranges if $contig =~ /^chr(\S+)/;
+ − 781 }
+ − 782 }
+ − 783
+ − 784 # Calculate the cDNA position of the leftmost (reference strand) base for each exon
+ − 785 for my $contig (keys %feature_range){
+ − 786 for my $parent (keys %{$feature_range{$contig}}){
+ − 787 my @feature_ranges = @{$feature_range{$contig}->{$parent}};
+ − 788 if($feature_strand{$parent} eq "-"){
+ − 789 # set up utr offset for correct CDS coordinates
+ − 790 my $feature_offset = 0;
+ − 791 for(my $i = $#feature_ranges; $i >= 0; $i--){
+ − 792 last if not $feature_cds_max{$parent};
+ − 793 # exon is completely 5' of the start
+ − 794 if($feature_ranges[$i]->[0] > $feature_cds_max{$parent}){
+ − 795 $feature_offset -= $feature_ranges[$i]->[1]-$feature_ranges[$i]->[0]+1;
+ − 796 }
+ − 797 # exon with the cds start
+ − 798 elsif($feature_ranges[$i]->[1] >= $feature_cds_max{$parent} and
+ − 799 $feature_ranges[$i]->[0] <= $feature_cds_max{$parent}){
+ − 800 $feature_offset += $feature_cds_max{$parent} - $feature_ranges[$i]->[1];
+ − 801 last;
+ − 802 }
+ − 803 else{
+ − 804 die "The CDS for $parent (on negative strand) ends downstream ",
+ − 805 "($feature_cds_max{$parent}) of the an exon",
+ − 806 " (", $feature_ranges[$i]->[0], "), which is illogical. Please revise the GFF file provided.\n";
+ − 807 }
+ − 808 }
+ − 809 for(my $i = $#feature_ranges; $i >= 0; $i--){
+ − 810 $feature_offset += $feature_ranges[$i]->[1]-$feature_ranges[$i]->[0]+1;
+ − 811 $feature_ranges[$i]->[2] = $feature_offset-1;
+ − 812 }
+ − 813 }
+ − 814 else{ # positive strand
+ − 815 # set up utr offset for correct CDS coordinates
+ − 816 my $feature_offset = 0;
+ − 817 for(my $i = 0; $i <= $#feature_ranges; $i++){
+ − 818 last if not $feature_cds_min{$parent};
+ − 819 # All 5' utr exon
+ − 820 if($feature_ranges[$i]->[1] < $feature_cds_min{$parent}){
+ − 821 $feature_offset -= $feature_ranges[$i]->[1]-$feature_ranges[$i]->[0]+1;
+ − 822 }
+ − 823 # exon with the cds start
+ − 824 elsif($feature_ranges[$i]->[1] >= $feature_cds_min{$parent} and
+ − 825 $feature_ranges[$i]->[0] <= $feature_cds_min{$parent}){
+ − 826 $feature_offset -= $feature_cds_min{$parent} - $feature_ranges[$i]->[0];
+ − 827 last;
+ − 828 }
+ − 829 else{
+ − 830 die "The CDS for $parent starts upstream ($feature_cds_min{$parent}) of the first exon",
+ − 831 " (", $feature_ranges[$i]->[0], "), which is illogical. Please revise the GFF file provided.\n";
+ − 832 }
+ − 833 }
+ − 834 # assign cDNA coords for each exon to the third array element
+ − 835 for(my $i = 0; $i <= $#feature_ranges; $i++){
+ − 836 $feature_ranges[$i]->[2] = $feature_offset;
+ − 837 $feature_offset += $feature_ranges[$i]->[1]-$feature_ranges[$i]->[0]+1;
+ − 838 }
+ − 839 }
+ − 840 }
+ − 841 }
+ − 842
+ − 843 print STDERR "Reading in gene name definitions...\n" unless $quiet;
+ − 844 die "Data file $genename_bed_file does not exist, aborting.\n" if not -e $genename_bed_file;
+ − 845 my %gene_ids;
+ − 846 open(TAB, $genename_bed_file)
+ − 847 or die "Cannot open gene name BED file $genename_bed_file for reading: $!\n";
+ − 848 while(<TAB>){
+ − 849 chomp;
+ − 850 # format should be "chr start stop gene_name ..."
+ − 851 my @fields = split /\t/, $_;
+ − 852 next if $#fields < 3;
+ − 853 my $c = $fields[0];
+ − 854 if(not exists $gene_ids{$c}){
+ − 855 $gene_ids{$c} = Set::IntervalTree->new();
+ − 856 }
+ − 857 $gene_ids{$c}->insert($fields[3], $fields[1], $fields[2]);
+ − 858 }
+ − 859
+ − 860 # Print output header
+ − 861 open(OUT, ">$output_file")
+ − 862 or die "Cannot open $output_file for writing: $!\n";
+ − 863
+ − 864 print OUT join("\t", "Feature type", "Transcript length", "Selected transcript", "Transcript HGVS", "Strand", "Chr", "DNA From", "DNA To", "Zygosity", "P-value", "Variant Reads", "Total Reads",
+ − 865 "Ref base", "Obs base", "Pop. freq. source", "Pop. freq.", "Variant DB ID"), "\t",
+ − 866 ($internal_snp ? "Internal pop. freq.\t" : ""),
+ − 867 join("\t", "Protein HGVS", "Closest exon junction (AA coding variants)", "Gene Name", "Caveats", "Phase", "Num rare variants in gene (MAF <= $rare_variant_prop)", "Num rare coding and splice site variants in gene (MAF <= $rare_variant_prop)"),"\n";
+ − 868
+ − 869 # If there is CNV data, load it.
+ − 870 # BED columns should be chr start stop caveats ploidy . ignored ignored r,g,b
+ − 871 # The dot means the strand doesn't matter.
+ − 872 # where the first five fields are required, others optional
+ − 873 # where r,g,b is overloaded with father,mother ploidies and "b" is integer representing affected status logical AND (father bit mask 1, mother bit mask 2)
+ − 874 if(defined $cnv_file){
+ − 875 print STDERR "Reading in CNV data...\n" unless $quiet;
+ − 876 open(CNV, $cnv_file)
+ − 877 or die "Cannot open $cnv_file for reading: $!\n";
+ − 878 while(<CNV>){
+ − 879 chomp;
+ − 880 my @F = split /\t/, $_, -1;
+ − 881 if(@F < 5){
+ − 882 print STDERR "Skipping unparseable line ($cnv_file #$.): $_\n";
+ − 883 next;
+ − 884 }
+ − 885 my $ploidy = $F[4];
+ − 886 my $cnv_chr = $F[0];
+ − 887 next if defined $which_chr and $cnv_chr ne $which_chr and "chr$cnv_chr" ne $which_chr and $cnv_chr ne "chr$which_chr";
+ − 888 my $cnv_start = $F[1];
+ − 889 my $cnv_end = $F[2];
+ − 890 my $p_value = "NA";
+ − 891 if($F[3] =~ s/p-value=(\S+?)(?:;|$)//){
+ − 892 $p_value = $1;
+ − 893 next if $min_pvalue < $p_value;
+ − 894 }
+ − 895
+ − 896 # Report a variant line for each gene that is found in this CNV
+ − 897 my $target_parents = $feature_intervaltree{$cnv_chr}->fetch($cnv_start, $cnv_end+1);
+ − 898
+ − 899 my $caveats = "";
+ − 900 if(@F == 9){
+ − 901 my @parents_ploidy = split /,/, $F[8];
+ − 902 if($parents_ploidy[2] == 0){ # neither parent affected
+ − 903 if($ploidy < $parents_ploidy[0] and $ploidy < $parents_ploidy[1]){
+ − 904 if($ploidy > 2){
+ − 905 $caveats = "Polyploidy is less severe than in either unaffected parents";
+ − 906 }
+ − 907 # else: no caveats, this offspring has fewer copies than normally observed, or in unaffected parents
+ − 908 elsif($ploidy < 2){
+ − 909 if($parents_ploidy[0] == 2 and $parents_ploidy[1] == 2){
+ − 910 $caveats = "De novo copy loss, unaffected parents are diploid";
+ − 911 }
+ − 912 else{
+ − 913 $caveats = "Copy loss is greater than in either unaffected parent";
+ − 914 }
+ − 915 }
+ − 916 }
+ − 917 elsif($ploidy >= $parents_ploidy[0] and $ploidy <= $parents_ploidy[1] or
+ − 918 $ploidy >= $parents_ploidy[1] and $ploidy <= $parents_ploidy[0]){
+ − 919 $caveats = "Aneuploidy likely inherited from an unaffected parent";
+ − 920 }
+ − 921 elsif($ploidy > $parents_ploidy[0] and $ploidy > $parents_ploidy[1]){
+ − 922 if($parents_ploidy[0] > 2){
+ − 923 if($parents_ploidy[1] > 2){
+ − 924 $caveats = "Lower polyploidy already exists in both unaffected parents";
+ − 925 }
+ − 926 else{
+ − 927 $caveats = "Lower polyploidy already exists in unaffected father";
+ − 928 }
+ − 929 }
+ − 930 else{
+ − 931 if($parents_ploidy[1] > 2){
+ − 932 $caveats = "Lower polyploidy already exists in unaffected mother";
+ − 933 }
+ − 934 # else no caveats, because both parents are "normal", yet we have polyploidy in the offspring
+ − 935 else{
+ − 936 $caveats = "De novo polyploidy, unaffected parents are diploid";
+ − 937 }
+ − 938 }
+ − 939 }
+ − 940 # else
+ − 941 else{
+ − 942 die "Oops! Error in program logic...how did we get here (unaffected parents)? $_";
+ − 943 }
+ − 944 }
+ − 945 elsif($parents_ploidy[2] == 1){ # father affected
+ − 946 if($ploidy == $parents_ploidy[1]){ # just like unaffected Mom
+ − 947 if($ploidy > 2){
+ − 948 if($ploidy == $parents_ploidy[0]){
+ − 949 $caveats = "Same polyploidy present in both affected and unaffected parents";
+ − 950 }
+ − 951 else{
+ − 952 $caveats = "Polyploidy inherited from unaffected mother";
+ − 953 }
+ − 954 }
+ − 955 elsif($ploidy < 2){
+ − 956 if($ploidy == $parents_ploidy[0]){
+ − 957 $caveats = "Same copy loss in both affected and unaffected parents";
+ − 958 }
+ − 959 else{
+ − 960 $caveats = "Copy loss is shared with unaffected mother";
+ − 961 }
+ − 962 }
+ − 963 else{
+ − 964 if($ploidy == $parents_ploidy[0]){
+ − 965 # Why was this even reported? parents and child have diploid status...
+ − 966 next;
+ − 967 }
+ − 968 $caveats = "Diploidy is shared with unaffected mother";
+ − 969 }
+ − 970 }
+ − 971 elsif($ploidy > 2){ # polyploidy
+ − 972 if($parents_ploidy[0] == 2){
+ − 973 if($parents_ploidy[1] > 2){
+ − 974 $caveats = "Unaffected mother has polyploidy (".$parents_ploidy[1]."x), but affected father is diploid";
+ − 975 }
+ − 976 elsif($parents_ploidy[1] == 2){
+ − 977 $caveats = "Both unaffected mother and affected father are diploid";
+ − 978 }
+ − 979 else{
+ − 980 $caveats = "Affected father is diploid, unaffected mother has copy loss (".$parents_ploidy[1]."x)";
+ − 981 }
+ − 982 }
+ − 983 elsif($parents_ploidy[0] < 2){
+ − 984 $caveats = "Polyploidy found, but affected father had copy loss (".$parents_ploidy[0]."x)";
+ − 985 }
+ − 986 elsif($ploidy < $parents_ploidy[1]){
+ − 987 $caveats = "Polyploidy is less severe than in unaffected mother (".$parents_ploidy[1]."x), or affected father (".$parents_ploidy[0]."x)";
+ − 988 }
+ − 989 # past here the ploidy is great than in the unaffected mother
+ − 990 elsif($parents_ploidy[1] < 2){
+ − 991 $caveats = "Polyploidy is also severe in affected father (".$parents_ploidy[0]."x), but unaffected mother actually had copy loss (". $parents_ploidy[1]. "x)";
+ − 992 }
+ − 993 elsif($parents_ploidy[1] == 2){
+ − 994 $caveats = "Polyploidy is also severe in affected father (".$parents_ploidy[0]."x), and mother is diploid";
+ − 995 }
+ − 996 elsif($ploidy < $parents_ploidy[0]){
+ − 997 $caveats = "Polyploidy is less severe than in affected father (".$parents_ploidy[0]."x), but more severe than unaffected mother (". $parents_ploidy[1]. "x)";
+ − 998 }
+ − 999 elsif($ploidy > $parents_ploidy[0]){
+ − 1000 $caveats = "Polyploidy is more severe than in affected father (".$parents_ploidy[0]."x)";
+ − 1001 }
+ − 1002 else{
+ − 1003 $caveats = "Polyploidy is as severe as in affected father";
+ − 1004 }
+ − 1005 }
+ − 1006 elsif($ploidy == 2){
+ − 1007 # Don't report diploid status, any funny recombination should show up in large indel analysis
+ − 1008 next;
+ − 1009 }
+ − 1010 else{ # copies < 2
+ − 1011 if($ploidy == $parents_ploidy[0]){
+ − 1012 if($ploidy > $parents_ploidy[1]){
+ − 1013 $caveats = "Copy loss is the same as affected father, but less than unaffected mother (". $parents_ploidy[1]. "x)";
+ − 1014 }
+ − 1015 else{
+ − 1016 $caveats = "Copy loss is as severe as in affected father";
+ − 1017 }
+ − 1018 }
+ − 1019 elsif($ploidy > $parents_ploidy[0]){
+ − 1020 if($ploidy > $parents_ploidy[1]){
+ − 1021 if($parents_ploidy[1] == 0 and $parents_ploidy[0] == 0){
+ − 1022 $caveats = "Poor mapping, or Mendelian inheritence violation is severe: no copies of region in either parent, but present in offspring";
+ − 1023 }
+ − 1024 elsif($ploidy == 2){
+ − 1025 next; # child got best of both parents, ignore from CNV standpoint (may still have SNPs of course, or translocation, etc.)
+ − 1026 }
+ − 1027 else{
+ − 1028 $caveats = "Copy loss is less severe than in unaffected mother (".$parents_ploidy[1]."x), or affected father (".$parents_ploidy[0]."x)";
+ − 1029 }
+ − 1030 }
+ − 1031 # else: child has less copies than unaffected mom, but more than affected Dad
+ − 1032 else{
+ − 1033 if($parents_ploidy[1] > 2){
+ − 1034 $caveats = "Copy loss was more severe in affected father (".$parents_ploidy[0]."x), but unaffected mother had polyploidy (".$parents_ploidy[1]."x)";
+ − 1035 }
+ − 1036 elsif($parents_ploidy[1] == 2){
+ − 1037 $caveats = "Copy loss was more severe in affected father (".$parents_ploidy[0]."x), but unaffected mother was diploid";
+ − 1038 }
+ − 1039 else{ # unaffected has loss
+ − 1040 $caveats = "Copy loss is more severe than unaffect mother (".$parents_ploidy[1]."x), but less severe than affected father (".$parents_ploidy[0]."x)";
+ − 1041 }
+ − 1042 }
+ − 1043 }
+ − 1044 # past here, ploidy is less than affected father
+ − 1045 elsif($parents_ploidy[1] > 2){
+ − 1046 $caveats = "Copy loss is more severe than affected father (".$parents_ploidy[0]."x), and unaffected mother had polyploidy (".$parents_ploidy[1]."x)";
+ − 1047 }
+ − 1048 elsif($parents_ploidy[1] == 2){
+ − 1049 $caveats = "Copy loss is more severe than in affected father (".$parents_ploidy[0]."x)";
+ − 1050 }
+ − 1051 else{
+ − 1052 $caveats = "Copy loss is more severe than in both unaffect mother (".$parents_ploidy[1]."x), and affected father (".$parents_ploidy[0]."x)";
+ − 1053 }
+ − 1054 }
+ − 1055 }
+ − 1056 elsif($parents_ploidy[2] == 2){ # mother affected
+ − 1057 if($ploidy == $parents_ploidy[0]){ # just like unaffected Dad
+ − 1058 if($ploidy > 2){
+ − 1059 if($ploidy == $parents_ploidy[1]){
+ − 1060 $caveats = "Same polyploidy present in both affected and unaffected parents";
+ − 1061 }
+ − 1062 else{
+ − 1063 $caveats = "Polyploidy inherited from unaffected father";
+ − 1064 }
+ − 1065 }
+ − 1066 elsif($ploidy < 2){
+ − 1067 if($ploidy == $parents_ploidy[1]){
+ − 1068 $caveats = "Same copy loss in both affected and unaffected parents";
+ − 1069 }
+ − 1070 else{
+ − 1071 $caveats = "Copy loss is shared with unaffected father";
+ − 1072 }
+ − 1073 }
+ − 1074 else{
+ − 1075 if($ploidy == $parents_ploidy[1]){
+ − 1076 # Why was this even reported? parents and child have diploid status...
+ − 1077 next;
+ − 1078 }
+ − 1079 $caveats = "Diploidy is shared with unaffected father";
+ − 1080 }
+ − 1081 }
+ − 1082 elsif($ploidy > 2){ # polyploidy
+ − 1083 if($parents_ploidy[1] == 2){
+ − 1084 if($parents_ploidy[0] > 2){
+ − 1085 $caveats = "Unaffected father has polyploidy (".$parents_ploidy[0]."x), but affected mother is diploid";
+ − 1086 }
+ − 1087 elsif($parents_ploidy[0] == 2){
+ − 1088 $caveats = "Both unaffected father and affected mother are diploid";
+ − 1089 }
+ − 1090 else{
+ − 1091 $caveats = "Affected mother is diploid, unaffected father has copy loss (".$parents_ploidy[1]."x)";
+ − 1092 }
+ − 1093 }
+ − 1094 elsif($parents_ploidy[1] < 2){
+ − 1095 $caveats = "Polyploidy found, but affected mother had copy loss (".$parents_ploidy[1]."x)";
+ − 1096 }
+ − 1097 elsif($ploidy < $parents_ploidy[0]){
+ − 1098 $caveats = "Polyploidy is less severe than in unaffected father (".$parents_ploidy[0]."x), or affected mother (".$parents_ploidy[1]."x)";
+ − 1099 }
+ − 1100 # past here the ploidy is great than in the unaffected father
+ − 1101 elsif($parents_ploidy[0] < 2){
+ − 1102 $caveats = "Polyploidy is also severe in affected mother (".$parents_ploidy[1]."x), but unaffected father actually had copy loss (". $parents_ploidy[0]. "x)";
+ − 1103 }
+ − 1104 elsif($parents_ploidy[0] == 2){
+ − 1105 $caveats = "Polyploidy is also severe in affected mother (".$parents_ploidy[1]."x), and unaffected father is diploid";
+ − 1106 }
+ − 1107 elsif($ploidy < $parents_ploidy[1]){
+ − 1108 $caveats = "Polyploidy is less severe than in affected mother (".$parents_ploidy[1]."x), but more severe than unaffected father (". $parents_ploidy[0]. "x)";
+ − 1109 }
+ − 1110 elsif($ploidy > $parents_ploidy[1]){
+ − 1111 $caveats = "Polyploidy is more severe than in affected mother (".$parents_ploidy[1]."x)";
+ − 1112 }
+ − 1113 else{
+ − 1114 $caveats = "Polyploidy is as severe as in affected mother";
+ − 1115 }
+ − 1116 }
+ − 1117 elsif($ploidy == 2){
+ − 1118 # Don't report diploid status, any funny recombination should show up in large indel analysis
+ − 1119 next;
+ − 1120 }
+ − 1121 else{ # copies < 2
+ − 1122 if($ploidy == $parents_ploidy[1]){
+ − 1123 if($ploidy > $parents_ploidy[0]){
+ − 1124 $caveats = "Copy loss is the same as affected mother, but less than unaffected father (". $parents_ploidy[0]. "x)";
+ − 1125 }
+ − 1126 else{
+ − 1127 $caveats = "Copy loss is as severe as in affected mother";
+ − 1128 }
+ − 1129 }
+ − 1130 elsif($ploidy > $parents_ploidy[1]){
+ − 1131 if($ploidy > $parents_ploidy[0]){
+ − 1132 if($parents_ploidy[1] == 0 and $parents_ploidy[0] == 0){
+ − 1133 $caveats = "Poor mapping, or Mendelian inheritence violation is severe: no copies of region in either parent, but present in offspring";
+ − 1134 }
+ − 1135 elsif($ploidy == 2){
+ − 1136 next; # child got best of both parents, ignore from CNV standpoint (may still have SNPs of course, or translocation, etc.)
+ − 1137 }
+ − 1138 else{
+ − 1139 $caveats = "Copy loss is less severe than in unaffected father (".$parents_ploidy[0]."x), or affected mother (".$parents_ploidy[1]."x)";
+ − 1140 }
+ − 1141 }
+ − 1142 # else: child has less copies than unaffected Dad, but more than affected Mom
+ − 1143 else{
+ − 1144 if($parents_ploidy[0] > 2){
+ − 1145 $caveats = "Copy loss was more severe in affected mother (".$parents_ploidy[1]."x), but unaffected father had polyploidy (".$parents_ploidy[0]."x)";
+ − 1146 }
+ − 1147 elsif($parents_ploidy[0] == 2){
+ − 1148 $caveats = "Copy loss was more severe in affected mother (".$parents_ploidy[1]."x), but unaffected father was diploid";
+ − 1149 }
+ − 1150 else{ # unaffected has loss
+ − 1151 $caveats = "Copy loss is more severe than unaffect father (".$parents_ploidy[0]."x), but less severe than affected mother (".$parents_ploidy[1]."x)";
+ − 1152 }
+ − 1153 }
+ − 1154 }
+ − 1155 # past here, ploidy is less than affected mother
+ − 1156 elsif($parents_ploidy[0] > 2){
+ − 1157 $caveats = "Copy loss is more severe than affected mother (".$parents_ploidy[1]."x), and unaffected father had polyploidy (".$parents_ploidy[0]."x)";
+ − 1158 }
+ − 1159 elsif($parents_ploidy[0] == 2){
+ − 1160 $caveats = "Copy loss is more severe than in affected mother (".$parents_ploidy[1]."x)";
+ − 1161 }
+ − 1162 else{
+ − 1163 $caveats = "Copy loss is more severe than in both unaffect father (".$parents_ploidy[0]."x), and affected mother (".$parents_ploidy[1]."x)";
+ − 1164 }
+ − 1165 }
+ − 1166 }
+ − 1167
+ − 1168 }
+ − 1169 if($F[3] and $F[3] ne "-"){ # prexisting caveat from CNV caller
+ − 1170 if(defined $caveats){
+ − 1171 $caveats .= "; $F[3]" unless $caveats =~ /\b$F[3]\b/;
+ − 1172 }
+ − 1173 else{
+ − 1174 $caveats = $F[3];
+ − 1175 }
+ − 1176 }
+ − 1177
+ − 1178 # Sort by start for consistency
+ − 1179 my @target_parents = sort {$feature_range{$cnv_chr}->{$a}->[0]->[0] <=> $feature_range{$cnv_chr}->{$b}->[0]->[0]} @$target_parents;
+ − 1180
+ − 1181 for my $target_parent (@target_parents){
+ − 1182 my $target_caveats = $caveats;
+ − 1183 my $strand = $feature_strand{$target_parent};
+ − 1184 # report the gain/loss of each gene separately, for simplicity in downstream analysis
+ − 1185 my $cnv_exon_start = 10000000000; # genomic coords
+ − 1186 my $cnv_exon_end = 0;
+ − 1187 my $cnv_cdna_start = 0; # cDNA coords
+ − 1188 my $cnv_cdna_end = 0;
+ − 1189 my $off5 = 0; # border outside the exon?
+ − 1190 my $off3 = 0;
+ − 1191 my @feature_ranges = @{$feature_range{$cnv_chr}->{$target_parent}};
+ − 1192 # find the first and last exons in the gene that are inside the CNV
+ − 1193 for my $subregion (@feature_ranges){
+ − 1194 # exon overlaps CNV?
+ − 1195 if($subregion->[0] <= $cnv_end and $subregion->[1] >= $cnv_start){
+ − 1196 if($cnv_exon_start > $subregion->[0]){
+ − 1197 if($cnv_start < $subregion->[0]){
+ − 1198 $cnv_exon_start = $subregion->[0]; $off5 = 1;
+ − 1199 $cnv_cdna_start = $subregion->[2];
+ − 1200 }
+ − 1201 else{
+ − 1202 $cnv_exon_start = $cnv_start; $off5 = 0;
+ − 1203 $cnv_cdna_start = $subregion->[2]+($strand eq "-" ? $subregion->[0]-$cnv_start: $cnv_start-$subregion->[0]);
+ − 1204 }
+ − 1205 }
+ − 1206 if($cnv_exon_end < $subregion->[1]){
+ − 1207 if($cnv_end > $subregion->[1]){
+ − 1208 $cnv_exon_end = $subregion->[1]; $off3 = 1;
+ − 1209 $cnv_cdna_end = $subregion->[2]+($strand eq "-" ? $subregion->[0]-$subregion->[1] : $subregion->[1]-$subregion->[0]);
+ − 1210 }
+ − 1211 else{
+ − 1212 $cnv_exon_end = $cnv_end; $off3 = 0;
+ − 1213 $cnv_cdna_end = $subregion->[2]+($strand eq "-" ? $subregion->[0]-$cnv_end : $cnv_end-$subregion->[0]);
+ − 1214 }
+ − 1215 }
+ − 1216 }
+ − 1217 }
+ − 1218
+ − 1219 my $ends_internally = 0;
+ − 1220 if($cnv_exon_end == 0){ # ends inside the exon
+ − 1221 $cnv_exon_end = $cnv_end;
+ − 1222 $ends_internally = 1;
+ − 1223 }
+ − 1224 # See if it's in the structural variant database
+ − 1225 my @gain_coverage; $#gain_coverage = $cnv_exon_end-$cnv_exon_start; # preallocate blanks
+ − 1226 my @loss_coverage; $#loss_coverage = $cnv_exon_end-$cnv_exon_start; # preallocate blanks
+ − 1227 my $dgv_loss_id; # report the DGV entry that covers most of the observed structural variant
+ − 1228 my $dgv_loss_length = 0; # report the DGV entry that covers most of the observed structural variant
+ − 1229 my $dgv_gain_id; # report the DGV entry that covers most of the observed structural variant
+ − 1230 my $dgv_gain_length = 0; # report the DGV entry that covers most of the observed structural variant
+ − 1231 my $gains;
+ − 1232 my $losses;
+ − 1233 my $dgv_chr = $cnv_chr;
+ − 1234 $dgv_chr =~ s/^chr//; # no prefix in DGV
+ − 1235 #open(DGV, "tabix $dgv_file $dgv_chr:$cnv_exon_start-$cnv_exon_end |") # check out CNV in this gene model region
+ − 1236 # or die "Cannot run tabix: $!\n";
+ − 1237 open(DGV, "/dev/null");
+ − 1238 while(<DGV>){
+ − 1239 my @C = split /\t/, $_;
+ − 1240 next if $C[4] ne "CNV"; # todo: handle indels?
+ − 1241 my $dgv_start = $C[2];
+ − 1242 my $dgv_end = $C[3];
+ − 1243 my $dgv_direction = $C[5];
+ − 1244 my $gain_cov_count = 0;
+ − 1245 my $loss_cov_count = 0;
+ − 1246 if($dgv_direction eq "Gain"){
+ − 1247 for(my $i = ($dgv_start < $cnv_exon_start ? $cnv_exon_start : $dgv_start); $i <= $dgv_end and $i <= $cnv_exon_end; $i++){
+ − 1248 $gain_coverage[$i-$cnv_exon_start] = 1 unless defined $gain_coverage[$i-$cnv_exon_start];
+ − 1249 $gain_cov_count++;
+ − 1250 }
+ − 1251 }
+ − 1252 elsif($dgv_direction eq "Loss"){
+ − 1253 for(my $i = ($dgv_start < $cnv_exon_start ? $cnv_exon_start : $dgv_start); $i <= $dgv_end and $i <= $cnv_exon_end; $i++){
+ − 1254 $loss_coverage[$i-$cnv_exon_start] = 1 unless defined $loss_coverage[$i-$cnv_exon_start];
+ − 1255 $loss_cov_count++;
+ − 1256 }
+ − 1257 }
+ − 1258 if($dgv_direction eq "Gain" and $gain_cov_count > $dgv_gain_length){
+ − 1259 $dgv_gain_id = $C[0];
+ − 1260 $dgv_gain_length = $gain_cov_count;
+ − 1261 }
+ − 1262 if($dgv_direction eq "Loss" and $loss_cov_count > $dgv_loss_length){
+ − 1263 $dgv_loss_id = $C[0];
+ − 1264 $dgv_loss_length = $loss_cov_count;
+ − 1265 }
+ − 1266 }
+ − 1267 close(DGV);
+ − 1268
+ − 1269 my $gain_coverage = 0;
+ − 1270 for my $count (@gain_coverage){
+ − 1271 $gain_coverage++ if defined $count;
+ − 1272 }
+ − 1273 $gain_coverage = sprintf "%.3f", $gain_coverage/($cnv_exon_end-$cnv_exon_start+1); # make it a proportion
+ − 1274 my $loss_coverage = 0;
+ − 1275 for my $count (@loss_coverage){
+ − 1276 $loss_coverage++ if defined $count;
+ − 1277 }
+ − 1278 $loss_coverage = sprintf "%.3f", $loss_coverage/($cnv_exon_end-$cnv_exon_start+1); # make it a proportion
+ − 1279
+ − 1280 my $src = "DGV";
+ − 1281 my $dgv_id = "NA";
+ − 1282 my $dgv_caveat;
+ − 1283 my $dgv_coverage;
+ − 1284 if($ploidy > 2){
+ − 1285 if(not defined $dgv_gain_id){
+ − 1286 if(defined $dgv_loss_id){
+ − 1287 $dgv_id = sprintf "%s/%.3f", $dgv_loss_id, $dgv_loss_length/($cnv_exon_end-$cnv_exon_start+1);
+ − 1288 $dgv_caveat = "; No gains are known in healthy controls, the DGV2 ID reported is for a loss in the same area";
+ − 1289 $dgv_coverage = $loss_coverage;
+ − 1290 }
+ − 1291 else{
+ − 1292 $dgv_id = "novel";
+ − 1293 $dgv_coverage = "NA";
+ − 1294 $src = "NA";
+ − 1295 }
+ − 1296 }
+ − 1297 else{
+ − 1298 $dgv_id = sprintf "%s/%.3f", $dgv_gain_id, $dgv_gain_length/($cnv_exon_end-$cnv_exon_start+1);
+ − 1299 $dgv_coverage = $gain_coverage;
+ − 1300 }
+ − 1301 }
+ − 1302 elsif($ploidy < 2){
+ − 1303 if(not defined $dgv_loss_id){
+ − 1304 if(defined $dgv_gain_id){
+ − 1305 $dgv_id = sprintf "%s/%.3f", $dgv_gain_id, $dgv_gain_length/($cnv_exon_end-$cnv_exon_start+1);
+ − 1306 $dgv_caveat = "; No losses are known in healthy controls, the DGV2 ID reported is for a gain in the same area";
+ − 1307 $dgv_coverage = $gain_coverage;
+ − 1308 }
+ − 1309 else{
+ − 1310 $dgv_id = "novel";
+ − 1311 $dgv_coverage = "NA";
+ − 1312 $src = "NA";
+ − 1313 }
+ − 1314 }
+ − 1315 else{
+ − 1316 $dgv_id = sprintf "%s/%.3f", $dgv_loss_id, $dgv_loss_length/($cnv_exon_end-$cnv_exon_start+1);
+ − 1317 $dgv_coverage = $loss_coverage;
+ − 1318 }
+ − 1319 }
+ − 1320
+ − 1321 my $non_coding = 0;
+ − 1322 if(not exists $feature_cds_max{$target_parent} or not defined $feature_cds_max{$target_parent} or $feature_cds_max{$target_parent} eq ""){
+ − 1323 $non_coding = 1;
+ − 1324 }
+ − 1325 $target_caveats .= $dgv_caveat if defined $dgv_caveat and $dgv_id ne "novel" and $target_caveats !~ /\Q$dgv_caveat\E/;
+ − 1326 #print "Recorded $cnv_chr:$cnv_start caveat $caveats\n";
+ − 1327 # if it doesn't overlap an exon, we need to find out which two exons it's between
+ − 1328 if($ends_internally){
+ − 1329 my $intron_found = 0;
+ − 1330 for(my $i = 0; $i < $#feature_ranges; $i++){
+ − 1331 if($feature_ranges[$i]->[1] < $cnv_start and $feature_ranges[$i+1]->[0] > $cnv_end){
+ − 1332 if($ploidy > 2){ # gain
+ − 1333 if($strand eq "-"){
+ − 1334 record_snv("$target_parent\t",
+ − 1335 ($non_coding ? "g.$cnv_start\_$cnv_end" :
+ − 1336 "c.".($feature_ranges[$i+1]->[2])."+".($feature_ranges[$i+1]->[0]-$cnv_end)."_".($feature_ranges[$i+1]->[2]+1)."-".($cnv_start-$feature_ranges[$i]->[1])),
+ − 1337 # pos Zygosity P-value Variant Reads Total Reads Ref Bases Var Bases Population Frequency Source Pop. freq. or DGV2 gain/loss coverage dbSNP or DGV2 ID
+ − 1338 "[".($ploidy-1)."]\t$strand\t$cnv_chr\t$cnv_start\t$cnv_end\tNA\t$p_value\tNA\tNA\t",
+ − 1339 "NA\tNA\t$src\t$dgv_coverage\t$dgv_id\tNA\tNA\t".range2genes($cnv_chr,$cnv_start,$cnv_end)."\t$target_caveats\t\n");
+ − 1340 }
+ − 1341 else{
+ − 1342 record_snv("$target_parent\t",
+ − 1343 ($non_coding ? "g.$cnv_start\_$cnv_end" :
+ − 1344 "c.".($feature_ranges[$i+1]->[2]-1)."+".($cnv_start-$feature_ranges[$i]->[1])."_".$feature_ranges[$i+1]->[2]."-".($feature_ranges[$i+1]->[0]-$cnv_end)),
+ − 1345 "[".($ploidy-1)."]\t$strand\t$cnv_chr\t$cnv_start\t$cnv_end\tNA\t$p_value\tNA\tNA\t",
+ − 1346 "NA\tNA\t$src\t$dgv_coverage\t$dgv_id\tNA\tNA\t".range2genes($cnv_chr,$cnv_start,$cnv_end)."\t$target_caveats\t\n");
+ − 1347 }
+ − 1348 }
+ − 1349 else{ # loss
+ − 1350 if($strand eq "-"){
+ − 1351 record_snv("$target_parent\t",
+ − 1352 ($non_coding ? "g.$cnv_start\_$cnv_end" :
+ − 1353 "c.".($feature_ranges[$i+1]->[2])."+".($feature_ranges[$i+1]->[0]-$cnv_end)."_".($feature_ranges[$i+1]->[2]+1)."-".($cnv_start-$feature_ranges[$i]->[1])),
+ − 1354 "del\t$strand\t$cnv_chr\t$cnv_start\t$cnv_end\t", ($ploidy == 1 ? "heterozygote" : "homozygote"), "\t$p_value\tNA\tNA\t",
+ − 1355 "NA\tNA\t$src\t$dgv_coverage\t$dgv_id\tNA\tNA\t".range2genes($cnv_chr,$cnv_start,$cnv_end)."\t$target_caveats\t\n");
+ − 1356 }
+ − 1357 else{
+ − 1358 record_snv("$target_parent\t",
+ − 1359 ($non_coding ? "g.$cnv_start\_$cnv_end" :
+ − 1360 "c.".($feature_ranges[$i+1]->[2]-1)."+".($cnv_start-$feature_ranges[$i]->[1])."_".$feature_ranges[$i+1]->[2]."-".($feature_ranges[$i+1]->[0]-$cnv_end)),
+ − 1361 "del\t$strand\t$cnv_chr\t$cnv_start\t$cnv_end\t", ($ploidy == 1 ? "heterozygote" : "homozygote"), "\t$p_value\tNA\tNA\t",
+ − 1362 "NA\tNA\t$src\t$dgv_coverage\t$dgv_id\tNA\tNA\t".range2genes($cnv_chr,$cnv_start,$cnv_end)."\t$target_caveats\t\n");
+ − 1363 }
+ − 1364 }
+ − 1365 $intron_found = 1; last;
+ − 1366 }
+ − 1367 }
+ − 1368 warn "Logic error: CNV overlaps a gene ($target_parent), but is neither intronic nor exonic. Offending CNV was $_\n" unless $intron_found;
+ − 1369 next;
+ − 1370 }
+ − 1371 if($strand eq "-"){
+ − 1372 my $tmp = $cnv_cdna_start;
+ − 1373 $cnv_cdna_start = $cnv_cdna_end;
+ − 1374 $cnv_cdna_end = $tmp;
+ − 1375 }
+ − 1376 # Make the location label pretty descriptive
+ − 1377 my $cnv_phase = "";
+ − 1378 if($cnv_exon_start > $cnv_start or $cnv_exon_end < $cnv_end){
+ − 1379 $cnv_phase = "CNV-$cnv_chr:$cnv_start-$cnv_end"; # Use phase to note that it's part of a bigger CNV than just the range of this feature
+ − 1380 }
+ − 1381 # if we get here, we're in a gained/deleted exon category
+ − 1382 # CNVs are fuzzy-edged (as opposed to large indels), so produce HGVS syntax that reflect this
+ − 1383 if($ploidy > 2){ # gain
+ − 1384 # find the exons encompassed by the CNV. NOTE that we assume that polyploidies are proximal
+ − 1385 record_snv("$target_parent\t",
+ − 1386 ($non_coding ? "g.".($cnv_exon_start > $cnv_start ? "$cnv_exon_start-?" : $cnv_start)."_".($cnv_exon_end < $cnv_end ? "$cnv_exon_end+?" : $cnv_end) :
+ − 1387 "c.$cnv_cdna_start".($off5?"-?":"")."_$cnv_cdna_end".($off3?"+?":"")),
+ − 1388 "[".($ploidy-1)."]\t$strand\t$cnv_chr\t$cnv_exon_start\t$cnv_exon_end\tNA\t$p_value\tNA\tNA\t",
+ − 1389 "NA\tNA\t$src\t$dgv_coverage\t$dgv_id\tNA\tNA\t".range2genes($cnv_chr,$cnv_start,$cnv_end)."\t$target_caveats\t$cnv_phase\n");
+ − 1390 }
+ − 1391 else{ # loss
+ − 1392 #translate genome coordinates into cDNA coordinates
+ − 1393 record_snv("$target_parent\t",
+ − 1394 ($non_coding ? "g.".($cnv_exon_start > $cnv_start ? "$cnv_exon_start-?" : $cnv_start)."_".($cnv_exon_end < $cnv_end ? "$cnv_exon_end+?" : $cnv_end) :
+ − 1395 "c.$cnv_cdna_start".($off5?"-?":"")."_$cnv_cdna_end".($off3?"+?":"")),
+ − 1396 "del\t$strand\t$cnv_chr\t$cnv_exon_start\t$cnv_exon_end\t", ($ploidy == 1 ? "heterozygote" : "homozygote"), "\t$p_value\tNA\tNA\t",
+ − 1397 "NA\tNA\t$src\t$dgv_coverage\t$dgv_id\tNA\tNA\t".range2genes($cnv_chr,$cnv_start,$cnv_end)."\t$target_caveats\t$cnv_phase\n");
+ − 1398 }
+ − 1399 }
+ − 1400 }
+ − 1401 close(CNV);
+ − 1402
+ − 1403 }
+ − 1404
+ − 1405
+ − 1406 #sort genes by start, then longest if tied
+ − 1407 my %rc = qw(A T T A G C C G N N S W W S K M M K Y R R Y X X);
+ − 1408 print STDERR "Processing variant calls..." unless $quiet;
+ − 1409 %chr_read = ();
+ − 1410 open(VCFIN, $input_file)
+ − 1411 or die "Cannot open $input_file for reading: $!\n";
+ − 1412 while(<VCFIN>){
+ − 1413 if(/^\s*(?:#|$)/){ # blank or hash comment lines
+ − 1414 next;
+ − 1415 }
+ − 1416 chomp;
+ − 1417 my @fields = split /\t/, $_;
+ − 1418
+ − 1419 next unless exists $feature_range{$fields[0]};
+ − 1420 if(not $quiet and not exists $chr_read{$fields[0]}){
+ − 1421 print STDERR " $fields[0]";
+ − 1422 $chr_read{$fields[0]} = 1;
+ − 1423 #print STDERR "(not in reference file!)" unless exists $feature_range{$fields[0]};
+ − 1424 }
+ − 1425
+ − 1426 next if $fields[4] eq "<NON_REF>"; # GVCF background stuff
+ − 1427 next if $fields[9] eq "./." or $fields[9] eq "."; # no call
+ − 1428 my $chr = $fields[0];
+ − 1429 next if defined $which_chr and $chr ne $which_chr and $chr ne "chr$which_chr" and "chr$chr" ne $which_chr;
+ − 1430 print STDERR "passes chr and field # test" if grep /dataset_7684.dat/, @ARGV;
+ − 1431 $chr = "chr$chr" if $chr !~ /^chr/;
+ − 1432 $fields[8] =~ s/\s+$//;
+ − 1433 my @values = split /:/, $fields[9];
+ − 1434 #print STDERR join(" / ", @values), "\n" if $. == 84;
+ − 1435 my @keys = split /:/, $fields[8];
+ − 1436 my $zygosity = "n/a";
+ − 1437 my $quality = "n/a";
+ − 1438 my $read_depth = "n/a";
+ − 1439 my $variant_depths = "n/a";
+ − 1440 for(my $i = 0; $i <= $#keys and $i <= $#values; $i++){ # values max index check because some genotypers are nasty and don't provide as many fields as they say they will
+ − 1441 if($keys[$i] eq "GT"){
+ − 1442 if($values[$i] =~ /\./ or $values[$i] =~ /0\/0/){ # one genotype not described
+ − 1443 $zygosity = "none";
+ − 1444 last;
+ − 1445 }
+ − 1446 else{ # genotypes described
+ − 1447 $zygosity = $values[$i] =~ /[02]/ ? "heterozygote" : "homozygote";
+ − 1448 }
+ − 1449 }
+ − 1450 elsif($keys[$i] eq "DNM_CONFIG" and $zygosity eq "n/a"){ # hack
+ − 1451 $zygosity = $values[$i] =~ /^(.)\1/ ? "homozygote" : "heterozygote";
+ − 1452 }
+ − 1453 elsif($keys[$i] eq "GQ" and $values[$i] ne "."){
+ − 1454 #print "Checking GQ (index $i) $values[$i] gq2p\n" if $. == 84;
+ − 1455 $quality = gq2p($values[$i]);
+ − 1456 }
+ − 1457 elsif($keys[$i] eq "RD"){ # from some tools like denovo variant finders
+ − 1458 $read_depth = $values[$i];
+ − 1459 }
+ − 1460 elsif($keys[$i] eq "DP"){
+ − 1461 $read_depth = $values[$i];
+ − 1462 }
+ − 1463 # the frequency of the variant can go by many names, here we have freebayes (A* are new and old versions) and atlas2_indel
+ − 1464 elsif($keys[$i] eq "AA" or $keys[$i] eq "VR" or $keys[$i] eq "AO"){
+ − 1465 $variant_depths = $values[$i];
+ − 1466 }
+ − 1467 # here we have GATK variant freq of form ref#,var#
+ − 1468 elsif($keys[$i] eq "AD"){
+ − 1469 $variant_depths = $values[$i];
+ − 1470 $variant_depths =~ s/^\d+,//;
+ − 1471 }
+ − 1472 else{
+ − 1473 #print STDERR "Ignoring field $keys[$i]\n";
+ − 1474 }
+ − 1475 }
+ − 1476 next if $zygosity eq "none"; # GVCF non-ref for example or when multiple samples are reported simultaneously
+ − 1477 $quality = z2p($1) if $fields[7] =~ /Z=(\d+\.\d+)/;
+ − 1478 if($quality eq "n/a" and $fields[5] ne "."){
+ − 1479 $quality = gq2p($fields[5]);
+ − 1480 }
+ − 1481 if($fields[5] eq "0" and $fields[6] eq "PASS"){ # not qual and a PASS in the filter column
+ − 1482 $quality = 1;
+ − 1483 }
+ − 1484 elsif($quality ne "n/a" and $quality > $min_pvalue){ # p-value cutoff
+ − 1485 #print "Checking call quality $fields[5] gq2p\n" if $. == 84;
+ − 1486 next unless gq2p($fields[5]) <= $min_pvalue; # in some cases, programs like FreeBayes give low genotype quality such as when a call is borderline homo/het
+ − 1487 # in these cases it would be silly to reject the call if their are many supporting reads.
+ − 1488 }
+ − 1489
+ − 1490 # Some tools like GATK don't report number of variant reads...infer from other data if possible
+ − 1491 if($variant_depths eq "n/a"){
+ − 1492 my @key_value_pairs = split /;/, $fields[7];
+ − 1493 for my $key_value_pair (@key_value_pairs){
+ − 1494 if($key_value_pair !~ /^(.*?)=(.*)$/){
+ − 1495 next;
+ − 1496 #next if $key_value_pair eq "INDEL"; # samtools peculiarity
+ − 1497 #die "Key-value pair field (column #8) does not have the format key=value for entry $key_value_pair (line #$. of ), please fix the VCF file\n";
+ − 1498 }
+ − 1499 if($1 eq "AB"){ # GATK: for het calls, AB is ref/(ref+var), only one variant reported per line
+ − 1500 $variant_depths = "";
+ − 1501 for my $ab (split /,/, $2){
+ − 1502 $variant_depths .= int((1-$ab)*$read_depth).",";
+ − 1503 }
+ − 1504 chop $variant_depths;
+ − 1505 }
+ − 1506 elsif($1 eq "MLEAC"){
+ − 1507 }
+ − 1508 elsif($1 eq "DP4"){ # samtools: high-quality ref-forward bases, ref-reverse, alt-forward and alt-reverse bases
+ − 1509 my @rds = split /,/, $2;
+ − 1510 $variant_depths = $rds[2]+$rds[3];
+ − 1511 $read_depth = $rds[0]+$rds[1]+$variant_depths;
+ − 1512 if($fields[4] =~ /,/){ # samtools doesn't break down compound het calls into individual frequencies
+ − 1513 my $num_alt_genotypes = $fields[4] =~ tr/,/,/;
+ − 1514 $num_alt_genotypes++;
+ − 1515 my $even_prop = sprintf "%.2f", $variant_depths/$read_depth/$num_alt_genotypes;
+ − 1516 $variant_depths = join(",", ($even_prop) x $num_alt_genotypes);
+ − 1517 if(not exists $chr2caveats{"$chr:$fields[1]"}){
+ − 1518 $chr2caveats{"$chr:$fields[1]"} = "compound het var freq n/a in orig VCF file, auto set to $even_prop";
+ − 1519 }
+ − 1520 else{
+ − 1521 $chr2caveats{"$chr:$fields[1]"} .= "; compound het var freq n/a in orig VCF file, auto set to $even_prop";
+ − 1522 }
+ − 1523 }
+ − 1524 }
+ − 1525 }
+ − 1526 }
+ − 1527 if($variant_depths eq "n/a"){ # usually homo var calls, can only assume all reads are variant
+ − 1528 if($zygosity eq "homozygote"){
+ − 1529 $variant_depths = $read_depth;
+ − 1530 if($read_depth ne "n/a"){
+ − 1531 if(not exists $chr2caveats{"$chr:$fields[1]"}){
+ − 1532 $chr2caveats{"$chr:$fields[1]"} = "homo var freq n/a in orig VCF file, auto set to 1.0";
+ − 1533 }
+ − 1534 else{
+ − 1535 $chr2caveats{"$chr:$fields[1]"} = "; homo var freq n/a in orig VCF file, auto set to 1.0";
+ − 1536 }
+ − 1537 }
+ − 1538 }
+ − 1539 else{
+ − 1540 if($read_depth ne "n/a"){
+ − 1541 $variant_depths = int($read_depth/2);
+ − 1542 if(not exists $chr2caveats{"$chr:$fields[1]"}){
+ − 1543 $chr2caveats{"$chr:$fields[1]"} = "het var freq n/a in orig VCF file, auto set to 0.5";
+ − 1544 }
+ − 1545 else{
+ − 1546 $chr2caveats{"$chr:$fields[1]"} = "; het var freq n/a in orig VCF file, auto set to 0.5";
+ − 1547 }
+ − 1548 }
+ − 1549 else{
+ − 1550 $variant_depths = $read_depth;
+ − 1551 }
+ − 1552 }
+ − 1553 }
+ − 1554
+ − 1555 my $target_parents = $feature_intervaltree{$chr}->fetch($fields[1]-$flanking_bases, $fields[1]+length($fields[3])+$flanking_bases);
+ − 1556 # Last ditch, if not in a gene model, is it at least in an enrichment region?
+ − 1557 if(not @$target_parents){
+ − 1558 next if not exists $enrichment_regions{$chr};
+ − 1559 my $regions_ref = $enrichment_regions{$chr};
+ − 1560 my $location = $fields[1];
+ − 1561 my $strand = "+"; # for lack of a better choice
+ − 1562 for(my $i = find_earliest_index($location-$flanking_bases, $regions_ref);
+ − 1563 $i < $#$regions_ref and $location <= $regions_ref->[$i]->[1]+$flanking_bases;
+ − 1564 $i++){
+ − 1565 next unless $regions_ref->[$i]->[0]-$flanking_bases <= $location and $regions_ref->[$i]->[1]+$flanking_bases >= $location;
+ − 1566
+ − 1567 my $feature_name = "enrichment_target_$chr:".$regions_ref->[$i]->[0]."-".$regions_ref->[$i]->[1];
+ − 1568 $feature_type{$feature_name} = "misc_enrichment_kit_target";
+ − 1569 $feature_length{$feature_name} = $regions_ref->[$i]->[1]-$regions_ref->[$i]->[0]+1;
+ − 1570 my @variants = split /,/, $fields[4];
+ − 1571 my @variant_depths = split /,/, $variant_depths;
+ − 1572 my $ref = uc($fields[3]);
+ − 1573 for(my $j = 0; $j <= $#variants; $j++){
+ − 1574 my $variant = $variants[$j];
+ − 1575 next if $variant eq "<NON_REF>" or $variant_depths[$j] eq "0"; # GVCF stuff
+ − 1576 my $variant_depth = $variant_depths[$j];
+ − 1577 if($min_prop){
+ − 1578 next unless $variant_depth >= $min_depth and $read_depth ne "n/a" and $variant_depth/$read_depth >= $min_prop;
+ − 1579 }
+ − 1580 if(length($ref) == 1 and length($variant) == 1){ # SNP
+ − 1581 record_snv("$feature_name\tg.$location",
+ − 1582 "$ref>$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 1583 join("\t",prop_info_key($chr,$location,$ref,$variant)),"\tNA\n");
+ − 1584 }
+ − 1585 elsif(length($ref) == 1 and length($variant) > 1){ # insertion
+ − 1586 record_snv("$feature_name\tg.$location\_",($location+1),
+ − 1587 "ins",substr($variant, 1),"\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 1588 join("\t",prop_info_key($chr,$location,$ref,$variant)),"\tNA\n");
+ − 1589 }
+ − 1590 elsif(length($variant) == 1 and length($ref) > 1){ # deletion
+ − 1591 record_snv("$feature_name\tg.$location\_",($location+length($ref)-1),
+ − 1592 "del",substr($ref, 1),"\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 1593 join("\t",prop_info_key($chr,$location,$ref,$variant)),"\tNA\n");
+ − 1594 }
+ − 1595 else{ # indel
+ − 1596 record_snv("$feature_name\tg.$location\_",($location+length($ref)-1),
+ − 1597 "delins$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 1598 join("\t",prop_info_key($chr,$location,$ref,$variant)),"\tNA\n");
+ − 1599 }
+ − 1600 } # end for variants
+ − 1601 next; # process next record, we've done all we can with a non-coding-gene SNP
+ − 1602 }
+ − 1603 }
+ − 1604
+ − 1605 for my $target_parent (@$target_parents){
+ − 1606
+ − 1607 print STDERR "Checking parent $target_parent for on $chr:$fields[1] $fields[3] -> $fields[4]\n" if grep /dataset_7684.dat/, @ARGV;
+ − 1608 my @feature_ranges = @{$feature_range{$chr}->{$target_parent}};
+ − 1609 # Calculate the position of the change within the feature range position
+ − 1610 my $strand = $feature_strand{$target_parent};
+ − 1611 my $trans_table = exists $feature_transl_table{$target_parent} ? $feature_transl_table{$target_parent} : $default_transl_table;
+ − 1612 $fields[4]=~tr/"//d; # sometime strangely surroundsed by quotes
+ − 1613 my @variants = split /,/, $fields[4];
+ − 1614 my @variant_depths = split /,/, $variant_depths;
+ − 1615 my @refs = (uc($fields[3])) x scalar(@variants);
+ − 1616 my @locations = ($fields[1]) x scalar(@variants);
+ − 1617
+ − 1618 for(my $j = 0; $j <= $#variants; $j++){
+ − 1619 my $ref = $refs[$j];
+ − 1620 my $location = $locations[$j];
+ − 1621 my $feature_offset = 0;
+ − 1622 my $feature_num = 0;
+ − 1623 my $variant = uc($variants[$j]);
+ − 1624 next if $variant eq "<NON_REF>" or $variant_depths[$j] eq "0"; # GVCF stuff
+ − 1625 my $variant_depth = $variant_depths[$j] || "n/a";
+ − 1626 #print STDERR "Evaluating target parent $target_parent for $chr:$fields[1]-".($fields[1]+length($fields[3]))." -> ",join("/", @$target_parents) , "\n" if $fields[1] == 982994;
+ − 1627
+ − 1628 # Break down MNPs into individual SNPs that are phased (TODO: skip if it's an inversion? would require amalgamating SNPs for tools that call them individually, phased :-P)
+ − 1629 if(length($variant) > 1 and length($variant) == length($ref)){
+ − 1630 my @subvariants;
+ − 1631 my @subrefs;
+ − 1632 my @sublocations;
+ − 1633 my $phase_range = $location."-".($location+length($ref)-1);
+ − 1634 for(my $k = 0; $k < length($variant); $k++){
+ − 1635 my $r = substr($ref, $k, 1);
+ − 1636 my $v = substr($variant, $k, 1);
+ − 1637 if($r ne $v){
+ − 1638 push @subvariants, $v;
+ − 1639 push @subrefs, $r;
+ − 1640 push @sublocations, $location+$k;
+ − 1641 }
+ − 1642 elsif(@variants == 1){
+ − 1643 next; # homo ref call
+ − 1644 }
+ − 1645 if($zygosity eq "heterozygote"){
+ − 1646 # need to ignore case where a homozygous call (variant or ref) is included in a double non-ref het MNP
+ − 1647 if(@variants > 1){
+ − 1648 my $homo = 1;
+ − 1649 for(my $l = 0; $l <= $#variants; $l++){ # using loop in case we handle oligoploid genomes in the future
+ − 1650 if(length($variants[$l]) <= $k or substr($variants[$l], $k, 1) ne $v){
+ − 1651 $homo = 0;
+ − 1652 last;
+ − 1653 }
+ − 1654 }
+ − 1655 next if $homo;
+ − 1656 }
+ − 1657 my $phase_key = "$chr:".($location+$k).":$v";
+ − 1658 my $phase_label = "M$j-$chr:$phase_range";
+ − 1659 if(exists $chr2phase{$phase_key}){
+ − 1660 if($chr2phase{$phase_key} !~ /$phase_label/){
+ − 1661 $chr2phase{$phase_key} .= ",$phase_label";
+ − 1662 }
+ − 1663 }
+ − 1664 else{
+ − 1665 $chr2phase{$phase_key} = $phase_label;
+ − 1666 }
+ − 1667 }
+ − 1668 }
+ − 1669 # recycle this MNP variant loop to start processing the individual SNPs
+ − 1670 splice(@refs, $j, 1, @subrefs);
+ − 1671 splice(@variants, $j, 1, @subvariants);
+ − 1672 splice(@locations, $j, 1, @sublocations);
+ − 1673 splice(@variant_depths, $j, 1, ($variant_depth) x scalar(@subvariants));
+ − 1674 $j--;
+ − 1675 next;
+ − 1676 }
+ − 1677
+ − 1678 if($min_prop != 0 and $variant_depth eq "n/a" or $variant_depth eq "."){
+ − 1679 print STDERR "Could not parse variant depth from $_\n" unless $quiet;
+ − 1680 next;
+ − 1681 }
+ − 1682 next unless $min_prop == 0 or $min_prop and $variant_depth >= $min_depth and $read_depth ne "n/a" and $variant_depth/$read_depth >= $min_prop;
+ − 1683 if($zygosity eq "NA"){
+ − 1684 # make the call ourselves
+ − 1685 $zygosity = $variant_depths/$read_depth > 1-$min_prop ? "homozygote" : "heterozygote";
+ − 1686 }
+ − 1687 # e.g. chr22 47857671 . CAAAG AAGAT,AAAAG 29.04 .
+ − 1688 if(length($variant) > 1 and length($variant) == length($ref)){
+ − 1689 for(my $k = 0; $k < length($variant); $k++){
+ − 1690 my $fixed_variant = $variant;
+ − 1691 substr($fixed_variant, $k, 1) = substr($ref, $k, 1);
+ − 1692 if($fixed_variant eq $ref){ # single base difference at base k between the two
+ − 1693 $ref = substr($ref, $k, 1);
+ − 1694 $variant = substr($variant, $k, 1);
+ − 1695 $location += $k;
+ − 1696 last;
+ − 1697 }
+ − 1698 }
+ − 1699 }
+ − 1700
+ − 1701 # samtools reports indels with long common tails, causing non-canonical HGVS reporting and a problem when looking up the variant in dbSNP etc.
+ − 1702 # remove common tails to variant calls in order to try to rectify this
+ − 1703 else{
+ − 1704 while(length($ref) > 1 and length($variant) > 1 and substr($ref, length($ref)-1) eq substr($variant, length($variant)-1)){
+ − 1705 chop $ref; chop $variant;
+ − 1706 }
+ − 1707 }
+ − 1708
+ − 1709 # See if a caveat should be added because the indel was in a polyhomomer region
+ − 1710 if(length($ref) > length($variant) and index($ref, $variant) == 0 and $fasta_index->fetch("$chr:".($location+1)."-".($location+length($ref)+1)) =~ /^([ACGT])\1+$/i){
+ − 1711 if(not exists $chr2caveats{"$chr:$location"}){
+ − 1712 $chr2caveats{"$chr:$location"} = "poly".uc($1)." region deletion";
+ − 1713 }
+ − 1714 elsif($chr2caveats{"$chr:$location"} !~ /poly/){
+ − 1715 $chr2caveats{"$chr:$location"} .= "; poly".uc($1)." region deletion";
+ − 1716 }
+ − 1717 }
+ − 1718 elsif(length($ref) < length($variant) and index($variant, $ref) == 0 and substr($variant, 1) =~ /^([ACGT])\1+$/i){
+ − 1719 if(not exists $chr2caveats{"$chr:$location"}){
+ − 1720 $chr2caveats{"$chr:$location"} .= "poly".uc($1)." region insertion";
+ − 1721 }
+ − 1722 elsif($chr2caveats{"$chr:$location"} !~ /poly/){
+ − 1723 $chr2caveats{"$chr:$location"} .= "; poly".uc($1)." region insertion";
+ − 1724 }
+ − 1725 }
+ − 1726
+ − 1727 # Not a protein-coding gene? Report genomic cooordinates for HGVS
+ − 1728 if(not exists $feature_cds_max{$target_parent} or not defined $feature_cds_max{$target_parent} or $feature_cds_max{$target_parent} eq ""){
+ − 1729 if(length($ref) == 1 and length($variant) == 1){ # SNP
+ − 1730 record_snv("$target_parent\tg.$location",
+ − 1731 "$ref>$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 1732 join("\t",prop_info_key($chr,$location,$ref,$variant)),"\tNA\n");
+ − 1733 }
+ − 1734 elsif(length($ref) == 1 and length($variant) > 1){ # insertion
+ − 1735 record_snv("$target_parent\tg.$location\_",($location+1),
+ − 1736 "ins",substr($variant, 1),"\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 1737 join("\t",prop_info_key($chr,$location,$ref,$variant)),"\tNA\n");
+ − 1738 }
+ − 1739 elsif(length($variant) == 1 and length($ref) > 1){ # deletion
+ − 1740 record_snv("$target_parent\tg.$location\_",($location+length($ref)-1),
+ − 1741 "del",substr($ref, 1),"\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 1742 join("\t",prop_info_key($chr,$location,$ref,$variant)),"\tNA\n");
+ − 1743 }
+ − 1744 else{ # indel
+ − 1745 record_snv("$target_parent\tg.$location\_",($location+length($ref)-1),
+ − 1746 "delins$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 1747 join("\t",prop_info_key($chr,$location,$ref,$variant)),"\tNA\n");
+ − 1748 }
+ − 1749 next; # process next record, we've done all we can with a non-coding-gene SNP
+ − 1750 }
+ − 1751
+ − 1752 if($strand eq "-"){
+ − 1753 # set up utr offset for correct CDS coordinates
+ − 1754 for(my $i = $#feature_ranges; $i >= 0; $i--){
+ − 1755 # exon is completely 5' of the start
+ − 1756 if($feature_ranges[$i]->[0] > $feature_cds_max{$target_parent}){
+ − 1757 #print STDERR "Whole 5' UTR exon $i: ",$feature_ranges[$i]->[1]-$feature_ranges[$i]->[0]+1,"\n";
+ − 1758 $feature_offset -= $feature_ranges[$i]->[1]-$feature_ranges[$i]->[0]+1;
+ − 1759 }
+ − 1760 # exon with the cds start
+ − 1761 elsif($feature_ranges[$i]->[1] >= $feature_cds_max{$target_parent} and
+ − 1762 $feature_ranges[$i]->[0] <= $feature_cds_max{$target_parent}){
+ − 1763 #print STDERR "Start codon in exon $i: ", $feature_cds_max{$target_parent} - $feature_ranges[$i]->[1],"\n";
+ − 1764 $feature_offset += $feature_cds_max{$target_parent} - $feature_ranges[$i]->[1];
+ − 1765 last;
+ − 1766 }
+ − 1767 else{
+ − 1768 die "The CDS for $target_parent (on negative strand) ends downstream ",
+ − 1769 "($feature_cds_max{$target_parent}) of the an exon",
+ − 1770 " (", $feature_ranges[$i]->[0], "), which is illogical. Please revise the GFF file provided.\n";
+ − 1771 }
+ − 1772 }
+ − 1773 for(my $i = $#feature_ranges; $i >= 0; $i--){
+ − 1774 my $feature = $feature_ranges[$i];
+ − 1775 # in the 3' UTR region of the gene
+ − 1776 if($location < $feature_cds_min{$target_parent}){
+ − 1777 my $feature_exon = 0;
+ − 1778 $feature = $feature_ranges[$feature_exon];
+ − 1779 while($location > $feature->[1]+$flanking_bases and
+ − 1780 $feature_exon < $#feature_ranges){
+ − 1781 $feature = $feature_ranges[++$feature_exon]; # find the 3' utr exon in which the variant is located
+ − 1782 }
+ − 1783 # utr exons passed entirely
+ − 1784 my $stop_exon = $feature_exon;
+ − 1785 while($feature_ranges[$stop_exon]->[1] < $feature_cds_min{$target_parent}){
+ − 1786 $stop_exon++;
+ − 1787 }
+ − 1788 my $post_offset = $feature_cds_min{$target_parent}-$feature_ranges[$stop_exon]->[0]; # offset from the stop codon in the final coding exon
+ − 1789 while($stop_exon > $feature_exon){
+ − 1790 $post_offset += $feature_ranges[$stop_exon]->[1]-$feature_ranges[$stop_exon]->[0]+1;
+ − 1791 $stop_exon--;
+ − 1792 }
+ − 1793
+ − 1794 my $pos = $feature->[1]-$location+1+$post_offset;
+ − 1795 my $junction_dist;
+ − 1796 if($location < $feature->[0]){ # after a UTR containing exon? set as .*DD+DD
+ − 1797 $junction_dist = ($feature->[0]-$location);
+ − 1798 $pos = ($post_offset+$feature->[1]-$feature->[0]+1)."+".$junction_dist;
+ − 1799 }
+ − 1800 elsif($location > $feature->[1]){ # before a total UTR exon? set as .*DD-DD
+ − 1801 $junction_dist = -($location-$feature->[1]);
+ − 1802 $pos = ($post_offset+1).$junction_dist;
+ − 1803 }
+ − 1804 else{ # in the UTR exon
+ − 1805 if($location - $feature->[0] < $feature->[1] - $location){ # location is closer to exon donor site
+ − 1806 $junction_dist = -($location - $feature->[0]+1); # +1 for HGVS syntax compatibility (there is no position 0, direct skip from -1 to +1)
+ − 1807 }
+ − 1808 else{
+ − 1809 $junction_dist = $feature->[1] - $location +1;
+ − 1810 }
+ − 1811 }
+ − 1812 if(length($ref) == 1 and length($variant) == 1){
+ − 1813 my $rc = join("",map({$rc{$_}} split(//,reverse($variant))));
+ − 1814 # 3' UTR SNP
+ − 1815 record_snv("$target_parent\tc.*$pos",
+ − 1816 "$rc{$ref}>$rc\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 1817 #"$ref>$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 1818 join("\t",prop_info_key($chr,$location,$ref,$variant,$junction_dist)),"\tNA\n");
+ − 1819 }
+ − 1820 elsif(length($ref) == 1 and length($variant) > 1 and substr($variant, 0, 1) eq $ref){
+ − 1821 my $rc = join("",map({$rc{$_}} split(//,reverse(substr($variant,1)))));
+ − 1822 # 3' UTR insertion
+ − 1823 record_snv("$target_parent\tc.*",
+ − 1824 hgvs_plus($pos,-1),"_*",$pos,
+ − 1825 "ins$rc\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 1826 #"ins",substr($variant,1),"\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 1827 join("\t", prop_info_key($chr,$location,$ref,$variant,$junction_dist)),"\tNA\n");
+ − 1828 }
+ − 1829 elsif(length($ref) > 1 and length($variant) == 1 and substr($ref, 0, 1) eq $variant){
+ − 1830 my $rc = join("",map({$rc{$_}} split(//,reverse($ref))));
+ − 1831 my $delBases = substr($rc,0,length($rc)-1);
+ − 1832 if(length($ref) == 2){
+ − 1833 # 3' UTR single base deletion
+ − 1834 record_snv("$target_parent\tc.*",hgvs_plus($pos,-1),
+ − 1835 "del$delBases\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 1836 join("\t", prop_info_key($chr,$location,$ref,$variant,$junction_dist)),"\tNA\n");
+ − 1837 }
+ − 1838 else{
+ − 1839 # longer 3' UTR deletion
+ − 1840 record_snv("$target_parent\tc.*",
+ − 1841 hgvs_plus($pos,-length($ref)+1),"_*",hgvs_plus($pos, -1),
+ − 1842 "del$delBases\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 1843 join("\t", prop_info_key($chr,$location,$ref,$variant,$junction_dist)),"\tNA\n");
+ − 1844 }
+ − 1845 }
+ − 1846 else{
+ − 1847 my $rc = join("",map({$rc{$_}} split(//,reverse($variant))));
+ − 1848 if($rc eq $ref and length($variant) > 3){
+ − 1849 # 3' UTR inversion
+ − 1850 record_snv("$target_parent\tc.*",
+ − 1851 hgvs_plus($pos,-length($ref)+1),"_*",$pos,
+ − 1852 "inv\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 1853 join("\t", prop_info_key($chr,$location,$ref,$variant,$junction_dist)),"\tNA\n");
+ − 1854 last;
+ − 1855 }
+ − 1856
+ − 1857 # complex substitution in 3' UTR
+ − 1858 # Will break if stop codon is modified
+ − 1859 record_snv("$target_parent\tc.*",
+ − 1860 hgvs_plus($pos,-length($ref)+1),"_*", $pos,
+ − 1861 "delins$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 1862 join("\t", prop_info_key($chr,$location,$ref,$variant,$junction_dist)),"\tNA\n");
+ − 1863 }
+ − 1864 last;
+ − 1865 }
+ − 1866 # in the feature
+ − 1867 elsif($location >= $feature->[0] and $location <= $feature->[1]){
+ − 1868 my $pos = $feature_offset+$feature->[1]-$location+1;
+ − 1869 if($location > $feature_cds_max{$target_parent}){ #since there is no position 0, the pos is in UTR, subtract one
+ − 1870 $pos = hgvs_plus($pos, -1);
+ − 1871 }
+ − 1872 my $first_exon_base = $feature_offset+1;
+ − 1873 my $exon_edge_dist = $feature->[1]-$location+1; # equiv of HGVS +X or -X intron syntax, but for exons
+ − 1874 $exon_edge_dist = $feature->[0]-$location-1 if abs($feature->[0]-$location-1) < $exon_edge_dist; # pick closer of donor and acceptor sites
+ − 1875 #print STDERR "Got ", $feature->[1]-$location+1, "vs. ", $feature->[0]-$location-1, ": chose $exon_edge_dist\n";
+ − 1876 if(length($ref) == 1 and length($variant) == 1){
+ − 1877 # SNP
+ − 1878 record_snv("$target_parent\tc.",
+ − 1879 $pos,
+ − 1880 "$rc{$ref}>$rc{$variant}\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 1881 join("\t", prop_info_key($chr,$location,$ref,$variant,$exon_edge_dist)),"\t",
+ − 1882 ($pos < 1 ? "NA" : hgvs_protein($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+ − 1883 #($pos < 1 ? "NA" : hgvs_protein_key($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+ − 1884 }
+ − 1885 elsif(length($ref) == 1 and length($variant) > 1 and substr($variant, 0, 1) eq $ref){
+ − 1886 my $rc = join("",map({$rc{$_}} split(//,reverse($variant))));
+ − 1887 my $insBases = substr($rc,1);
+ − 1888 # insertion
+ − 1889 record_snv("$target_parent\tc.",
+ − 1890 hgvs_plus_exon($pos, -1, $first_exon_base),"_",$pos,"ins$insBases",
+ − 1891 "\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 1892 join("\t", prop_info_key($chr,$location,$ref,$variant,$exon_edge_dist)),"\t",
+ − 1893 ($pos < 1 ? "NA" : hgvs_protein($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+ − 1894 #($pos < 1 ? "NA" : hgvs_protein_key($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+ − 1895 }
+ − 1896 elsif(length($ref) > 1 and length($variant) == 1 and substr($ref, 0, 1) eq $variant){
+ − 1897 my $rc = join("",map({$rc{$_}} split(//,reverse($ref))));
+ − 1898 my $delBases = substr($rc,0,length($rc)-1);
+ − 1899 # single nucleotide deletion
+ − 1900 if(length($ref) == 2){
+ − 1901 record_snv("$target_parent\tc.",
+ − 1902 hgvs_plus_exon($pos, -1, $first_exon_base),
+ − 1903 "del$delBases\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 1904 join("\t", prop_info_key($chr,$location,$ref,$variant,$exon_edge_dist)),"\t",
+ − 1905 ($pos-1 < 1 ? "NA" : $pos-1 < $first_exon_base ? "p.?" : hgvs_protein($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+ − 1906 #($pos-1 < 1 ? "NA" : $pos-1 < $first_exon_base ? "p.?" : hgvs_protein_key($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+ − 1907 }
+ − 1908 # longer deletion
+ − 1909 else{
+ − 1910 $exon_edge_dist = $feature->[1]-$location-length($ref)+1 if abs($feature->[1]-$location-length($ref)+1) < $exon_edge_dist;
+ − 1911 record_snv("$target_parent\tc.",
+ − 1912 hgvs_plus_exon($pos, -length($ref)+1, $first_exon_base),"_",
+ − 1913 hgvs_plus_exon($pos, -1, $first_exon_base),
+ − 1914 "del$delBases\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 1915 join("\t", prop_info_key($chr,$location,$ref,$variant,$exon_edge_dist)),"\t",
+ − 1916 ($pos-1 < 1 ? "NA" : $pos-length($ref)+1 < $first_exon_base ? "p.?" : hgvs_protein($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+ − 1917 #($pos-1 < 1 ? "NA" : $pos-length($ref)+1 < $first_exon_base ? "p.?" : hgvs_protein_key($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+ − 1918 }
+ − 1919 }
+ − 1920 else{
+ − 1921 # complex substitution
+ − 1922 $exon_edge_dist = $feature->[1]-$location-length($ref)+1 if abs($feature->[1]-$location-length($ref)+1) < $exon_edge_dist;
+ − 1923 my $rc = join("",map({$rc{$_}} split(//,reverse($variant))));
+ − 1924 if($rc eq $variant and length($variant) > 3){
+ − 1925 # inversion
+ − 1926 record_snv("$target_parent\tc.",
+ − 1927 hgvs_plus_exon($pos,-length($ref)+1,$first_exon_base),"_",
+ − 1928 $pos,
+ − 1929 "inv",
+ − 1930 "\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 1931 join("\t", prop_info_key($chr,$location,$ref,$variant,$exon_edge_dist)),"\t",
+ − 1932 ($pos < 1 ? "NA" : $pos-length($ref)+1 < $first_exon_base ? "p.?" : hgvs_protein($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+ − 1933
+ − 1934 last;
+ − 1935 }
+ − 1936 record_snv("$target_parent\tc.",
+ − 1937 hgvs_plus_exon($pos,-length($ref)+1,$first_exon_base),"_",
+ − 1938 $pos,
+ − 1939 "delins$rc",
+ − 1940 "\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 1941 join("\t", prop_info_key($chr,$location,$ref,$variant,$exon_edge_dist)),"\t",
+ − 1942 ($pos < 1 ? "NA" : $pos-length($ref)+1 < $first_exon_base ? "p.?" : hgvs_protein($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+ − 1943 #($pos < 1 ? "NA" : $pos-length($ref)+1 < $first_exon_base ? "p.?" : hgvs_protein_key($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+ − 1944 }
+ − 1945 last;
+ − 1946 }
+ − 1947 # 5' of feature (on negative strand)
+ − 1948 elsif($location > $feature->[1]){
+ − 1949 if(length($ref) == 1 and length($variant) == 1){
+ − 1950 # intronic SNP
+ − 1951 if($i == $#feature_ranges or $feature->[1]-$location >= -1*$flanking_bases){
+ − 1952 # Closer to acceptor site
+ − 1953 record_snv("$target_parent\tc.",$feature_offset+1,
+ − 1954 ($feature->[1]-$location),
+ − 1955 "$rc{$ref}>$rc{$variant}\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 1956 #"$ref>$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 1957 join("\t", prop_info_key($chr,$location,$ref,$variant, $feature->[1]-$location)),"\tNA\n");
+ − 1958 }
+ − 1959 else{
+ − 1960 # Closer to donor site
+ − 1961 record_snv("$target_parent\tc.",$feature_offset,"+",
+ − 1962 ($feature_ranges[$i+1]->[0]-$location),
+ − 1963 "$rc{$ref}>$rc{$variant}\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 1964 #"$ref>$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 1965 join("\t", prop_info_key($chr,$location,$ref,$variant, $feature_ranges[$i+1]->[0]-$location)),"\tNA\n");
+ − 1966 }
+ − 1967 }
+ − 1968 elsif(length($ref) == 1 and length($variant) > 1 and substr($variant, 0, 1) eq $ref){
+ − 1969 my $rc = join("",map({$rc{$_}} split(//,reverse(substr($variant,1)))));
+ − 1970 if($i == $#feature_ranges or $feature->[1]-$location >= -1*$flanking_bases){
+ − 1971 # intronic insertion near acceptor
+ − 1972 my $pos = ($feature_offset+1).($feature->[1]-$location);
+ − 1973 record_snv("$target_parent\tc.",
+ − 1974 hgvs_plus($pos,-1),"_",$pos,
+ − 1975 "ins",
+ − 1976 $rc,"\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 1977 #substr($variant, 1),"\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 1978 join("\t", prop_info_key($chr,$location,$ref,$variant,$feature->[1]-$location-1)),"\tNA\n");
+ − 1979 }
+ − 1980 else{
+ − 1981 # intronic insertion near donor
+ − 1982 my $pos = $feature_offset."+".($feature_ranges[$i+1]->[0]-$location);
+ − 1983 record_snv("$target_parent\tc.",
+ − 1984 hgvs_plus($pos,-1),"_",$pos,
+ − 1985 "ins",
+ − 1986 $rc,"\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 1987 #substr($variant, 1),"\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 1988 join("\t", prop_info_key($chr,$location,$ref,$variant,$feature_ranges[$i+1]->[0]-$location+1)),"\tNA\n");
+ − 1989 }
+ − 1990 }
+ − 1991 elsif(length($ref) > 1 and length($variant) == 1 and substr($ref, 0, 1) eq $variant){
+ − 1992 # intronic deletion
+ − 1993 # single nucleotide deletion
+ − 1994 my $rc = reverse($ref);
+ − 1995 $rc=~tr/ACGT/TGCA/;
+ − 1996 my $delBases = substr($rc, 0, length($rc)-1);
+ − 1997 if(length($ref) == 2){
+ − 1998 # single intronic deletion near acceptor
+ − 1999 if($i == $#feature_ranges or $feature->[1]-$location >= -1*$flanking_bases){
+ − 2000 my $off = $feature->[1]-$location-1;
+ − 2001 record_snv("$target_parent\tc.",
+ − 2002 ($feature_offset+1),$off,
+ − 2003 "del$delBases\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 2004 join("\t", prop_info_key($chr,$location,$ref,$variant,$off)),"\t",($off >= -2 ? "p.?" : "NA"),"\n");
+ − 2005 }
+ − 2006 # single intronic deletion near donor
+ − 2007 else{
+ − 2008 my $pos = $feature_offset;
+ − 2009 my $off = $feature_ranges[$i+1]->[0]-$location+1;
+ − 2010 record_snv("$target_parent\tc.",
+ − 2011 hgvs_plus_exon($pos, $off, $pos),
+ − 2012 "del$delBases\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 2013 join("\t", prop_info_key($chr,$location,$ref,$variant,$off)),"\t",($off <= 2 ? "p.?" : "NA"),"\n");
+ − 2014 }
+ − 2015 }
+ − 2016 # longer deletion
+ − 2017 else{
+ − 2018 if($i == $#feature_ranges or $feature->[1]-$location >= -1*$flanking_bases){
+ − 2019 # long intronic deletion near acceptor
+ − 2020 my $off = $feature->[1]-$location-1;
+ − 2021 my $pos = ($feature_offset+1).$off;
+ − 2022 record_snv("$target_parent\tc.",
+ − 2023 hgvs_plus($pos,-length($ref)+2),"_",$pos,
+ − 2024 "del$delBases\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 2025 join("\t", prop_info_key($chr,$location,$ref,$variant,$off)),"\t",($off >= -2 ? "p.?" : "NA"),"\n");
+ − 2026 }
+ − 2027 else{
+ − 2028 # long intronic deletion near donor
+ − 2029 my $off = $feature_ranges[$i+1]->[0]-$location+1;
+ − 2030 my $pos = ($feature_offset)."+".$off;
+ − 2031 record_snv("$target_parent\tc.",
+ − 2032 $pos,"_",hgvs_plus($pos,-length($ref)-1),
+ − 2033 "del$delBases\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 2034 join("\t", prop_info_key($chr,$location,$ref,$variant,$off)),"\t",($off-length($ref)+1 <= 2 ? "p.?" : "NA"),"\n");
+ − 2035 }
+ − 2036 last;
+ − 2037 }
+ − 2038 }
+ − 2039 else{
+ − 2040 my $rc = reverse($ref);
+ − 2041 $rc=~tr/ACGT/TGCA/;
+ − 2042 if($rc eq $variant and length($variant) > 3){
+ − 2043 # intronic inversion near acceptor site
+ − 2044 if($i == $#feature_ranges or $feature->[1]-$location >= -1*$flanking_bases){
+ − 2045 my $pos = ($feature_offset+1).($feature->[1]-$location);
+ − 2046 record_snv("$target_parent\tc.",
+ − 2047 hgvs_plus($pos,-length($ref)+1),"_",$pos,
+ − 2048 "inv\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 2049 join("\t", prop_info_key($chr,$location,$ref,$variant,$feature->[1]-$location)),"\tNA\n");
+ − 2050 }
+ − 2051 else{
+ − 2052 my $pos = ($feature_offset)."+".($feature_ranges[$i+1]->[0]-$location);
+ − 2053 record_snv("$target_parent\tc.",
+ − 2054 $pos,"_",hgvs_plus($pos, length($ref)-1),
+ − 2055 "inv\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 2056 join("\t", prop_info_key($chr,$location,$ref,$variant,$feature_ranges[$i+1]->[0]-$location)),"\tNA\n");
+ − 2057 }
+ − 2058 last;
+ − 2059 }
+ − 2060 $rc = reverse($variant);
+ − 2061 $rc=~tr/ACGT/TGCA/;
+ − 2062 # Intronic complex substitution
+ − 2063 if($i == $#feature_ranges or $feature->[1]-$location >= -1*$flanking_bases){
+ − 2064 # complex intronic substitution near acceptor site
+ − 2065 my $pos = ($feature_offset+1).($feature->[1]-$location);
+ − 2066 record_snv("$target_parent\tc.",
+ − 2067 hgvs_plus($pos, -length($ref)+1),"_",$pos,
+ − 2068 "delins$rc\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 2069 #"delins$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 2070 join("\t", prop_info_key($chr,$location,$ref,$variant,$feature->[1]-$location)),"\tNA\n");
+ − 2071 }
+ − 2072 else{
+ − 2073 # complex intronic substitution near donor site
+ − 2074 my $pos = $feature_offset."+".($feature_ranges[$i+1]->[0]-$location);
+ − 2075 record_snv("$target_parent\tc.",
+ − 2076 $pos,"_",hgvs_plus($pos, length($ref)-1),
+ − 2077 "delins$rc\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 2078 #"delins$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 2079 join("\t", prop_info_key($chr,$location,$ref,$variant,$feature_ranges[$i+1]->[0]-$location)),"\tNA\n");
+ − 2080 }
+ − 2081 }
+ − 2082 last;
+ − 2083 }
+ − 2084 else{
+ − 2085 #print STDERR "Offset going from ", $feature_offset, " to ", $feature_offset+$feature->[1]-$feature->[0]+1,"\n";
+ − 2086 $feature_offset += $feature->[1]-$feature->[0]+1;
+ − 2087 #print STDERR "Set feature offset to $feature_offset by adding ",$feature->[1],"-", $feature->[0],"+1\n";
+ − 2088 }
+ − 2089 }
+ − 2090 }
+ − 2091 else{
+ − 2092 # forward strand
+ − 2093
+ − 2094 # set up utr offset for correct CDS coordinates
+ − 2095 for(my $i = 0; $i <= $#feature_ranges; $i++){
+ − 2096 # All 5' utr exon
+ − 2097 if($feature_ranges[$i]->[1] < $feature_cds_min{$target_parent}){
+ − 2098 $feature_offset -= $feature_ranges[$i]->[1]-$feature_ranges[$i]->[0]+1;
+ − 2099 }
+ − 2100 # exon with the cds start
+ − 2101 elsif($feature_ranges[$i]->[1] >= $feature_cds_min{$target_parent} and
+ − 2102 $feature_ranges[$i]->[0] <= $feature_cds_min{$target_parent}){
+ − 2103 $feature_offset -= $feature_cds_min{$target_parent} - $feature_ranges[$i]->[0];
+ − 2104 last;
+ − 2105 }
+ − 2106 else{
+ − 2107 die "The CDS for $target_parent starts upstream ($feature_cds_max{$target_parent}) of the first exon",
+ − 2108 " (", $feature_ranges[$i]->[0], "), which is illogical. Please revise the GFF file provided.\n";
+ − 2109 }
+ − 2110 }
+ − 2111 for(my $i = 0; $i <= $#feature_ranges; $i++){
+ − 2112 my $feature = $feature_ranges[$i];
+ − 2113 # 3' of last coding position
+ − 2114 if($location > $feature_cds_max{$target_parent}){
+ − 2115 # find the exon with the stop codon
+ − 2116 while($feature->[1] < $feature_cds_max{$target_parent}){
+ − 2117 $feature = $feature_ranges[++$i];
+ − 2118 }
+ − 2119 my $post_offset = $feature->[0]-$feature_cds_max{$target_parent};
+ − 2120 while($location > $feature->[1]+$flanking_bases and
+ − 2121 $i < $#feature_ranges){
+ − 2122 $post_offset += $feature->[1]-$feature->[0]+1;
+ − 2123 $feature = $feature_ranges[++$i]; # find the 3' utr exon in which the variant is located
+ − 2124 }
+ − 2125 my $pos = $location-$feature->[0]+$post_offset;
+ − 2126 #print STDERR "Positive strand: got $pos for $location, exon #$i of $#feature_ranges, post_offset is $post_offset\n" if $location-$feature->[1] > $flanking_bases;
+ − 2127 my $off;
+ − 2128 if($location > $feature->[1]){ # after a UTR containing exon? set as .*DD+DD
+ − 2129 $off = $location-$feature->[1];
+ − 2130 $pos = ($post_offset+$feature->[1]-$feature->[0]+1)."+".$off;
+ − 2131 }
+ − 2132 elsif($location < $feature->[0]){ # before a total UTR exon? set as .*DD-DD
+ − 2133 $off = -($feature->[0]-$location);
+ − 2134 $pos = ($post_offset+1).$off;
+ − 2135 }
+ − 2136 else{
+ − 2137 if($location-$feature->[0] < $feature->[1]-$location){
+ − 2138 $off = $location-$feature->[0]+1; # +1 since HGVS skips right from -1 to +1 at exon boundary
+ − 2139 }
+ − 2140 else{
+ − 2141 $off = $location-$feature->[1]-1; # will be negative
+ − 2142 }
+ − 2143 }
+ − 2144 if(length($ref) == 1 and length($variant) == 1){
+ − 2145 # 3' UTR SNP
+ − 2146 record_snv("$target_parent\tc.*$pos",
+ − 2147 "$ref>$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 2148 join("\t", prop_info_key($chr,$location,$ref,$variant, $off)),"\tNA\n");
+ − 2149 }
+ − 2150 elsif(length($ref) == 1 and length($variant) > 1 and substr($variant, 0, 1) eq $ref){
+ − 2151 # 3' UTR insertion
+ − 2152 record_snv("$target_parent\tc.*",
+ − 2153 hgvs_plus($pos,1),"_*",hgvs_plus($pos,2),
+ − 2154 "ins",substr($variant,1),"\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 2155 join("\t", prop_info_key($chr,$location,$ref,$variant, $off)),"\tNA\n");
+ − 2156 }
+ − 2157 elsif(length($ref) > 1 and length($variant) == 1 and substr($ref, 0, 1) eq $variant){
+ − 2158 my $delBases = substr($ref, 1);
+ − 2159 if(length($ref) == 2){
+ − 2160 # 3' UTR single base deletion
+ − 2161 record_snv("$target_parent\tc.*",hgvs_plus($pos,1),
+ − 2162 "del$delBases\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 2163 join("\t", prop_info_key($chr,$location,$ref,$variant,$off)),"\tNA\n");
+ − 2164 }
+ − 2165 else{
+ − 2166 # longer 3' UTR deletion
+ − 2167 record_snv("$target_parent\tc.*",
+ − 2168 hgvs_plus($pos,1),"_*",hgvs_plus($pos,length($ref)-1),
+ − 2169 "del$delBases\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 2170 join("\t", prop_info_key($chr,$location,$ref,$variant, $off)),"\tNA\n");
+ − 2171 }
+ − 2172 }
+ − 2173 else{
+ − 2174 my $rc = reverse($ref);
+ − 2175 $rc=~tr/ACGT/TGCA/;
+ − 2176 if($rc eq $variant and length($variant) > 3){
+ − 2177 # 3' UTR inversion
+ − 2178 record_snv("$target_parent\tc.*$pos",
+ − 2179 "_*",hgvs_plus($pos,length($ref)-1),
+ − 2180 "inv\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 2181 join("\t", prop_info_key($chr,$location,$ref,$variant,$off)),"\tNA\n");
+ − 2182 last;
+ − 2183 }
+ − 2184 # complex substitution in 3' UTR
+ − 2185 record_snv("$target_parent\tc.*$pos",
+ − 2186 "_*",hgvs_plus($pos,length($ref)-1),
+ − 2187 "delins$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 2188 join("\t", prop_info_key($chr,$location,$ref,$variant,$off)),"\tNA\n");
+ − 2189 }
+ − 2190 last;
+ − 2191 }
+ − 2192 # in the exon
+ − 2193 elsif($location >= $feature->[0] and $location <= $feature->[1]){
+ − 2194 my $pos = $feature_offset+$location-$feature->[0]+1;
+ − 2195 my $last_exon_base = $feature_offset+$feature->[1]-$feature->[0]+1;
+ − 2196 my $exon_edge_dist = $location-$feature->[0]+1; # equiv of HGVS +X or -X intron syntax, but for exons
+ − 2197 $exon_edge_dist = $location-$feature->[1]-1 if abs($location-$feature->[1]-1) < $exon_edge_dist; # pick closer of donor and acceptor sites
+ − 2198 #print STDERR "Got ", $location-$feature->[0]+1, "vs. ", $location-$feature->[1]-1, ": chose $exon_edge_dist\n";
+ − 2199 if($location < $feature_cds_min{$target_parent}){ #since there is no position 0, the pos is in UTR, subtract one
+ − 2200 $pos = hgvs_plus($pos, -1);
+ − 2201 }
+ − 2202 if(length($ref) == 1 and length($variant) == 1){
+ − 2203 # SNP
+ − 2204 record_snv("$target_parent\tc.$pos",
+ − 2205 "$ref>$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 2206 join("\t", prop_info_key($chr,$location,$ref,$variant,$exon_edge_dist)),"\t",
+ − 2207 ($pos < 1 ? "NA" : hgvs_protein($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+ − 2208 #($pos < 1 ? "NA" : hgvs_protein_key($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+ − 2209 }
+ − 2210 elsif(length($ref) == 1 and length($variant) > 1 and substr($variant, 0, 1) eq $ref){
+ − 2211 # insertion
+ − 2212 record_snv("$target_parent\tc.$pos",
+ − 2213 "_",hgvs_plus_exon($pos,1,$last_exon_base),"ins",
+ − 2214 substr($variant, 1),"\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 2215 join("\t", prop_info_key($chr,$location,$ref,$variant,$exon_edge_dist)),"\t",
+ − 2216 ($pos < 1 ? "NA" : hgvs_protein($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+ − 2217 #($pos < 1 ? "NA" : hgvs_protein_key($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+ − 2218 }
+ − 2219 elsif(length($ref) > 1 and length($variant) == 1 and substr($ref, 0, 1) eq $variant){
+ − 2220 my $delBases = substr($ref, 1);
+ − 2221 # single nucleotide deletion
+ − 2222 if(length($delBases) == 1){
+ − 2223 record_snv("$target_parent\tc.",
+ − 2224 hgvs_plus_exon($pos,1,$last_exon_base),
+ − 2225 "del$delBases\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 2226 join("\t", prop_info_key($chr,$location,$ref,$variant,$exon_edge_dist)),"\t",
+ − 2227 ($pos < 1 or $pos > $last_exon_base ? "NA" : hgvs_protein($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+ − 2228 #($pos < 1 or $pos > $last_exon_base ? "NA" : hgvs_protein_key($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+ − 2229 }
+ − 2230 # longer deletion
+ − 2231 else{
+ − 2232 $exon_edge_dist = $feature->[1]-$location-length($ref)-1 if abs($feature->[1]-$location-length($ref)-1) < $exon_edge_dist;
+ − 2233 record_snv("$target_parent\tc.",
+ − 2234 hgvs_plus_exon($pos,1,$last_exon_base),"_",
+ − 2235 hgvs_plus_exon($pos,length($ref)-1,$last_exon_base),
+ − 2236 "del$delBases\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 2237 join("\t", prop_info_key($chr,$location,$ref,$variant,$exon_edge_dist)),"\t",
+ − 2238 ($pos+length($ref)-1 > $last_exon_base ? "p.?" : $pos < 1 or $pos > $last_exon_base ? "NA" : hgvs_protein($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+ − 2239 #($pos+length($ref)-1 > $last_exon_base ? "p.?" : $pos < 1 or $pos > $last_exon_base ? "NA" : hgvs_protein_key($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+ − 2240 }
+ − 2241 }
+ − 2242 else{
+ − 2243 $exon_edge_dist = $feature->[1]-$location-length($ref)-1 if abs($feature->[1]-$location-length($ref)-1) < $exon_edge_dist;
+ − 2244 my $rc = reverse($ref);
+ − 2245 $rc=~tr/ACGT/TGCA/;
+ − 2246 if($rc eq $variant and length($variant) > 3){
+ − 2247 # inversion
+ − 2248 record_snv("$target_parent\tc.$pos",
+ − 2249 "_",hgvs_plus_exon($pos,length($ref)-1, $last_exon_base),
+ − 2250 "inv\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 2251 join("\t", prop_info_key($chr,$location,$ref,$variant,$exon_edge_dist)),"\t",
+ − 2252 ($pos+length($ref)-1 > $last_exon_base ? "p.?" : $pos < 1 or $pos > $last_exon_base ? "NA" : hgvs_protein($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+ − 2253 #($pos+length($ref)-1 > $last_exon_base ? "p.?" : $pos < 1 or $pos > $last_exon_base ? "NA" : hgvs_protein_key($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+ − 2254 last;
+ − 2255 }
+ − 2256 # complex substitution
+ − 2257 record_snv("$target_parent\tc.$pos",
+ − 2258 "_",hgvs_plus_exon($pos, length($ref)-1, $last_exon_base),
+ − 2259 "delins$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 2260 join("\t", prop_info_key($chr,$location,$ref,$variant,$exon_edge_dist)),"\t",
+ − 2261 ($pos+length($ref)-1 > $last_exon_base ? "p.?" : $pos < 1 or $pos > $last_exon_base ? "NA" : hgvs_protein($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+ − 2262 #($pos+length($ref)-1 > $last_exon_base ? "p.?" : $pos < 1 or $pos > $last_exon_base ? "NA" : hgvs_protein_key($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+ − 2263 }
+ − 2264 last;
+ − 2265 }
+ − 2266 # 5' of feature
+ − 2267 elsif($location < $feature->[0]){
+ − 2268 if(length($ref) == 1 and length($variant) == 1){
+ − 2269 # intronic SNP
+ − 2270 if($i != 0 and $location-$feature->[0] < -1*$flanking_bases){
+ − 2271 # Closer to donor site
+ − 2272 record_snv("$target_parent\tc.",$feature_offset,"+",
+ − 2273 ($location-$feature_ranges[$i-1]->[1]),
+ − 2274 "$ref>$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 2275 join("\t", prop_info_key($chr,$location,$ref,$variant,$location-$feature_ranges[$i-1]->[1])),"\tNA\n");
+ − 2276 }
+ − 2277 else{
+ − 2278 # Closer to acceptor site
+ − 2279 record_snv("$target_parent\tc.",$feature_offset+1,
+ − 2280 ($location-$feature->[0]),
+ − 2281 "$ref>$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 2282 join("\t", prop_info_key($chr,$location,$ref,$variant,$location-$feature->[0])),"\tNA\n");
+ − 2283 }
+ − 2284 }
+ − 2285 elsif(length($ref) == 1 and length($variant) > 1 and substr($variant, 0, 1) eq $ref){
+ − 2286 if($i != 0 and $location-$feature->[0] < -1*$flanking_bases){
+ − 2287 # intronic insertion near donor
+ − 2288 my $pos = $feature_offset."+".($location-$feature_ranges[$i-1]->[1]);
+ − 2289 record_snv("$target_parent\tc.",
+ − 2290 $pos,"_",hgvs_plus($pos,1),
+ − 2291 "ins",
+ − 2292 substr($variant, 1),"\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 2293 join("\t", prop_info_key($chr,$location,$ref,$variant,$location-$feature_ranges[$i-1]->[1])),"\tNA\n");
+ − 2294 }
+ − 2295 else{
+ − 2296 # intronic insertion near acceptor
+ − 2297 my $pos = ($feature_offset+1).($location-$feature->[0]);
+ − 2298 record_snv("$target_parent\tc.",
+ − 2299 $pos,"_",hgvs_plus($pos,1),
+ − 2300 "ins",
+ − 2301 substr($variant, 1),"\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 2302 join("\t", prop_info_key($chr,$location,$ref,$variant,$location-$feature->[0])),"\tNA\n");
+ − 2303 }
+ − 2304 }
+ − 2305 elsif(length($ref) > 1 and length($variant) == 1 and substr($ref, 0, 1) eq $variant){
+ − 2306 # intronic deletion
+ − 2307 # single nucleotide deletion
+ − 2308 my $delBases = substr($ref, 1);
+ − 2309 if(length($ref) == 2){
+ − 2310 # single intronic deletion near donor
+ − 2311 if($i != 0 and $location-$feature->[0] < -1*$flanking_bases){
+ − 2312 my $off = $location-$feature_ranges[$i-1]->[1]+1;
+ − 2313 record_snv("$target_parent\tc.",
+ − 2314 $feature_offset,"+",$off,
+ − 2315 "del$delBases\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 2316 join("\t", prop_info_key($chr,$location,$ref,$variant,$off)),"\t",($off <= 2 ? "p.?" : "NA"),"\n");
+ − 2317 }
+ − 2318 # single intronic deletion near acceptor
+ − 2319 else{
+ − 2320 my $pos = ($feature_offset+1);
+ − 2321 my $off = $location-$feature->[0];
+ − 2322 record_snv("$target_parent\tc.",
+ − 2323 hgvs_plus_exon($pos, $off, $pos),
+ − 2324 "del$delBases\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 2325 join("\t", prop_info_key($chr,$location,$ref,$variant,$off)),"\t",($off >= -2 ? "p.?" : "NA"),"\n");
+ − 2326 }
+ − 2327 }
+ − 2328 # longer deletion
+ − 2329 else{
+ − 2330 if($i != 0 and $location-$feature->[0] < -1*$flanking_bases){
+ − 2331 # long intronic deletion near donor
+ − 2332 my $off = $location-$feature_ranges[$i-1]->[1]+1;
+ − 2333 my $pos = $feature_offset."+".$off;
+ − 2334 record_snv("$target_parent\tc.",
+ − 2335 $pos,"_",hgvs_plus($pos,length($ref)-2),
+ − 2336 "del$delBases\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 2337 join("\t", prop_info_key($chr,$location,$ref,$variant,$off)),"\t",($off <= 2 ? "p.?" : "NA"),"\n");
+ − 2338 }
+ − 2339 else{
+ − 2340 # long intronic deletion near acceptor
+ − 2341 my $off = $location-$feature->[0]+1;
+ − 2342 my $pos = ($feature_offset+1).$off;
+ − 2343 record_snv("$target_parent\tc.",
+ − 2344 $pos,"_",hgvs_plus($pos,length($ref)-2),
+ − 2345 "del$delBases\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 2346 join("\t", prop_info_key($chr,$location,$ref,$variant,$off)),"\t",($off+length($ref)-2 >= -2 ? "p.?" : "NA"),"\n");
+ − 2347 }
+ − 2348 }
+ − 2349 }
+ − 2350 else{
+ − 2351 my $rc = reverse($ref);
+ − 2352 $rc=~tr/ACGT/TGCA/;
+ − 2353 if($rc eq $variant and length($variant) > 3){
+ − 2354 # intronic inversion near donor site
+ − 2355 if($i != 0 and $location-$feature->[0] < -1*$flanking_bases){
+ − 2356 my $pos = $feature_offset."+".($location-$feature_ranges[$i-1]->[1]);
+ − 2357 record_snv("$target_parent\tc.",
+ − 2358 $pos,"_",hgvs_plus($pos,length($ref)-1),
+ − 2359 "inv\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 2360 join("\t", prop_info_key($chr,$location,$ref,$variant,$location-$feature_ranges[$i-1]->[1])),"\tNA\n");
+ − 2361 }
+ − 2362 else{
+ − 2363 my $pos = ($feature_offset+1).($location-$feature->[0]);
+ − 2364 record_snv("$target_parent\tc.",
+ − 2365 $pos,"_",hgvs_plus($pos, length($ref)-1),
+ − 2366 "inv\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 2367 join("\t", prop_info_key($chr,$location,$ref,$variant,$location-$feature->[0])),"\tNA\n");
+ − 2368 }
+ − 2369 last;
+ − 2370 }
+ − 2371 # Intronic complex substitution
+ − 2372 # Note: sub maybe have comma in it to denote two possible variants
+ − 2373 if($i != 0 and $location-$feature->[0] < -1*$flanking_bases){
+ − 2374 # complex intronic substitution near donor site
+ − 2375 my $pos = $feature_offset."+".($location-$feature_ranges[$i-1]->[1]);
+ − 2376 record_snv("$target_parent\tc.",
+ − 2377 $pos,"_",hgvs_plus($pos, length($ref)-1),
+ − 2378 "delins$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 2379 join("\t", prop_info_key($chr,$location,$ref,$variant,$location-$feature_ranges[$i-1]->[1])),"\tNA\n");
+ − 2380 }
+ − 2381 else{
+ − 2382 # complex intronic substitution near acceptor site
+ − 2383 my $pos = ($feature_offset+1).($location-$feature->[0]);
+ − 2384 record_snv("$target_parent\tc.",
+ − 2385 $pos,"_",hgvs_plus($pos, length($ref)-1),
+ − 2386 "delins$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+ − 2387 join("\t", prop_info_key($chr,$location,$ref,$variant,$location-$feature->[0])),"\tNA\n");
+ − 2388 }
+ − 2389 }
+ − 2390 last;
+ − 2391 }
+ − 2392 else{
+ − 2393 # feature is past this exon
+ − 2394 $feature_offset += $feature->[1]-$feature->[0]+1;
+ − 2395 }
+ − 2396 }
+ − 2397 }
+ − 2398 } # for each variant on the line
+ − 2399 } # for each gene overlapping the site the VCF describes
+ − 2400 } # for each VCF line
+ − 2401 print STDERR "\n" unless $quiet;
+ − 2402 close(VCFIN);
+ − 2403
+ − 2404 # Before we can start printing the variants, we need to look at the phase information and calculate the real haplotype HGVS changes
+ − 2405 #if(keys %chr2phase){
+ − 2406 # Note that we could have samtools read-backed haplotype info, MNPs in the VCF, and pre-existing haplotypes in the input VCF (e.g. imputed or based on Mendelian inheritance where trios exist)
+ − 2407 # We need to create new disjoint sets of phased blocks from the (consistent) union of these data.
+ − 2408 # my $chr2phase2variants = combine_phase_data(\%chr2phase);
+ − 2409
+ − 2410 # TODO: Calculate protein HGVS syntax for each variant, now that all phase data has been incorporated
+ − 2411 #for my $chr (keys %$chr2phase2variants){
+ − 2412 # for my $phase (keys %{$chr2phase2variants{$chr}){
+ − 2413 # # apply all phased changes to the reference chromosomal seq
+ − 2414 # my $phased_seq = $seq{$chr}; #reference
+ − 2415 # # sort the variants from 3' to 5' so that edits after indels don't need adjustment in their ref coordinate
+ − 2416 # my @sorted_variants = sort {my($a_pos) = $a =~ /:(\d+):/; my($b_pos) = $b =~ /:(\d+):/; $b_pos <=> $a_pos} @{$chr2phase2variants{$chr}->{$phase}};
+ − 2417 # for my $variant (@sorted_variants){
+ − 2418 # }
+ − 2419 # }
+ − 2420 #}
+ − 2421 #}
+ − 2422
+ − 2423 # retrieve the MAF info en masse for each chromosome, as this is much more efficient
+ − 2424 my @out_lines;
+ − 2425 for my $snv (@snvs){
+ − 2426 chomp $snv;
+ − 2427 my @fields = split /\t/, $snv;
+ − 2428 # For CNVs, all the fields are already filled out
+ − 2429 if(@fields > 13){
+ − 2430 push @out_lines, join("\t", $feature_type{$fields[0]}, ($fields[0] =~ /\S/ ? $feature_length{$fields[0]} : "NA"), @fields);
+ − 2431 next;
+ − 2432 }
+ − 2433 my $variant_key = $fields[9];
+ − 2434 $fields[9] = join("\t", prop_info($dbsnp,$internal_snp,$variant_key));
+ − 2435 my $from = $fields[4];
+ − 2436 my $chr_pos_key = $fields[3].":".$from;
+ − 2437 my $to = $fields[4]; # true for SNPs and insertions
+ − 2438 my @variant_key = split /:/, $variant_key;
+ − 2439 # For deletions and complex variants, calculate the affected reference genome range and set the 'to'
+ − 2440 if(length($variant_key[2]) > 1){
+ − 2441 $to += length($variant_key[2])-1;
+ − 2442 }
+ − 2443 splice(@fields, 5, 0, $to);
+ − 2444
+ − 2445 # Otherwise expand the key into the relevant MAF values
+ − 2446 $fields[0] =~ s/\/chr.*$//; # some transcript ids are repeated... we expect "id/chr#" in the GTF file to distinguish these, but should get rid of them at reporting time
+ − 2447 # the offset from the nearest exon border if coding
+ − 2448 push @fields, ($#variant_key > 3 ? $variant_key[4] : "");
+ − 2449 # add gene name(s)
+ − 2450 push @fields, range2genes($fields[3], $from, $to+1);
+ − 2451 # add caveats
+ − 2452 my $c = $fields[3];
+ − 2453 if(not exists $chr2mappability{$c}){
+ − 2454 if($c =~ s/^chr//){
+ − 2455 # nothing more
+ − 2456 }
+ − 2457 else{
+ − 2458 $c = "chr$c";
+ − 2459 }
+ − 2460 }
+ − 2461 my $mappability_caveats = exists $chr2mappability{$c} ? $chr2mappability{$c}->fetch($fields[4], $fields[4]+1) : [];
+ − 2462 if(ref $mappability_caveats eq "ARRAY" and @$mappability_caveats){
+ − 2463 my %h;
+ − 2464 if(exists $chr2caveats{$chr_pos_key}){
+ − 2465 if($chr2caveats{$chr_pos_key} !~ /non-unique/){
+ − 2466 $chr2caveats{$chr_pos_key} = join("; ", grep {not $h{$_}++} @$mappability_caveats)."; ".$chr2caveats{$chr_pos_key};
+ − 2467 }
+ − 2468 }
+ − 2469 else{
+ − 2470 $chr2caveats{$chr_pos_key} = join("; ", grep {not $h{$_}++} @$mappability_caveats);
+ − 2471 }
+ − 2472 }
+ − 2473 push @fields, (exists $chr2caveats{$chr_pos_key} ? $chr2caveats{$chr_pos_key} : "");
+ − 2474 # add phase
+ − 2475 push @fields, find_phase($variant_key);
+ − 2476 push @out_lines, join("\t", $feature_type{$fields[0]}, $feature_length{$fields[0]}, @fields);
+ − 2477 }
+ − 2478
+ − 2479 # Now tabulate the rare variant numbers
+ − 2480 my %gene2rares;
+ − 2481 my %gene2aa_rares;
+ − 2482 for my $line (@out_lines){
+ − 2483 my @F = split /\t/, $line, -1;
+ − 2484 if($F[15] eq "NA" or $F[15] < $rare_variant_prop and (!$internal_snp or $F[17] < $rare_variant_prop)){
+ − 2485 my $gene_list = $internal_snp ? $F[20] : $F[19];
+ − 2486 next unless defined $gene_list;
+ − 2487 for my $g (split /; /, $gene_list){
+ − 2488 $gene2rares{$g}++;
+ − 2489 # Check the cDNA HGVS syntax for relevance
+ − 2490 if($F[3] =~ /c.\d+/ or # coding
+ − 2491 $F[3] =~ /c.\d+.*-[12]/ or # splicing acceptor
+ − 2492 $F[3] =~ /c.\d+\+[12345]/){ # splicing donor
+ − 2493 $gene2aa_rares{$g}++;
+ − 2494 }
+ − 2495 }
+ − 2496 }
+ − 2497 }
+ − 2498
+ − 2499 # Print the results
+ − 2500 for my $line (@out_lines){
+ − 2501 my @F = split /\t/, $line, -1;
+ − 2502 my $gene_list = $internal_snp ? $F[20] : $F[19];
+ − 2503 if(not defined $gene_list){
+ − 2504 print OUT join("\t", @F, "", ""), "\n"; next;
+ − 2505 }
+ − 2506
+ − 2507 my $max_gene_rare = 0;
+ − 2508 my $max_gene_aa_rare = 0;
+ − 2509 for my $g (split /; /, $gene_list){
+ − 2510 next unless exists $gene2rares{$g};
+ − 2511 if($gene2rares{$g} > $max_gene_rare){
+ − 2512 $max_gene_rare = $gene2rares{$g};
+ − 2513 }
+ − 2514 next unless exists $gene2aa_rares{$g};
+ − 2515 if($gene2aa_rares{$g} > $max_gene_aa_rare){
+ − 2516 $max_gene_aa_rare = $gene2aa_rares{$g};
+ − 2517 }
+ − 2518 }
+ − 2519 print OUT join("\t", @F, $max_gene_rare, $max_gene_aa_rare), "\n";
+ − 2520 }
+ − 2521 close(OUT);
+ − 2522
+ − 2523 sub range2genes{
+ − 2524 my ($c, $from, $to) = @_;
+ − 2525 if(not exists $gene_ids{$c}){
+ − 2526 if($c =~ s/^chr//){
+ − 2527 # nothing more
+ − 2528 }
+ − 2529 else{
+ − 2530 $c = "chr$c";
+ − 2531 }
+ − 2532 }
+ − 2533 if(exists $gene_ids{$c}){
+ − 2534 my %have;
+ − 2535 return join("; ", grep {not $have{$_}++} @{$gene_ids{$c}->fetch($from, $to+1)});
+ − 2536 }
+ − 2537 else{
+ − 2538 return "";
+ − 2539 }
+ − 2540 }
+ − 2541 sub combine_phase_data{
+ − 2542 my ($phases) = @_; # map from variant to its phase data
+ − 2543
+ − 2544 # Create a reverse mapping from phase regions to their variants
+ − 2545 my %chr2phase_region2variants;
+ − 2546 my @variants = keys %$phases;
+ − 2547 for my $variant (@variants){
+ − 2548 my ($chr) = $variant =~ /^\S+?-(\S+):/;
+ − 2549 $chr2phase_region2variants{$chr} = {} if not exists $chr2phase_region2variants{$chr};
+ − 2550 for my $phase_region (split /,/, $phases->{$variant}){
+ − 2551 $chr2phase_region2variants{$chr}->{$phase_region} = [] if not exists $chr2phase_region2variants{$chr}->{$phase_region};
+ − 2552 push @{$chr2phase_region2variants{$phase_region}}, $variant;
+ − 2553 }
+ − 2554 }
+ − 2555
+ − 2556 # Now for each phased block known so far, see if any variant in it is also part of another block
+ − 2557 # If so, do a union since phasing is both transitive and commutative.
+ − 2558 # The quickest way to do this is to check for overlapping intervals, then check for common members amongst those that do overlap
+ − 2559 for my $chr (keys %chr2phase_region2variants){
+ − 2560 my @ordered_phase_regions = sort {my($a_pos) = $a =~ /:(\d+)/; my($b_pos) = $b =~ /:(\d+)/; $a_pos <=> $b_pos} keys %{$chr2phase_region2variants{$chr}};
+ − 2561 my $sets = new DisjointSets(scalar(@ordered_phase_regions));
+ − 2562
+ − 2563 for (my $i = 0; $i < $#ordered_phase_regions; $i++){
+ − 2564 my ($start, $stop, $variant) = $ordered_phase_regions[$i];
+ − 2565 for (my $j = $i+1; $j <= $#ordered_phase_regions; $j++){
+ − 2566 my ($start2, $stop2, $variant2) = $ordered_phase_regions[$j];
+ − 2567 if($start2 > $stop){ # won't overlap any regions after this in the sorted list
+ − 2568 last;
+ − 2569 }
+ − 2570 # If we get here, it is implicit that $stop >= $start2 and $start < $stop2, i.e. there is overlap
+ − 2571 # Now check if there is a shared variant (otherwise we might erroneously join blocks from different physical chromosomal arms)
+ − 2572 my $have_shared_variant = 0;
+ − 2573 my $overlapping_phase_region = $ordered_phase_regions[$j];
+ − 2574 for my $variant (@{$chr2phase_region2variants{$chr}->{$ordered_phase_regions[$i]}}){
+ − 2575 if($phases->{$variant} =~ /\b$overlapping_phase_region\b/){
+ − 2576 $have_shared_variant = 1; last;
+ − 2577 }
+ − 2578 }
+ − 2579 # sanity check that there aren't conflicting variants in the new block (i.e. two different variants in the same position)
+ − 2580 my %pos2base;
+ − 2581 my $have_conflicting_variant = 0;
+ − 2582 for my $variant (@{$chr2phase_region2variants{$chr}->{$ordered_phase_regions[$i]}}, @{$chr2phase_region2variants{$chr}->{$ordered_phase_regions[$j]}}){
+ − 2583 my ($pos, $base) = $variant =~ /(\d+):(.+?)$/;
+ − 2584 if(exists $pos2base{$pos} and $pos2base{$pos} ne $base){
+ − 2585 # conflict, note with a caveat
+ − 2586 if(exists $chr2caveats{"$chr:$pos"}){
+ − 2587 $chr2caveats{"$chr:$pos"} .= "; inconsistent haplotype phasing" unless $chr2caveats{"$chr:$pos"} =~ /inconsistent haplotype phasing/;
+ − 2588 }
+ − 2589 else{
+ − 2590 $chr2caveats{"$chr:$pos"} = "inconsistent haplotype phasing";
+ − 2591 }
+ − 2592 $have_conflicting_variant ||= 1;
+ − 2593 }
+ − 2594 elsif(not exists $pos2base{$pos}){
+ − 2595 $pos2base{$pos} = $base;
+ − 2596 }
+ − 2597 }
+ − 2598
+ − 2599 $sets->union($i+1, $j+1) if $have_shared_variant and not $have_conflicting_variant; # indexes are one-based for sets rather than 0-based
+ − 2600 }
+ − 2601 }
+ − 2602 my $phase_sets = $sets->sets; #disjoint haplotype sets
+ − 2603 my %region_counts;
+ − 2604 for my $phase_set (@$phase_sets){
+ − 2605 next if scalar(@$phase_set) == 1; # No change to existing phase region (is disjoint from all others)
+ − 2606 # Generate a new haploblock to replace the old ones that are being merged
+ − 2607 my $merged_start = 10000000000;
+ − 2608 my $merged_end = 0;
+ − 2609 for my $ordered_phase_region_index (@$phase_set){
+ − 2610 my ($start, $end) = $ordered_phase_regions[$ordered_phase_region_index-1] =~ /(\d+)-(\d+)$/;
+ − 2611 $merged_start = $start if $start < $merged_start;
+ − 2612 $merged_end = $end if $end > $merged_end;
+ − 2613 }
+ − 2614 # At the start of the region is a unique prefix so we can tell the arms apart if two haploblocks have the exact same boundary
+ − 2615 my $region_count = $region_counts{"$merged_start-$merged_end"}++;
+ − 2616 my $merged_haploblock_name = "Y$region_count-$chr:$merged_start-$merged_end";
+ − 2617 # Assign this new name to overwrite the premerge values for each variant contained within
+ − 2618 for my $ordered_phase_region_index (@$phase_set){
+ − 2619 for my $variant (@{$chr2phase_region2variants{$chr}->{$ordered_phase_regions[$ordered_phase_region_index-1]}}){ # incl. one-based set correction in 0-based array index
+ − 2620 print STDERR "Merging $variant from ", $phases->{$variant}, " into new block $merged_haploblock_name\n";
+ − 2621 $phases->{$variant} = $merged_haploblock_name;
+ − 2622 }
+ − 2623 }
+ − 2624 }
+ − 2625 # TODO: if there are overlapping phase blocks still, but with different variants in the same position, we can infer that they are on the opposite strands...
+ − 2626 }
+ − 2627 }
+ − 2628
+ − 2629 # Sees if the positions of the variant are in the range of a phased haplotype, returning which allele it belongs to
+ − 2630 sub find_phase{
+ − 2631 my ($chr,$pos,$ref,$variant) = split /:/, $_[0];
+ − 2632 return "" if length($ref) != length($variant); # Can only deal with SNPs (and broken down MNPs) for now
+ − 2633 for(my $i = 0; $i < length($ref); $i++){
+ − 2634 my $key = "$chr:".($pos+$i).":".substr($variant, $i, 1);
+ − 2635 #print STDERR "Checking phase for $key\n" if $pos == 12907379;
+ − 2636 if(exists $chr2phase{$key}){
+ − 2637 #print STDERR "returning phase data $chr2phase{$key}\n";
+ − 2638 return $chr2phase{$key};
+ − 2639 }
+ − 2640 elsif(exists $chr2phase{"chr".$key}){
+ − 2641 #print STDERR "returning phase data ", $chr2phase{"chr".$key}, "\n";
+ − 2642 return $chr2phase{"chr".$key};
+ − 2643 }
+ − 2644 }
+ − 2645 return "";
+ − 2646 }
+ − 2647
+ − 2648 sub find_earliest_index{
+ − 2649 # employs a binary search to find the smallest index that must be the starting point of a search of [start,end] elements sorted in an array by start
+ − 2650 my ($query, $array) = @_;
+ − 2651
+ − 2652 return 0 if $query < $array->[0]->[0];
+ − 2653
+ − 2654 my ($left, $right, $prevCenter) = (0, $#$array, -1);
+ − 2655
+ − 2656 while(1){
+ − 2657 my $center = int (($left + $right)/2);
+ − 2658
+ − 2659 my $cmp = $query <=> $array->[$center]->[0] || ($center == 0 || $query != $array->[$center-1]->[0] ? 0 : -1);
+ − 2660
+ − 2661 return $center if $cmp == 0;
+ − 2662 if ($center == $prevCenter) {
+ − 2663 return $left;
+ − 2664 }
+ − 2665 $right = $center if $cmp < 0;
+ − 2666 $left = $center if $cmp > 0;
+ − 2667 $prevCenter = $center;
+ − 2668 }
+ − 2669 }
+ − 2670
