From ac794b8af19421f1e5498471eef46817ea190334 Mon Sep 17 00:00:00 2001 From: Henrik Stranneheim Date: Tue, 27 Aug 2019 09:49:35 +0200 Subject: [PATCH 1/2] feat(subsample_mt): Fixed bug in downsampling of MT bam for IGV --- CHANGELOG.md | 3 + lib/MIP/Program/Alignment/Bedtools.pm | 59 ++++++++++++------- .../Recipes/Analysis/Samtools_subsample_mt.pm | 15 +++-- t/bedtools_genomecov.t | 42 ++++++------- 4 files changed, 68 insertions(+), 51 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 862a76741..517e58023 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,9 @@ This project adheres to [Semantic Versioning](http://semver.org/). ## [7.1.1] - Fixed bug when skipping evaluation in QC_collect +## [7.1.2] +- Update samtools_subsample_mt to fix bug in downsampling of MT bam + ## [7.1.0] - Updated TIDDIT to enable faster processing - Updated GATK for faster haplotypecalling diff --git a/lib/MIP/Program/Alignment/Bedtools.pm b/lib/MIP/Program/Alignment/Bedtools.pm index 321d48374..63de157d9 100644 --- a/lib/MIP/Program/Alignment/Bedtools.pm +++ b/lib/MIP/Program/Alignment/Bedtools.pm @@ -14,6 +14,7 @@ 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}; @@ -21,7 +22,7 @@ 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 index e3c9cb74c..6cee59c29 100644 --- a/lib/MIP/Recipes/Analysis/Samtools_subsample_mt.pm +++ b/lib/MIP/Recipes/Analysis/Samtools_subsample_mt.pm @@ -26,7 +26,7 @@ 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 }; @@ -34,7 +34,6 @@ BEGIN { } ## 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::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 index 1d9190326..d830837f9 100644 --- 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, } ); } From cba70121553ab35abdf4f90f31d5c4fb5ff5be7f Mon Sep 17 00:00:00 2001 From: Henrik Stranneheim Date: Tue, 27 Aug 2019 09:55:25 +0200 Subject: [PATCH 2/2] feat(subsample_mt): Bumped version and fixed import --- CHANGELOG.md | 5 +++-- lib/MIP/Constants.pm | 2 +- lib/MIP/Recipes/Analysis/Samtools_subsample_mt.pm | 2 +- 3 files changed, 5 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 517e58023..9d42e8ef5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,12 +1,13 @@ # 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.1] -- Fixed bug when skipping evaluation in QC_collect ## [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 + ## [7.1.0] - Updated TIDDIT to enable faster processing - Updated GATK for faster haplotypecalling diff --git a/lib/MIP/Constants.pm b/lib/MIP/Constants.pm index f022da80a..f4112d284 100644 --- 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/Recipes/Analysis/Samtools_subsample_mt.pm b/lib/MIP/Recipes/Analysis/Samtools_subsample_mt.pm index 6cee59c29..950e9ebb3 100644 --- a/lib/MIP/Recipes/Analysis/Samtools_subsample_mt.pm +++ b/lib/MIP/Recipes/Analysis/Samtools_subsample_mt.pm @@ -142,7 +142,7 @@ sub analysis_samtools_subsample_mt { use MIP::Language::Awk qw{ awk }; use MIP::Parse::File qw{ parse_io_outfiles }; use MIP::Program::Alignment::Samtools qw{ samtools_index samtools_view }; - use MIP::Program::Bedtools qw{ bedtools_genomecov }; + 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 };