--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/annotate_low_coverage	Wed Mar 25 15:49:28 2015 -0600
@@ -0,0 +1,194 @@
+#!/usr/bin/env perl
+
+use strict;
+use warnings;
+use File::Basename;
+#@ARGV >= 4 and @ARGV <= 6 or die "Usage: $0 <capture kit.bed> <poor regions input.bed> <annotated output.hgvs.txt> [target gene list.txt] [gene name pattern for detail reporting]\n";
+
+my $dirname = dirname(__FILE__);
+my $tool_dir = shift @ARGV;
+
+#read config file
+if(not -e "$tool_dir/report_poor_converage" ){
+  system("cp $dirname/tool-data/report_poor_coverage $tool_dir/report_poor_coverage");
+}
+my %configuration_file;
+open FILE, '<', "$tool_dir/report_poor_coverage"; 
+while(<FILE>){
+	(my $key, my $value) = split(/\s+/,$_);
+	$configuration_file{$key} = $value;
+}
+close FILE;
+my $ref_fasta = $configuration_file{"ref_fasta"};
+my $ref_flat = $configuration_file{"ref_flat"};
+my $ref_gencode = $configuration_file{"ref_gencode"};
+my $clinVar = $configuration_file{"clinVar"};
+my $dgv2 = $configuration_file{"dgv2"};
+
+my $capture_kit_bed = $configuration_file{"capturekits_directory"} . (shift @ARGV) . ".bed" ;
+my $low_input_bed = shift @ARGV;
+my $low_output_hgvs = shift @ARGV;
+my $target_list = @ARGV ? shift @ARGV : undef;
+my $gene_name_regex = @ARGV ? shift @ARGV : undef;
+my $tmp_hgvs = "$$.tmp.hgvs.txt";
+
+my $tmp_bed;
+if(defined $target_list and $target_list ne "/dev/null"){
+  $tmp_bed = "$$.targeted.tmp.bed";
+  #print STDERR "filter_by_list false $low_input_bed $target_list $tmp_bed 1 3\n";
+  system "$dirname/add_names_to_bed $low_input_bed $ref_flat - /dev/null | $dirname/filter_by_list false - $target_list $tmp_bed 1 3"; # 3 indicates 4th column of the bed file, which should have the gene name
+  $low_input_bed = $tmp_bed;
+  open(TARGETS, "$dirname/add_names_to_bed $capture_kit_bed $ref_flat - /dev/null | $dirname/filter_by_list false - $target_list - 1 3 |")
+    or die "Cannot run filter_by_list, aborting: $!\n";
+}
+else{
+  open(TARGETS, "$dirname/add_names_to_bed $capture_kit_bed $ref_flat - /dev/null |")
+    or die "Cannot run add_names_to_bed on $capture_kit_bed: $!\n";
+}
+
+#print STDERR "/export/achri_galaxy/galaxy-dist/tools/achri/vcf2hgvs_table -q -p 0.05 -c $low_input_bed -g /export/achri_galaxy/dbs/DGV2/hg19.2012-03-29.txt.gz -o - -r /export/achri_galaxy/dbs/hg19.fa -e /export/achri_galaxy/dbs/hg19_refGene_gencode_ultrasensitive.gtf -i /dev/null -b /export/achri_galaxy/dbs/hg19_refFlat_2014-07-10.bed | /export/achri_galaxy/galaxy-dist/tools/achri/filter_by_index_gamma /export/achri_galaxy/dbs/ClinVar/ClinVarFullRelease_2014-03.xml. ClinVar - $tmp_hgvs \"\"";
+system "$dirname/vcf2hgvs_table -q -p 0.05 -c $low_input_bed -g $dgv2 -o - -r $ref_fasta -e $ref_gencode -i /dev/null -b $ref_flat | $dirname/filter_by_index_gamma $clinVar. ClinVar - $tmp_hgvs \"\"";
+my $tmp_output = "$$.tmp.annotated.hgvs.txt";
+#system "/export/achri_galaxy/galaxy-dist/tools/achri/hgvs_table_annotate /export/achri_galaxy/dbs/sift/hg19 /export/achri_galaxy/dbs/polyphen2/hg19.txt.gz /export/achri_galaxy/dbs/gerp/hg19 /export/achri_galaxy/dbs/TissueDistributionDBs/human.v2009-07-30.tab /export/achri_galaxy/dbs/pathways/KEGG.human.2012-09-25.txt /export/achri_galaxy/dbs/interpro_supermatch_hg19.bed $tmp_hgvs $tmp_output /export/achri_galaxy/dbs/hg19.fa";
+open(OUT, ">$low_output_hgvs")
+  or die "Cannot open $low_output_hgvs for writing: $!\n";
+open(COLLAPSE, "$dirname/hgvs_collapse_transcripts $tmp_hgvs - 1 |")
+  or die "cannot run hgvs_collapse_transcripts: $!\n";
+my $header = <COLLAPSE>;
+chomp $header;
+my @headers = split /\t/, $header;
+my ($chr_column, $from_column, $to_column, $ftype_column, $clinvar_column, $loc_column, $maf_column, $gene_column, $gene_desc_column, $cdna_hgvs_column, $transcript_column, $zygosity_column, $phase_column);
+my %req_columns = (
+        "Chr" => \$chr_column,
+        "DNA From" => \$from_column,
+        "DNA To" => \$to_column,
+        "Feature type" => \$ftype_column,
+#        "Variant context" => \$loc_column,
+        "Pop. freq." => \$maf_column,
+        "Gene Name" => \$gene_column,
+#        "Gene Function" => \$gene_desc_column,
+        "ClinVar Text Matches" => \$clinvar_column,
+        "Transcript HGVS" => \$cdna_hgvs_column,
+        "Selected transcript" => \$transcript_column,
+        "Zygosity" => \$zygosity_column,
+        "Phase" => \$phase_column);
+&load_header_columns(\%req_columns, \@headers);
+my %col2name = reverse %req_columns;
+
+my %target_sizes;
+<TARGETS>; # header
+while(<TARGETS>){
+  chomp;
+  my @F = split /\t/, $_;
+  next unless $#F >= 3; # nameless records
+  for my $gene_name (split /;\s*/, $F[3]){
+    #print STDERR "Processing $_\n";
+    $target_sizes{$gene_name} += $F[2]-$F[1]+1; # assumes non-overlapping targets in BED file
+  }
+}
+close(TARGETS);
+my $target_total = 0;
+for my $gene_name (keys %target_sizes){
+  #print STDERR "Adding size of $gene_name ($target_sizes{$gene_name})\n";
+  $target_total += $target_sizes{$gene_name};
+}
+
+my %missing_homo;
+my %missing_coding_homo;
+my %missing_het;
+my %missing_coding_het;
+my @detail_lines;
+while(<COLLAPSE>){
+  chomp;
+  my @F = split /\t/, $_, -1;
+ 
+  for my $gene_name (split /;\s*/, uc($F[$gene_column])){
+    next if not exists $target_sizes{$gene_name}; # not a gene of interest
+    #print STDERR "Processing $_\n";
+    next if defined $gene_name_regex and $F[$transcript_column] !~ /$gene_name_regex/o;
+    #print STDERR "Passes filters\n";
+    if(not exists $missing_homo{$gene_name}){
+      $missing_homo{$gene_name} = {};
+      $missing_coding_homo{$gene_name} = {};
+      $missing_het{$gene_name} = {};
+      $missing_coding_het{$gene_name} = {};
+    }
+ 
+    if($F[$zygosity_column] =~ /homo/){
+      for my $pos ($F[$from_column]..$F[$to_column]){
+        $missing_homo{$gene_name}->{$pos}++;
+        $missing_coding_homo{$gene_name}->{$pos}++ if $F[$ftype_column] =~ /protein_coding/;# and $F[$loc_column] =~ /exon/;
+      }
+    }
+    else{
+      for my $pos ($F[$from_column]..$F[$to_column]){
+        $missing_het{$gene_name}->{$pos}++;
+        $missing_coding_het{$gene_name}->{$pos}++ if $F[$ftype_column] =~ /protein_coding/;# and $F[$loc_column] =~ /exon/;
+        #print STDERR "Missing $gene_name het $pos ($F[$from_column]..$F[$to_column]) $F[$ftype_column] $F[$loc_column]\n";
+      }
+    }
+    next unless $F[$ftype_column] =~ /protein_coding/;# and $F[$loc_column] =~ /exon/;
+    #push @detail_lines, [$F[$chr_column], $F[$from_column], $F[$to_column], $F[$maf_column], $gene_name, $F[$cdna_hgvs_column], $F[$transcript_column], ($F[$zygosity_column] =~ /homo/ ? "<4x" : "<20x"), $F[$gene_desc_column], $F[$domain_column], $F[$gene_column]];
+    $F[$clinvar_column] = "" if not defined $F[$clinvar_column];
+    $F[$clinvar_column] =~ s/&gt;/>/;
+    push @detail_lines, [$F[$chr_column], $F[$from_column], $F[$to_column], $F[$maf_column], $gene_name, $F[$cdna_hgvs_column], $F[$transcript_column], ($F[$zygosity_column] =~ /homo/ ? "<4x" : "<20x"), $F[$gene_column], $F[$clinvar_column]];
+  }
+}
+my $missing_homo = 0;
+my $missing_coding_homo = 0;
+my $missing_het = 0;
+my $missing_coding_het = 0;
+for my $gene_name (keys %target_sizes){
+  $missing_homo{$gene_name} = defined $missing_homo{$gene_name} ? keys %{$missing_homo{$gene_name}} : 0;
+  $missing_homo += $missing_homo{$gene_name};
+  $missing_coding_homo{$gene_name} = defined $missing_coding_homo{$gene_name} ? keys %{$missing_coding_homo{$gene_name}} : 0;
+  $missing_coding_homo += $missing_coding_homo{$gene_name};
+  $missing_het{$gene_name} = defined $missing_het{$gene_name} ? keys %{$missing_het{$gene_name}} : 0;
+  $missing_het += $missing_het{$gene_name};
+  $missing_coding_het{$gene_name} = defined $missing_coding_het{$gene_name} ? keys %{$missing_coding_het{$gene_name}} : 0;
+  $missing_coding_het += $missing_coding_het{$gene_name};
+}
+if(defined $gene_name_regex){
+  print OUT "Note: only considering coding sequence target bases in gene models with IDs matching \"$gene_name_regex\"\n";
+}
+print OUT "\ttotal size\tcoding <4x coverage\t%\tcoding<20x coverage\t%\ttotal <4x coverage\t%\ttotal <20x coverage\n";
+printf OUT "Total target gene bases designed for capture\t%d\t%d\t%.1f\t%d\t%.1f\t%d\t%.1f\t%d\t%.1f\n", $target_total, $missing_coding_homo, $missing_coding_homo/$target_total*100, $missing_coding_het, $missing_coding_het/$target_total*100,
+            $missing_homo, $missing_homo/$target_total*100, $missing_het, $missing_het/$target_total*100;
+print OUT "Per gene missing coding sequence coverage breakdown:\n";
+for my $gene_name (sort keys %missing_homo){
+  printf OUT "%s\t%d\t%d\t%.1f\t%d\t%.1f\t%d\t%.1f\t%d\t%.1f\n", $gene_name, $target_sizes{$gene_name}, $missing_coding_homo{$gene_name},
+            $missing_coding_homo{$gene_name}/$target_sizes{$gene_name}*100, $missing_coding_het{$gene_name}, $missing_coding_het{$gene_name}/$target_sizes{$gene_name}*100,
+            $missing_homo{$gene_name}, $missing_homo{$gene_name}/$target_sizes{$gene_name}*100, 
+            $missing_het{$gene_name}, $missing_het{$gene_name}/$target_sizes{$gene_name}*100;
+}
+# Print revised header that's only a few of the fields
+@detail_lines = sort {$a->[4] cmp $b->[4] or $a->[1] <=> $b->[1]} @detail_lines;
+print OUT "\n\nCoding region low coverage details (alphabetical by gene name):\n";
+print OUT join("\t", $col2name{\$chr_column}, $col2name{\$from_column}, $col2name{\$to_column}, $col2name{\$maf_column}, $col2name{\$gene_column}, $col2name{\$cdna_hgvs_column}, 
+                     $col2name{\$transcript_column}, "Low coverage type", "Genes in this region (gene overlaps possible)", $col2name{\$clinvar_column}), "\n";
+print OUT join("\n", map {join("\t", @{$_})} @detail_lines), "\n";
+
+unlink($tmp_hgvs);
+unlink($tmp_output);
+unlink($tmp_bed) if defined $tmp_bed;
+
+sub load_header_columns{
+  my ($reqs_hash_ref, $headers_array_ref) = @_;
+  my %unfulfilled;
+  for my $field_name (keys %$reqs_hash_ref){
+    $unfulfilled{$field_name} = 1;
+  }
+  for(my $i = 0; $i <= $#{$headers_array_ref}; $i++){
+    for my $req_header_name (keys %$reqs_hash_ref){
+      if($req_header_name eq $headers_array_ref->[$i]){
+        ${$reqs_hash_ref->{$req_header_name}} = $i;
+        delete $unfulfilled{$req_header_name};
+        last;
+      }
+    }
+  }
+  if(keys %unfulfilled){
+    die "Aborting. Could not find headers in the input file for the following required fields: ", join(", ", sort keys %unfulfilled), "\n";
+  }
+}
+