--roi-file

--reference-sequence

--normal-tumor-bam-pair

--output-file

EOS ); }

sub _doc_authors { return " Cyriac Kandoth, Ph.D."; }

sub _doc_see_also { return <<EOS genome-music-bmr(1), genome-music(1), genome(1) EOS }

sub execute { my $self = shift; my $roi_file = $self->roi_file; my $ref_seq = $self->reference_sequence; my $tumor_bam = $self->tumor_bam; my $normal_bam = $self->normal_bam; my $output_file = $self->final_output_file; my $normal_min_depth = $self->normal_min_depth; my $tumor_min_depth = $self->tumor_min_depth; my $min_mapq = $self->min_mapq;

  # Check on all the input data before starting work
  print STDERR "ROI file not found or is empty: $roi_file\n" unless( -s $roi_file );
  print STDERR "Reference sequence file not found: $ref_seq\n" unless( -e $ref_seq );
  print STDERR "Normal BAM file not found or is empty: $normal_bam\n" unless( -s $normal_bam );
  print STDERR "Tumor BAM file not found or is empty: $tumor_bam\n" unless( -s $tumor_bam );
  return undef unless( -s $roi_file && -e $ref_seq && -s $normal_bam && -s $tumor_bam );

  # Check whether the annotated regions of interest are clumped together by chromosome
  my $roiFh = IO::File->new( $roi_file ) or die "ROI file could not be opened. $!\n";
  my @chroms = ( "" );
  while( my $line = $roiFh->getline ) # Emulate Unix's uniq command on the chromosome column
  {
    my ( $chrom ) = ( $line =~ m/^(\S+)/ );
    push( @chroms, $chrom ) if( $chrom ne $chroms[-1] );
  }
  $roiFh->close;
  my %chroms = map { $_ => 1 } @chroms; # Get the actual number of unique chromosomes
  if( scalar( @chroms ) != scalar( keys %chroms ))
  {
    print STDERR "ROIs from the same chromosome must be listed adjacent to each other in file. ";
    print STDERR "If in UNIX, try:\nsort -k 1,1 $roi_file\n";
    return undef;
  }

  # If the reference sequence FASTA file hasn't been indexed, do it
  my $ref_seq_idx = "$ref_seq.fai";
  system( "samtools faidx $ref_seq" ) unless( -e $ref_seq_idx );

  $normal_bam = '' unless( defined $normal_bam );
  $tumor_bam = '' unless( defined $tumor_bam );
  print STDERR "Normal BAM not found: \"$normal_bam\"\n" unless( -e $normal_bam );
  print STDERR "Tumor BAM not found: \"$tumor_bam\"\n" unless( -e $tumor_bam );
  next unless( -e $normal_bam && -e $tumor_bam );

  # Construct the command that calculates coverage per ROI
  my $calcRoiCovg_cmd = "calcRoiCovg $normal_bam $tumor_bam $roi_file $ref_seq $output_file $normal_min_depth $tumor_min_depth $min_mapq";

  # If the calcRoiCovg output was already generated, then don't rerun it
  if( -s $output_file )
  {
    print "Output file $output_file found. Skipping re-calculation.\n";
  }
  # Run the calcRoiCovg command on this tumor-normal pair. This could take a while
  elsif( system( "$calcRoiCovg_cmd" ) != 0 )
  {
    print STDERR "Failed to execute: $calcRoiCovg_cmd\n";
    return;
  }
  else
  {
    print "$output_file generated and stored.\n";
    return 1;
  }

}

1;