Metabat override fix

src/bin/coverm.rs

@@ -1679,16 +1679,32 @@ impl FilterParameters {
 
     pub fn add_metabat_filtering_if_required(&mut self, m: &clap::ArgMatches) {
         if doing_metabat(m) {
+            // A threshold given with --min-read-percent-identity is kept; the
+            // 97% default only replaces an unset one.
+            if m.contains_id("min-read-percent-identity") {
+                info!(
+                    "Using the given single read percent identity threshold of {} for \
+                     MetaBAT adjusted coverage",
+                    self.min_percent_identity_single
+                );
+            } else {
+                // we use >= where metabat uses >. Gah.
+                self.min_percent_identity_single = 0.97001;
+                info!(
+                    "Setting single read percent identity threshold at 0.97 for \
+                     MetaBAT adjusted coverage"
+                );
+            }
+            // Improper pair and supplementary alignments are already included
+            // unless the user asked otherwise with --proper-pairs-only or
+            // --exclude-supplementary, which are respected. Secondary
+            // alignments are excluded by default, so include them here.
+            self.flag_filters.include_secondary = true;
             info!(
-                "Setting single read percent identity threshold at 0.97 for \
-                 MetaBAT adjusted coverage, and not filtering out supplementary, \
-                 secondary and improper pair alignments"
+                "For MetaBAT adjusted coverage, including secondary alignments; \
+                 supplementary alignments included: {}; improper pair alignments included: {}",
+                self.flag_filters.include_supplementary, self.flag_filters.include_improper_pairs
             );
-            // we use >= where metabat uses >. Gah.
-            self.min_percent_identity_single = 0.97001;
-            self.flag_filters.include_improper_pairs = true;
-            self.flag_filters.include_supplementary = true;
-            self.flag_filters.include_secondary = true;
         }
     }
 

src/cli.rs

@@ -208,8 +208,15 @@ fn add_thresholding_options(manual: Manual) -> Manual {
                     .long("--min-read-percent-identity")
                     .help(&format!(
                         "Exclude reads by overall percent \
-        identity e.g. 95 for 95%. {}",
-                        default_roff("0")
+        identity e.g. 95 for 95%. Reads at exactly the threshold \
+        are kept ({}). {} With {} in {} the default is {} \
+        (see the {} method); a given value is used as is.",
+                        &monospace_roff(">="),
+                        default_roff("0"),
+                        &monospace_roff("--methods metabat"),
+                        &monospace_roff("coverm contig"),
+                        &monospace_roff("97"),
+                        &monospace_roff("metabat")
                     )),
             )
             .option(
@@ -965,7 +972,13 @@ pub fn contig_full_help() -> Manual {
                     &[&monospace_roff("variance"), "Variance of coverage depths"],
                     &[&monospace_roff("length"), "Length of each contig in base pairs"],
                     &[&monospace_roff("count"), "Number of reads aligned to each contig. Note that supplementary alignments are not counted."],
-                    &[&monospace_roff("metabat"), "(\"MetaBAT adjusted coverage\") Coverage as defined in Kang et al 2015 https://doi.org/10.7717/peerj.1165"],
+                    &[&monospace_roff("metabat"), &format!("(\"MetaBAT adjusted coverage\") Coverage as defined in Kang et al 2015 https://doi.org/10.7717/peerj.1165. By default, reads with identity not strictly greater than 97% are excluded, matching MetaBAT's {}. A given {} uses CoverM's {} instead, so {} is not identical to the default. This method includes secondary alignments, and by default improper pair and supplementary alignments too; use {} or {} to exclude the latter two.",
+                        &monospace_roff(">"),
+                        &monospace_roff("--min-read-percent-identity"),
+                        &monospace_roff(">="),
+                        &monospace_roff("97"),
+                        &monospace_roff("--proper-pairs-only"),
+                        &monospace_roff("--exclude-supplementary"))],
                     &[&monospace_roff("reads_per_base"), "Number of reads aligned divided by the length of the contig"],
                     &[&monospace_roff("anir"), "Average BLAST-like identity of mapped reads"],
                     &[&monospace_roff("rpkm"), "Reads mapped per kilobase of contig, per million mapped reads"],

tests/test_cmdline.rs

@@ -1956,6 +1956,52 @@ k141_2005182	225	1.9333	1.9333	0.0631",
             .unwrap();
     }
 
+    #[test]
+    fn test_metabat_given_percent_identity_is_kept() {
+        // Reads below 97% identity, excluded by the metabat default, pass a
+        // given 90% threshold.
+        Assert::main_binary()
+            .with_args(&[
+                "contig",
+                "-m",
+                "metabat",
+                "--min-read-percent-identity",
+                "90",
+                "-b",
+                "tests/data/k141_2005182.head11.bam",
+            ])
+            .succeeds()
+            .stdout()
+            .contains(
+                "contigName	contigLen	totalAvgDepth	k141_2005182.head11.bam	k141_2005182.head11.bam-var
+k141_2005182	225	6.04	6.04	1.12",
+            )
+            .unwrap();
+    }
+
+    #[test]
+    fn test_metabat_given_97_uses_greater_or_equal() {
+        // A given 97 uses CoverM's >=, unlike the metabat default which
+        // approximates MetaBAT's > via 0.97001. Same BAM is 1.9333 by default.
+        Assert::main_binary()
+            .with_args(&[
+                "contig",
+                "-m",
+                "metabat",
+                "--min-read-percent-identity",
+                "97",
+                "-b",
+                "tests/data/k141_2005182.head11.bam",
+            ])
+            .succeeds()
+            .stdout()
+            .contains(
+                "contigName	contigLen	totalAvgDepth	k141_2005182.head11.bam	k141_2005182.head11.bam-var
+k141_2005182	225	2.9333	2.9333	0.0631",
+            )
+            .unwrap();
+    }
+
     #[test]
     fn test_deletions_count_towards_perc_id() {
         Assert::main_binary()