poor_gene_coverage: vcf2hgvs

comparison vcf2hgvs_table @ 0:7cdd13ff182a default tip

initial commit

author	Yusuf Ali <ali@yusuf.email>
date	Wed, 25 Mar 2015 15:49:28 -0600
parents
children

comparison

equal deleted inserted replaced

--1:000000000000
+:7cdd13ff182a
+#!/usr/bin/env perl
+BEGIN{
+my $prog_dir = `dirname $0`;
+chomp $prog_dir;
+push @INC, $prog_dir; # so DisjointSets.pm can be found no matter the working directory
+}
+use DisjointSets; # homebrew module
+use Bio::DB::Sam; # for FastA reference pulls
+use Bio::SeqUtils;
+use Bio::Tools::CodonTable;
+use Statistics::Zed;
+use Getopt::Long;
+use Set::IntervalTree;
+use strict;
+use warnings;
+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);
+if(@ARGV == 1 and $ARGV[0] eq "-v"){
+print "Version 1.0\n";
+exit;
+}
+#$max_args = `getconf ARG_MAX`; # largest number of args you can send to a system command (enviroment included, see limits.h)
+#chomp $max_args;
+$max_args = 4096; # if not defined $max_args or $max_args < 1; # the minimum since System V
+$max_args -= 50;
+# find out if a variant appears in the user provided data
+sub internal_prop($$$$){
+my ($chr,$pos,$ref,$variant) = @_;
+my $key = "$chr:$pos:$ref:$variant";
+if(exists $internal_prop{$key}){
+return $internal_prop{$key};
+}
+#print STDERR "Checking if internal_prop for $key exists: ";
+if(exists $chr2internal_vcf_lines{$chr}->{$pos}){
+for(@{$chr2internal_vcf_lines{$chr}->{$pos}}){
+my @fields = split /\t/, $_;
+if($pos == $fields[1] and length($fields[3]) == length($ref) and $fields[4] eq $variant){
+#print STDERR "yes\n";
+if(/MAF=(\d\.\d+)/){
+$internal_prop{$key} = $1; # change from percent to proportion
+return $1;
+}
+}
+}
+}
+else{
+#print STDERR "no\n";
+}
+$internal_prop{$key} = "NA";
+return "NA";
+}
+# find out if a variant appears in the NCBI's dbSNP
+sub dbsnp_info($$$$){
+my ($chr,$pos,$ref,$variant) = @_;
+my $key = "$chr:$pos:$ref:$variant";
+if(exists $dbsnp_info{$key}){
+return @{$dbsnp_info{$key}};
+}
+if(exists $chr2dbsnp_vcf_lines{$chr}->{$pos}){
+#print STDERR "Checking existing SNP data for $chr:$pos -> ", join("\n", @{$chr2dbsnp_vcf_lines{$chr}->{$pos}}), "\n";
+for(@{$chr2dbsnp_vcf_lines{$chr}->{$pos}}){
+	my @fields = split /\t/, $_;
+	for my $var (split /,/, $fields[4]){
+# Allows for different reference seqs between dbSNP and input, assuming patches only
+	    if(length($fields[3]) == length($ref) and ($var eq $variant or $ref eq $var and $variant eq $fields[3])){
+		my ($freq, $subpop) = ("","");
+		$freq = $1 if $fields[7] =~ /(?:\A|;)MMAF=(0\.\d+)(?:;|\Z)/;
+		$subpop = $1 if $fields[7] =~ /(?:\A|;)MMAF_SRC=(\S+?)(?:;|\Z)/;
+$dbsnp_info{$key} = [$subpop, $freq || "NA", $fields[2]];
+		return @{$dbsnp_info{$key}};
+	    }
+	}
+}
+}
+$dbsnp_info{$key} = ["novel", "NA", "NA"];
+return @{$dbsnp_info{$key}};
+}
+sub record_snv{
+my $line = join("", @_);
+push @snvs, $line;
+my @fields = split /\t/, $line;
+my $prop_info_key = $fields[9];
+my ($chr,$pos,$ref,$variant) = split /:/, $prop_info_key;
+$chr2variant_locs{$chr} = {} unless exists $chr2variant_locs{$chr};
+return unless $ref; # ref not defined for CNVs
+# Need to grab whole range for MNPs
+for(my $i = 0; $i < length($ref); $i++){
+$chr2variant_locs{$chr}->{$pos+$i} = 1;
+}
+}
+sub retrieve_vcf_lines($$$){
+my ($dbsnp_file, $internal_snp_file, $chr) = @_;
+my (%dbsnp_lines, %internal_snp_lines);
+if(not defined $dbsnp_file or not exists $chr2variant_locs{$chr}){
+return ({}, {}, {}, {}); # no data requested for this chromosome
+}
+# build up the request
+my @tabix_regions;
+my @var_locs = keys %{$chr2variant_locs{$chr}};
+# sort by variant start location
+for my $var_loc (sort {$a <=> $b} @var_locs){
+push @tabix_regions, $chr.":".$var_loc."-".$var_loc;
+}
+for(my $i = 0; $i <= $#tabix_regions; $i += $max_args){ # chunkify tabix request if too many for the system to handle
+my $end = $i + $max_args > $#tabix_regions ? $#tabix_regions : $i + $max_args;
+my $regions = "'".join("' '", @tabix_regions[$i..$end])."'";
+# From file is very slow for some reason
+#my $regions_file =  "/tmp/vcf2hgvs_$$.bed";
+#open(REQ_BED, ">$regions_file")
+#  or die "Cannot open $regions_file for writing: $!\n";
+#print REQ_BED join("\n", @tabix_regions), "\n";
+#close(REQ_BED);
+# retrieve the data
+die "Cannot find dbSNP VCF file $dbsnp_file\n" if not -e $dbsnp_file;
+open(VCF, "tabix $dbsnp_file $regions |")
+or die "Cannot run tabix on $dbsnp_file (args ".substr($regions, 0, length($regions)>100? 100 : length($regions))."): $!\n";
+while(<VCF>){
+#if(/^(\S+\t(\d+)(?:\t\S+){6})/ and grep {$_ eq $2} @var_locs){ # take only main columns to save room, if possible
+if(/^(\S+\t(\d+)(?:\t\S+){6})/ and exists $chr2variant_locs{$chr}->{$2}){ # take only main columns to save room, if possible
+$dbsnp_lines{$2} = [] unless exists $dbsnp_lines{$2};
+push @{$dbsnp_lines{$2}}, $1;
+}
+}
+close(VCF);
+if($internal_snp_file){
+die "Cannot find internal VCF file $internal_snp_file\n" if not -e $internal_snp_file;
+open(VCF, "tabix $internal_snp_file $regions |")
+or die "Cannot run tabix on $internal_snp_file: $!\n";
+while(<VCF>){
+#if(/^(\S+\t(\d+)(?:\t\S+){6})/ and grep {$_ eq $2} @var_locs){ # take only main columns to save room, if possible
+if(/^(\S+\t(\d+)(?:\t\S+){5})/ and exists $chr2variant_locs{$chr}->{$2}){ # take only main columns to save room, if possible
+$internal_snp_lines{$2} = [] unless exists $internal_snp_lines{$2};
+push @{$internal_snp_lines{$2}}, $1;
+}
+}
+close(VCF);
+}
+}
+#unlink $regions_file;
+return (\%dbsnp_lines, \%internal_snp_lines);
+}
+sub prop_info_key{
+my($chr,$pos,$ref,$variant,$exon_edge_dist) = @_;
+$chr =~ s/^chr//;
+if($chr eq "M"){
+$chr = "MT"; # NCBI uses different name for mitochondrial chromosome
+$pos-- if $pos >= 3107;  # also, doesn't keep the old positioning (historical)
+}
+return join(":", $chr,$pos,$ref,$variant, ($exon_edge_dist ? $exon_edge_dist : ""));
+}
+sub prop_info($$$){
+my($snpfile,$internal_snps_file,$prop_info_key) = @_;
+my($chr,$pos,$ref,$variant) = split /:/, $prop_info_key;
+# is this the first call for this chromosome? If so, retrieve the VCF lines for it en masse
+if(not exists $chr2dbsnp_vcf_lines{$chr}){
+($chr2dbsnp_vcf_lines{$chr}, $chr2internal_vcf_lines{$chr}) = retrieve_vcf_lines($snpfile,$internal_snps_file,$chr);
+}
+my $internal_maf = 0;
+if($internal_snps_file){
+$internal_maf = internal_prop($chr,$pos,$ref,$variant);
+$internal_maf = 0 if $internal_maf eq "NA";
+}
+my @results = dbsnp_info($chr,$pos,$ref,$variant);
+# Not all entries have a proportion in dbSNP
+return $internal_snps_file ? ($ref, $variant, @results, $internal_maf) : ($ref, $variant, @results);
+}
+#offset a given HGVS nomenclature position (single position only) by a given number of bases
+sub hgvs_plus($$){
+my ($hgvs, $offset) = @_;
+if($hgvs =~ /^(\S+)(-\d+)(.*)/){
+	# all negative
+	if($2+$offset<0){
+	    return $1.($2+$offset).$3;
+	}
+	# switches to positive, need to mod
+	else{
+	    return $1+($2+$offset);
+	}
+}
+elsif($hgvs =~ /^(\S+)\+(\d+)(.*)/){
+	# all positive
+	if($2+$offset>0){
+	    return $1."+".($2+$offset).$3;
+	}
+	# switches to negative, need to mod
+	else{
+	    return $1+($2+$offset);
+	}
+}
+elsif($hgvs =~ /^(-?\d+)(.*)/){
+	# special case if offset spans -/+ since there is no position 0
+	if($1 < 0 and $1+$offset >= 0){
+	    $offset++;
+	}
+	elsif($1 > 0 and $1+$offset <= 0){
+	    $offset--;
+	}
+	return ($1+$offset).$2;
+}
+else{
+	die "Cannot convert $hgvs to a new offset ($offset), only single base position nomenclature is currently supported\n";
+}
+}
+# offset a given position by a given number of bases,
+# taking into account that if the new offset crosses the threshold in the last argument,
+# HGVS boundary nomenclature has to be introduced
+sub hgvs_plus_exon($$$){
+my ($pos, $offset, $boundary) = @_;
+# special case if offset spans -/+ since there is no position 0
+if($pos =~ /^(-?\d+)(.*)/){
+if($1 < 0 and $1+$offset >= 0){
+$offset++;
+}
+elsif($1 > 0 and $1+$offset <= 0){
+$offset--;
+}
+}
+my $new_pos = $pos + $offset;
+if($new_pos > $boundary and $pos <= $boundary){
+	# just moved into an intron 3'
+	$new_pos = $boundary."+".($new_pos-$boundary);
+}
+elsif($new_pos < $boundary and $pos >= $boundary){
+	# just moved into an intron 5'
+	$new_pos = $boundary.($new_pos-$boundary);
+}
+return $new_pos;
+}
+# given a nucleotide position, calculates the AA there (assumes coding region)
+sub getCodonFromSeq($$$$){
+my ($chr_ref, $location, $frame_offset, $strand) = @_;
+my $codon;
+if($strand eq "+"){
+$codon = substr($$chr_ref, $location-1-$frame_offset, 3);
+}
+else{
+$codon = substr($$chr_ref, $location-3+$frame_offset, 3);
+$codon = reverse($codon);
+$codon =~ tr/ACGTacgt/TGCAtgca/;
+}
+return $codon;
+}
+sub getCodonFromSeqIndex($$$$){
+my ($chr, $location, $frame_offset, $strand) = @_;
+my $codon;
+if($strand eq "+"){
+$codon = $fasta_index->fetch($chr.":".($location-$frame_offset)."-".($location-$frame_offset+2));
+}
+else{
+$codon = $fasta_index->fetch($chr.":".($location-2+$frame_offset)."-".($location+$frame_offset));
+$codon = reverse($codon);
+$codon =~ tr/ACGTacgt/TGCAtgca/;
+}
+return $codon;
+}
+sub getAAFromSeq($$$$$){
+return $_[4]->translate(getCodonFromSeq($_[0], $_[1], $_[2], $_[3]));
+}
+sub getAAFromSeqIndex($$$$$){
+# convert codon to AA
+if(exists $transl_except{"$_[0]:$_[1]"}){
+return $transl_except{"$_[0]:$_[1]"};
+}
+else{
+return $_[4]->translate(getCodonFromSeqIndex($_[0], $_[1], $_[2], $_[3]));
+}
+}
+sub hgvs_protein{
+my ($chr, $location, $ref, $variant, $cdna_pos, $strand, $transl_table) = @_;
+if(substr($ref,0,1) eq substr($variant,0,1)){
+substr($ref,0,1) = "";
+substr($variant,0,1) = "";
+$location++;
+if($strand eq "-"){
+$cdna_pos--;
+}
+else{
+$cdna_pos++;
+}
+}
+if($cdna_pos !~ /^\d+/){
+die "Aborting: got illegal cDNA position ($cdna_pos) for protein HGVS conversion of position ",
+"$location, ref $ref, variant $variant. Please correct the program code.\n";
+}
+# Get the correct frame for the protein translation, to know what codons are affected
+my $aapos = int(($cdna_pos-1)/3)+1;
+# does it destroy the start codon?
+if($cdna_pos < 4){ # assumes animal codon usage
+return "p.0?"; # indicates start codon missing, unsure of effect
+}
+my $table = $transl_table ne $default_transl_table ?  # non standard translation table requested
+Bio::Tools::CodonTable->new(-id=>$transl_table) : $codonTable;
+my $frame_offset = ($cdna_pos-1)%3;
+my $origAA = getAAFromSeqIndex($chr, $location, $frame_offset, $strand, $table);
+# take 100000 bp on either side for translation context of variant seq
+my $five_prime_buffer = $location < 10000 ? $location-1 : 10000;
+my $mutSeq = $fasta_index->fetch($chr.":".($location-$five_prime_buffer)."-".($location+10000));
+# substitute all of the immediately adjacent variants in phase with this one to get the correct local effect
+substr($mutSeq, $five_prime_buffer, length($ref)) = $variant;
+# does it cause a frameshift?
+my $length_diff = length($variant)-length($ref);
+if($length_diff%3){ # insertion or deletion not a multiple of three
+my $fs_codon = getCodonFromSeq(\$mutSeq, $five_prime_buffer+1, $frame_offset, $strand);
+my $ext = 0;
+my $newAA;
+do{
+$ext++;
+# The "NA"s below make it so that we don't pick up any translation exceptions from the original reference annotation
+if($strand eq "+"){
+$newAA = getAAFromSeq(\$mutSeq, $five_prime_buffer+1+$ext*3, $frame_offset, $strand, $table);
+}
+else{
+$newAA = getAAFromSeq(\$mutSeq, $five_prime_buffer+1-$ext*3, $frame_offset, $strand, $table);
+}
+} while($newAA ne "*");
+return "p.".$origAA.$aapos.$table->translate($fs_codon)."fs*$ext";
+}
+# does it cause a stop codon to be lost?
+if($origAA eq "*"){
+my $stopChangeCodon = getCodonFromSeq(\$mutSeq, $five_prime_buffer+1, $frame_offset, $strand);
+# still a stop after the mutation (ignore translation exceptions)
+if($table->is_ter_codon($stopChangeCodon)){
+return "p.*$aapos=";
+}
+# calculate the new stop, assuming there aren't mutations downstream in candidate stop codons
+my $ext = 0;
+my $newCodon;
+do{
+if($strand eq "+"){
+$newCodon = getCodonFromSeq(\$mutSeq, $five_prime_buffer+1+(++$ext*3), $frame_offset, $strand);
+}
+else{
+$newCodon = getCodonFromSeq(\$mutSeq, $five_prime_buffer+1-(++$ext*3), $frame_offset, $strand);
+}
+} while(not $table->is_ter_codon($newCodon));
+return "p.*".$aapos.$table->translate($stopChangeCodon)."ext*".$ext;
+}
+# if we get this far, it's a "regular" AA level change
+my $origAAs = "";
+for(my $i = 0; $i < length($ref)+$frame_offset; $i+=3){
+my $oldAA = getAAFromSeqIndex($chr, $location+$i, $frame_offset, $strand, $table);
+if($strand eq "+"){
+$origAAs .= $oldAA;
+}
+else{
+$origAAs = $oldAA . $origAAs;
+}
+}
+my $newAAs = "";
+for(my $i = 0; $i < length($variant)+$frame_offset; $i+=3){
+# NA means we don't take translation exceptions from the original
+my $newAA = getAAFromSeq(\$mutSeq, $five_prime_buffer+1+$i, $frame_offset, $strand, $table);
+if($strand eq "+"){
+$newAAs .= $newAA;
+}
+else{
+$newAAs = $newAA . $newAAs;
+}
+}
+# silent
+if($origAAs eq $newAAs){
+return "p.".$origAAs.$aapos."=";
+}
+# minimize the difference if there are leading or trailing AAs the same
+my $delLength = length($ref);
+while(substr($newAAs, 0, 1) eq substr($origAAs, 0, 1)){
+$newAAs = substr($newAAs, 1);
+$origAAs = substr($origAAs, 1);
+$location+=3;
+$delLength-=3;
+$aapos++;
+}
+while(substr($newAAs, -1) eq substr($origAAs, -1)){
+$newAAs = substr($newAAs, 0, length($newAAs)-1);
+$origAAs = substr($origAAs, 0, length($origAAs)-1);
+}
+# insertion
+if(length($origAAs) == 0){
+my $insAAs = getAAFromSeqIndex($chr,$location-3,$frame_offset,$strand,$table).($aapos-1)."_".
+getAAFromSeqIndex($chr,$location,$frame_offset,$strand,$table);
+return "p.".$insAAs.$aapos."ins".$newAAs;
+}
+# deletion
+elsif(length($newAAs) == 0){
+my $delAAs;
+if(length($origAAs) == 1){
+$delAAs = getAAFromSeqIndex($chr,$location,$frame_offset,$strand,$table).$aapos; # single AA deletion
+}
+else{ # deleting a stretch
+if($strand eq "+"){
+my $endPoint = $location+$delLength-1;
+$delAAs = getAAFromSeqIndex($chr,$location,$frame_offset,$strand,$table).$aapos."_".
+getAAFromSeqIndex($chr,$endPoint,$frame_offset,$strand,$table).($aapos+int(($delLength-1)/3));
+}
+else{
+my $endPoint = $location-$delLength+1;
+$delAAs = getAAFromSeqIndex($chr,$endPoint,$frame_offset,$strand,$table).($aapos-int(($delLength-1)/3))."_".
+getAAFromSeqIndex($chr,$location,$frame_offset,$strand,$table).$aapos;
+}
+}
+return "p.".$delAAs."del";
+}
+else{
+# substitution
+if(length($origAAs) == 1 and length($newAAs) == 1){
+return "p.".$origAAs.$aapos.$newAAs;
+}
+# indel
+elsif(length($origAAs) != 1){
+# convert ref stretch into range syntax
+if($strand eq "+"){
+$origAAs = substr($origAAs, 0, 1).$aapos."_".substr($origAAs, -1).($aapos+length($origAAs)-1);
+}
+else{
+$origAAs = substr($origAAs, 0, 1).($aapos-length($origAAs)+1)."_".substr($origAAs, -1).$aapos;
+}
+}
+return "p.".$origAAs."delins".$newAAs;
+}
+return ("NA", "");
+}
+sub z2p{
+if(not defined $zed){
+$zed = new Statistics::Zed;
+}
+my $p = $zed->z2p(value => $_[0]);
+return $p < 0.0000000001 ? 0 : $p;
+}
+sub gq2p{
+return $_[0] > 200 ? 0 : 10**($_[0]/-10);
+}
+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);
+$multi_phased = 0;
+$min_depth = 2;
+$flanking_bases = 30;
+$min_pvalue = 0.01;
+$min_prop = 0.14;
+$rare_variant_prop = 0.05;
+$input_file = "-"; # STDIN by default
+$output_file = "-"; # STDOUT by default
+$default_transl_table = "1"; # assumes NCBI 'Standard' table, unless it is an argument to the script...
+&GetOptions("d=i"       => \$min_depth,
+"f=i"       => \$flanking_bases,
+"s=s"       => \$dbsnp,
+"t=s"       => \$dir_1000G,
+"n=s"       => \$dir_esp6500,
+"u=s"       => \$internal_snp,
+"q"         => \$quiet,
+"p=f"       => \$min_pvalue,
+"h=f"       => \$min_prop,
+"m=s"       => \$mappability_file,
+"r=s"       => \$reference_file,
+"z=s"       => \$samtools_phasing_file,
+"e=s"       => \$exons_file,
+"i=s"       => \$input_file,
+"c=s"       => \$cnv_file,
+"g=s"       => \$dgv_file,
+"b=s"       => \$genename_bed_file,
+"w=s"       => \$which_chr,
+"o=s"       => \$output_file,
+"a=i"       => \$default_transl_table,
+"v=f"       => \$rare_variant_prop,
+"x=s"       => \$enrichment_regions_file); # if enrichment regions are specified, variants without a transcript model but in these ranges will be reported
+if(($input_file ne "/dev/null" and not defined $reference_file) or
+not defined $exons_file or
+(defined $cnv_file and not defined $dgv_file)){
+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",
+"  [-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",
+"  [-t(housand) genomes_integrated_vcfs_gz_dir] [-n ESP6500_dir] [-u(ser) specified_population.vcf.gz] [-m(appability) crg_file.bed]\n",
+"  [-x enrichment_regions_file.bed] [-a(mino) acid translation table number from NCBI]\n",
+"  [-i(nput) genotypes.vcf <-r(eference) sequence_file.fasta>] [-o(utput) hgvs_file.tsv] [-s(np) database_from_ncbi.vcf.gz]\n",
+"  <-b(ed) file of named gene regions.bed> <-e(xons) file.gtf>\n\n",
+"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";
+}
+print STDERR "Considering $flanking_bases flanking bases for variants as well\n" unless $quiet;
+$codonTable = new Bio::Tools::CodonTable(id => $default_transl_table);
+my %enrichment_regions;
+# Note, we assume the regions are non-overlapping
+if(defined $enrichment_regions_file){
+print STDERR "Loading enrichment regions...\n" unless $quiet;
+open(BED, $enrichment_regions_file)
+or die "Cannot open $enrichment_regions_file for reading: $!\n";
+while(<BED>){
+chomp;
+my @F = split /\t/, $_;
+$enrichment_regions{$F[0]} = [] if not exists $enrichment_regions{$F[0]};
+push @{$enrichment_regions{$F[0]}}, [$F[1], $F[2]];
+}
+close(BED);
+}
+for my $chr (keys %enrichment_regions){ # sort by start
+$enrichment_regions{$chr} = [sort {$a->[0] <=> $b->[0]} @{$enrichment_regions{$chr}}];
+}
+if(defined $reference_file){
+print STDERR "Scanning reference FastA info\n" unless $quiet;
+if(not -e $reference_file){
+die "Reference FastA file ($reference_file) does not exist.\n";
+}
+if(not -e $reference_file.".fai" and not -w dirname($reference_file)){
+die "Reference FastA file ($reference_file) is not indexed, and the directory is not writable.\n";
+}
+$fasta_index = Bio::DB::Sam::Fai->load($reference_file);
+}
+my %chr2mappability;
+if(defined $mappability_file){
+print STDERR "Reading in mappability data\n" unless $quiet;
+my ($nmer) = $mappability_file =~ /(\d+).*?$/;
+die "Cannot determine nmer from nmer file name $mappability_file, aborting\n" unless $nmer;
+open(MAP, $mappability_file)
+or die "Cannot open mappability data file $mappability_file for reading: $!\n";
+<MAP>; # header
+while(<MAP>){
+next if /^#/;
+chomp;
+my @F = split /\t/, $_;
+my $x = int(1/$F[3]+0.5);
+$chr2mappability{$F[0]} = Set::IntervalTree->new() if not exists $chr2mappability{$F[0]};
+$chr2mappability{$F[0]}->insert("non-unique mapping region (x$x)", $F[1], $F[2]+$nmer-1);
+}
+close(MAP);
+}
+# Is phasing data provided?
+if(defined $samtools_phasing_file){
+print STDERR "Reading in phasing data\n" unless $quiet;
+open(PHASE, $samtools_phasing_file)
+or die "Cannot open phasing data file $samtools_phasing_file for reading: $!\n";
+my $phase_range;
+while(<PHASE>){
+if(/^PS/){
+chomp;
+my @F = split /\t/, $_;
+$phase_range = "$F[2]-$F[3]";
+}
+if(/^M[12]/){
+chomp;
+my @F = split /\t/, $_;
+#ignore strange cases where haplotype reference has no cases (weird samtools call)
+next if $F[9] == 0 or $F[7] == 0;
+my $chr = $F[1];
+next if defined $which_chr and not $chr eq $which_chr;
+my $pos = $F[3];
+#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;
+if(($F[10]+$F[8])/($F[9]+$F[7]) >= $min_prop){ # error meets reporting threshold
+$chr2caveats{"$chr:$pos"} .= "; " if exists $chr2caveats{"$chr:$pos"};
+$chr2caveats{"$chr:$pos"} .= "inconsistent haplotype phasing";
+}
+else{ # appears to be a genuine phasing
+$chr2phase{"$chr:$pos:$F[4]"} = "A-$chr:$phase_range"; # grouping for haplotype
+$chr2phase{"$chr:$pos:$F[5]"} = "B-$chr:$phase_range"; # grouping for haplotype
+}
+}
+}
+close(PHASE);
+}
+# Check the VCF file to see if contains phase data
+open(VCFIN, $input_file)
+or die "Cannot open $input_file for reading: $!\n";
+my $phase_chr = "";
+my @phase_dataA;
+my @phase_dataB;
+while(<VCFIN>){
+if(/^\s*(?:#|$)/){ # blank or hash comment
+next;
+}
+my @F = split /\t/, $_;
+next if exists $chr2caveats{"$F[0]:$F[1]"} and $chr2caveats{"$F[0]:$F[1]"} =~ /inconsistent haplotype phasing/;
+# | indicates phased
+if($F[8] =~ m(^(\d+)\|(\d+):)){
+next if $1 eq $2;  # not useful to us (actually would mess up phase combining later on), but is provided sometimes
+# start of a phasing block
+if($phase_chr eq ""){
+$phase_chr = $F[0];
+}
+my @vars = split /,/, $F[4];
+if($1 > @vars){
+die "Invalid VCF file (line #$.): First haplotype listed as $1, but only ", scalar(@vars), " variants were provided (", join(",", @vars), "\n";
+}
+if($2 > @vars){
+die "Invalid VCF file (line #$.): Second haplotype listed as $1, but only ", scalar(@vars), " variants were provided (", join(",", @vars), "\n";
+}
+unshift @vars, $F[3];
+push @phase_dataA, [$F[1], $vars[$1]];
+push @phase_dataB, [$F[1], $vars[$2]];
+}
+# non phased het call, ends any phasing block there might be
+elsif($F[8] =~ m(^0/1)){
+# Did we just finish a phased block? If so, output it.
+if(@phase_dataA > 1){
+my $phase_def = "G-$phase_chr:".$phase_dataA[0]->[0]."-".$phase_dataA[$#phase_dataA]->[0];
+for my $d (@phase_dataA){
+my ($p, $v) = @$d;
+if(exists $chr2phase{"$phase_chr:$p:$v"}){
+$chr2phase{"$phase_chr:$p:$v"} .= ",$phase_def";
+$multi_phased ||= 1;
+}
+else{
+$chr2phase{"$phase_chr:$p:$v"} = $phase_def;
+}
+}
+$phase_def = "H-$phase_chr:".$phase_dataB[0]->[0]."-".$phase_dataB[$#phase_dataB]->[0];
+for my $d (@phase_dataB){
+my ($p, $v) = @$d;
+if(exists $chr2phase{"$phase_chr:$p:$v"}){
+$chr2phase{"$phase_chr:$p:$v"} = ",$phase_def";
+$multi_phased ||= 1;
+}
+else{
+$chr2phase{"$phase_chr:$p:$v"} = $phase_def;
+}
+}
+}
+if($phase_chr ne ""){
+$phase_chr = "";
+@phase_dataA = ();
+@phase_dataB = ();
+}
+}
+}
+print STDERR "Reading in feature GTF data..." unless $quiet;
+my %feature_range; # chr => transcript_id => [[genomic_exon_start,genomic_exon_end,cdna_start_pos],...]
+my %feature_intervaltree; # chr => transcript_id => [[genomic_exon_start,genomic_exon_end,cdna_start_pos],...]
+my %feature_strand; # transcript_id => +|-
+my $feature_count = 0;
+my %feature_min;
+my %feature_max;
+my %feature_cds_min;
+my %feature_cds_max;
+my %feature_contig;
+my %feature_length;
+my %feature_type;
+my %feature_transl_table; # note alternate translation table usage
+my %chr_read;
+open(GTF, $exons_file)
+or die "Cannot open $exons_file for reading: $!\n";
+while(<GTF>){
+next if /^\s*#/;
+my @fields = split /\t/, $_;
+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";
+if($fields[2] eq "exon" or $fields[2] eq "CDS"){
+	next unless $fields[$#fields] =~ /transcript_id \"(.*?)\"/o;
+	my $parent = $1;
+	if(not $quiet and not exists $chr_read{$fields[0]}){
+	    print STDERR " $fields[0]";
+	    $chr_read{$fields[0]} = 1;
+	}
+	if(not exists $feature_strand{$parent}){
+	    $feature_strand{$parent} = $fields[6];
+	    $feature_contig{$parent} = $fields[0];
+if($fields[$#fields] =~ /transcript_type \"(.*?)\"/){
+$feature_type{$parent} = $1;
+}
+else{
+$feature_type{$parent} = "NA";
+}
+	}
+	if($fields[2] eq "CDS"){
+	    #print STDERR "CDS value for $parent is $fields[2]..$fields[3]\n";
+	    if(not exists $feature_cds_min{$parent} or $fields[3] < $feature_cds_min{$parent}){
+		$feature_cds_min{$parent} = $fields[3];
+	    }
+	    if(not exists $feature_cds_max{$parent} or $fields[4] > $feature_cds_max{$parent}){
+		$feature_cds_max{$parent} = $fields[4];
+	    }
+if($fields[$#fields] =~ /transl_table \"(\d+)\"/){
+$feature_transl_table{$parent} = $1; #assume one translation table per CDS, which should be reasonable
+}
+while($fields[$#fields] =~ /transl_except \"pos:(\S+?),aa:(\S+?)\"/g){
+my $pos = $1;
+my $new_aa = $2; # needs to change from three letter code to 1
+if($new_aa =~ /^ter/i){ # can be funny so have special case (allows TERM, etc.)
+$new_aa = "*";
+}
+elsif(length($new_aa) == 3){
+$new_aa = Bio::SeqUtils->new()->seq3in($new_aa);
+}
+if($pos =~ /^(\d+)\.\.(\d+)/){
+for my $p ($1..$2){
+$transl_except{"$fields[0]:$p"} = $new_aa;
+}
+}
+else{
+$transl_except{"$fields[0]:$pos"} = $new_aa;
+}
+}
+	    next;
+	}
+	if(not exists $feature_min{$parent} or $fields[3] < $feature_min{$parent}){
+	    $feature_min{$parent} = $fields[3];
+	}
+	if(not exists $feature_max{$parent} or $fields[4] > $feature_max{$parent}){
+	    $feature_max{$parent} = $fields[4];
+	}
+	$feature_count++;
+	if(not exists $feature_range{$fields[0]}){
+	    $feature_range{$fields[0]} = {}; # Chr => {parentID => [start,stop]}
+$feature_intervaltree{$fields[0]} = Set::IntervalTree->new();
+	}
+	if(not exists $feature_range{$fields[0]}->{$parent}){
+	    $feature_range{$fields[0]}->{$parent} = [];
+	}
+push @{$feature_range{$fields[0]}->{$parent}}, [$fields[3],$fields[4]];
+$feature_intervaltree{$fields[0]}->insert($parent, $fields[3], $fields[4]+1); # ranges need to have positive length for module to work properly
+$feature_length{$parent} += $fields[4]-$fields[3]+1;
+}
+}
+close(GTF);
+print STDERR "\nFound $feature_count exons on ", scalar(keys %feature_range), " contigs in the GTF file\n" unless $quiet;
+for my $contig (keys %feature_range){
+for my $parent (keys %{$feature_range{$contig}}){
+	# sort by subrange start
+	my @feature_ranges = sort {$a->[0] <=> $b->[0]} @{$feature_range{$contig}->{$parent}};
+	$feature_range{$contig}->{$parent} = \@feature_ranges;
+	$feature_range{"chr".$contig}->{$parent} = \@feature_ranges if not $contig =~ /^chr/;
+	$feature_range{$1}->{$parent} = \@feature_ranges if $contig =~ /^chr(\S+)/;
+}
+}
+# Calculate the cDNA position of the leftmost (reference strand) base for each exon
+for my $contig (keys %feature_range){
+for my $parent (keys %{$feature_range{$contig}}){
+my @feature_ranges = @{$feature_range{$contig}->{$parent}};
+if($feature_strand{$parent} eq "-"){
+# set up utr offset for correct CDS coordinates
+my $feature_offset = 0;
+for(my $i = $#feature_ranges; $i >= 0; $i--){
+last if not $feature_cds_max{$parent};
+# exon is completely 5' of the start
+if($feature_ranges[$i]->[0] > $feature_cds_max{$parent}){
+$feature_offset -= $feature_ranges[$i]->[1]-$feature_ranges[$i]->[0]+1;
+}
+# exon with the cds start
+elsif($feature_ranges[$i]->[1] >= $feature_cds_max{$parent} and
+$feature_ranges[$i]->[0] <= $feature_cds_max{$parent}){
+$feature_offset += $feature_cds_max{$parent} - $feature_ranges[$i]->[1];
+last;
+}
+else{
+die "The CDS for $parent (on negative strand) ends downstream ",
+"($feature_cds_max{$parent}) of the an exon",
+" (", $feature_ranges[$i]->[0], "), which is illogical.  Please revise the GFF file provided.\n";
+}
+}
+for(my $i = $#feature_ranges; $i >= 0; $i--){
+$feature_offset += $feature_ranges[$i]->[1]-$feature_ranges[$i]->[0]+1;
+$feature_ranges[$i]->[2] = $feature_offset-1;
+}
+}
+else{ # positive strand
+# set up utr offset for correct CDS coordinates
+my $feature_offset = 0;
+for(my $i = 0; $i <= $#feature_ranges; $i++){
+last if not $feature_cds_min{$parent};
+# All 5' utr exon
+if($feature_ranges[$i]->[1] < $feature_cds_min{$parent}){
+$feature_offset -= $feature_ranges[$i]->[1]-$feature_ranges[$i]->[0]+1;
+}
+# exon with the cds start
+elsif($feature_ranges[$i]->[1] >= $feature_cds_min{$parent} and
+$feature_ranges[$i]->[0] <= $feature_cds_min{$parent}){
+$feature_offset -= $feature_cds_min{$parent} - $feature_ranges[$i]->[0];
+last;
+}
+else{
+die "The CDS for $parent starts upstream ($feature_cds_min{$parent}) of the first exon",
+" (", $feature_ranges[$i]->[0], "), which is illogical.  Please revise the GFF file provided.\n";
+}
+}
+# assign cDNA coords for each exon to the third array element
+for(my $i = 0; $i <= $#feature_ranges; $i++){
+$feature_ranges[$i]->[2] = $feature_offset;
+$feature_offset += $feature_ranges[$i]->[1]-$feature_ranges[$i]->[0]+1;
+}
+}
+}
+}
+print STDERR "Reading in gene name definitions...\n" unless $quiet;
+die "Data file $genename_bed_file does not exist, aborting.\n" if not -e $genename_bed_file;
+my %gene_ids;
+open(TAB, $genename_bed_file)
+or die "Cannot open gene name BED file $genename_bed_file for reading: $!\n";
+while(<TAB>){
+chomp;
+# format should be "chr start stop gene_name ..."
+my @fields = split /\t/, $_;
+next if $#fields < 3;
+my $c = $fields[0];
+if(not exists $gene_ids{$c}){
+$gene_ids{$c} = Set::IntervalTree->new();
+}
+$gene_ids{$c}->insert($fields[3], $fields[1], $fields[2]);
+}
+# Print output header
+open(OUT, ">$output_file")
+or die "Cannot open $output_file for writing: $!\n";
+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",
+	       "Ref base", "Obs base", "Pop. freq. source", "Pop. freq.", "Variant DB ID"), "\t",
+($internal_snp ? "Internal pop. freq.\t" : ""),
+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";
+# If there is CNV data, load it.
+# BED columns should be chr start stop caveats ploidy . ignored ignored r,g,b
+# The dot means the strand doesn't matter.
+# where the first five fields are required, others optional
+# 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)
+if(defined $cnv_file){
+print STDERR "Reading in CNV data...\n" unless $quiet;
+open(CNV, $cnv_file)
+or die "Cannot open $cnv_file for reading: $!\n";
+while(<CNV>){
+chomp;
+my @F = split /\t/,  $_, -1;
+if(@F < 5){
+print STDERR "Skipping unparseable line ($cnv_file #$.): $_\n";
+next;
+}
+my $ploidy = $F[4];
+my $cnv_chr = $F[0];
+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";
+my $cnv_start = $F[1];
+my $cnv_end = $F[2];
+my $p_value = "NA";
+if($F[3] =~ s/p-value=(\S+?)(?:;|$)//){
+$p_value = $1;
+next if $min_pvalue < $p_value;
+}
+# Report a variant line for each gene that is found in this CNV
+my $target_parents = $feature_intervaltree{$cnv_chr}->fetch($cnv_start, $cnv_end+1);
+my $caveats = "";
+if(@F == 9){
+my @parents_ploidy = split /,/, $F[8];
+if($parents_ploidy[2] == 0){ # neither parent affected
+if($ploidy < $parents_ploidy[0] and $ploidy < $parents_ploidy[1]){
+if($ploidy > 2){
+$caveats = "Polyploidy is less severe than in either unaffected parents";
+}
+# else: no caveats, this offspring has fewer copies than normally observed, or in unaffected parents
+elsif($ploidy < 2){
+if($parents_ploidy[0] == 2 and $parents_ploidy[1] == 2){
+$caveats = "De novo copy loss, unaffected parents are diploid";
+}
+else{
+$caveats = "Copy loss is greater than in either unaffected parent";
+}
+}
+}
+elsif($ploidy >= $parents_ploidy[0] and $ploidy <= $parents_ploidy[1] or
+$ploidy >= $parents_ploidy[1] and $ploidy <= $parents_ploidy[0]){
+$caveats = "Aneuploidy likely inherited from an unaffected parent";
+}
+elsif($ploidy > $parents_ploidy[0] and $ploidy > $parents_ploidy[1]){
+if($parents_ploidy[0] > 2){
+if($parents_ploidy[1] > 2){
+$caveats = "Lower polyploidy already exists in both unaffected parents";
+}
+else{
+$caveats = "Lower polyploidy already exists in unaffected father";
+}
+}
+else{
+if($parents_ploidy[1] > 2){
+$caveats = "Lower polyploidy already exists in unaffected mother";
+}
+# else no caveats, because both parents are "normal", yet we have polyploidy in the offspring
+else{
+$caveats = "De novo polyploidy, unaffected parents are diploid";
+}
+}
+}
+# else
+else{
+die "Oops! Error in program logic...how did we get here (unaffected parents)? $_";
+}
+}
+elsif($parents_ploidy[2] == 1){ # father affected
+if($ploidy == $parents_ploidy[1]){ # just like unaffected Mom
+if($ploidy > 2){
+if($ploidy == $parents_ploidy[0]){
+$caveats = "Same polyploidy present in both affected and unaffected parents";
+}
+else{
+$caveats = "Polyploidy inherited from unaffected mother";
+}
+}
+elsif($ploidy < 2){
+if($ploidy == $parents_ploidy[0]){
+$caveats = "Same copy loss in both affected and unaffected parents";
+}
+else{
+$caveats = "Copy loss is shared with unaffected mother";
+}
+}
+else{
+if($ploidy == $parents_ploidy[0]){
+# Why was this even reported? parents and child have diploid status...
+next;
+}
+$caveats = "Diploidy is shared with unaffected mother";
+}
+}
+elsif($ploidy > 2){ # polyploidy
+if($parents_ploidy[0] == 2){
+if($parents_ploidy[1] > 2){
+$caveats = "Unaffected mother has polyploidy (".$parents_ploidy[1]."x), but affected father is diploid";
+}
+elsif($parents_ploidy[1] == 2){
+$caveats = "Both unaffected mother and affected father are diploid";
+}
+else{
+$caveats = "Affected father is diploid, unaffected mother has copy loss (".$parents_ploidy[1]."x)";
+}
+}
+elsif($parents_ploidy[0] < 2){
+$caveats = "Polyploidy found, but affected father had copy loss (".$parents_ploidy[0]."x)";
+}
+elsif($ploidy < $parents_ploidy[1]){
+$caveats = "Polyploidy is less severe than in unaffected mother (".$parents_ploidy[1]."x), or affected father (".$parents_ploidy[0]."x)";
+}
+# past here the ploidy is great than in the unaffected mother
+elsif($parents_ploidy[1] < 2){
+$caveats = "Polyploidy is also severe in affected father (".$parents_ploidy[0]."x), but unaffected mother actually had copy loss (". $parents_ploidy[1]. "x)";
+}
+elsif($parents_ploidy[1] == 2){
+$caveats = "Polyploidy is also severe in affected father (".$parents_ploidy[0]."x), and mother is diploid";
+}
+elsif($ploidy < $parents_ploidy[0]){
+$caveats = "Polyploidy is less severe than in affected father (".$parents_ploidy[0]."x), but more severe than unaffected mother (". $parents_ploidy[1]. "x)";
+}
+elsif($ploidy > $parents_ploidy[0]){
+$caveats = "Polyploidy is more severe than in affected father (".$parents_ploidy[0]."x)";
+}
+else{
+$caveats = "Polyploidy is as severe as in affected father";
+}
+}
+elsif($ploidy == 2){
+# Don't report diploid status, any funny recombination should show up in large indel analysis
+next;
+}
+else{ # copies < 2
+if($ploidy == $parents_ploidy[0]){
+if($ploidy > $parents_ploidy[1]){
+$caveats = "Copy loss is the same as affected father, but less than unaffected mother (". $parents_ploidy[1]. "x)";
+}
+else{
+$caveats = "Copy loss is as severe as in affected father";
+}
+}
+elsif($ploidy > $parents_ploidy[0]){
+if($ploidy > $parents_ploidy[1]){
+if($parents_ploidy[1] == 0 and $parents_ploidy[0] == 0){
+$caveats = "Poor mapping, or Mendelian inheritence violation is severe: no copies of region in either parent, but present in offspring";
+}
+elsif($ploidy == 2){
+next; # child got best of both parents, ignore from CNV standpoint (may still have SNPs of course, or translocation, etc.)
+}
+else{
+$caveats = "Copy loss is less severe than in unaffected mother (".$parents_ploidy[1]."x), or affected father (".$parents_ploidy[0]."x)";
+}
+}
+# else: child has less copies than unaffected mom, but more than affected Dad
+else{
+if($parents_ploidy[1] > 2){
+$caveats = "Copy loss was more severe in affected father (".$parents_ploidy[0]."x), but unaffected mother had polyploidy (".$parents_ploidy[1]."x)";
+}
+elsif($parents_ploidy[1] == 2){
+$caveats = "Copy loss was more severe in affected father (".$parents_ploidy[0]."x), but unaffected mother was diploid";
+}
+else{ # unaffected has loss
+$caveats = "Copy loss is more severe than unaffect mother (".$parents_ploidy[1]."x), but less severe than affected father (".$parents_ploidy[0]."x)";
+}
+}
+}
+# past here, ploidy is less than affected father
+elsif($parents_ploidy[1] > 2){
+$caveats = "Copy loss is more severe than affected father (".$parents_ploidy[0]."x), and unaffected mother had polyploidy (".$parents_ploidy[1]."x)";
+}
+elsif($parents_ploidy[1] == 2){
+$caveats = "Copy loss is more severe than in affected father (".$parents_ploidy[0]."x)";
+}
+else{
+$caveats = "Copy loss is more severe than in both unaffect mother (".$parents_ploidy[1]."x), and affected father (".$parents_ploidy[0]."x)";
+}
+}
+}
+elsif($parents_ploidy[2] == 2){ # mother affected
+if($ploidy == $parents_ploidy[0]){ # just like unaffected Dad
+if($ploidy > 2){
+if($ploidy == $parents_ploidy[1]){
+$caveats = "Same polyploidy present in both affected and unaffected parents";
+}
+else{
+$caveats = "Polyploidy inherited from unaffected father";
+}
+}
+elsif($ploidy < 2){
+if($ploidy == $parents_ploidy[1]){
+$caveats = "Same copy loss in both affected and unaffected parents";
+}
+else{
+$caveats = "Copy loss is shared with unaffected father";
+}
+}
+else{
+if($ploidy == $parents_ploidy[1]){
+# Why was this even reported? parents and child have diploid status...
+next;
+}
+$caveats = "Diploidy is shared with unaffected father";
+}
+}
+elsif($ploidy > 2){ # polyploidy
+if($parents_ploidy[1] == 2){
+if($parents_ploidy[0] > 2){
+$caveats = "Unaffected father has polyploidy (".$parents_ploidy[0]."x), but affected mother is diploid";
+}
+elsif($parents_ploidy[0] == 2){
+$caveats = "Both unaffected father and affected mother are diploid";
+}
+else{
+$caveats = "Affected mother is diploid, unaffected father has copy loss (".$parents_ploidy[1]."x)";
+}
+}
+elsif($parents_ploidy[1] < 2){
+$caveats = "Polyploidy found, but affected mother had copy loss (".$parents_ploidy[1]."x)";
+}
+elsif($ploidy < $parents_ploidy[0]){
+$caveats = "Polyploidy is less severe than in unaffected father (".$parents_ploidy[0]."x), or affected mother (".$parents_ploidy[1]."x)";
+}
+# past here the ploidy is great than in the unaffected father
+elsif($parents_ploidy[0] < 2){
+$caveats = "Polyploidy is also severe in affected mother (".$parents_ploidy[1]."x), but unaffected father actually had copy loss (". $parents_ploidy[0]. "x)";
+}
+elsif($parents_ploidy[0] == 2){
+$caveats = "Polyploidy is also severe in affected mother (".$parents_ploidy[1]."x), and unaffected father is diploid";
+}
+elsif($ploidy < $parents_ploidy[1]){
+$caveats = "Polyploidy is less severe than in affected mother (".$parents_ploidy[1]."x), but more severe than unaffected father (". $parents_ploidy[0]. "x)";
+}
+elsif($ploidy > $parents_ploidy[1]){
+$caveats = "Polyploidy is more severe than in affected mother (".$parents_ploidy[1]."x)";
+}
+else{
+$caveats = "Polyploidy is as severe as in affected mother";
+}
+}
+elsif($ploidy == 2){
+# Don't report diploid status, any funny recombination should show up in large indel analysis
+next;
+}
+else{ # copies < 2
+if($ploidy == $parents_ploidy[1]){
+if($ploidy > $parents_ploidy[0]){
+$caveats = "Copy loss is the same as affected mother, but less than unaffected father (". $parents_ploidy[0]. "x)";
+}
+else{
+$caveats = "Copy loss is as severe as in affected mother";
+}
+}
+elsif($ploidy > $parents_ploidy[1]){
+if($ploidy > $parents_ploidy[0]){
+if($parents_ploidy[1] == 0 and $parents_ploidy[0] == 0){
+$caveats = "Poor mapping, or Mendelian inheritence violation is severe: no copies of region in either parent, but present in offspring";
+}
+elsif($ploidy == 2){
+next; # child got best of both parents, ignore from CNV standpoint (may still have SNPs of course, or translocation, etc.)
+}
+else{
+$caveats = "Copy loss is less severe than in unaffected father (".$parents_ploidy[0]."x), or affected mother (".$parents_ploidy[1]."x)";
+}
+}
+# else: child has less copies than unaffected Dad, but more than affected Mom
+else{
+if($parents_ploidy[0] > 2){
+$caveats = "Copy loss was more severe in affected mother (".$parents_ploidy[1]."x), but unaffected father had polyploidy (".$parents_ploidy[0]."x)";
+}
+elsif($parents_ploidy[0] == 2){
+$caveats = "Copy loss was more severe in affected mother (".$parents_ploidy[1]."x), but unaffected father was diploid";
+}
+else{ # unaffected has loss
+$caveats = "Copy loss is more severe than unaffect father (".$parents_ploidy[0]."x), but less severe than affected mother (".$parents_ploidy[1]."x)";
+}
+}
+}
+# past here, ploidy is less than affected mother
+elsif($parents_ploidy[0] > 2){
+$caveats = "Copy loss is more severe than affected mother (".$parents_ploidy[1]."x), and unaffected father had polyploidy (".$parents_ploidy[0]."x)";
+}
+elsif($parents_ploidy[0] == 2){
+$caveats = "Copy loss is more severe than in affected mother (".$parents_ploidy[1]."x)";
+}
+else{
+$caveats = "Copy loss is more severe than in both unaffect father (".$parents_ploidy[0]."x), and affected mother (".$parents_ploidy[1]."x)";
+}
+}
+}
+}
+if($F[3] and $F[3] ne "-"){ # prexisting caveat from CNV caller
+if(defined $caveats){
+$caveats .= "; $F[3]" unless $caveats =~ /\b$F[3]\b/;
+}
+else{
+$caveats  = $F[3];
+}
+}
+# Sort by start for consistency
+my @target_parents = sort {$feature_range{$cnv_chr}->{$a}->[0]->[0] <=> $feature_range{$cnv_chr}->{$b}->[0]->[0]} @$target_parents;
+for my $target_parent (@target_parents){
+my $target_caveats = $caveats;
+my $strand = $feature_strand{$target_parent};
+# report the gain/loss of each gene separately, for simplicity in downstream analysis
+my $cnv_exon_start = 10000000000; # genomic coords
+my $cnv_exon_end = 0;
+my $cnv_cdna_start = 0; # cDNA coords
+my $cnv_cdna_end = 0;
+my $off5 = 0; # border outside the exon?
+my $off3 = 0;
+my @feature_ranges = @{$feature_range{$cnv_chr}->{$target_parent}};
+# find the first and last exons in the gene that are inside the CNV
+for my $subregion (@feature_ranges){
+# exon overlaps CNV?
+if($subregion->[0] <= $cnv_end and $subregion->[1] >= $cnv_start){
+if($cnv_exon_start > $subregion->[0]){
+if($cnv_start < $subregion->[0]){
+$cnv_exon_start = $subregion->[0]; $off5 = 1;
+$cnv_cdna_start = $subregion->[2];
+}
+else{
+$cnv_exon_start = $cnv_start; $off5 = 0;
+$cnv_cdna_start = $subregion->[2]+($strand eq "-" ? $subregion->[0]-$cnv_start: $cnv_start-$subregion->[0]);
+}
+}
+if($cnv_exon_end < $subregion->[1]){
+if($cnv_end > $subregion->[1]){
+$cnv_exon_end = $subregion->[1]; $off3 = 1;
+$cnv_cdna_end = $subregion->[2]+($strand eq "-" ? $subregion->[0]-$subregion->[1] : $subregion->[1]-$subregion->[0]);
+}
+else{
+$cnv_exon_end = $cnv_end; $off3 = 0;
+$cnv_cdna_end = $subregion->[2]+($strand eq "-" ? $subregion->[0]-$cnv_end : $cnv_end-$subregion->[0]);
+}
+}
+}
+}
+my $ends_internally = 0;
+if($cnv_exon_end == 0){ # ends inside the exon
+$cnv_exon_end = $cnv_end;
+$ends_internally = 1;
+}
+# See if it's in the structural variant database
+my @gain_coverage; $#gain_coverage = $cnv_exon_end-$cnv_exon_start; # preallocate blanks
+my @loss_coverage; $#loss_coverage = $cnv_exon_end-$cnv_exon_start; # preallocate blanks
+my $dgv_loss_id; # report the DGV entry that covers most of the observed structural variant
+my $dgv_loss_length = 0; # report the DGV entry that covers most of the observed structural variant
+my $dgv_gain_id; # report the DGV entry that covers most of the observed structural variant
+my $dgv_gain_length = 0; # report the DGV entry that covers most of the observed structural variant
+my $gains;
+my $losses;
+my $dgv_chr = $cnv_chr;
+$dgv_chr =~ s/^chr//; # no prefix in DGV
+#open(DGV, "tabix $dgv_file $dgv_chr:$cnv_exon_start-$cnv_exon_end |") # check out CNV in this gene model region
+#  or die "Cannot run tabix: $!\n";
+open(DGV, "/dev/null");
+while(<DGV>){
+my @C = split /\t/, $_;
+next if $C[4] ne "CNV"; # todo: handle indels?
+my $dgv_start = $C[2];
+my $dgv_end = $C[3];
+my $dgv_direction = $C[5];
+my $gain_cov_count = 0;
+my $loss_cov_count = 0;
+if($dgv_direction eq "Gain"){
+for(my $i = ($dgv_start < $cnv_exon_start ? $cnv_exon_start : $dgv_start); $i <= $dgv_end and $i <= $cnv_exon_end; $i++){
+$gain_coverage[$i-$cnv_exon_start] = 1 unless defined $gain_coverage[$i-$cnv_exon_start];
+$gain_cov_count++;
+}
+}
+elsif($dgv_direction eq "Loss"){
+for(my $i = ($dgv_start < $cnv_exon_start ? $cnv_exon_start : $dgv_start); $i <= $dgv_end and $i <= $cnv_exon_end; $i++){
+$loss_coverage[$i-$cnv_exon_start] = 1 unless defined $loss_coverage[$i-$cnv_exon_start];
+$loss_cov_count++;
+}
+}
+if($dgv_direction eq "Gain" and $gain_cov_count > $dgv_gain_length){
+$dgv_gain_id = $C[0];
+$dgv_gain_length = $gain_cov_count;
+}
+if($dgv_direction eq "Loss" and $loss_cov_count > $dgv_loss_length){
+$dgv_loss_id = $C[0];
+$dgv_loss_length = $loss_cov_count;
+}
+}
+close(DGV);
+my $gain_coverage = 0;
+for my $count (@gain_coverage){
+$gain_coverage++ if defined $count;
+}
+$gain_coverage = sprintf "%.3f", $gain_coverage/($cnv_exon_end-$cnv_exon_start+1); # make it a proportion
+my $loss_coverage = 0;
+for my $count (@loss_coverage){
+$loss_coverage++ if defined $count;
+}
+$loss_coverage = sprintf "%.3f", $loss_coverage/($cnv_exon_end-$cnv_exon_start+1); # make it a proportion
+my $src = "DGV";
+my $dgv_id = "NA";
+my $dgv_caveat;
+my $dgv_coverage;
+if($ploidy > 2){
+if(not defined $dgv_gain_id){
+if(defined $dgv_loss_id){
+$dgv_id = sprintf "%s/%.3f", $dgv_loss_id, $dgv_loss_length/($cnv_exon_end-$cnv_exon_start+1);
+$dgv_caveat = "; No gains are known in healthy controls, the DGV2 ID reported is for a loss in the same area";
+$dgv_coverage = $loss_coverage;
+}
+else{
+$dgv_id = "novel";
+$dgv_coverage = "NA";
+$src = "NA";
+}
+}
+else{
+$dgv_id = sprintf "%s/%.3f", $dgv_gain_id, $dgv_gain_length/($cnv_exon_end-$cnv_exon_start+1);
+$dgv_coverage = $gain_coverage;
+}
+}
+elsif($ploidy < 2){
+if(not defined $dgv_loss_id){
+if(defined $dgv_gain_id){
+$dgv_id = sprintf "%s/%.3f", $dgv_gain_id, $dgv_gain_length/($cnv_exon_end-$cnv_exon_start+1);
+$dgv_caveat = "; No losses are known in healthy controls, the DGV2 ID reported is for a gain in the same area";
+$dgv_coverage = $gain_coverage;
+}
+else{
+$dgv_id = "novel";
+$dgv_coverage = "NA";
+$src = "NA";
+}
+}
+else{
+$dgv_id = sprintf "%s/%.3f", $dgv_loss_id, $dgv_loss_length/($cnv_exon_end-$cnv_exon_start+1);
+$dgv_coverage = $loss_coverage;
+}
+}
+my $non_coding = 0;
+if(not exists $feature_cds_max{$target_parent} or not defined $feature_cds_max{$target_parent} or $feature_cds_max{$target_parent} eq ""){
+$non_coding = 1;
+}
+$target_caveats .= $dgv_caveat if defined $dgv_caveat and $dgv_id ne "novel" and $target_caveats !~ /\Q$dgv_caveat\E/;
+#print "Recorded $cnv_chr:$cnv_start caveat $caveats\n";
+# if it doesn't overlap an exon, we need to find out which two exons it's between
+if($ends_internally){
+my $intron_found = 0;
+for(my $i = 0; $i < $#feature_ranges; $i++){
+if($feature_ranges[$i]->[1] < $cnv_start and $feature_ranges[$i+1]->[0] > $cnv_end){
+if($ploidy > 2){ # gain
+if($strand eq "-"){
+record_snv("$target_parent\t",
+($non_coding ? "g.$cnv_start\_$cnv_end" :
+"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])),
+# 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
+"[".($ploidy-1)."]\t$strand\t$cnv_chr\t$cnv_start\t$cnv_end\tNA\t$p_value\tNA\tNA\t",
+"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");
+}
+else{
+record_snv("$target_parent\t",
+($non_coding ? "g.$cnv_start\_$cnv_end" :
+"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)),
+"[".($ploidy-1)."]\t$strand\t$cnv_chr\t$cnv_start\t$cnv_end\tNA\t$p_value\tNA\tNA\t",
+"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");
+}
+}
+else{ # loss
+if($strand eq "-"){
+record_snv("$target_parent\t",
+($non_coding ? "g.$cnv_start\_$cnv_end" :
+"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])),
+"del\t$strand\t$cnv_chr\t$cnv_start\t$cnv_end\t", ($ploidy == 1 ? "heterozygote" : "homozygote"), "\t$p_value\tNA\tNA\t",
+"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");
+}
+else{
+record_snv("$target_parent\t",
+($non_coding ? "g.$cnv_start\_$cnv_end" :
+"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)),
+"del\t$strand\t$cnv_chr\t$cnv_start\t$cnv_end\t", ($ploidy == 1 ? "heterozygote" : "homozygote"), "\t$p_value\tNA\tNA\t",
+"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");
+}
+}
+$intron_found = 1; last;
+}
+}
+warn "Logic error: CNV overlaps a gene ($target_parent), but is neither intronic nor exonic. Offending CNV was $_\n" unless $intron_found;
+next;
+}
+if($strand eq "-"){
+my $tmp = $cnv_cdna_start;
+$cnv_cdna_start = $cnv_cdna_end;
+$cnv_cdna_end = $tmp;
+}
+# Make the location label pretty descriptive
+my $cnv_phase = "";
+if($cnv_exon_start > $cnv_start or $cnv_exon_end < $cnv_end){
+$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
+}
+# if we get here, we're in a gained/deleted exon category
+# CNVs are fuzzy-edged (as opposed to large indels), so produce HGVS syntax that reflect this
+if($ploidy > 2){ # gain
+# find the exons encompassed by the CNV. NOTE that we assume that polyploidies are proximal
+record_snv("$target_parent\t",
+($non_coding ? "g.".($cnv_exon_start > $cnv_start ? "$cnv_exon_start-?" : $cnv_start)."_".($cnv_exon_end < $cnv_end ? "$cnv_exon_end+?" : $cnv_end) :
+"c.$cnv_cdna_start".($off5?"-?":"")."_$cnv_cdna_end".($off3?"+?":"")),
+"[".($ploidy-1)."]\t$strand\t$cnv_chr\t$cnv_exon_start\t$cnv_exon_end\tNA\t$p_value\tNA\tNA\t",
+"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");
+}
+else{ # loss
+#translate genome coordinates into cDNA coordinates
+record_snv("$target_parent\t",
+($non_coding ? "g.".($cnv_exon_start > $cnv_start ? "$cnv_exon_start-?" : $cnv_start)."_".($cnv_exon_end < $cnv_end ? "$cnv_exon_end+?" : $cnv_end) :
+"c.$cnv_cdna_start".($off5?"-?":"")."_$cnv_cdna_end".($off3?"+?":"")),
+"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",
+"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");
+}
+}
+}
+close(CNV);
+}
+#sort genes by start, then longest if tied
+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);
+print STDERR "Processing variant calls..." unless $quiet;
+%chr_read = ();
+open(VCFIN, $input_file)
+or die "Cannot open $input_file for reading: $!\n";
+while(<VCFIN>){
+if(/^\s*(?:#|$)/){ # blank or hash comment lines
+next;
+}
+chomp;
+my @fields = split /\t/, $_;
+next unless exists $feature_range{$fields[0]};
+if(not $quiet and not exists $chr_read{$fields[0]}){
+	print STDERR " $fields[0]";
+	$chr_read{$fields[0]} = 1;
+	#print STDERR "(not in reference file!)" unless exists $feature_range{$fields[0]};
+}
+next if $fields[4] eq "<NON_REF>"; # GVCF background stuff
+next if $fields[9] eq "./." or $fields[9] eq "."; # no call
+my $chr = $fields[0];
+next if defined $which_chr and $chr ne $which_chr and $chr ne "chr$which_chr" and "chr$chr" ne $which_chr;
+print STDERR "passes chr and field # test" if grep /dataset_7684.dat/, @ARGV;
+$chr = "chr$chr" if $chr !~ /^chr/;
+$fields[8] =~ s/\s+$//;
+my @values = split /:/, $fields[9];
+#print STDERR join(" / ", @values), "\n" if $. == 84;
+my @keys = split /:/, $fields[8];
+my $zygosity = "n/a";
+my $quality = "n/a";
+my $read_depth = "n/a";
+my $variant_depths = "n/a";
+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
+	if($keys[$i] eq "GT"){
+if($values[$i] =~ /\./ or $values[$i] =~ /0\/0/){ # one genotype not described
+$zygosity = "none";
+last;
+}
+	    else{ # genotypes described
+$zygosity = $values[$i] =~ /[02]/ ? "heterozygote" : "homozygote";
+}
+	}
+elsif($keys[$i] eq "DNM_CONFIG" and $zygosity eq "n/a"){ # hack
+$zygosity = $values[$i] =~ /^(.)\1/ ? "homozygote" : "heterozygote";
+}
+	elsif($keys[$i] eq "GQ" and $values[$i] ne "."){
+#print "Checking GQ (index $i) $values[$i] gq2p\n" if $. == 84;
+	    $quality = gq2p($values[$i]);
+	}
+	elsif($keys[$i] eq "RD"){ # from some tools like denovo variant finders
+	    $read_depth = $values[$i];
+	}
+	elsif($keys[$i] eq "DP"){
+	    $read_depth = $values[$i];
+	}
+# the frequency of the variant can go by many names, here we have freebayes (A* are new and old versions) and atlas2_indel
+	elsif($keys[$i] eq "AA" or $keys[$i] eq "VR" or $keys[$i] eq "AO"){
+	    $variant_depths = $values[$i];
+	}
+# here we have GATK variant freq of form ref#,var#
+	elsif($keys[$i] eq "AD"){
+	    $variant_depths = $values[$i];
+$variant_depths =~ s/^\d+,//;
+	}
+else{
+#print STDERR "Ignoring field $keys[$i]\n";
+}
+}
+next if $zygosity eq "none"; # GVCF non-ref for example or when multiple samples are reported simultaneously
+$quality = z2p($1) if $fields[7] =~ /Z=(\d+\.\d+)/;
+if($quality eq "n/a" and $fields[5] ne "."){
+$quality = gq2p($fields[5]);
+}
+if($fields[5] eq "0" and $fields[6] eq "PASS"){ # not qual and a PASS in the filter column
+$quality = 1;
+}
+elsif($quality ne "n/a" and $quality > $min_pvalue){ # p-value cutoff
+#print "Checking call quality $fields[5]  gq2p\n" if $. == 84;
+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
+# in these cases it would be silly to reject the call if their are many supporting reads.
+}
+# Some tools like GATK don't report number of variant reads...infer from other data if possible
+if($variant_depths eq "n/a"){
+my @key_value_pairs = split /;/, $fields[7];
+for my $key_value_pair (@key_value_pairs){
+if($key_value_pair !~ /^(.*?)=(.*)$/){
+next;
+#next if $key_value_pair eq "INDEL"; # samtools peculiarity
+#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";
+}
+if($1 eq "AB"){ # GATK: for het calls, AB is ref/(ref+var), only one variant reported per line
+$variant_depths = "";
+for my $ab (split /,/, $2){
+$variant_depths .= int((1-$ab)*$read_depth).",";
+}
+chop $variant_depths;
+}
+elsif($1 eq "MLEAC"){
+}
+elsif($1 eq "DP4"){ # samtools: high-quality ref-forward bases, ref-reverse, alt-forward and alt-reverse bases
+my @rds = split /,/, $2;
+$variant_depths = $rds[2]+$rds[3];
+$read_depth = $rds[0]+$rds[1]+$variant_depths;
+if($fields[4] =~ /,/){ # samtools doesn't break down compound het calls into individual frequencies
+my $num_alt_genotypes = $fields[4] =~ tr/,/,/;
+$num_alt_genotypes++;
+my $even_prop = sprintf "%.2f", $variant_depths/$read_depth/$num_alt_genotypes;
+$variant_depths = join(",", ($even_prop) x $num_alt_genotypes);
+if(not exists $chr2caveats{"$chr:$fields[1]"}){
+$chr2caveats{"$chr:$fields[1]"} = "compound het var freq n/a in orig VCF file, auto set to $even_prop";
+}
+else{
+$chr2caveats{"$chr:$fields[1]"} .= "; compound het var freq n/a in orig VCF file, auto set to $even_prop";
+}
+}
+}
+}
+}
+if($variant_depths eq "n/a"){ # usually homo var calls, can only assume all reads are variant
+if($zygosity eq "homozygote"){
+$variant_depths = $read_depth;
+if($read_depth ne "n/a"){
+if(not exists $chr2caveats{"$chr:$fields[1]"}){
+$chr2caveats{"$chr:$fields[1]"} = "homo var freq n/a in orig VCF file, auto set to 1.0";
+}
+else{
+$chr2caveats{"$chr:$fields[1]"} = "; homo var freq n/a in orig VCF file, auto set to 1.0";
+}
+}
+}
+else{
+if($read_depth ne "n/a"){
+$variant_depths = int($read_depth/2);
+if(not exists $chr2caveats{"$chr:$fields[1]"}){
+$chr2caveats{"$chr:$fields[1]"} = "het var freq n/a in orig VCF file, auto set to 0.5";
+}
+else{
+$chr2caveats{"$chr:$fields[1]"} = "; het var freq n/a in orig VCF file, auto set to 0.5";
+}
+}
+else{
+$variant_depths = $read_depth;
+}
+}
+}
+my $target_parents = $feature_intervaltree{$chr}->fetch($fields[1]-$flanking_bases, $fields[1]+length($fields[3])+$flanking_bases);
+# Last ditch, if not in a gene model, is it at least in an enrichment region?
+if(not @$target_parents){
+next if not exists $enrichment_regions{$chr};
+my $regions_ref = $enrichment_regions{$chr};
+my $location = $fields[1];
+my $strand = "+"; # for lack of a better choice
+for(my $i = find_earliest_index($location-$flanking_bases, $regions_ref);
+$i < $#$regions_ref and $location <= $regions_ref->[$i]->[1]+$flanking_bases;
+$i++){
+next unless $regions_ref->[$i]->[0]-$flanking_bases <= $location and $regions_ref->[$i]->[1]+$flanking_bases >= $location;
+my $feature_name = "enrichment_target_$chr:".$regions_ref->[$i]->[0]."-".$regions_ref->[$i]->[1];
+$feature_type{$feature_name} = "misc_enrichment_kit_target";
+$feature_length{$feature_name} = $regions_ref->[$i]->[1]-$regions_ref->[$i]->[0]+1;
+my @variants = split /,/, $fields[4];
+my @variant_depths = split /,/, $variant_depths;
+my $ref = uc($fields[3]);
+for(my $j = 0; $j <= $#variants; $j++){
+my $variant = $variants[$j];
+next if $variant eq "<NON_REF>" or $variant_depths[$j] eq "0"; # GVCF stuff
+my $variant_depth = $variant_depths[$j];
+if($min_prop){
+next unless $variant_depth >= $min_depth and $read_depth ne "n/a" and $variant_depth/$read_depth >= $min_prop;
+}
+if(length($ref) == 1 and length($variant) == 1){ # SNP
+record_snv("$feature_name\tg.$location",
+"$ref>$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+join("\t",prop_info_key($chr,$location,$ref,$variant)),"\tNA\n");
+}
+elsif(length($ref) == 1 and length($variant) > 1){ # insertion
+record_snv("$feature_name\tg.$location\_",($location+1),
+"ins",substr($variant, 1),"\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+join("\t",prop_info_key($chr,$location,$ref,$variant)),"\tNA\n");
+}
+elsif(length($variant) == 1 and length($ref) > 1){ # deletion
+record_snv("$feature_name\tg.$location\_",($location+length($ref)-1),
+"del",substr($ref, 1),"\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+join("\t",prop_info_key($chr,$location,$ref,$variant)),"\tNA\n");
+}
+else{ # indel
+record_snv("$feature_name\tg.$location\_",($location+length($ref)-1),
+"delins$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+join("\t",prop_info_key($chr,$location,$ref,$variant)),"\tNA\n");
+}
+} # end for variants
+next; # process next record, we've done all we can with a non-coding-gene SNP
+}
+}
+for my $target_parent (@$target_parents){
+print STDERR "Checking parent $target_parent for on $chr:$fields[1] $fields[3] -> $fields[4]\n" if grep /dataset_7684.dat/, @ARGV;
+my @feature_ranges = @{$feature_range{$chr}->{$target_parent}};
+# Calculate the position of the change within the feature range position
+my $strand = $feature_strand{$target_parent};
+my $trans_table = exists $feature_transl_table{$target_parent} ? $feature_transl_table{$target_parent} : $default_transl_table;
+$fields[4]=~tr/"//d; # sometime strangely surroundsed by quotes
+my @variants = split /,/, $fields[4];
+my @variant_depths = split /,/, $variant_depths;
+my @refs = (uc($fields[3])) x scalar(@variants);
+my @locations = ($fields[1]) x scalar(@variants);
+for(my $j = 0; $j <= $#variants; $j++){
+	my $ref = $refs[$j];
+	my $location = $locations[$j];
+	my $feature_offset = 0;
+	my $feature_num = 0;
+	my $variant = uc($variants[$j]);
+next if $variant eq "<NON_REF>" or $variant_depths[$j] eq "0"; # GVCF stuff
+	my $variant_depth = $variant_depths[$j] || "n/a";
+#print STDERR "Evaluating target parent $target_parent for $chr:$fields[1]-".($fields[1]+length($fields[3]))." -> ",join("/", @$target_parents) , "\n" if $fields[1] == 982994;
+# 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)
+if(length($variant) > 1 and length($variant) == length($ref)){
+my @subvariants;
+my @subrefs;
+my @sublocations;
+my $phase_range = $location."-".($location+length($ref)-1);
+for(my $k = 0; $k < length($variant); $k++){
+my $r = substr($ref, $k, 1);
+my $v = substr($variant, $k, 1);
+if($r ne $v){
+push @subvariants, $v;
+push @subrefs, $r;
+push @sublocations, $location+$k;
+}
+elsif(@variants == 1){
+next; # homo ref call
+}
+if($zygosity eq "heterozygote"){
+# need to ignore case where a homozygous call (variant or ref) is included in a double non-ref het MNP
+if(@variants > 1){
+my $homo = 1;
+for(my $l = 0; $l <= $#variants; $l++){ # using loop in case we handle oligoploid genomes in the future
+if(length($variants[$l]) <= $k or substr($variants[$l], $k, 1) ne $v){
+$homo = 0;
+last;
+}
+}
+next if $homo;
+}
+my $phase_key = "$chr:".($location+$k).":$v";
+my $phase_label = "M$j-$chr:$phase_range";
+if(exists $chr2phase{$phase_key}){
+if($chr2phase{$phase_key} !~ /$phase_label/){
+$chr2phase{$phase_key} .= ",$phase_label";
+}
+}
+else{
+$chr2phase{$phase_key} = $phase_label;
+}
+}
+}
+# recycle this MNP variant loop to start processing the individual SNPs
+splice(@refs, $j, 1, @subrefs);
+splice(@variants, $j, 1, @subvariants);
+splice(@locations, $j, 1, @sublocations);
+splice(@variant_depths, $j, 1, ($variant_depth) x scalar(@subvariants));
+$j--;
+next;
+}
+if($min_prop != 0 and $variant_depth eq "n/a" or $variant_depth eq "."){
+print STDERR "Could not parse variant depth from $_\n" unless $quiet;
+next;
+}
+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;
+if($zygosity eq "NA"){
+# make the call ourselves
+$zygosity = $variant_depths/$read_depth > 1-$min_prop ? "homozygote" : "heterozygote";
+}
+	# e.g. chr22   47857671   . CAAAG   AAGAT,AAAAG    29.04  .
+	if(length($variant) > 1 and length($variant) == length($ref)){
+	    for(my $k = 0; $k < length($variant); $k++){
+		my $fixed_variant = $variant;
+		substr($fixed_variant, $k, 1) = substr($ref, $k, 1);
+		if($fixed_variant eq $ref){ # single base difference at base k between the two
+		    $ref = substr($ref, $k, 1);
+		    $variant = substr($variant, $k, 1);
+		    $location += $k;
+		    last;
+		}
+	    }
+	}
+# samtools reports indels with long common tails, causing non-canonical HGVS reporting and a problem when looking up the variant in dbSNP etc.
+# remove common tails to variant calls in order to try to rectify this
+else{
+while(length($ref) > 1 and length($variant) > 1 and substr($ref, length($ref)-1) eq substr($variant, length($variant)-1)){
+chop $ref; chop $variant;
+}
+}
+# See if a caveat should be added because the indel was in a polyhomomer region
+if(length($ref) > length($variant) and index($ref, $variant) == 0 and $fasta_index->fetch("$chr:".($location+1)."-".($location+length($ref)+1)) =~ /^([ACGT])\1+$/i){
+if(not exists $chr2caveats{"$chr:$location"}){
+$chr2caveats{"$chr:$location"} = "poly".uc($1)." region deletion";
+}
+elsif($chr2caveats{"$chr:$location"} !~ /poly/){
+$chr2caveats{"$chr:$location"} .= "; poly".uc($1)." region deletion";
+}
+}
+elsif(length($ref) < length($variant) and index($variant, $ref) == 0 and substr($variant, 1) =~ /^([ACGT])\1+$/i){
+if(not exists $chr2caveats{"$chr:$location"}){
+$chr2caveats{"$chr:$location"} .= "poly".uc($1)." region insertion";
+}
+elsif($chr2caveats{"$chr:$location"} !~ /poly/){
+$chr2caveats{"$chr:$location"} .= "; poly".uc($1)." region insertion";
+}
+}
+# Not a protein-coding gene?  Report genomic cooordinates for HGVS
+if(not exists $feature_cds_max{$target_parent} or not defined $feature_cds_max{$target_parent} or $feature_cds_max{$target_parent} eq ""){
+if(length($ref) == 1 and length($variant) == 1){ # SNP
+record_snv("$target_parent\tg.$location",
+"$ref>$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+join("\t",prop_info_key($chr,$location,$ref,$variant)),"\tNA\n");
+}
+elsif(length($ref) == 1 and length($variant) > 1){ # insertion
+record_snv("$target_parent\tg.$location\_",($location+1),
+"ins",substr($variant, 1),"\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+join("\t",prop_info_key($chr,$location,$ref,$variant)),"\tNA\n");
+}
+elsif(length($variant) == 1 and length($ref) > 1){ # deletion
+record_snv("$target_parent\tg.$location\_",($location+length($ref)-1),
+"del",substr($ref, 1),"\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+join("\t",prop_info_key($chr,$location,$ref,$variant)),"\tNA\n");
+}
+else{ # indel
+record_snv("$target_parent\tg.$location\_",($location+length($ref)-1),
+"delins$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+join("\t",prop_info_key($chr,$location,$ref,$variant)),"\tNA\n");
+}
+next; # process next record, we've done all we can with a non-coding-gene SNP
+}
+	if($strand eq "-"){
+	    # set up utr offset for correct CDS coordinates
+	    for(my $i = $#feature_ranges; $i >= 0; $i--){
+		# exon is completely 5' of the start
+		if($feature_ranges[$i]->[0] > $feature_cds_max{$target_parent}){
+		    #print STDERR "Whole 5' UTR exon $i: ",$feature_ranges[$i]->[1]-$feature_ranges[$i]->[0]+1,"\n";
+		    $feature_offset -= $feature_ranges[$i]->[1]-$feature_ranges[$i]->[0]+1;
+		}
+		# exon with the cds start
+		elsif($feature_ranges[$i]->[1] >= $feature_cds_max{$target_parent} and
+		      $feature_ranges[$i]->[0] <= $feature_cds_max{$target_parent}){
+		    #print STDERR "Start codon in exon $i: ", $feature_cds_max{$target_parent} - $feature_ranges[$i]->[1],"\n";
+		    $feature_offset += $feature_cds_max{$target_parent} - $feature_ranges[$i]->[1];
+		    last;
+		}
+		else{
+		    die "The CDS for $target_parent (on negative strand) ends downstream ",
+		    "($feature_cds_max{$target_parent}) of the an exon",
+		    " (", $feature_ranges[$i]->[0], "), which is illogical.  Please revise the GFF file provided.\n";
+		}
+	    }
+	    for(my $i = $#feature_ranges; $i >= 0; $i--){
+		my $feature = $feature_ranges[$i];
+		# in the 3' UTR region of the gene
+		if($location < $feature_cds_min{$target_parent}){
+		    my $feature_exon = 0;
+		    $feature = $feature_ranges[$feature_exon];
+		    while($location > $feature->[1]+$flanking_bases and
+			  $feature_exon < $#feature_ranges){
+			$feature = $feature_ranges[++$feature_exon]; # find the 3' utr exon in which the variant is located
+		    }
+		    # utr exons passed entirely
+		    my $stop_exon = $feature_exon;
+		    while($feature_ranges[$stop_exon]->[1] < $feature_cds_min{$target_parent}){
+			$stop_exon++;
+		    }
+		    my $post_offset = $feature_cds_min{$target_parent}-$feature_ranges[$stop_exon]->[0]; # offset from the stop codon in the final coding exon
+		    while($stop_exon > $feature_exon){
+			$post_offset += $feature_ranges[$stop_exon]->[1]-$feature_ranges[$stop_exon]->[0]+1;
+			$stop_exon--;
+		    }
+		    my $pos = $feature->[1]-$location+1+$post_offset;
+my $junction_dist;
+		    if($location < $feature->[0]){ # after a UTR containing exon? set as .*DD+DD
+$junction_dist = ($feature->[0]-$location);
+			$pos = ($post_offset+$feature->[1]-$feature->[0]+1)."+".$junction_dist;
+		    }
+		    elsif($location > $feature->[1]){ # before a total UTR exon? set as .*DD-DD
+$junction_dist = -($location-$feature->[1]);
+			$pos = ($post_offset+1).$junction_dist;
+		    }
+else{ # in the UTR exon
+if($location - $feature->[0] < $feature->[1] - $location){ # location is closer to exon donor site
+$junction_dist = -($location - $feature->[0]+1); # +1 for HGVS syntax compatibility (there is no position 0, direct skip from -1 to +1)
+}
+else{
+$junction_dist = $feature->[1] - $location +1;
+}
+}
+		    if(length($ref) == 1 and length($variant) == 1){
+my $rc = join("",map({$rc{$_}} split(//,reverse($variant))));
+			# 3' UTR SNP
+			record_snv("$target_parent\tc.*$pos",
+			"$rc{$ref}>$rc\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			#"$ref>$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			join("\t",prop_info_key($chr,$location,$ref,$variant,$junction_dist)),"\tNA\n");
+		    }
+		    elsif(length($ref) == 1 and length($variant) > 1 and substr($variant, 0, 1) eq $ref){
+my $rc = join("",map({$rc{$_}} split(//,reverse(substr($variant,1)))));
+			# 3' UTR insertion
+			record_snv("$target_parent\tc.*",
+			hgvs_plus($pos,-1),"_*",$pos,
+			"ins$rc\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			#"ins",substr($variant,1),"\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			join("\t", prop_info_key($chr,$location,$ref,$variant,$junction_dist)),"\tNA\n");
+		    }
+		    elsif(length($ref) > 1 and length($variant) == 1 and substr($ref, 0, 1) eq $variant){
+my $rc = join("",map({$rc{$_}} split(//,reverse($ref))));
+my $delBases = substr($rc,0,length($rc)-1);
+			if(length($ref) == 2){
+			    # 3' UTR single base deletion
+			    record_snv("$target_parent\tc.*",hgvs_plus($pos,-1),
+			    "del$delBases\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			    join("\t", prop_info_key($chr,$location,$ref,$variant,$junction_dist)),"\tNA\n");
+			}
+			else{
+			    # longer 3' UTR deletion
+			    record_snv("$target_parent\tc.*",
+			    hgvs_plus($pos,-length($ref)+1),"_*",hgvs_plus($pos, -1),
+			    "del$delBases\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			    join("\t", prop_info_key($chr,$location,$ref,$variant,$junction_dist)),"\tNA\n");
+			}
+		    }
+		    else{
+			my $rc = join("",map({$rc{$_}} split(//,reverse($variant))));
+			if($rc eq $ref and length($variant) > 3){
+			    # 3' UTR inversion
+			    record_snv("$target_parent\tc.*",
+			    hgvs_plus($pos,-length($ref)+1),"_*",$pos,
+			    "inv\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			    join("\t", prop_info_key($chr,$location,$ref,$variant,$junction_dist)),"\tNA\n");
+			    last;
+			}
+			# complex substitution in 3' UTR
+			# Will break if stop codon is modified
+			record_snv("$target_parent\tc.*",
+			hgvs_plus($pos,-length($ref)+1),"_*", $pos,
+			"delins$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			join("\t", prop_info_key($chr,$location,$ref,$variant,$junction_dist)),"\tNA\n");
+		    }
+		    last;
+		}
+		# in the feature
+		elsif($location >= $feature->[0] and $location <= $feature->[1]){
+		    my $pos = $feature_offset+$feature->[1]-$location+1;
+		    if($location > $feature_cds_max{$target_parent}){ #since there is no position 0, the pos is in UTR, subtract one
+			$pos = hgvs_plus($pos, -1);
+		    }
+		    my $first_exon_base = $feature_offset+1;
+my $exon_edge_dist = $feature->[1]-$location+1; # equiv of HGVS +X or -X intron syntax, but for exons
+$exon_edge_dist = $feature->[0]-$location-1 if abs($feature->[0]-$location-1) < $exon_edge_dist; # pick closer of donor and acceptor sites
+#print STDERR "Got ", $feature->[1]-$location+1, "vs. ", $feature->[0]-$location-1, ": chose $exon_edge_dist\n";
+		    if(length($ref) == 1 and length($variant) == 1){
+			# SNP
+			record_snv("$target_parent\tc.",
+			$pos,
+			"$rc{$ref}>$rc{$variant}\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			join("\t", prop_info_key($chr,$location,$ref,$variant,$exon_edge_dist)),"\t",
+($pos < 1 ? "NA" : hgvs_protein($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+#($pos < 1 ? "NA" : hgvs_protein_key($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+		    }
+		    elsif(length($ref) == 1 and length($variant) > 1 and substr($variant, 0, 1) eq $ref){
+			my $rc = join("",map({$rc{$_}} split(//,reverse($variant))));
+my $insBases = substr($rc,1);
+			# insertion
+			record_snv("$target_parent\tc.",
+			hgvs_plus_exon($pos, -1, $first_exon_base),"_",$pos,"ins$insBases",
+			"\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			join("\t", prop_info_key($chr,$location,$ref,$variant,$exon_edge_dist)),"\t",
+($pos < 1 ? "NA" : hgvs_protein($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+#($pos < 1 ? "NA" : hgvs_protein_key($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+		    }
+		    elsif(length($ref) > 1 and length($variant) == 1 and substr($ref, 0, 1) eq $variant){
+			my $rc = join("",map({$rc{$_}} split(//,reverse($ref))));
+my $delBases = substr($rc,0,length($rc)-1);
+			# single nucleotide deletion
+			if(length($ref) == 2){
+			    record_snv("$target_parent\tc.",
+			    hgvs_plus_exon($pos, -1, $first_exon_base),
+			    "del$delBases\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			    join("\t", prop_info_key($chr,$location,$ref,$variant,$exon_edge_dist)),"\t",
+($pos-1 < 1 ? "NA" : $pos-1 < $first_exon_base ? "p.?" : hgvs_protein($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+#($pos-1 < 1 ? "NA" : $pos-1 < $first_exon_base ? "p.?" : hgvs_protein_key($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+			}
+			# longer deletion
+			else{
+$exon_edge_dist = $feature->[1]-$location-length($ref)+1 if abs($feature->[1]-$location-length($ref)+1) < $exon_edge_dist;
+			    record_snv("$target_parent\tc.",
+			    hgvs_plus_exon($pos, -length($ref)+1, $first_exon_base),"_",
+			    hgvs_plus_exon($pos, -1, $first_exon_base),
+			    "del$delBases\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			    join("\t", prop_info_key($chr,$location,$ref,$variant,$exon_edge_dist)),"\t",
+($pos-1 < 1 ? "NA" : $pos-length($ref)+1 < $first_exon_base ? "p.?" : hgvs_protein($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+#($pos-1 < 1 ? "NA" : $pos-length($ref)+1 < $first_exon_base ? "p.?" : hgvs_protein_key($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+			}
+		    }
+		    else{
+			# complex substitution
+$exon_edge_dist = $feature->[1]-$location-length($ref)+1 if abs($feature->[1]-$location-length($ref)+1) < $exon_edge_dist;
+			my $rc = join("",map({$rc{$_}} split(//,reverse($variant))));
+if($rc eq $variant and length($variant) > 3){
+# inversion
+record_snv("$target_parent\tc.",
+hgvs_plus_exon($pos,-length($ref)+1,$first_exon_base),"_",
+$pos,
+"inv",
+"\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+join("\t", prop_info_key($chr,$location,$ref,$variant,$exon_edge_dist)),"\t",
+($pos < 1 ? "NA" : $pos-length($ref)+1 < $first_exon_base ? "p.?" : hgvs_protein($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+last;
+}
+			record_snv("$target_parent\tc.",
+			hgvs_plus_exon($pos,-length($ref)+1,$first_exon_base),"_",
+			$pos,
+			"delins$rc",
+			"\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			join("\t", prop_info_key($chr,$location,$ref,$variant,$exon_edge_dist)),"\t",
+($pos < 1 ? "NA" : $pos-length($ref)+1 < $first_exon_base ? "p.?" : hgvs_protein($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+#($pos < 1 ? "NA" : $pos-length($ref)+1 < $first_exon_base ? "p.?" : hgvs_protein_key($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+		    }
+		    last;
+		}
+		# 5' of feature (on negative strand)
+		elsif($location > $feature->[1]){
+		    if(length($ref) == 1 and length($variant) == 1){
+			# intronic SNP
+			if($i == $#feature_ranges or $feature->[1]-$location >= -1*$flanking_bases){
+			    # Closer to acceptor site
+			    record_snv("$target_parent\tc.",$feature_offset+1,
+			    ($feature->[1]-$location),
+			    "$rc{$ref}>$rc{$variant}\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			    #"$ref>$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			    join("\t", prop_info_key($chr,$location,$ref,$variant, $feature->[1]-$location)),"\tNA\n");
+			}
+			else{
+			    # Closer to donor site
+			    record_snv("$target_parent\tc.",$feature_offset,"+",
+			    ($feature_ranges[$i+1]->[0]-$location),
+			    "$rc{$ref}>$rc{$variant}\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			    #"$ref>$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			    join("\t", prop_info_key($chr,$location,$ref,$variant, $feature_ranges[$i+1]->[0]-$location)),"\tNA\n");
+			}
+		    }
+		    elsif(length($ref) == 1 and length($variant) > 1 and substr($variant, 0, 1) eq $ref){
+my $rc = join("",map({$rc{$_}} split(//,reverse(substr($variant,1)))));
+			if($i == $#feature_ranges or $feature->[1]-$location >= -1*$flanking_bases){
+			    # intronic insertion near acceptor
+			    my $pos = ($feature_offset+1).($feature->[1]-$location);
+			    record_snv("$target_parent\tc.",
+			    hgvs_plus($pos,-1),"_",$pos,
+			    "ins",
+			    $rc,"\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			    #substr($variant, 1),"\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			    join("\t", prop_info_key($chr,$location,$ref,$variant,$feature->[1]-$location-1)),"\tNA\n");
+			}
+			else{
+			    # intronic insertion near donor
+			    my $pos = $feature_offset."+".($feature_ranges[$i+1]->[0]-$location);
+			    record_snv("$target_parent\tc.",
+			    hgvs_plus($pos,-1),"_",$pos,
+			    "ins",
+			    $rc,"\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			    #substr($variant, 1),"\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			    join("\t", prop_info_key($chr,$location,$ref,$variant,$feature_ranges[$i+1]->[0]-$location+1)),"\tNA\n");
+			}
+		    }
+		    elsif(length($ref) > 1 and length($variant) == 1 and substr($ref, 0, 1) eq $variant){
+			# intronic deletion
+			# single nucleotide deletion
+			my $rc = reverse($ref);
+			$rc=~tr/ACGT/TGCA/;
+my $delBases = substr($rc, 0, length($rc)-1);
+			if(length($ref) == 2){
+			    # single intronic deletion near acceptor
+			    if($i == $#feature_ranges or $feature->[1]-$location >= -1*$flanking_bases){
+my $off = $feature->[1]-$location-1;
+				record_snv("$target_parent\tc.",
+				($feature_offset+1),$off,
+				"del$delBases\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+				join("\t", prop_info_key($chr,$location,$ref,$variant,$off)),"\t",($off >= -2 ? "p.?" : "NA"),"\n");
+			    }
+			    # single intronic deletion near donor
+			    else{
+				my $pos = $feature_offset;
+my $off = $feature_ranges[$i+1]->[0]-$location+1;
+				record_snv("$target_parent\tc.",
+				hgvs_plus_exon($pos, $off, $pos),
+				"del$delBases\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+				join("\t", prop_info_key($chr,$location,$ref,$variant,$off)),"\t",($off <= 2 ? "p.?" : "NA"),"\n");
+			    }
+			}
+			# longer deletion
+			else{
+			    if($i == $#feature_ranges or $feature->[1]-$location >= -1*$flanking_bases){
+				# long intronic deletion near acceptor
+my $off = $feature->[1]-$location-1;
+				my $pos = ($feature_offset+1).$off;
+				record_snv("$target_parent\tc.",
+				hgvs_plus($pos,-length($ref)+2),"_",$pos,
+				"del$delBases\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+				join("\t", prop_info_key($chr,$location,$ref,$variant,$off)),"\t",($off >= -2 ? "p.?" : "NA"),"\n");
+			    }
+			    else{
+				# long intronic deletion near donor
+my $off = $feature_ranges[$i+1]->[0]-$location+1;
+				my $pos = ($feature_offset)."+".$off;
+				record_snv("$target_parent\tc.",
+				$pos,"_",hgvs_plus($pos,-length($ref)-1),
+				"del$delBases\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+				join("\t", prop_info_key($chr,$location,$ref,$variant,$off)),"\t",($off-length($ref)+1 <= 2 ? "p.?" : "NA"),"\n");
+			    }
+			    last;
+			}
+		    }
+		    else{
+			my $rc = reverse($ref);
+			$rc=~tr/ACGT/TGCA/;
+			if($rc eq $variant and length($variant) > 3){
+			    # intronic inversion near acceptor site
+			    if($i == $#feature_ranges or $feature->[1]-$location >= -1*$flanking_bases){
+				my $pos = ($feature_offset+1).($feature->[1]-$location);
+				record_snv("$target_parent\tc.",
+				hgvs_plus($pos,-length($ref)+1),"_",$pos,
+				"inv\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+				join("\t", prop_info_key($chr,$location,$ref,$variant,$feature->[1]-$location)),"\tNA\n");
+			    }
+			    else{
+				my $pos = ($feature_offset)."+".($feature_ranges[$i+1]->[0]-$location);
+				record_snv("$target_parent\tc.",
+				$pos,"_",hgvs_plus($pos, length($ref)-1),
+				"inv\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+				join("\t", prop_info_key($chr,$location,$ref,$variant,$feature_ranges[$i+1]->[0]-$location)),"\tNA\n");
+			    }
+			    last;
+			}
+$rc = reverse($variant);
+$rc=~tr/ACGT/TGCA/;
+			# Intronic complex substitution
+			if($i == $#feature_ranges or $feature->[1]-$location >= -1*$flanking_bases){
+			    # complex intronic substitution near acceptor site
+			    my $pos = ($feature_offset+1).($feature->[1]-$location);
+			    record_snv("$target_parent\tc.",
+			    hgvs_plus($pos, -length($ref)+1),"_",$pos,
+			    "delins$rc\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			    #"delins$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			    join("\t", prop_info_key($chr,$location,$ref,$variant,$feature->[1]-$location)),"\tNA\n");
+			}
+			else{
+			    # complex intronic substitution near donor site
+			    my $pos = $feature_offset."+".($feature_ranges[$i+1]->[0]-$location);
+			    record_snv("$target_parent\tc.",
+			    $pos,"_",hgvs_plus($pos, length($ref)-1),
+			    "delins$rc\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			    #"delins$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			    join("\t", prop_info_key($chr,$location,$ref,$variant,$feature_ranges[$i+1]->[0]-$location)),"\tNA\n");
+			}
+		    }
+		    last;
+		}
+		else{
+		    #print STDERR "Offset going from ", $feature_offset, " to ", $feature_offset+$feature->[1]-$feature->[0]+1,"\n";
+		    $feature_offset += $feature->[1]-$feature->[0]+1;
+		    #print STDERR "Set feature offset to  $feature_offset by adding ",$feature->[1],"-", $feature->[0],"+1\n";
+		}
+	    }
+	}
+	else{
+	    # forward strand
+	    # set up utr offset for correct CDS coordinates
+	    for(my $i = 0; $i <= $#feature_ranges; $i++){
+		# All 5' utr exon
+		if($feature_ranges[$i]->[1] < $feature_cds_min{$target_parent}){
+		    $feature_offset -= $feature_ranges[$i]->[1]-$feature_ranges[$i]->[0]+1;
+		}
+		# exon with the cds start
+		elsif($feature_ranges[$i]->[1] >= $feature_cds_min{$target_parent} and
+		      $feature_ranges[$i]->[0] <= $feature_cds_min{$target_parent}){
+		    $feature_offset -= $feature_cds_min{$target_parent} - $feature_ranges[$i]->[0];
+		    last;
+		}
+		else{
+		    die "The CDS for $target_parent starts upstream ($feature_cds_max{$target_parent}) of the first exon",
+		    " (", $feature_ranges[$i]->[0], "), which is illogical.  Please revise the GFF file provided.\n";
+		}
+	    }
+	    for(my $i = 0; $i <= $#feature_ranges; $i++){
+		my $feature = $feature_ranges[$i];
+		# 3' of last coding position
+		if($location > $feature_cds_max{$target_parent}){
+		    # find the exon with the stop codon
+		    while($feature->[1] < $feature_cds_max{$target_parent}){
+			$feature = $feature_ranges[++$i];
+		    }
+		    my $post_offset = $feature->[0]-$feature_cds_max{$target_parent};
+		    while($location > $feature->[1]+$flanking_bases and
+			  $i < $#feature_ranges){
+			$post_offset += $feature->[1]-$feature->[0]+1;
+			$feature = $feature_ranges[++$i]; # find the 3' utr exon in which the variant is located
+		    }
+		    my $pos = $location-$feature->[0]+$post_offset;
+		    #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;
+my $off;
+		    if($location > $feature->[1]){ # after a UTR containing exon? set as .*DD+DD
+$off = $location-$feature->[1];
+			$pos = ($post_offset+$feature->[1]-$feature->[0]+1)."+".$off;
+		    }
+		    elsif($location < $feature->[0]){ # before a total UTR exon? set as .*DD-DD
+$off = -($feature->[0]-$location);
+			$pos = ($post_offset+1).$off;
+		    }
+else{
+if($location-$feature->[0] < $feature->[1]-$location){
+$off = $location-$feature->[0]+1; # +1 since HGVS skips right from -1 to +1 at exon boundary
+}
+else{
+$off = $location-$feature->[1]-1; # will be negative
+}
+}
+		    if(length($ref) == 1 and length($variant) == 1){
+			# 3' UTR SNP
+			record_snv("$target_parent\tc.*$pos",
+			"$ref>$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			join("\t", prop_info_key($chr,$location,$ref,$variant, $off)),"\tNA\n");
+		    }
+		    elsif(length($ref) == 1 and length($variant) > 1 and substr($variant, 0, 1) eq $ref){
+			# 3' UTR insertion
+			record_snv("$target_parent\tc.*",
+			hgvs_plus($pos,1),"_*",hgvs_plus($pos,2),
+			"ins",substr($variant,1),"\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			join("\t", prop_info_key($chr,$location,$ref,$variant, $off)),"\tNA\n");
+		    }
+		    elsif(length($ref) > 1 and length($variant) == 1 and substr($ref, 0, 1) eq $variant){
+my $delBases = substr($ref, 1);
+			if(length($ref) == 2){
+			    # 3' UTR single base deletion
+			    record_snv("$target_parent\tc.*",hgvs_plus($pos,1),
+			    "del$delBases\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			    join("\t", prop_info_key($chr,$location,$ref,$variant,$off)),"\tNA\n");
+			}
+			else{
+			    # longer 3' UTR deletion
+			    record_snv("$target_parent\tc.*",
+			    hgvs_plus($pos,1),"_*",hgvs_plus($pos,length($ref)-1),
+			    "del$delBases\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			    join("\t", prop_info_key($chr,$location,$ref,$variant, $off)),"\tNA\n");
+			}
+		    }
+		    else{
+			my $rc = reverse($ref);
+			$rc=~tr/ACGT/TGCA/;
+			if($rc eq $variant and length($variant) > 3){
+			    # 3' UTR inversion
+			    record_snv("$target_parent\tc.*$pos",
+			    "_*",hgvs_plus($pos,length($ref)-1),
+			    "inv\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			    join("\t", prop_info_key($chr,$location,$ref,$variant,$off)),"\tNA\n");
+			    last;
+			}
+			# complex substitution in 3' UTR
+			record_snv("$target_parent\tc.*$pos",
+			"_*",hgvs_plus($pos,length($ref)-1),
+			"delins$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			join("\t", prop_info_key($chr,$location,$ref,$variant,$off)),"\tNA\n");
+		    }
+		    last;
+		}
+		# in the exon
+		elsif($location >= $feature->[0] and $location <= $feature->[1]){
+		    my $pos = $feature_offset+$location-$feature->[0]+1;
+		    my $last_exon_base = $feature_offset+$feature->[1]-$feature->[0]+1;
+my $exon_edge_dist = $location-$feature->[0]+1; # equiv of HGVS +X or -X intron syntax, but for exons
+$exon_edge_dist = $location-$feature->[1]-1 if abs($location-$feature->[1]-1) < $exon_edge_dist; # pick closer of donor and acceptor sites
+#print STDERR "Got ", $location-$feature->[0]+1, "vs. ", $location-$feature->[1]-1, ": chose $exon_edge_dist\n";
+		    if($location < $feature_cds_min{$target_parent}){ #since there is no position 0, the pos is in UTR, subtract one
+			$pos = hgvs_plus($pos, -1);
+		    }
+		    if(length($ref) == 1 and length($variant) == 1){
+			# SNP
+			record_snv("$target_parent\tc.$pos",
+			"$ref>$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			join("\t", prop_info_key($chr,$location,$ref,$variant,$exon_edge_dist)),"\t",
+($pos < 1 ? "NA" : hgvs_protein($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+#($pos < 1 ? "NA" : hgvs_protein_key($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+		    }
+		    elsif(length($ref) == 1 and length($variant) > 1 and substr($variant, 0, 1) eq $ref){
+			# insertion
+			record_snv("$target_parent\tc.$pos",
+			"_",hgvs_plus_exon($pos,1,$last_exon_base),"ins",
+			substr($variant, 1),"\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			join("\t", prop_info_key($chr,$location,$ref,$variant,$exon_edge_dist)),"\t",
+($pos < 1 ? "NA" : hgvs_protein($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+#($pos < 1 ? "NA" : hgvs_protein_key($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+		    }
+		    elsif(length($ref) > 1 and length($variant) == 1 and substr($ref, 0, 1) eq $variant){
+my $delBases = substr($ref, 1);
+			# single nucleotide deletion
+			if(length($delBases) == 1){
+			    record_snv("$target_parent\tc.",
+			    hgvs_plus_exon($pos,1,$last_exon_base),
+			    "del$delBases\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			    join("\t", prop_info_key($chr,$location,$ref,$variant,$exon_edge_dist)),"\t",
+($pos < 1 or $pos > $last_exon_base ? "NA" : hgvs_protein($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+#($pos < 1 or $pos > $last_exon_base ? "NA" : hgvs_protein_key($chr,$location,$ref,$variant,$pos,$strand,$trans_table)),"\n");
+			}
+			# longer deletion
+			else{
+$exon_edge_dist = $feature->[1]-$location-length($ref)-1 if abs($feature->[1]-$location-length($ref)-1) < $exon_edge_dist;
+			    record_snv("$target_parent\tc.",
+			    hgvs_plus_exon($pos,1,$last_exon_base),"_",
+			    hgvs_plus_exon($pos,length($ref)-1,$last_exon_base),
+			    "del$delBases\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			    join("\t", prop_info_key($chr,$location,$ref,$variant,$exon_edge_dist)),"\t",
+($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");
+#($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");
+			}
+		    }
+		    else{
+$exon_edge_dist = $feature->[1]-$location-length($ref)-1 if abs($feature->[1]-$location-length($ref)-1) < $exon_edge_dist;
+			my $rc = reverse($ref);
+			$rc=~tr/ACGT/TGCA/;
+			if($rc eq $variant and length($variant) > 3){
+			    # inversion
+			    record_snv("$target_parent\tc.$pos",
+			    "_",hgvs_plus_exon($pos,length($ref)-1, $last_exon_base),
+			    "inv\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			    join("\t", prop_info_key($chr,$location,$ref,$variant,$exon_edge_dist)),"\t",
+($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");
+#($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");
+			    last;
+			}
+			# complex substitution
+			record_snv("$target_parent\tc.$pos",
+			"_",hgvs_plus_exon($pos, length($ref)-1, $last_exon_base),
+			"delins$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			join("\t", prop_info_key($chr,$location,$ref,$variant,$exon_edge_dist)),"\t",
+($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");
+#($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");
+		    }
+		    last;
+		}
+		# 5' of feature
+		elsif($location < $feature->[0]){
+		    if(length($ref) == 1 and length($variant) == 1){
+			# intronic SNP
+			if($i != 0 and $location-$feature->[0] < -1*$flanking_bases){
+			    # Closer to donor site
+			    record_snv("$target_parent\tc.",$feature_offset,"+",
+			    ($location-$feature_ranges[$i-1]->[1]),
+			    "$ref>$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			    join("\t", prop_info_key($chr,$location,$ref,$variant,$location-$feature_ranges[$i-1]->[1])),"\tNA\n");
+			}
+			else{
+			    # Closer to acceptor site
+			    record_snv("$target_parent\tc.",$feature_offset+1,
+			    ($location-$feature->[0]),
+			    "$ref>$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			    join("\t", prop_info_key($chr,$location,$ref,$variant,$location-$feature->[0])),"\tNA\n");
+			}
+		    }
+		    elsif(length($ref) == 1 and length($variant) > 1 and substr($variant, 0, 1) eq $ref){
+			if($i != 0 and $location-$feature->[0] < -1*$flanking_bases){
+			    # intronic insertion near donor
+			    my $pos = $feature_offset."+".($location-$feature_ranges[$i-1]->[1]);
+			    record_snv("$target_parent\tc.",
+			    $pos,"_",hgvs_plus($pos,1),
+			    "ins",
+			    substr($variant, 1),"\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			    join("\t", prop_info_key($chr,$location,$ref,$variant,$location-$feature_ranges[$i-1]->[1])),"\tNA\n");
+			}
+			else{
+			    # intronic insertion near acceptor
+			    my $pos = ($feature_offset+1).($location-$feature->[0]);
+			    record_snv("$target_parent\tc.",
+			    $pos,"_",hgvs_plus($pos,1),
+			    "ins",
+			    substr($variant, 1),"\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			    join("\t", prop_info_key($chr,$location,$ref,$variant,$location-$feature->[0])),"\tNA\n");
+			}
+		    }
+		    elsif(length($ref) > 1 and length($variant) == 1 and substr($ref, 0, 1) eq $variant){
+			# intronic deletion
+			# single nucleotide deletion
+my $delBases = substr($ref, 1);
+			if(length($ref) == 2){
+			    # single intronic deletion near donor
+			    if($i != 0 and $location-$feature->[0] < -1*$flanking_bases){
+my $off = $location-$feature_ranges[$i-1]->[1]+1;
+				record_snv("$target_parent\tc.",
+				$feature_offset,"+",$off,
+				"del$delBases\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+				join("\t", prop_info_key($chr,$location,$ref,$variant,$off)),"\t",($off <= 2 ? "p.?" : "NA"),"\n");
+			    }
+			    # single intronic deletion near acceptor
+			    else{
+				my $pos = ($feature_offset+1);
+my $off = $location-$feature->[0];
+				record_snv("$target_parent\tc.",
+				hgvs_plus_exon($pos, $off, $pos),
+				"del$delBases\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+				join("\t", prop_info_key($chr,$location,$ref,$variant,$off)),"\t",($off >= -2 ? "p.?" : "NA"),"\n");
+			    }
+			}
+			# longer deletion
+			else{
+			    if($i != 0 and $location-$feature->[0] < -1*$flanking_bases){
+				# long intronic deletion near donor
+my $off = $location-$feature_ranges[$i-1]->[1]+1;
+				my $pos = $feature_offset."+".$off;
+				record_snv("$target_parent\tc.",
+				$pos,"_",hgvs_plus($pos,length($ref)-2),
+				"del$delBases\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+				join("\t", prop_info_key($chr,$location,$ref,$variant,$off)),"\t",($off <= 2 ? "p.?" : "NA"),"\n");
+			    }
+			    else{
+				# long intronic deletion near acceptor
+my $off = $location-$feature->[0]+1;
+				my $pos = ($feature_offset+1).$off;
+				record_snv("$target_parent\tc.",
+				$pos,"_",hgvs_plus($pos,length($ref)-2),
+				"del$delBases\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+				join("\t", prop_info_key($chr,$location,$ref,$variant,$off)),"\t",($off+length($ref)-2 >= -2 ? "p.?" : "NA"),"\n");
+			    }
+			}
+		    }
+		    else{
+			my $rc = reverse($ref);
+			$rc=~tr/ACGT/TGCA/;
+			if($rc eq $variant and length($variant) > 3){
+			    # intronic inversion near donor site
+			    if($i != 0 and $location-$feature->[0] < -1*$flanking_bases){
+				my $pos = $feature_offset."+".($location-$feature_ranges[$i-1]->[1]);
+				record_snv("$target_parent\tc.",
+				$pos,"_",hgvs_plus($pos,length($ref)-1),
+				"inv\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+				join("\t", prop_info_key($chr,$location,$ref,$variant,$location-$feature_ranges[$i-1]->[1])),"\tNA\n");
+			    }
+			    else{
+				my $pos = ($feature_offset+1).($location-$feature->[0]);
+				record_snv("$target_parent\tc.",
+				$pos,"_",hgvs_plus($pos, length($ref)-1),
+				"inv\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+				join("\t", prop_info_key($chr,$location,$ref,$variant,$location-$feature->[0])),"\tNA\n");
+			    }
+			    last;
+			}
+			# Intronic complex substitution
+			# Note: sub maybe have  comma in it to denote two possible variants
+			if($i != 0 and $location-$feature->[0] < -1*$flanking_bases){
+			    # complex intronic substitution near donor site
+			    my $pos = $feature_offset."+".($location-$feature_ranges[$i-1]->[1]);
+			    record_snv("$target_parent\tc.",
+			    $pos,"_",hgvs_plus($pos, length($ref)-1),
+			    "delins$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			    join("\t", prop_info_key($chr,$location,$ref,$variant,$location-$feature_ranges[$i-1]->[1])),"\tNA\n");
+			}
+			else{
+			    # complex intronic substitution near acceptor site
+			    my $pos = ($feature_offset+1).($location-$feature->[0]);
+			    record_snv("$target_parent\tc.",
+			    $pos,"_",hgvs_plus($pos, length($ref)-1),
+			    "delins$variant\t$strand\t$chr\t$location\t$zygosity\t$quality\t$variant_depth\t$read_depth\t",
+			    join("\t", prop_info_key($chr,$location,$ref,$variant,$location-$feature->[0])),"\tNA\n");
+			}
+		    }
+		    last;
+		}
+		else{
+		    # feature is past this exon
+		    $feature_offset += $feature->[1]-$feature->[0]+1;
+		}
+	    }
+	}
+}  # for each variant on the line
+} # for each  gene overlapping the site the VCF describes
+} # for each VCF line
+print STDERR "\n" unless $quiet;
+close(VCFIN);
+# Before we can start printing the variants, we need to look at the phase information and calculate the real haplotype HGVS changes
+#if(keys %chr2phase){
+# 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)
+# We need to create new disjoint sets of phased blocks from the (consistent) union of these data.
+#  my $chr2phase2variants = combine_phase_data(\%chr2phase);
+# TODO: Calculate protein HGVS syntax for each variant, now that all phase data has been incorporated
+#for my $chr (keys %$chr2phase2variants){
+#  for my $phase (keys %{$chr2phase2variants{$chr}){
+#    # apply all phased changes to the reference chromosomal seq
+#    my $phased_seq = $seq{$chr}; #reference
+#    # sort the variants from 3' to 5' so that edits after indels don't need adjustment in their ref coordinate
+#    my @sorted_variants = sort {my($a_pos) = $a =~ /:(\d+):/; my($b_pos) = $b =~ /:(\d+):/; $b_pos <=> $a_pos} @{$chr2phase2variants{$chr}->{$phase}};
+#    for my $variant (@sorted_variants){
+#    }
+#  }
+#}
+#}
+# retrieve the MAF info en masse for each chromosome, as this is much more efficient
+my @out_lines;
+for my $snv (@snvs){
+chomp $snv;
+my @fields = split /\t/, $snv;
+# For CNVs, all the fields are already filled out
+if(@fields > 13){
+push @out_lines, join("\t", $feature_type{$fields[0]}, ($fields[0] =~ /\S/ ? $feature_length{$fields[0]} : "NA"), @fields);
+next;
+}
+my $variant_key = $fields[9];
+$fields[9] = join("\t", prop_info($dbsnp,$internal_snp,$variant_key));
+my $from = $fields[4];
+my $chr_pos_key = $fields[3].":".$from;
+my $to = $fields[4]; # true for SNPs and insertions
+my @variant_key = split /:/, $variant_key;
+# For deletions and complex variants, calculate the affected reference genome range and set the 'to'
+if(length($variant_key[2]) > 1){
+$to += length($variant_key[2])-1;
+}
+splice(@fields, 5, 0, $to);
+# Otherwise expand the key into the relevant MAF values
+$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
+# the offset from the nearest exon border if coding
+push @fields, ($#variant_key > 3 ? $variant_key[4] : "");
+# add gene name(s)
+push @fields, range2genes($fields[3], $from, $to+1);
+# add caveats
+my $c = $fields[3];
+if(not exists $chr2mappability{$c}){
+if($c =~ s/^chr//){
+# nothing more
+}
+else{
+$c = "chr$c";
+}
+}
+my $mappability_caveats = exists $chr2mappability{$c} ? $chr2mappability{$c}->fetch($fields[4], $fields[4]+1) : [];
+if(ref $mappability_caveats eq "ARRAY" and @$mappability_caveats){
+my %h;
+if(exists $chr2caveats{$chr_pos_key}){
+if($chr2caveats{$chr_pos_key} !~ /non-unique/){
+$chr2caveats{$chr_pos_key} = join("; ", grep {not $h{$_}++} @$mappability_caveats)."; ".$chr2caveats{$chr_pos_key};
+}
+}
+else{
+$chr2caveats{$chr_pos_key} = join("; ", grep {not $h{$_}++} @$mappability_caveats);
+}
+}
+push @fields, (exists $chr2caveats{$chr_pos_key} ? $chr2caveats{$chr_pos_key} : "");
+# add phase
+push @fields, find_phase($variant_key);
+push @out_lines, join("\t", $feature_type{$fields[0]}, $feature_length{$fields[0]}, @fields);
+}
+# Now tabulate the rare variant numbers
+my %gene2rares;
+my %gene2aa_rares;
+for my $line (@out_lines){
+my @F = split /\t/, $line, -1;
+if($F[15] eq "NA" or $F[15] < $rare_variant_prop and (!$internal_snp or $F[17] < $rare_variant_prop)){
+my $gene_list = $internal_snp ? $F[20] : $F[19];
+next unless defined $gene_list;
+for my $g (split /; /, $gene_list){
+$gene2rares{$g}++;
+# Check the cDNA HGVS syntax for relevance
+if($F[3] =~ /c.\d+/ or # coding
+$F[3] =~ /c.\d+.*-[12]/ or # splicing acceptor
+$F[3] =~ /c.\d+\+[12345]/){ # splicing donor
+$gene2aa_rares{$g}++;
+}
+}
+}
+}
+# Print the results
+for my $line (@out_lines){
+my @F = split /\t/, $line, -1;
+my $gene_list = $internal_snp ? $F[20] : $F[19];
+if(not defined $gene_list){
+print OUT join("\t", @F, "", ""), "\n"; next;
+}
+my $max_gene_rare = 0;
+my $max_gene_aa_rare = 0;
+for my $g (split /; /, $gene_list){
+next unless exists $gene2rares{$g};
+if($gene2rares{$g} > $max_gene_rare){
+$max_gene_rare = $gene2rares{$g};
+}
+next unless exists $gene2aa_rares{$g};
+if($gene2aa_rares{$g} > $max_gene_aa_rare){
+$max_gene_aa_rare = $gene2aa_rares{$g};
+}
+}
+print OUT join("\t", @F, $max_gene_rare, $max_gene_aa_rare), "\n";
+}
+close(OUT);
+sub range2genes{
+my ($c, $from, $to) = @_;
+if(not exists $gene_ids{$c}){
+if($c =~ s/^chr//){
+# nothing more
+}
+else{
+$c = "chr$c";
+}
+}
+if(exists $gene_ids{$c}){
+my %have;
+return join("; ", grep {not $have{$_}++} @{$gene_ids{$c}->fetch($from, $to+1)});
+}
+else{
+return "";
+}
+}
+sub combine_phase_data{
+my ($phases) = @_; # map from variant to its phase data
+# Create a reverse mapping from phase regions to their variants
+my %chr2phase_region2variants;
+my @variants = keys %$phases;
+for my $variant (@variants){
+my ($chr) = $variant =~ /^\S+?-(\S+):/;
+$chr2phase_region2variants{$chr} = {} if not exists $chr2phase_region2variants{$chr};
+for my $phase_region (split /,/, $phases->{$variant}){
+$chr2phase_region2variants{$chr}->{$phase_region} = [] if not exists $chr2phase_region2variants{$chr}->{$phase_region};
+push @{$chr2phase_region2variants{$phase_region}}, $variant;
+}
+}
+# Now for each phased block known so far, see if any variant in it is also part of another block
+# If so, do a union since phasing is both transitive and commutative.
+# The quickest way to do this is to check for overlapping intervals, then check for common members amongst those that do overlap
+for my $chr (keys %chr2phase_region2variants){
+my @ordered_phase_regions = sort {my($a_pos) = $a =~ /:(\d+)/; my($b_pos) = $b =~ /:(\d+)/; $a_pos <=> $b_pos} keys %{$chr2phase_region2variants{$chr}};
+my $sets = new DisjointSets(scalar(@ordered_phase_regions));
+for (my $i = 0; $i < $#ordered_phase_regions; $i++){
+my ($start, $stop, $variant) = $ordered_phase_regions[$i];
+for (my $j = $i+1; $j <= $#ordered_phase_regions; $j++){
+my ($start2, $stop2, $variant2) = $ordered_phase_regions[$j];
+if($start2 > $stop){ # won't overlap any regions after this in the sorted list
+last;
+}
+# If we get here, it is implicit that $stop >= $start2 and $start < $stop2, i.e. there is overlap
+# Now check if there is a shared variant (otherwise we might erroneously join blocks from different physical chromosomal arms)
+my $have_shared_variant = 0;
+my $overlapping_phase_region = $ordered_phase_regions[$j];
+for my $variant (@{$chr2phase_region2variants{$chr}->{$ordered_phase_regions[$i]}}){
+if($phases->{$variant} =~ /\b$overlapping_phase_region\b/){
+$have_shared_variant = 1; last;
+}
+}
+# sanity check that there aren't conflicting variants in the new block (i.e. two different variants in the same position)
+my %pos2base;
+my $have_conflicting_variant = 0;
+for my $variant (@{$chr2phase_region2variants{$chr}->{$ordered_phase_regions[$i]}}, @{$chr2phase_region2variants{$chr}->{$ordered_phase_regions[$j]}}){
+my ($pos, $base) = $variant =~ /(\d+):(.+?)$/;
+if(exists $pos2base{$pos} and $pos2base{$pos} ne $base){
+# conflict, note with a caveat
+if(exists $chr2caveats{"$chr:$pos"}){
+$chr2caveats{"$chr:$pos"} .= "; inconsistent haplotype phasing" unless $chr2caveats{"$chr:$pos"} =~ /inconsistent haplotype phasing/;
+}
+else{
+$chr2caveats{"$chr:$pos"} = "inconsistent haplotype phasing";
+}
+$have_conflicting_variant ||= 1;
+}
+elsif(not exists $pos2base{$pos}){
+$pos2base{$pos} = $base;
+}
+}
+$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
+}
+}
+my $phase_sets = $sets->sets; #disjoint haplotype sets
+my %region_counts;
+for my $phase_set (@$phase_sets){
+next if scalar(@$phase_set) == 1; # No change to existing phase region (is disjoint from all others)
+# Generate a new haploblock to replace the old ones that are being merged
+my $merged_start = 10000000000;
+my $merged_end = 0;
+for my $ordered_phase_region_index (@$phase_set){
+my ($start, $end) = $ordered_phase_regions[$ordered_phase_region_index-1] =~ /(\d+)-(\d+)$/;
+$merged_start = $start if $start < $merged_start;
+$merged_end = $end if $end > $merged_end;
+}
+# 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
+my $region_count = $region_counts{"$merged_start-$merged_end"}++;
+my $merged_haploblock_name = "Y$region_count-$chr:$merged_start-$merged_end";
+# Assign this new name to overwrite the premerge values for each variant contained within
+for my $ordered_phase_region_index (@$phase_set){
+for my $variant (@{$chr2phase_region2variants{$chr}->{$ordered_phase_regions[$ordered_phase_region_index-1]}}){ # incl. one-based set correction in 0-based array index
+print STDERR "Merging $variant from ", $phases->{$variant}, " into new block $merged_haploblock_name\n";
+$phases->{$variant} = $merged_haploblock_name;
+}
+}
+}
+# 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...
+}
+}
+# Sees if the positions of the variant are in the range of a phased haplotype, returning which allele it belongs to
+sub find_phase{
+my ($chr,$pos,$ref,$variant) = split /:/, $_[0];
+return "" if length($ref) != length($variant); # Can only deal with SNPs (and broken down MNPs) for now
+for(my $i = 0; $i < length($ref); $i++){
+my $key = "$chr:".($pos+$i).":".substr($variant, $i, 1);
+#print STDERR "Checking phase for $key\n" if $pos == 12907379;
+if(exists $chr2phase{$key}){
+#print STDERR "returning phase data $chr2phase{$key}\n";
+return $chr2phase{$key};
+}
+elsif(exists $chr2phase{"chr".$key}){
+#print STDERR "returning phase data ", $chr2phase{"chr".$key}, "\n";
+return $chr2phase{"chr".$key};
+}
+}
+return "";
+}
+sub find_earliest_index{
+# 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
+my ($query, $array) = @_;
+return 0 if $query < $array->[0]->[0];
+my ($left, $right, $prevCenter) = (0, $#$array, -1);
+while(1){
+my $center = int (($left + $right)/2);
+my $cmp = $query <=> $array->[$center]->[0] || ($center == 0 || $query != $array->[$center-1]->[0] ? 0 : -1);
+return $center if $cmp == 0;
+if ($center == $prevCenter) {
+return $left;
+}
+$right = $center if $cmp < 0;
+$left  = $center if $cmp > 0;
+$prevCenter = $center;
+}
+}

Mercurial > repos > yusuf > poor_gene_coverage

comparison vcf2hgvs_table @ 0:7cdd13ff182a default tip