Add "full-bamout" option to output reference reads

src/config/option_collation.cpp

@@ -2061,6 +2061,11 @@ boost::optional<fs::path> split_bamout_request(const OptionMap& options)
     return boost::none;
 }
 
+bool full_bamouts_requested(const OptionMap& options)
+{
+    return options.at("full-bamout").as<bool>();
+}
+
 unsigned max_open_read_files(const OptionMap& options)
 {
     return 2 * std::min(as_unsigned("max-open-read-files", options), count_read_paths(options));

src/config/option_collation.hpp

@@ -88,6 +88,7 @@ bool annotate_filter_output(const OptionMap& options);
 
 boost::optional<fs::path> bamout_request(const OptionMap& options);
 boost::optional<fs::path> split_bamout_request(const OptionMap& options);
+bool full_bamouts_requested(const OptionMap& options);
 
 unsigned estimate_max_open_files(const OptionMap& options);
 

src/config/option_parser.cpp

@@ -184,13 +184,17 @@ OptionMap parse_options(const int argc, const char** argv)
      po::value<fs::path>(),
      "Output realigned BAM files")
 
-     ("split-bamout",
-      po::value<fs::path>(),
-      "Output split realigned BAM files")
+    ("split-bamout",
+     po::value<fs::path>(),
+     "Output split realigned BAM files")
 
-     ("data-profile",
-      po::value<fs::path>(),
-      "Output a profile of polymorphisms and errors found in the data")
+    ("full-bamout",
+     po::bool_switch()->default_value(false),
+     "Output all reads when producing realigned bam outputs rather than just variant read minibams")
+
+    ("data-profile",
+     po::value<fs::path>(),
+     "Output a profile of polymorphisms and errors found in the data")
     ;
 
     po::options_description transforms("Read transformations");

src/core/calling_components.cpp

@@ -168,6 +168,11 @@ boost::optional<GenomeCallingComponents::Path> GenomeCallingComponents::split_ba
     return components_.split_bamout;
 }
 
+bool GenomeCallingComponents::write_full_bamouts() const noexcept
+{
+    return components_.write_full_bamouts;
+}
+
 boost::optional<GenomeCallingComponents::Path> GenomeCallingComponents::data_profile() const
 {
     return components_.data_profile;
@@ -501,6 +506,7 @@ GenomeCallingComponents::Components::Components(ReferenceGenome&& reference, Rea
 , filter_request {}
 , bamout {options::bamout_request(options)}
 , split_bamout {options::split_bamout_request(options)}
+, write_full_bamouts {options::full_bamouts_requested(options)}
 , data_profile {options::data_profile_request(options)}
 {
     drop_unused_samples(this->samples, this->read_manager);

src/core/calling_components.hpp

@@ -72,6 +72,7 @@ class GenomeCallingComponents
     boost::optional<Path> filter_request() const;
     boost::optional<Path> bamout() const;
     boost::optional<Path> split_bamout() const;
+    bool write_full_bamouts() const noexcept;
     boost::optional<Path> data_profile() const;
 
 private:
@@ -108,7 +109,9 @@ class GenomeCallingComponents
         bool sites_only;
         boost::optional<Path> legacy;
         boost::optional<Path> filter_request;
-        boost::optional<Path> bamout, split_bamout, data_profile;
+        boost::optional<Path> bamout, split_bamout;
+        bool write_full_bamouts;
+        boost::optional<Path> data_profile;
         // Components that require temporary directory during construction appear last to make
         // exception handling easier.
         boost::optional<Path> temp_directory;

src/core/octopus.cpp

@@ -1551,9 +1551,11 @@ void run_bam_realign(GenomeCallingComponents& components)
     if (is_bam_realignment_requested(components)) {
         if (check_bam_realign(components)) {
             components.read_manager().close();
+            BAMRealigner::Config config {};
+            config.copy_hom_ref_reads = components.write_full_bamouts();
             if (components.read_manager().paths().size() == 1) {
                 realign(components.read_manager().paths().front(), get_bam_realignment_vcf(components),
-                        *components.bamout(), components.reference());
+                        *components.bamout(), components.reference(), config);
             } else {
                 namespace fs = boost::filesystem;
                 const auto bamout_directory = *components.bamout();
@@ -1574,7 +1576,7 @@ void run_bam_realign(GenomeCallingComponents& components)
                     auto bamout_path = bamout_directory;
                     bamout_path /= bamin_path.filename();
                     if (bamin_path != bamout_path) {
-                        realign(bamin_path, get_bam_realignment_vcf(components), bamout_path, components.reference());
+                        realign(bamin_path, get_bam_realignment_vcf(components), bamout_path, components.reference(), config);
                     } else {
                         logging::WarningLogger warn_log {};
                         stream(warn_log) << "Cannot make evidence bam " << bamout_path << " as it is an input bam";
@@ -1586,10 +1588,12 @@ void run_bam_realign(GenomeCallingComponents& components)
     if (is_split_bam_realignment_requested(components)) {
         if (check_bam_realign(components)) {
             components.read_manager().close();
+            BAMRealigner::Config config {};
+            config.copy_hom_ref_reads = components.write_full_bamouts();
             if (components.read_manager().paths().size() == 1) {
                 auto out_paths = get_haplotype_bam_paths(*components.split_bamout(), get_max_ploidy(components));
                 realign(components.read_manager().paths().front(), get_bam_realignment_vcf(components),
-                        std::move(out_paths), components.reference());
+                        std::move(out_paths), components.reference(), config);
             } else {
                 namespace fs = boost::filesystem;
                 const auto bamout_directory = *components.split_bamout();
@@ -1612,7 +1616,7 @@ void run_bam_realign(GenomeCallingComponents& components)
                     bamout_prefix /= bamin_path.filename().stem();
                     const auto max_ploidy = get_max_called_ploidy(realignment_vcf, bamin_path);
                     auto bamout_paths = get_haplotype_bam_paths(bamout_prefix, max_ploidy);
-                    realign(bamin_path, realignment_vcf, std::move(bamout_paths), components.reference());
+                    realign(bamin_path, realignment_vcf, std::move(bamout_paths), components.reference(), config);
                 }
             }
         }