Merge pull request #45 from bigbio/fix/sdrf-condition-mapping-and-directlfq-norm-warning

Fix features2proteins LFQ correctness: per-run intensity routing, SDRF condition mapping, DirectLFQ norm warning

Author: ypriverol

python/mokume/io/feature.py

@@ -215,22 +215,47 @@ def _detect_qpx_format(self) -> None:
         self._has_is_decoy = "is_decoy" in cols
         self._has_anchor_protein = "anchor_protein" in cols
 
+        # New-QPX LFQ vs TMT discriminator. In LFQ the ``intensities.label`` names
+        # the run each (possibly MBR-transferred) intensity belongs to, so it — not
+        # the row's anchor ``run_file_name`` — is the per-run/per-sample key. In TMT
+        # the label is a reporter channel and ``run_file_name`` is the run. Detect by
+        # checking whether every distinct label is a member of the run_file_name
+        # domain (true for LFQ, never for TMT channels).
+        self._label_is_run = False
+        if self._is_new_qpx:
+            try:
+                n_off = self.parquet_db.execute(
+                    "SELECT count(*) FROM ("
+                    " SELECT DISTINCT unnest.label AS lbl"
+                    " FROM parquet_db_raw, UNNEST(intensities) AS unnest"
+                    ") WHERE lbl NOT IN"
+                    " (SELECT DISTINCT run_file_name FROM parquet_db_raw)"
+                ).fetchone()[0]
+                self._label_is_run = n_off == 0
+            except duckdb.Error as exc:
+                logger.debug("label-is-run detection failed: %s", exc)
+
     def _create_unnest_view(self) -> None:
         """Create the long-format DuckDB view by unnesting intensities."""
         if self._is_new_qpx:
+            # LFQ: the intensity's owning run is its label (keep run_file_name only
+            # as provenance). TMT: label is a channel, run stays run_file_name.
+            run_key = "unnest.label" if self._label_is_run else "run_file_name"
             unnest_sql = (
-                "run_file_name as sample_accession,\n"
+                f"{run_key} as sample_accession,\n"
                 "                    unnest.label as channel,\n"
                 "                    unnest.intensity"
             )
-            sa_default = "run_file_name"
+            sa_default = run_key
+            run_expr = run_key
         else:
             unnest_sql = (
                 "unnest.sample_accession,\n"
                 "                    unnest.channel,\n"
                 "                    unnest.intensity"
             )
             sa_default = "unnest.sample_accession"
+            run_expr = self._run_col
 
         charge_col, run_col = self._charge_col, self._run_col
         # Normalize pg_accessions: extract accession strings from struct if needed
@@ -265,7 +290,7 @@ def _create_unnest_view(self) -> None:
                     unnest_sql,
                     ",",
                     " ",
-                    run_col,
+                    run_expr,
                     " as run,",
                     " ",
                     sa_default,

rust/crates/mokume-io/src/qpx.rs

@@ -140,20 +140,42 @@ pub fn flatten_qpx_batch(batch: &RecordBatch) -> Result<Vec<QpxFeatureRecord>> {
             .transpose()?
             .flatten();
 
-        for entry in intensity_entries(intensities, row)? {
+        let entries = intensity_entries(intensities, row)?;
+        // In quantms.io LFQ the per-run intensities are labeled by run file (the
+        // row's anchor run appears among them); in isobaric (TMT/iTRAQ) the labels
+        // are reporter channels and never equal a run file. When the labels are
+        // runs, the run/sample for each intensity is its own `label`, not the row's
+        // anchor `run_file_name` -- otherwise every run collapses onto the anchor
+        // sample and all other runs are dropped (see mokume-io qpx tests).
+        let row_key = crate::sdrf::normalize_file_key(&run_file_name);
+        let labels_are_runs = entries.iter().any(|entry| {
+            entry
+                .label
+                .as_deref()
+                .is_some_and(|label| crate::sdrf::normalize_file_key(label) == row_key)
+        });
+        for entry in entries {
             if entry.intensity.is_finite() {
+                let (entry_run_file_name, entry_label) = if labels_are_runs {
+                    (
+                        entry.label.clone().unwrap_or_else(|| run_file_name.clone()),
+                        None,
+                    )
+                } else {
+                    (run_file_name.clone(), entry.label)
+                };
                 records.push(QpxFeatureRecord {
                     sequence: sequence.clone(),
                     peptidoform: peptidoform.clone(),
                     charge,
-                    run_file_name: run_file_name.clone(),
+                    run_file_name: entry_run_file_name,
                     sample_accession: entry.sample_accession,
                     protein_accessions: protein_accessions.clone(),
                     anchor_protein: anchor_protein.clone(),
                     unique,
                     is_decoy,
                     pg_global_qvalue,
-                    label: entry.label,
+                    label: entry_label,
                     intensity: entry.intensity,
                 });
             }
@@ -507,6 +529,50 @@ mod tests {
         Ok(())
     }
 
+    #[test]
+    fn routes_lfq_run_labels_to_distinct_runs() -> Result<(), Box<dyn std::error::Error>> {
+        // LFQ: the intensity labels are run files (the row's anchor run appears
+        // among them), so each intensity must map to its own run/sample instead of
+        // collapsing onto the anchor `run_file_name`.
+        let intensities = lfq_run_intensities()?;
+        let proteins = new_qpx_proteins()?;
+        let schema = Arc::new(Schema::new(vec![
+            Field::new("sequence", DataType::Utf8, false),
+            Field::new("peptidoform", DataType::Utf8, false),
+            Field::new("charge", DataType::Int16, false),
+            Field::new("run_file_name", DataType::Utf8, false),
+            Field::new("anchor_protein", DataType::Utf8, false),
+            Field::new("unique", DataType::Boolean, false),
+            Field::new("is_decoy", DataType::Boolean, false),
+            Field::new("intensities", intensities.data_type().clone(), true),
+            Field::new("pg_accessions", proteins.data_type().clone(), true),
+        ]));
+        let batch = arrow::record_batch::RecordBatch::try_new(
+            schema,
+            vec![
+                Arc::new(StringArray::from(vec!["PEPTIDEK"])) as ArrayRef,
+                Arc::new(StringArray::from(vec!["PEPTIDEK"])) as ArrayRef,
+                Arc::new(Int16Array::from(vec![2])) as ArrayRef,
+                Arc::new(StringArray::from(vec!["Run_A_01.mzML"])) as ArrayRef,
+                Arc::new(StringArray::from(vec!["sp|P12345|PROT_HUMAN"])) as ArrayRef,
+                Arc::new(BooleanArray::from(vec![true])) as ArrayRef,
+                Arc::new(BooleanArray::from(vec![false])) as ArrayRef,
+                intensities,
+                proteins,
+            ],
+        )?;
+
+        let records = flatten_qpx_batch(&batch)?;
+
+        assert_eq!(records.len(), 3);
+        let runs: Vec<_> = records.iter().map(|r| r.run_file_name.as_str()).collect();
+        assert_eq!(runs, ["Run_A_01.mzML", "Run_B_01.mzML", "Run_B_02.mzML"]);
+        // LFQ has no reporter channel; run-file labels must not leak into `label`.
+        assert!(records.iter().all(|record| record.label.is_none()));
+        assert_eq!(records[1].intensity, 20.0);
+        Ok(())
+    }
+
     #[test]
     fn reads_integer_unique_and_decoy_flags() -> Result<(), Box<dyn std::error::Error>> {
         let intensities = new_qpx_intensities()?;
@@ -597,6 +663,33 @@ mod tests {
         Ok(Arc::new(builder.finish()) as ArrayRef)
     }
 
+    fn lfq_run_intensities() -> Result<ArrayRef, Box<dyn std::error::Error>> {
+        // Labels are run files; the anchor run (`Run_A_01.mzML`) appears among them.
+        let fields = Fields::from(vec![
+            Field::new("label", DataType::Utf8, false),
+            Field::new("intensity", DataType::Float32, false),
+        ]);
+        let struct_builder = StructBuilder::new(
+            fields,
+            vec![
+                Box::new(StringBuilder::new()),
+                Box::new(Float32Builder::new()),
+            ],
+        );
+        let mut builder = ListBuilder::new(struct_builder);
+        for (label, value) in [
+            ("Run_A_01.mzML", 10.0f32),
+            ("Run_B_01.mzML", 20.0),
+            ("Run_B_02.mzML", 30.0),
+        ] {
+            append_string_field(builder.values(), 0, label)?;
+            append_f32_field(builder.values(), 1, value)?;
+            builder.values().append(true);
+        }
+        builder.append(true);
+        Ok(Arc::new(builder.finish()) as ArrayRef)
+    }
+
     fn new_qpx_proteins() -> Result<ArrayRef, Box<dyn std::error::Error>> {
         let fields = Fields::from(vec![
             Field::new("accession", DataType::Utf8, false),

rust/crates/mokume-io/src/sdrf.rs

@@ -305,10 +305,20 @@ fn push_key(keys: &mut Vec<String>, value: &str) {
 }
 
 pub fn normalize_file_key(value: &str) -> String {
-    basename(value)
+    let key = basename(value)
         .trim()
         .trim_start_matches("file://")
-        .to_ascii_lowercase()
+        .to_ascii_lowercase();
+    // Strip a trailing run extension so the SDRF `comment[data file]` key matches
+    // the QPX `run_file_name`, which is stored extension-less. Without this the
+    // per-feature SDRF lookup misses and the Condition column falls back to the
+    // run filename. Mirrors `mokume-pipeline` de.rs::strip_run_extension.
+    for ext in [".raw", ".mzml", ".d", ".wiff"] {
+        if let Some(stem) = key.strip_suffix(ext) {
+            return stem.to_owned();
+        }
+    }
+    key
 }
 
 pub fn normalize_label_key(value: &str) -> String {
@@ -452,4 +462,43 @@ mod tests {
         assert_eq!(other.pooled_sample.as_deref(), Some("not pooled"));
         Ok(())
     }
+
+    #[test]
+    fn normalize_file_key_strips_run_extensions() {
+        // The QPX `run_file_name` is stored extension-less while the SDRF
+        // `comment[data file]` carries a run extension; both must normalize to
+        // the same key so the per-feature lookup matches.
+        for name in [
+            "Run_Condition_A_01.raw",
+            "Run_Condition_A_01.mzML",
+            "Run_Condition_A_01.d",
+            "Run_Condition_A_01.wiff",
+            "/data/Run_Condition_A_01.raw",
+            "Run_Condition_A_01",
+        ] {
+            assert_eq!(
+                normalize_file_key(name),
+                "run_condition_a_01",
+                "key for {name}"
+            );
+        }
+    }
+
+    #[test]
+    fn lookup_matches_extensionless_run_key() -> Result<()> {
+        // Regression: SDRF data files carry `.raw`, but the QPX run key is
+        // extension-less. Without extension stripping the lookup misses and the
+        // Condition column falls back to the run filename.
+        let input = concat!(
+            "source name\tcomment[data file]\tfactor value[group]\n",
+            "A1\tRun_Condition_A_01.raw\tA\n",
+            "B1\tRun_Condition_B_01.raw\tB\n",
+        );
+        let table = SdrfTable::from_reader(input.as_bytes())?;
+        let a = table
+            .lookup("Run_Condition_A_01", None)
+            .ok_or_else(|| invalid_input("missing extension-less lookup"))?;
+        assert_eq!(a.condition.as_deref(), Some("A"));
+        Ok(())
+    }
 }

rust/crates/mokume-pipeline/src/lib.rs

@@ -4642,6 +4642,16 @@ pub fn run_features_to_proteins(config: &FeatureToProteinsConfig) -> Result<()>
 
     validate_features_to_proteins(config)?;
     validate_implemented_subset(config)?;
+    if let Some(method) = ignored_directlfq_sample_normalization(config) {
+        warn!(
+            sample_normalization = method,
+            quant_method = %config.quantification,
+            "`--sample-normalization {method}` has no effect on the DirectLFQ route \
+             (used by `directlfq` and, unless `maxlfq.force_builtin` is set, `maxlfq`): \
+             DirectLFQ performs its own internal sample normalization, so mokume does \
+             not apply another. The protein matrix is unchanged by this option."
+        );
+    }
     configure_thread_pool(config.runtime.threads.or(config.directlfq.cores));
 
     let sdrf = config
@@ -5283,6 +5293,34 @@ fn is_cell_based_linear_quant(method: QuantMethod) -> bool {
 /// they center the summed-canonical-peptide matrix in log2 space, NOT the raw
 /// per-feature intensities, so they belong on the cell path rather than the
 /// ingest-time factor path.
+/// Returns the requested `--sample-normalization` method name when it is a
+/// dataset-level method that will be silently ignored because the effective
+/// quantification route is DirectLFQ (which performs its own internal sample
+/// normalization; see `Aggregation::applies_dataset_normalization`). Used to warn
+/// the user rather than dropping the option silently.
+fn ignored_directlfq_sample_normalization(config: &FeatureToProteinsConfig) -> Option<&str> {
+    let routes_to_directlfq = config.quantification == QuantMethod::DirectLfq
+        || (config.quantification == QuantMethod::MaxLfq && !config.maxlfq.force_builtin);
+    if !routes_to_directlfq {
+        return None;
+    }
+    match parse_sample_normalization_method(&config.normalization.sample_method)
+        .ok()
+        .flatten()
+    {
+        Some(
+            SampleNormalizationMethod::Quantile
+            | SampleNormalizationMethod::Rlr
+            | SampleNormalizationMethod::Loess
+            | SampleNormalizationMethod::Hierarchical
+            | SampleNormalizationMethod::MedianCenter
+            | SampleNormalizationMethod::MeanCenter
+            | SampleNormalizationMethod::Tmm,
+        ) => Some(config.normalization.sample_method.as_str()),
+        _ => None,
+    }
+}
+
 fn dataset_sample_normalization_method(
     config: &FeatureToProteinsConfig,
 ) -> Result<Option<SampleNormalizationMethod>> {