Add Sage converter for label-free quantification

DESCRIPTION

@@ -33,6 +33,7 @@ Suggests:
     rmarkdown
 LinkingTo: Rcpp
 Collate: 
+    'clean_Sage.R'
     'clean_MZMine.R'
     'clean_ProteinProspector.R'
     'clean_Metamorpheus.R'
@@ -64,6 +65,7 @@ Collate:
     'converters_PhilosophertoMSstatsTMTFormat.R'
     'converters_ProgenesistoMSstatsFormat.R'
     'converters_ProteinProspectortoMSstatsTMTFormat.R'
+    'converters_SagetoMSstatsFormat.R'
     'converters_SkylinetoMSstatsFormat.R'
     'converters_SpectroMinetoMSstatsTMTFormat.R'
     'converters_SpectronauttoMSstatsFormat.R'

NAMESPACE

@@ -26,6 +26,7 @@ export(PDtoMSstatsTMTFormat)
 export(PhilosophertoMSstatsTMTFormat)
 export(ProgenesistoMSstatsFormat)
 export(ProteinProspectortoMSstatsTMTFormat)
+export(SagetoMSstatsFormat)
 export(SkylinetoMSstatsFormat)
 export(SpectroMinetoMSstatsTMTFormat)
 export(SpectronauttoMSstatsFormat)

R/MSstatsConvert_core_functions.R

@@ -75,6 +75,10 @@ setClass("MSstatsProteinProspectorFiles", contains = "MSstatsInputFiles")
 #' @rdname MSstatsInputFiles
 #' @keywords internal
 setClass("MSstatsMZMineFiles", contains = "MSstatsInputFiles")
+#' MSstatsSageFiles: class for Sage files.
+#' @rdname MSstatsInputFiles
+#' @keywords internal
+setClass("MSstatsSageFiles", contains = "MSstatsInputFiles")
 
 
 #' Get one of files contained in an instance of `MSstatsInputFiles` class.
@@ -304,6 +308,13 @@ setMethod("MSstatsClean", signature = "MSstatsProteinProspectorFiles",
 #' @return data.table
 setMethod("MSstatsClean", signature = "MSstatsMZMineFiles",
           .cleanRawMZMine)
+#' Clean Sage files
+#' @include clean_Sage.R
+#' @rdname MSstatsClean
+#' @inheritParams .cleanRawSage
+#' @return data.table
+setMethod("MSstatsClean", signature = "MSstatsSageFiles",
+          .cleanRawSage)
 
 
 #' Preprocess outputs from MS signal processing tools for analysis with MSstats

R/clean_Sage.R

@@ -0,0 +1,70 @@
+#' Clean raw Sage LFQ output
+#'
+#' Operates on Sage's `lfq.tsv` report (produced when `quant.lfq: true`). This
+#' is a wide-format table with fixed columns `peptide`, `charge`, `proteins`,
+#' `q_value`, `score`, `spectral_angle`, followed by one intensity column per
+#' input mzML, each headed by the run's file name. Intensity columns are
+#' identified as every column that is not one of the six fixed columns (matching
+#' is by name, not by file extension, so renamed files are handled). The table is
+#' melted to long format, columns are renamed to the MSstats standard, and zero
+#' intensities -- which Sage writes for precursors it did not quantify in a run --
+#' are converted to `NA`.
+#'
+#' @param msstats_object an object of class `MSstatsSageFiles`.
+#' @return data.table
+#' @keywords internal
+.cleanRawSage = function(msstats_object) {
+    Intensity = NULL
+
+    sage_input = getInputFile(msstats_object, "input")
+    sage_input = data.table::as.data.table(sage_input)
+
+    fixed_columns = c("peptide", "charge", "proteins", "q_value",
+                      "score", "spectral_angle")
+    required_columns = c("peptide", "charge", "proteins", "q_value")
+    missing_columns = setdiff(required_columns, colnames(sage_input))
+    if (length(missing_columns) > 0) {
+        msg = paste("The following required columns are missing from the Sage",
+                    "input:", paste(missing_columns, sep = ", ", collapse = ", "))
+        getOption("MSstatsLog")("ERROR", msg)
+        stop(msg)
+    }
+
+    intensity_columns = setdiff(colnames(sage_input), fixed_columns)
+    if (length(intensity_columns) == 0) {
+        msg = paste("No intensity columns found in the Sage input. Expected at",
+                    "least one per-run intensity column in addition to the fixed",
+                    "columns:", paste(fixed_columns, sep = ", ", collapse = ", "))
+        getOption("MSstatsLog")("ERROR", msg)
+        stop(msg)
+    }
+
+    id_columns = intersect(c("proteins", "peptide", "charge", "q_value"),
+                           colnames(sage_input))
+    sage_input = sage_input[, c(id_columns, intensity_columns), with = FALSE]
+
+    long = data.table::melt(sage_input,
+                            id.vars = id_columns,
+                            measure.vars = intensity_columns,
+                            variable.name = "Run",
+                            value.name = "Intensity",
+                            variable.factor = FALSE)
+
+    data.table::setnames(long,
+                         c("proteins", "peptide", "charge"),
+                         c("ProteinName", "PeptideSequence", "PrecursorCharge"))
+
+    long[, Intensity := as.numeric(Intensity)]
+    long[Intensity == 0, Intensity := NA_real_]
+
+    if (all(long$PrecursorCharge == -1)) {
+        msg = paste("** All PrecursorCharge values are -1: Sage combined charge",
+                    "states (combine_charge_states = true), so the feature key is",
+                    "effectively the peptide sequence alone.")
+        getOption("MSstatsLog")("INFO", msg)
+        getOption("MSstatsMsg")("INFO", msg)
+    }
+
+    .logSuccess("Sage", "clean")
+    long
+}

R/converters_SagetoMSstatsFormat.R

@@ -0,0 +1,176 @@
+#' Import Sage LFQ files
+#'
+#' Converts the label-free quantification report (`lfq.tsv`) produced by the Sage
+#' search engine into a `data.frame` in the format required by MSstats. The input
+#' is wide (one row per precursor, one intensity column per run); it is reshaped
+#' to long format, q-value filtered, and returned ready for `dataProcess`.
+#'
+#' @inheritParams .sharedParametersAmongConverters
+#' @param input Sage `lfq.tsv` report, as a `data.frame`/`data.table` or a path.
+#'   Wide format with fixed columns `peptide`, `charge`, `proteins`, `q_value`,
+#'   `score`, `spectral_angle`, followed by one intensity column per input mzML,
+#'   each headed by the run's file name.
+#' @param annotation `data.frame` with `Run`, `Condition` and `BioReplicate`
+#'   columns (a `Fraction` column may also be supplied). This argument is
+#'   **required**: Sage's `lfq.tsv` carries no experimental design, so condition
+#'   and replicate information must be provided separately.
+#' @param qvalue_cutoff Cutoff for the `q_value` column. Default is 0.01.
+#' @param filter_with_Qvalue TRUE (default) replaces intensities whose `q_value`
+#'   exceeds `qvalue_cutoff` with `NA` (treated as censored missing downstream);
+#'   FALSE performs no q-value filtering. See the "FDR filtering" section for
+#'   why this matters.
+#'
+#' @return `data.frame` in the MSstats required format.
+#'
+#' @section Input file:
+#' Use `lfq.tsv`, which Sage writes only when `quant.lfq: true` is set in the
+#' search configuration. Do **not** use `results.sage.tsv`: its `ms2_intensity`
+#' column is the summed intensity of matched b/y fragment ions -- a PSM score
+#' feature -- and is not a quantitative measure of precursor abundance.
+#'
+#' @section Acquisition mode:
+#' Sage's LFQ quantification is MS1-based regardless of how the data was
+#' acquired; the `wide_window` setting changes only how spectra are searched,
+#' not how they are quantified. For DDA, MS1 precursor areas are the standard
+#' quantification approach, so this converter is appropriate. For DIA, MS1
+#' quantification is not recommended: the MS1 signal carries substantial
+#' interference from co-eluting precursors, which is why DIA workflows normally
+#' quantify at the fragment level. Sage does write per-fragment MS2 intensities
+#' to `matched_fragments.sage.tsv` when `annotate_matches` is enabled, but those
+#' are matched peaks from individual spectra and Sage does not trace them across
+#' runs, so there is no fragment-level quantification available to convert. In
+#' short, this converter is recommended for DDA experiments; DIA users should
+#' expect MS1-based quantification with the limitations above.
+#'
+#' @section FDR filtering:
+#' `lfq.tsv` is not FDR-filtered. Sage writes every quantified peptide and
+#' charge row regardless of its q-value, leaving the filtering choice to
+#' downstream tools. The `lfq_settings.peptide_q_value` setting in the Sage
+#' configuration is an internal threshold used when building the LFQ traces;
+#' it does not filter what is written to the file, so a report produced with
+#' `peptide_q_value` set to `0.01` will still contain rows well above `0.01`.
+#' This converter applies the filter: `filter_with_Qvalue` defaults to `TRUE`
+#' and `qvalue_cutoff` defaults to `0.01`, so a default call returns
+#' FDR-filtered output. The filter is load-bearing on real data -- in two
+#' files from the issue author, 61 percent of rows in an eight-run file and
+#' 23 percent in a single-run file were above `0.01`. Set
+#' `filter_with_Qvalue = FALSE` to return unfiltered data.
+#'
+#' @section Charge states:
+#' Sage's `combine_charge_states` option (default `true`) sums charge states and
+#' writes `charge` as `-1` for every row, so `PrecursorCharge` will be `-1`
+#' throughout and the feature key reduces to the peptide sequence. Setting
+#' `combine_charge_states: false` reports real precursor charges, but is
+#' considerably slower across multiple files.
+#'
+#' @section Run name matching:
+#' MSstatsConvert standardizes column names by removing spaces and dots (`.`)
+#' while preserving hyphens and underscores. The melted `Run` values (the
+#' intensity column headers) and the annotation `Run` values are both passed
+#' through this same standardization before merging, so they match automatically.
+#' For example, a run named `B.naive_01steady-state.mzML.gz` in the annotation
+#' becomes `Bnaive_01steady-statemzMLgz`; supply the raw file name in the
+#' annotation and the merge resolves it. Note that the `Run` values in the
+#' returned table are the standardized form.
+#'
+#' @section Shared peptides:
+#' Sage pre-joins shared proteins into a single semicolon-delimited `proteins`
+#' value (e.g. `sp|A|X;sp|B|Y`). Because that is one `ProteinName` string rather
+#' than several, MSstats' shared-peptide removal sees a single protein per
+#' peptide and does not treat these rows as shared. Consequently
+#' `useUniquePeptide = TRUE` has no effect on peptides that Sage reports against a
+#' shared (semicolon-joined) protein group.
+#'
+#' @export
+#'
+#' @examples
+#' sage_lfq = system.file("tinytest/raw_data/Sage/lfq.tsv",
+#'                        package = "MSstatsConvert")
+#' annot_path = system.file("tinytest/raw_data/Sage/annotation.csv",
+#'                          package = "MSstatsConvert")
+#' if (nzchar(sage_lfq) && nzchar(annot_path)) {
+#'     sage_input = data.table::fread(sage_lfq)
+#'     annotation = read.csv(annot_path)
+#'     sage_imported = SagetoMSstatsFormat(sage_input, annotation,
+#'                                         use_log_file = FALSE)
+#'     head(sage_imported)
+#' }
+#'
+SagetoMSstatsFormat = function(
+        input, annotation, useUniquePeptide = TRUE,
+        removeFewMeasurements = TRUE, removeProtein_with1Peptide = FALSE,
+        qvalue_cutoff = 0.01, filter_with_Qvalue = TRUE,
+        use_log_file = TRUE, append = FALSE, verbose = TRUE, log_file_path = NULL,
+        ...
+) {
+    IsotopeLabelType = NULL
+
+    validation_config = list(
+        input = input,
+        annotation = annotation,
+        filter_with_Qvalue = filter_with_Qvalue,
+        qvalue_cutoff = qvalue_cutoff,
+        useUniquePeptide = useUniquePeptide,
+        removeFewMeasurements = removeFewMeasurements,
+        removeProtein_with1Feature = removeProtein_with1Peptide,
+        use_log_file = use_log_file,
+        append = append,
+        verbose = verbose,
+        log_file_path = log_file_path
+    )
+    .validateMSstatsConverterParameters(validation_config)
+
+    MSstatsConvert::MSstatsLogsSettings(use_log_file, append, verbose,
+                                        log_file_path)
+
+    input = MSstatsConvert::MSstatsImport(list(input = input),
+                                          "MSstats", "Sage", ...)
+    input = MSstatsConvert::MSstatsClean(input)
+
+    if (inherits(annotation, "data.frame") &&
+        is.element("IsotopeLabelType", colnames(annotation))) {
+        annotation = data.table::as.data.table(annotation)
+        annotation[, IsotopeLabelType := NULL]
+        msg = paste("** An IsotopeLabelType column was found in the annotation",
+                    "and has been dropped. Sage LFQ is label-free;",
+                    "IsotopeLabelType is set to 'L' for all rows.")
+        getOption("MSstatsLog")("INFO", msg)
+        getOption("MSstatsMsg")("INFO", msg)
+    }
+    annotation = MSstatsConvert::MSstatsMakeAnnotation(input, annotation)
+
+    qval_filter = list(score_column = "q_value",
+                       score_threshold = qvalue_cutoff,
+                       direction = "smaller",
+                       behavior = "fill",
+                       handle_na = "keep",
+                       fill_value = NA_real_,
+                       filter = filter_with_Qvalue,
+                       drop_column = TRUE)
+
+    feature_columns = c("PeptideSequence", "PrecursorCharge")
+    input = MSstatsConvert::MSstatsPreprocess(
+        input,
+        annotation,
+        feature_columns,
+        remove_shared_peptides = useUniquePeptide,
+        remove_single_feature_proteins = removeProtein_with1Peptide,
+        feature_cleaning = list(
+            remove_features_with_few_measurements = removeFewMeasurements,
+            summarize_multiple_psms = function(x, na.rm = TRUE) {
+                if (all(is.na(x))) NA_real_ else max(x, na.rm = na.rm)
+            }),
+        score_filtering = list(qvalue = qval_filter),
+        columns_to_fill = list("FragmentIon" = NA,
+                               "ProductCharge" = NA,
+                               "IsotopeLabelType" = "L"))
+    input = MSstatsConvert::MSstatsBalancedDesign(input, feature_columns,
+                                                  remove_few = removeFewMeasurements)
+
+    msg_final = paste("** Finished preprocessing. The dataset is ready",
+                      "to be processed by the dataProcess function.")
+    getOption("MSstatsLog")("INFO", msg_final)
+    getOption("MSstatsMsg")("INFO", msg_final)
+    getOption("MSstatsLog")("INFO", "\n")
+    input
+}

inst/tinytest/raw_data/Sage/annotation.csv

@@ -0,0 +1,9 @@
+Run,Condition,BioReplicate,IsotopeLabelType
+B.naive_01steady-state.mzML.gz,B.naive,B.naive_1,L
+B.naive_02steady-state.mzML.gz,B.naive,B.naive_2,L
+B.naive_03steady-state.mzML.gz,B.naive,B.naive_3,L
+B.naive_04steady-state.mzML.gz,B.naive,B.naive_4,L
+T4.naive_01steady-state.mzML.gz,T4.naive,T4.naive_1,L
+T4.naive_02steady-state.mzML.gz,T4.naive,T4.naive_2,L
+T4.naive_03steady-state.mzML.gz,T4.naive,T4.naive_3,L
+T4.naive_04steady-state.mzML.gz,T4.naive,T4.naive_4,L

inst/tinytest/raw_data/Sage/annotation_charge_resolved.csv

@@ -0,0 +1,2 @@
+Run,Condition,BioReplicate
+B.naive_01steady-state.mzML.gz,B.naive,B.naive_1