Merge pull request #934 from Clinical-Genomics/subsample_mt_mast

Subsample mt mast
Clinical-Genomics · Aug 27, 2019 · 8b04cfa · 8b04cfa
2 parents 0194e1e + cba7012
commit 8b04cfa
Show file tree

Hide file tree

Showing 5 changed files with 70 additions and 52 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -1,6 +1,10 @@
 # Change Log
 All notable changes to this project will be documented in this file.
 This project adheres to [Semantic Versioning](http://semver.org/).
+
+## [7.1.2]
+- Update samtools_subsample_mt to fix bug in downsampling of MT bam
+
 ## [7.1.1]
 - Fixed bug when skipping evaluation in QC_collect
 

diff --git a/lib/MIP/Constants.pm b/lib/MIP/Constants.pm
@@ -60,7 +60,7 @@ BEGIN {
 ## Constants
 ## Set MIP version
 ## Constants
-Readonly our $MIP_VERSION => q{v7.1.1};
+Readonly our $MIP_VERSION => q{v7.1.2};
 
 ## Cli
 Readonly our $MOOSEX_APP_SCEEN_WIDTH => 160;

diff --git a/lib/MIP/Program/Alignment/Bedtools.pm b/lib/MIP/Program/Alignment/Bedtools.pm
@@ -14,14 +14,15 @@ use warnings qw{ FATAL utf8 };
 use Readonly;
 
 ## MIPs lib/
+use MIP::Constants qw{ $SPACE };
 use MIP::Unix::Standard_streams qw{unix_standard_streams};
 use MIP::Unix::Write_to_file qw{unix_write_to_file};
 
 BEGIN {
     require Exporter;
 
     # Set the version for version checking
-    our $VERSION = 1.02;
+    our $VERSION = 1.03;
 
     # Inherit from Exporter to export functions and variables
     use base qw { Exporter };
@@ -31,14 +32,13 @@ BEGIN {
 
 }
 
-## Constants
-Readonly my $SPACE => q{ };
-
 sub bedtools_genomecov {
 
 ## Function : Perl wrapper for writing bedtools genomecov recipe to $FILEHANDLE. Based on bedtools 2.26.0.
 ## Returns  : "@commands"
-## Arguments: $FILEHANDLE             => Sbatch filehandle to write to
+## Arguments: $bam_infile_path        => BAM infile path
+##          : $depth_each_position    => Report the depth at each genome position (with zero-based coordinates)
+##          : $FILEHANDLE             => Sbatch filehandle to write to
 ##          : $infile_path            => Infile paths
 ##          : $max_coverage           => Combine all positions with a depth >= max into a single bin in the histogram
 ##          : $referencefile_path     => Genome reference file
@@ -49,56 +49,75 @@ sub bedtools_genomecov {
     my ($arg_href) = @_;
 
     ## Flatten argument(s)
+    my $bam_infile_path;
     my $FILEHANDLE;
     my $infile_path;
+    my $max_coverage;
     my $referencefile_path;
     my $stderrfile_path;
     my $stderrfile_path_append;
     my $stdoutfile_path;
 
     ## Default(s)
-    my $max_coverage;
+    my $depth_each_position;
 
     my $tmpl = {
+        bam_infile_path => {
+            store       => \$bam_infile_path,
+            strict_type => 1,
+        },
+        depth_each_position => {
+            allow       => qr/\A \d+ \z/sxm,
+            default     => 0,
+            store       => \$depth_each_position,
+            strict_type => 1,
+        },
         FILEHANDLE  => { store => \$FILEHANDLE },
         infile_path => {
-            required    => 1,
-            defined     => 1,
             strict_type => 1,
             store       => \$infile_path
         },
         max_coverage => {
-            allow       => qr/^\d+$/,
+            allow       => qr/\A \d+ \z/sxm,
             strict_type => 1,
             store       => \$max_coverage
         },
         referencefile_path => {
-            required    => 1,
-            defined     => 1,
             strict_type => 1,
             store       => \$referencefile_path
         },
-        stderrfile_path => { strict_type => 1, store => \$stderrfile_path },
-        stderrfile_path_append =>
-          { strict_type => 1, store => \$stderrfile_path_append },
-        stdoutfile_path => { strict_type => 1, store => \$stdoutfile_path },
+        stderrfile_path        => { strict_type => 1, store => \$stderrfile_path },
+        stderrfile_path_append => { strict_type => 1, store => \$stderrfile_path_append },
+        stdoutfile_path        => { strict_type => 1, store => \$stdoutfile_path },
     };
 
     check( $tmpl, $arg_href, 1 ) or croak qw{Could not parse arguments!};
 
     ## Array @commands stores commands depending on input parameters
-    my @commands = qw{bedtools genomecov};
+    my @commands = qw{ bedtools genomecov };
 
     ## Options
+    if ($bam_infile_path) {
+
+        push @commands, q{-ibam} . $SPACE . $bam_infile_path;
+    }
+    if ($infile_path) {
+
+        push @commands, q{-i} . $SPACE . $infile_path;
+    }
+
     if ($max_coverage) {
 
         push @commands, q{-max} . $SPACE . $max_coverage;
     }
+    if ($depth_each_position) {
 
-    ## Infile
-    push @commands, q{-ibam} . $SPACE . $infile_path;
+        push @commands, q{-dz};
+    }
+    if ($referencefile_path) {
 
-    push @commands, q{-g} . $SPACE . $referencefile_path;
+        push @commands, q{-g} . $SPACE . $referencefile_path;
+    }
 
     # Redirect stderr output to program specific stderr file
     push @commands,
@@ -113,8 +132,8 @@ sub bedtools_genomecov {
     unix_write_to_file(
         {
             commands_ref => \@commands,
-            separator    => $SPACE,
             FILEHANDLE   => $FILEHANDLE,
+            separator    => $SPACE,
         }
     );
     return @commands;

diff --git a/lib/MIP/Recipes/Analysis/Samtools_subsample_mt.pm b/lib/MIP/Recipes/Analysis/Samtools_subsample_mt.pm
@@ -26,15 +26,14 @@ BEGIN {
     use base qw{ Exporter };
 
     # Set the version for version checking
-    our $VERSION = 1.09;
+    our $VERSION = 1.10;
 
     # Functions and variables which can be optionally exported
     our @EXPORT_OK = qw{ analysis_samtools_subsample_mt };
 
 }
 
 ## Constants
-Readonly my $MAX_DEPTH_TRESHOLD          => 500_000;
 Readonly my $MAX_LIMIT_SEED              => 100;
 Readonly my $SAMTOOLS_UNMAPPED_READ_FLAG => 4;
 
@@ -142,8 +141,8 @@ sub analysis_samtools_subsample_mt {
     use MIP::Get::Parameter qw{ get_recipe_attributes get_recipe_resources };
     use MIP::Language::Awk qw{ awk };
     use MIP::Parse::File qw{ parse_io_outfiles };
-    use MIP::Program::Alignment::Samtools
-      qw{ samtools_depth samtools_index samtools_view };
+    use MIP::Program::Alignment::Samtools qw{ samtools_index samtools_view };
+    use MIP::Program::Alignment::Bedtools qw{ bedtools_genomecov };
     use MIP::Processmanagement::Processes qw{ submit_recipe };
     use MIP::Sample_info qw{ set_recipe_outfile_in_sample_info };
     use MIP::Script::Setup_script qw{ setup_script };
@@ -237,11 +236,11 @@ sub analysis_samtools_subsample_mt {
     print {$FILEHANDLE} q{MT_COVERAGE=} . $BACKTICK;
 
     # Get depth per base
-    samtools_depth(
+    bedtools_genomecov(
         {
-            FILEHANDLE         => $FILEHANDLE,
-            infile_path        => $infile_path,
-            max_depth_treshold => $MAX_DEPTH_TRESHOLD,
+            depth_each_position => 1,
+            FILEHANDLE          => $FILEHANDLE,
+            bam_infile_path     => $infile_path,
         }
     );
 

diff --git a/t/bedtools_genomecov.t b/t/bedtools_genomecov.t
@@ -46,11 +46,7 @@ GetOptions(
     # Display version number
     q{v|version} => sub {
         done_testing();
-        say {*STDOUT} $NEWLINE
-          . basename($PROGRAM_NAME)
-          . $SPACE
-          . $VERSION
-          . $NEWLINE;
+        say {*STDOUT} $NEWLINE . basename($PROGRAM_NAME) . $SPACE . $VERSION . $NEWLINE;
         exit;
     },
     q{vb|verbose} => $VERBOSE,
@@ -99,7 +95,7 @@ diag(   q{Test bedtools_genomecov from Bedtools.pm v}
       . $EXECUTABLE_NAME );
 
 ## Base arguments
-my @function_base_commands = qw{ bedtools };
+my @function_base_commands = qw{ bedtools genomecov };
 
 my %base_argument = (
     FILEHANDLE => {
@@ -126,30 +122,29 @@ my %required_argument = (
         input           => undef,
         expected_output => \@function_base_commands,
     },
-    infile_path => {
-        input           => q{infile.test},
-        expected_output => q{-ibam infile.test},
-    },
-    referencefile_path => {
-        input           => q{pathToRef.test},
-        expected_output => q{-g pathToRef.test},
-
-    },
 );
 
 ## Specific arguments
 my %specific_argument = (
+    bam_infile_path => {
+        input           => q{infile.bam},
+        expected_output => q{-ibam infile.bam},
+    },
+    depth_each_position => {
+        input           => 1,
+        expected_output => q{-dz},
+    },
     infile_path => {
-        input           => q{infile.test},
-        expected_output => q{-ibam infile.test},
+        input           => q{infile.bed},
+        expected_output => q{-i infile.bed},
     },
     max_coverage => {
-        input           => q{500},
-        expected_output => q{-max 500},
+        input           => 2,
+        expected_output => q{-max 2},
     },
     referencefile_path => {
-        input           => q{pathToRef.test},
-        expected_output => q{-g pathToRef.test},
+        input           => q{path_to_ref.test},
+        expected_output => q{-g path_to_ref.test},
 
     },
 );
@@ -165,9 +160,10 @@ foreach my $argument_href (@arguments) {
     my @commands = test_function(
         {
             argument_href              => $argument_href,
-            required_argument_href     => \%required_argument,
-            module_function_cref       => $module_function_cref,
+            do_test_base_command       => 1,
             function_base_commands_ref => \@function_base_commands,
+            module_function_cref       => $module_function_cref,
+            required_argument_href     => \%required_argument,
         }
     );
 }