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,162 @@
+#' 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 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

inst/tinytest/raw_data/Sage/lfq.tsv

@@ -0,0 +1,27 @@
+peptide	charge	proteins	q_value	score	spectral_angle	B.naive_01steady-state.mzML.gz	B.naive_02steady-state.mzML.gz	B.naive_03steady-state.mzML.gz	B.naive_04steady-state.mzML.gz	T4.naive_01steady-state.mzML.gz	T4.naive_02steady-state.mzML.gz	T4.naive_03steady-state.mzML.gz	T4.naive_04steady-state.mzML.gz
+SEGFDTYR	-1	sp|P12004|PCNA_HUMAN	0.0017778707	0.8596611544837863	0.9529590363688335	77817957.28898548	54730854.00169097	44278836.27082016	79736868.87154807	0.0	166636508.30663136	90507.02519561221	0.0
+VLYLGAASGTTVSHVSDIVGPDGLVYAVEFSHR	-1	sp|P22087|FBRL_HUMAN	0.0059178164	0.7689931383115772	0.9400131025469994	8167823.2774881115	6714142.398753918	5391514.995625099	11163861.941858275	220498241.90016308	1369872189.8564844	3424035.1178108207	823966765.049154
+QINEDNERWETNR	-1	sp|Q92620|PRP16_HUMAN	0.01138836	0.7042599144076336	0.8916808811404437	2338583.2452762416	11908007.34856142	7027569.145615996	9766125.904102754	38625778.3048101	38440497.17166741	18126309.721357092	43445868.098005034
+YLLQPVLAPLHR	-1	sp|Q9C0B7|TNG6_HUMAN	0.0063552475	0.7646499537315222	0.9209709867931564	170933.8067667412	313176.7755232925	564752.8904693746	230549.4767265296	23015919.078579333	0.0	7619412.940611036	21620470.4114524
+KIEDLIK	-1	sp|O75935|DCTN3_HUMAN	0.0008149959	0.9260075351426054	0.9790878726135238	119660.06945591315	889238.6806568115	5469118.4792797025	849206.4551180338	157658514.34854704	276963873.44115895	59520615.210684635	289946595.26979446
+AQGEPVAGHESPK	-1	sp|O94979|SC31A_HUMAN	0.05990019	0.41293955380751407	0.7497528020043859	586045.345304421	25312606.74640692	620830.1387602019	76141.46130667519	13626364.145430267	6143.119122577197	1792599.2700471443	41490673.92943229
+STAGDTHLGGEDFDNR	-1	sp|P11142|HSP7C_HUMAN;sp|P54652|HSP72_HUMAN	0.006671817	0.7610814290609372	0.9306927272260735	4277476762.8728995	5401129431.49641	2988697335.9601717	3879914287.010699	11661596334.61851	9799036.062493572	1377037237.292962	12075911760.157255
+DLLDQILM[+15.994915]LDPAK	-1	sp|Q13523|PRP4B_HUMAN	0.002228933	0.8362442785549401	0.9421304772494757	743508.1825100655	23031808.928539995	24198397.73073512	26660489.78436039	327182.4742234778	401353.99213329196	236647.34163009128	107380.94719313568
+IDSILEVVQTGR	-1	sp|Q9NZW5|PALS2_HUMAN	0.064401604	0.39847785284792436	0.8592752928588256	40024289.28258917	34939304.71416704	42476703.561981075	48743102.68411822	11593139.282803845	24797539.813975845	7476026.899138306	22065720.749139514
+ADNTWDPEIPVC[+57.021465]EK	-1	sp|P20023|CR2_HUMAN	0.0015804922	0.8651594046154879	0.9634082723710683	40624720.612872265	77054506.0736038	57327535.089999124	49428406.766203314	0.0	0.0	0.0	283069.43880510033
+AIQSLKK	-1	sp|O75964|ATP5L_HUMAN;sp|Q7Z4Y8|AT5L2_HUMAN	0.0008149959	0.9591453463988604	0.9906303720856269	0.0	7328662.742700506	0.0	47651.952119284484	0.0	0.0	18509.13971851515	0.0
+GFGFVTFDDHDPVDK	-1	sp|P22626|ROA2_HUMAN	0.028423676	0.5610882303897902	0.9524681927608495	9442272946.076956	745595963.1325712	7636709076.084447	8307305230.95171	6142772.814928512	7125065.134709419	12563652.814297214	17889814.35587386
+LTALER	-1	sp|Q5VZP5|STYL2_HUMAN;sp|Q8WUW1|BRK1_HUMAN	0.0031524044	0.8179238374286326	0.9468500231177203	137914663.99828738	66697028.02119711	24960682.03221422	26534231.722668294	1146742.627962545	13938566.088161986	4607497.752185434	6851775.428993263
+GTAVAIC[+57.021465]R	-1	sp|P16615|AT2A2_HUMAN;sp|Q93084|AT2A3_HUMAN	0.0013819312	0.8859071754424623	0.9604234154507731	5510079.084221685	220288626.77732605	170159172.58043033	100883408.81250122	7925901.40969537	15466045.49329941	18407395.445723712	41532517.26298703
+DLSLEEIQK	-1	sp|P16949|STMN1_HUMAN;sp|Q93045|STMN2_HUMAN	0.007042969	0.7560161124390581	0.9238835969893973	248555277.27070376	907025.1900574858	159649919.04960185	633606.0077070697	27517169.758237816	340975012.55704695	6202235.105255563	24661025.83516855
+SLTAEIDR	-1	sp|P54136|SYRC_HUMAN	0.0040523116	0.797296911999848	0.9272710349927167	167213584.75358346	35879524.214564264	13938517.723104605	31234697.89723105	29362277.405670784	240963953.93132198	132886029.7659108	278062917.3039626
+VGILTLSDYVLELGHPYLWVQK	-1	sp|Q13769|THOC5_HUMAN	0.024243616	0.5921044903672563	0.8397166894002543	0.0	0.0	0.0	0.0	67403257.38741247	123701495.03355601	16797087.397493932	75340695.68453664
+NPSAMAVESFMATAPFVQIGR	-1	sp|Q9BTV4|TMM43_HUMAN	0.0027505157	0.8264128592192239	0.9384238181947318	0.0	312188.15119036747	0.0	32542.996896974815	13209000.056127746	90622863.80383575	503115809.00002736	672286616.1006802
+AVDAALK	-1	sp|P01042|KNG1_HUMAN	0.0063552475	0.7649703077288139	0.925923633524725	2154744.172844004	14859919.302213334	11997683.193815425	11404039.256254723	2412446.902966986	20637635.162529614	13810458.152812112	13461619.338058997
+SC[+57.021465]STFEQWFNAPFAMTGEK	-1	sp|P51532|SMCA4_HUMAN	0.0039590127	0.7987676059851891	0.9406410329468484	118742.59420617236	87040.81437970039	54636.499999999985	54159.21173634676	27972753.755372617	764578.9053270838	16220841.498979999	25091889.82963039
+DMNHTNSYGIIR	-1	sp|Q6P2Q9|PRP8_HUMAN	0.0025925206	0.828113734753338	0.9507899034568713	34151412.76742061	33592535.707059756	8289100.97522942	25535702.62727023	378105665.2935424	20824276.437241968	231793023.67943525	382648534.6238318
+YRQFPQLTR	-1	sp|O95178|NDUB2_HUMAN	0.0012557555	0.8935027614581135	0.9631604323972088	1724858.2180483385	611073.1034039721	459408.461838078	585393.118838761	290079691.5456341	613227349.5816588	250345906.41322634	596961620.8246347
+HIM[+15.994915]GQNVADYM[+15.994915]R	-1	sp|P46777|RL5_HUMAN	0.0028177642	0.8246557875916052	0.9486897923731843	566052.2905316218	157034960.93516892	108829606.00649722	3901780.2195805656	21775041.622718714	10384437.815386051	6238520.984804116	17239062.113721438
+TFVNLTK	-1	sp|Q6ZN55|ZN574_HUMAN	0.007119441	0.7516215588107297	0.9092146179459012	111238.54414194662	142111.64229268857	0.0	0.0	42254615.26886397	39490841.18555253	51858573.53238733	80533227.82593004
+FGSPKPPVAVKPSSEEKPDK	-1	sp|O15117|FYB1_HUMAN	0.009162225	0.730180940894752	0.9045383725134525	19487077.996781565	1509622.1153875482	641658.3577688276	1007946.277595538	318465415.7884141	611465083.5754277	177489610.00745344	502316059.11888725
+DELLSYIK	-1	sp|Q96Q11|TRNT1_HUMAN	0.0013819312	0.8879309471956518	0.9700104282749398	4352138.113172804	0.0	185327.70449523636	383326.61100796773	87652167.32207315	415171884.8467694	105483562.54983357	74805899.42414689