Exclude sample groups from quantification when MapAlignerIdentification doesn't finds shared feature maps

CHANGELOG.md

@@ -13,6 +13,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### `Fixed`
 
 - Fixed `UNTAR`, `UNZIP` and `EPICORE` version `eval` crashing on parenthesis-containing tool banners [#455](https://github.com/nf-core/mhcquant/pull/455)
+- Fixed `FeatureFinderIdentification` (OpenMS 3.5.0) losing all peptide identifications on FAIMS data, which crashed `SUMMARIZE_RESULTS`, by disabling `faims:merge_features` [#474](https://github.com/nf-core/mhcquant/pull/474)
+- Fixed `MapAlignerIdentification` aborting with `no data points for 'linear' model` when a run shares no RT landmarks with its group; the whole sample/condition group is now excluded from quantification, exported with identifications only, and a warning is emitted [#474](https://github.com/nf-core/mhcquant/pull/474)
 
 ### `Changed`
 

conf/modules.config

@@ -342,6 +342,10 @@ process {
     }
 
     withName: 'OPENMS_MAPALIGNERIDENTIFICATION' {
+        // MapAlignerIdentification exits with 8 (OpenMS UNKNOWN_ERROR, "no data points for 'linear' model") when a run
+        // shares no RT landmarks with the other runs of its group. The code is generic, but no other known failure of
+        // this tool maps to it; ignore it and let MAP_ALIGNMENT exclude such groups from quantification with a warning.
+        errorStrategy = { task.exitStatus == 8 ? 'ignore' : task.exitStatus in ((130..145) + 104 + (175..177)) ? 'retry' : 'finish' }
         ext.args    = [
             "-model:type linear",
             "-algorithm:max_rt_shift ${params.max_rt_alignment_shift}"
@@ -442,7 +446,9 @@ process {
                 "-extract:rt_window ${params.quantification_rt_window}",
                 "-detect:mapping_tolerance ${params.quantification_mapping_tolerance}",
                 "-detect:peak_width ${params.quantification_peak_width}",
-                "-detect:min_peak_width ${params.quantification_min_peak_width}"
+                "-detect:min_peak_width ${params.quantification_min_peak_width}",
+                // OpenMS 3.5.0 drops all peptide identifications from features when merging FAIMS CV groups
+                "-faims:merge_features false"
             ].join(' ').trim()
         publishDir  = [
             path: {"${params.outdir}/intermediate_results/features"},

docs/output.md

@@ -78,7 +78,7 @@ This folder contains the intermediate results from various steps of the MHCquant
 <summary>Output files</summary>
 
 - `intermediate_results/`
-  - `alignment`: Contains the `trafoXML` files of each run that document the retention time shift after alignment in quantification mode.
+  - `alignment`: Contains the `trafoXML` files of each run that document the retention time shift after alignment in quantification mode. Sample/condition groups for which `MapAlignerIdentification` fails because a run shares no RT landmarks with the others have no `trafoXML` files; the sample is excluded from quantification and exported with identifications only, and a warning is printed.
   - `comet`: Contains pin files generated by comet after database search
   - `features`: Holds information of quantified features in `featureXML` files as a result of the [FeatureFinderIdentification](https://openms.de/doxygen/release/3.0.0/html/TOPP_FeatureFinderIdentification.html) in the quantification mode.
   - `ion_annotations`(if `--annotate_ions` is specified)

docs/usage.md

@@ -179,6 +179,8 @@ nextflow run nf-core/mhcquant \
 
 When quantification is enabled, the pipeline performs retention time alignment and feature processing as detailed in the README documentation.
 
+If the retention times of a sample/condition group cannot be aligned because one of its runs shares no peptide identifications with the others within `--max_rt_alignment_shift`, the pipeline prints a warning, excludes the whole group from quantification and exports it with identifications only. Check the alignment log of that group and consider splitting runs acquired with different LC gradients or instruments into separate conditions.
+
 The quantification workflow produces a ConsensusXML file containing integrated peak areas for identified peptides across all samples.
 
 ## Spectrum Library Generation

modules/local/openms/mapaligneridentification/main.nf

@@ -29,9 +29,9 @@ process OPENMS_MAPALIGNERIDENTIFICATION {
     """
 
     stub:
+    def out_names = idxmls.collect { it.baseName.replace('_fdr_filtered','')+'.trafoXML' }.join(' ')
 
     """
-    touch test1.consensusXML
-    touch test2.consensusXML
+    touch ${out_names}
     """
 }

subworkflows/local/map_alignment/main.nf

@@ -17,33 +17,44 @@ workflow MAP_ALIGNMENT {
         // Compute group-wise alignment rt transformation
         OPENMS_MAPALIGNERIDENTIFICATION( ch_runs_to_be_aligned )
 
-        // Join run specific trafoXMLs with meta information
-        merge_meta_map
-            .flatMap { group_meta, metas -> metas }
-            .map { meta -> [[spectra:meta.spectra], meta]}
-            .join( OPENMS_MAPALIGNERIDENTIFICATION.out.trafoxml
-                    .flatMap { group_meta, trafoxmls -> [trafoxmls].flatten().collect { trafoxml -> [[spectra: trafoxml.baseName], trafoxml] } })
-            .map { spectra, meta, trafoxml -> [meta, trafoxml] }
+        // Run-specific trafoXMLs: [[spectra], trafoxml]
+        OPENMS_MAPALIGNERIDENTIFICATION.out.trafoxml
+            .flatMap { group_meta, trafoxmls -> [trafoxmls].flatten().collect { trafoxml -> [[spectra: trafoxml.baseName], trafoxml] } }
             .set { ch_trafos }
 
-        // Align mzML files using trafoXMLs
-        ch_trafos_mzmls = ch_mzml.join(ch_trafos)
-        OPENMS_MAPRTTRANSFORMERMZML(ch_trafos_mzmls)
-
-        // Align idXMLfiles using trafoXMLs
+        // Runs with their meta, idXML and mzML: [[spectra], meta, idxml, mzml]
         ch_runs_to_be_aligned
             .flatMap { group_meta, idxmls -> [idxmls].flatten().collect { idxml -> [[spectra: idxml.baseName.replace("_fdr_filtered","")], idxml] } }
             .join( merge_meta_map
                     .flatMap { group_meta, metas -> metas }
                     .map { meta -> [[spectra:meta.spectra], meta]} )
-            .map { group_meta, idxml, meta -> [meta, idxml] }
-            .join( ch_trafos )
-            .set { ch_trafos_idxml }
-
-        OPENMS_MAPRTTRANSFORMERIDXML(ch_trafos_idxml)
+            .join( ch_mzml.map { meta, mzml -> [[spectra: meta.spectra], mzml] } )
+            // Groups whose alignment failed (MapAlignerIdentification exit 8 is ignored) have no trafoXMLs
+            .join( ch_trafos, remainder: true )
+            .map { spectra, idxml, meta, mzml, trafoxml -> [meta, idxml, mzml, trafoxml] }
+            .branch { meta, idxml, mzml, trafoxml ->
+                aligned: trafoxml
+                failed:  true
+            }
+            .set { ch_runs_by_alignment }
+
+        // Groups that could not be aligned are excluded from quantification: [[id: sample_condition]]
+        ch_runs_by_alignment.failed
+            .map { meta, idxml, mzml, trafoxml -> [id: "${meta.sample}_${meta.condition}"] }
+            .unique()
+            .set { ch_failed_groups }
+
+        ch_failed_groups.subscribe { group_meta ->
+            log.warn "RT alignment of sample '${group_meta.id}' failed: at least one run has no peptide IDs in common with the other runs within --max_rt_alignment_shift. Skipping quantification for this sample, only identifications are reported."
+        }
+
+        // Align mzML and idXML files using trafoXMLs
+        OPENMS_MAPRTTRANSFORMERMZML(ch_runs_by_alignment.aligned.map { meta, idxml, mzml, trafoxml -> [meta, mzml, trafoxml] })
+        OPENMS_MAPRTTRANSFORMERIDXML(ch_runs_by_alignment.aligned.map { meta, idxml, mzml, trafoxml -> [meta, idxml, trafoxml] })
 
     emit:
         aligned_idxml = OPENMS_MAPRTTRANSFORMERIDXML.out.aligned
-        aligned_mzml = OPENMS_MAPRTTRANSFORMERMZML.out.aligned
-        trafoxml = OPENMS_MAPALIGNERIDENTIFICATION.out.trafoxml
+        aligned_mzml  = OPENMS_MAPRTTRANSFORMERMZML.out.aligned
+        trafoxml      = OPENMS_MAPALIGNERIDENTIFICATION.out.trafoxml
+        failed_groups = ch_failed_groups
 }

subworkflows/local/quant/main.nf

@@ -85,6 +85,7 @@ workflow QUANT {
         OPENMS_MZTABEXPORTER(PROCESS_FEATURE.out.consensusxml)
 
     emit:
-        consensusxml = PROCESS_FEATURE.out.consensusxml
-        trafoxml = MAP_ALIGNMENT.out.trafoxml
+        consensusxml  = PROCESS_FEATURE.out.consensusxml
+        trafoxml      = MAP_ALIGNMENT.out.trafoxml
+        failed_groups = MAP_ALIGNMENT.out.failed_groups
 }

workflows/mhcquant.nf

@@ -176,7 +176,9 @@ workflow MHCQUANT {
     //
     if (params.quantify) {
         QUANT(merge_meta_map, RESCORE.out.rescored_runs, RESCORE.out.fdr_filtered, ch_clean_mzml_file)
-        ch_output = QUANT.out.consensusxml.mix(RESCORE.out.fdr_filtered_empty)
+        // Samples whose RT alignment failed are exported with identifications only, like empty samples
+        ch_alignment_failed = RESCORE.out.fdr_filtered.join(QUANT.out.failed_groups)
+        ch_output = QUANT.out.consensusxml.mix(RESCORE.out.fdr_filtered_empty, ch_alignment_failed)
     } else {
         ch_output = RESCORE.out.fdr_filtered.mix(RESCORE.out.fdr_filtered_empty)
     }