Add inverted soft clipped masking

Author: Daniel Cooke

src/config/option_collation.cpp

@@ -691,7 +691,7 @@ bool allow_assembler_generation(const OptionMap& options)
     return options.at("assembly-candidate-generator").as<bool>() && !is_fast_mode(options);
 }
 
-auto make_read_transformers(const OptionMap& options)
+auto make_read_transformers(const ReferenceGenome& reference, const OptionMap& options)
 {
     using namespace octopus::readpipe;
     ReadTransformer prefilter_transformer {}, postfilter_transformer {};
@@ -739,6 +739,9 @@ auto make_read_transformers(const OptionMap& options)
         if (options.at("overlap-masking").as<bool>()) {
             postfilter_transformer.add(MaskStrandOfDuplicatedBases {});
         }
+        if (options.at("mask-inverted-soft-clipping").as<bool>()) {
+            prefilter_transformer.add(MaskInvertedSoftClippedReadEnds {reference, 5, 300});
+        }
         prefilter_transformer.shrink_to_fit();
         postfilter_transformer.shrink_to_fit();
     }
@@ -834,9 +837,9 @@ boost::optional<readpipe::Downsampler> make_downsampler(const OptionMap& options
     return boost::none;
 }
 
-ReadPipe make_read_pipe(ReadManager& read_manager, std::vector<SampleName> samples, const OptionMap& options)
+ReadPipe make_read_pipe(ReadManager& read_manager, const ReferenceGenome& reference, std::vector<SampleName> samples, const OptionMap& options)
 {
-    auto transformers = make_read_transformers(options);
+    auto transformers = make_read_transformers(reference, options);
     if (transformers.second.num_transforms() > 0) {
         return ReadPipe {read_manager, std::move(transformers.first), make_read_filterer(options),
                          std::move(transformers.second), make_downsampler(options), std::move(samples)};
@@ -1970,10 +1973,10 @@ ReadPipe make_default_filter_read_pipe(ReadManager& read_manager, std::vector<Sa
     return ReadPipe {read_manager, std::move(transformer), std::move(filterer), boost::none, std::move(samples)};
 }
 
-ReadPipe make_call_filter_read_pipe(ReadManager& read_manager, std::vector<SampleName> samples, const OptionMap& options)
+ReadPipe make_call_filter_read_pipe(ReadManager& read_manager, const ReferenceGenome& reference, std::vector<SampleName> samples, const OptionMap& options)
 {
     if (use_calling_read_pipe_for_call_filtering(options)) {
-        return make_read_pipe(read_manager, std::move(samples), options);
+        return make_read_pipe(read_manager, reference, std::move(samples), options);
     } else {
         return make_default_filter_read_pipe(read_manager, std::move(samples));
     }

src/config/option_collation.hpp

@@ -51,7 +51,7 @@ boost::optional<std::vector<SampleName>> get_user_samples(const OptionMap& optio
 
 ReadManager make_read_manager(const OptionMap& options);
 
-ReadPipe make_read_pipe(ReadManager& read_manager, std::vector<SampleName> samples, const OptionMap& options);
+ReadPipe make_read_pipe(ReadManager& read_manager, const ReferenceGenome& reference, std::vector<SampleName> samples, const OptionMap& options);
 
 bool call_sites_only(const OptionMap& options);
 
@@ -73,7 +73,7 @@ bool use_calling_read_pipe_for_call_filtering(const OptionMap& options) noexcept
 
 bool keep_unfiltered_calls(const OptionMap& options) noexcept;
 
-ReadPipe make_call_filter_read_pipe(ReadManager& read_manager, std::vector<SampleName> samples, const OptionMap& options);
+ReadPipe make_call_filter_read_pipe(ReadManager& read_manager, const ReferenceGenome& reference, std::vector<SampleName> samples, const OptionMap& options);
 
 boost::optional<fs::path> get_output_path(const OptionMap& options);
 

src/config/option_parser.cpp

@@ -226,6 +226,10 @@ OptionMap parse_options(const int argc, const char** argv)
     ("overlap-masking",
      po::value<bool>()->default_value(true),
      "Enable read segment overlap masking")
+     
+    ("mask-inverted-soft-clipping",
+    po::value<bool>()->default_value(true),
+    "Enable read segment overlap masking")
     ;
 
     po::options_description filters("Read filtering");

src/core/calling_components.cpp

@@ -491,7 +491,7 @@ GenomeCallingComponents::Components::Components(ReferenceGenome&& reference, Rea
 , regions {get_search_regions(options, this->reference, this->read_manager)}
 , contigs {get_contigs(this->regions, this->reference, options::get_contig_output_order(options))}
 , reads_profile {profile_reads(this->samples, this->regions, this->read_manager)}
-, read_pipe {options::make_read_pipe(this->read_manager, this->samples, options)}
+, read_pipe {options::make_read_pipe(this->read_manager, this->reference, this->samples, options)}
 , caller_factory {options::make_caller_factory(this->reference, this->read_pipe, this->regions, options, this->reads_profile)}
 , filter_read_pipe {}
 , output {std::move(output)}
@@ -578,7 +578,7 @@ void GenomeCallingComponents::Components::setup_writers(const options::OptionMap
 void GenomeCallingComponents::Components::setup_filter_read_pipe(const options::OptionMap& options)
 {
     if (!options::use_calling_read_pipe_for_call_filtering(options)) {
-        filter_read_pipe = options::make_call_filter_read_pipe(read_manager, samples, options);
+        filter_read_pipe = options::make_call_filter_read_pipe(read_manager, reference, samples, options);
     }
 }
 

src/readpipe/transformers/read_transform.cpp

@@ -9,6 +9,7 @@
 #include <cassert>
 
 #include "utils/maths.hpp"
+#include "utils/sequence_utils.hpp"
 
 namespace octopus { namespace readpipe {
 
@@ -164,6 +165,17 @@ MaskLowAverageQualitySoftClippedTails::MaskLowAverageQualitySoftClippedTails(Bas
 
 namespace {
 
+auto get_soft_clip_head_size(const AlignedRead& read) noexcept
+{
+    CigarOperation::Size front_size, back_size;
+    std::tie(front_size, back_size) = get_soft_clipped_sizes(read);
+    if (read.is_marked_reverse_mapped()) {
+        return back_size;
+    } else {
+        return front_size;
+    }
+}
+
 auto get_soft_clip_tail_size(const AlignedRead& read) noexcept
 {
     CigarOperation::Size front_size, back_size;
@@ -192,6 +204,15 @@ auto mean_tail_quality(const AlignedRead& read, const std::size_t num_bases) noe
     }
 }
 
+void zero_head_base_qualities(AlignedRead& read, const std::size_t num_bases) noexcept
+{
+    if (read.is_marked_reverse_mapped()) {
+        zero_back_qualities(read, num_bases);
+    } else {
+        zero_front_qualities(read, num_bases);
+    }
+}
+
 void zero_tail_base_qualities(AlignedRead& read, const std::size_t num_bases) noexcept
 {
     if (read.is_marked_reverse_mapped()) {
@@ -214,6 +235,78 @@ void MaskLowAverageQualitySoftClippedTails::operator()(AlignedRead& read) const
     }
 }
 
+MaskInvertedSoftClippedReadEnds::MaskInvertedSoftClippedReadEnds(const ReferenceGenome& reference,
+                                                                 AlignedRead::NucleotideSequence::size_type min_clip_length,
+                                                                 GenomicRegion::Size max_flank_search)
+: reference_ {reference}
+, min_clip_length_ {min_clip_length}
+, max_flank_search_ {max_flank_search}
+{}
+
+namespace {
+
+auto copy_head_sequence(const AlignedRead& read, const AlignedRead::NucleotideSequence::size_type length)
+{
+    if (length >= sequence_size(read)) return read.sequence();
+    if (read.is_marked_reverse_mapped()) {
+        return AlignedRead::NucleotideSequence {std::prev(std::cend(read.sequence()), length), std::cend(read.sequence())};
+    } else {
+        return AlignedRead::NucleotideSequence {std::cbegin(read.sequence()), std::next(std::cbegin(read.sequence()), length)};
+    }
+}
+
+auto copy_tail_sequence(const AlignedRead& read, const AlignedRead::NucleotideSequence::size_type length)
+{
+    if (length >= sequence_size(read)) return read.sequence();
+    if (read.is_marked_reverse_mapped()) {
+        return AlignedRead::NucleotideSequence {std::cbegin(read.sequence()), std::next(std::cbegin(read.sequence()), length)};
+    } else {
+        return AlignedRead::NucleotideSequence {std::prev(std::cend(read.sequence()), length), std::cend(read.sequence())};
+    }
+}
+
+auto copy_soft_clipped_head_sequence(const AlignedRead& read)
+{
+    return copy_head_sequence(read, get_soft_clip_head_size(read));
+}
+
+auto copy_soft_clipped_tail_sequence(const AlignedRead& read)
+{
+    return copy_tail_sequence(read, get_soft_clip_tail_size(read));
+}
+
+template <typename Range1, typename Range2>
+bool includes(const Range1& target, const Range2& query)
+{
+    return std::search(std::cbegin(target), std::cend(target), std::cbegin(query), std::cend(query)) != std::cend(target);
+}
+
+} // namespace
+
+void MaskInvertedSoftClippedReadEnds::operator()(AlignedRead& read) const
+{
+    if (is_soft_clipped(read)) {
+        const auto soft_clipped_head_length = get_soft_clip_head_size(read);
+        if (soft_clipped_head_length >= min_clip_length_) {
+            auto query = copy_head_sequence(read, soft_clipped_head_length);
+            utils::reverse_complement(query);
+            auto target = reference_.get().fetch_sequence(expand(mapped_region(read), max_flank_search_));
+            if (includes(target, query)) {
+                zero_head_base_qualities(read, soft_clipped_head_length);
+            }
+        }
+        const auto soft_clipped_tail_length = get_soft_clip_tail_size(read);
+        if (soft_clipped_tail_length >= min_clip_length_) {
+            auto query = copy_tail_sequence(read, soft_clipped_tail_length);
+            utils::reverse_complement(query);
+            auto target = reference_.get().fetch_sequence(expand(mapped_region(read), max_flank_search_));
+            if (includes(target, query)) {
+                zero_tail_base_qualities(read, soft_clipped_tail_length);
+            }
+        }
+    }
+}
+
 // template transforms
 
 void mask_adapter_contamination(AlignedRead& forward, AlignedRead& reverse) noexcept

src/readpipe/transformers/read_transform.hpp

@@ -8,6 +8,7 @@
 #include <vector>
 
 #include "basics/aligned_read.hpp"
+#include "io/reference/reference_genome.hpp"
 
 namespace octopus { namespace readpipe {
 
@@ -133,6 +134,20 @@ struct MaskLowAverageQualitySoftClippedTails
     const Length min_tail_length_;
 };
 
+struct MaskInvertedSoftClippedReadEnds
+{
+    MaskInvertedSoftClippedReadEnds(const ReferenceGenome& reference,
+                                    AlignedRead::NucleotideSequence::size_type min_clip_length,
+                                    GenomicRegion::Size max_flank_search);
+    
+    void operator()(AlignedRead& read) const;
+    
+private:
+    std::reference_wrapper<const ReferenceGenome> reference_;
+    AlignedRead::NucleotideSequence::size_type min_clip_length_;
+    GenomicRegion::Size max_flank_search_;
+};
+
 using ReadReferenceVector = std::vector<std::reference_wrapper<AlignedRead>>;
 
 struct MaskTemplateAdapters