comparison annotate_low_coverage @ 0:7cdd13ff182a default tip

initial commit

0:7cdd13ff182a

#!/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";
  }
}