diff --git a/DESCRIPTION b/DESCRIPTION index d409af11..1df9759b 100644 --- a/DESCRIPTION +++ b/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' diff --git a/NAMESPACE b/NAMESPACE index 94f287cf..c15e35df 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -26,6 +26,7 @@ export(PDtoMSstatsTMTFormat) export(PhilosophertoMSstatsTMTFormat) export(ProgenesistoMSstatsFormat) export(ProteinProspectortoMSstatsTMTFormat) +export(SagetoMSstatsFormat) export(SkylinetoMSstatsFormat) export(SpectroMinetoMSstatsTMTFormat) export(SpectronauttoMSstatsFormat) diff --git a/R/MSstatsConvert_core_functions.R b/R/MSstatsConvert_core_functions.R index 7e4c77ae..e084ab00 100644 --- a/R/MSstatsConvert_core_functions.R +++ b/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 diff --git a/R/clean_Sage.R b/R/clean_Sage.R new file mode 100644 index 00000000..d74e103c --- /dev/null +++ b/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 +} diff --git a/R/converters_SagetoMSstatsFormat.R b/R/converters_SagetoMSstatsFormat.R new file mode 100644 index 00000000..ddbcda25 --- /dev/null +++ b/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 +} diff --git a/inst/tinytest/raw_data/Sage/annotation.csv b/inst/tinytest/raw_data/Sage/annotation.csv new file mode 100644 index 00000000..c5c0cebb --- /dev/null +++ b/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 diff --git a/inst/tinytest/raw_data/Sage/annotation_charge_resolved.csv b/inst/tinytest/raw_data/Sage/annotation_charge_resolved.csv new file mode 100644 index 00000000..95076395 --- /dev/null +++ b/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 diff --git a/inst/tinytest/raw_data/Sage/lfq.tsv b/inst/tinytest/raw_data/Sage/lfq.tsv new file mode 100644 index 00000000..1095697f --- /dev/null +++ b/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 diff --git a/inst/tinytest/raw_data/Sage/lfq_charge_resolved.tsv b/inst/tinytest/raw_data/Sage/lfq_charge_resolved.tsv new file mode 100644 index 00000000..c0c37a4a --- /dev/null +++ b/inst/tinytest/raw_data/Sage/lfq_charge_resolved.tsv @@ -0,0 +1,17 @@ +peptide charge proteins q_value score spectral_angle B.naive_01steady-state.mzML.gz +TIPWLEDRVPQK 2 sp|O43707|ACTN4_HUMAN 0.0026025237 0.8062962329450262 0.9559143416649949 71339186.66976815 +TIPWLEDRVPQK 3 sp|O43707|ACTN4_HUMAN 0.0026025237 0.8066638416368576 0.9566809449114555 203432318.96055964 +GWQDVTATSAYKK 2 sp|P55327|TPD52_HUMAN 0.0010507355 0.924242775957422 0.9762487405050265 55586267.52244711 +GWQDVTATSAYKK 3 sp|P55327|TPD52_HUMAN 0.0010507355 0.9096214284739153 0.9710733284139133 157347364.0188433 +HNYGVVESFTVQR 3 sp|P01911|DRB1_HUMAN 0.0016844276 0.8445341635955704 0.9516889139060459 56154030.8900844 +HNYGVVESFTVQR 2 sp|P01911|DRB1_HUMAN 0.004562829 0.7228290743061249 0.9035824913746858 25695257.588136908 +AVSDWIDEQEK 2 sp|Q12906|ILF3_HUMAN 0.0010507355 0.9332034028718135 0.9772194550760043 256956990.39009976 +AVSDWIDEQEK 3 sp|Q12906|ILF3_HUMAN 0.0016383065 0.849620908242692 0.9471273920960341 1141479.7267733715 +GVSQTGTPVC[+57.021465]EEDGDAGLGIR 3 sp|Q12888|TP53B_HUMAN 0.00061576354 0.9627451501906364 0.9889407762648862 9601575.47511771 +GVSQTGTPVC[+57.021465]EEDGDAGLGIR 2 sp|Q12888|TP53B_HUMAN 0.008568041 0.6432843525811024 0.8651660339502047 13369710.457350656 +LRDHDDAAESLIEQTTALNK 3 sp|Q9NVK5|FGOP2_HUMAN 0.0019852011 0.8262793758164204 0.9540490534935053 12669608.02702529 +FC[+57.021465]AEIK 2 sp|P49903|SPS1_HUMAN 0.008841005 0.639872471608157 0.8717617294315736 132139217.61067836 +QSGESIDIITR 2 sp|Q16531|DDB1_HUMAN 0.0010507355 0.9027512508686987 0.9708218239508123 177966814.78880996 +DIDAFWLQR 2 sp|O75643|U520_HUMAN 0.0013855582 0.8782748894145357 0.9712187573469061 49539732.8735502 +ELC[+57.021465]HVC[+57.021465]C[+57.021465]MK 2 sp|O14672|ADA10_HUMAN 0.011006618 0.6106634885078682 0.8522182216802359 4339921.698240255 +VSQAAADLK 3 sp|A0A804HLA8|GBG5B_HUMAN;sp|P63218|GBG5_HUMAN 0.048834465 0.3398385671534383 0.7429746919519163 44554.79772631668 diff --git a/inst/tinytest/test_clean_Sage.R b/inst/tinytest/test_clean_Sage.R new file mode 100644 index 00000000..7d812fdc --- /dev/null +++ b/inst/tinytest/test_clean_Sage.R @@ -0,0 +1,71 @@ +# Test .cleanRawSage on the Sage lfq.tsv fixture --------------------------- +sage_path = system.file("tinytest/raw_data/Sage/lfq.tsv", + package = "MSstatsConvert") +if (!nzchar(sage_path)) { + exit_file("Sage fixtures not present in inst/tinytest/raw_data/Sage/") +} +sage_raw = data.table::fread(sage_path) + +fixed_cols = c("peptide", "charge", "proteins", "q_value", "score", + "spectral_angle") +intensity_cols = setdiff(colnames(sage_raw), fixed_cols) + +msstats_input = MSstatsConvert::MSstatsImport( + list(input = sage_raw), "MSstats", "Sage") +cleaned = MSstatsConvert:::.cleanRawSage(msstats_input) +cleaned = data.table::as.data.table(cleaned) + +# Renamed to canonical MSstats names; wide reshaped to long +expect_true(all(c("ProteinName", "PeptideSequence", "PrecursorCharge", + "q_value", "Run", "Intensity") %in% colnames(cleaned))) +expect_false("proteins" %in% colnames(cleaned)) +expect_false("peptide" %in% colnames(cleaned)) +expect_false("charge" %in% colnames(cleaned)) +# score / spectral_angle are dropped in cleaning +expect_false("score" %in% colnames(cleaned)) +expect_false("spectral_angle" %in% colnames(cleaned)) + +# One Run per intensity column; every wide row melted across every run +expect_equal(data.table::uniqueN(cleaned$Run), length(intensity_cols)) +expect_equal(nrow(cleaned), nrow(sage_raw) * length(intensity_cols)) + +# Zero intensities converted to NA; NA count equals the number of zero cells +expect_false(any(cleaned$Intensity == 0, na.rm = TRUE)) +n_zero_cells = sum(as.matrix(sage_raw[, intensity_cols, with = FALSE]) == 0, + na.rm = TRUE) +expect_equal(sum(is.na(cleaned$Intensity)), n_zero_cells) + +# q_value is retained for the downstream converter-level filter +expect_true("q_value" %in% colnames(cleaned)) + +# Combined charge states -> PrecursorCharge is -1 throughout on this fixture +expect_true(all(cleaned$PrecursorCharge == -1)) + +# A missing required column is an error naming the offending column +bad_input = data.table::copy(sage_raw) +bad_input$q_value = NULL +msstats_bad = MSstatsConvert::MSstatsImport( + list(input = bad_input), "MSstats", "Sage") +expect_error(MSstatsConvert:::.cleanRawSage(msstats_bad), "q_value") + +# No intensity columns is an error +only_fixed = sage_raw[, intersect(fixed_cols, colnames(sage_raw)), with = FALSE] +msstats_only_fixed = MSstatsConvert::MSstatsImport( + list(input = only_fixed), "MSstats", "Sage") +expect_error(MSstatsConvert:::.cleanRawSage(msstats_only_fixed), "intensity") + + +# Test .cleanRawSage on the charge-resolved fixture ------------------------ +sage_cr_path = system.file("tinytest/raw_data/Sage/lfq_charge_resolved.tsv", + package = "MSstatsConvert") +if (nzchar(sage_cr_path)) { + sage_cr = data.table::fread(sage_cr_path) + msstats_cr = MSstatsConvert::MSstatsImport( + list(input = sage_cr), "MSstats", "Sage") + cleaned_cr = MSstatsConvert:::.cleanRawSage(msstats_cr) + cleaned_cr = data.table::as.data.table(cleaned_cr) + + # Real precursor charges survive cleaning (not all -1) + expect_true(all(c(2, 3) %in% unique(cleaned_cr$PrecursorCharge))) + expect_false(all(cleaned_cr$PrecursorCharge == -1)) +} diff --git a/inst/tinytest/test_converters_SagetoMSstatsFormat.R b/inst/tinytest/test_converters_SagetoMSstatsFormat.R new file mode 100644 index 00000000..610c0a6d --- /dev/null +++ b/inst/tinytest/test_converters_SagetoMSstatsFormat.R @@ -0,0 +1,130 @@ +# Test SagetoMSstatsFormat --------------------------- +sage_path = 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_path) || !nzchar(annot_path)) { + exit_file("Sage fixtures not present in inst/tinytest/raw_data/Sage/") +} + +sage_raw = data.table::fread(sage_path) +annotation = data.table::fread(annot_path) + +fixed_cols = c("peptide", "charge", "proteins", "q_value", "score", + "spectral_angle") +intensity_cols = setdiff(colnames(sage_raw), fixed_cols) + +output = SagetoMSstatsFormat(sage_raw, annotation, use_log_file = FALSE) +output = data.table::as.data.table(output) + +# The 11 expected MSstats columns exist +expected_cols = c("ProteinName", "PeptideSequence", "PrecursorCharge", + "FragmentIon", "ProductCharge", "IsotopeLabelType", + "Condition", "BioReplicate", "Run", "Fraction", "Intensity") +for (col in expected_cols) { + expect_true(col %in% colnames(output)) +} +expect_equal(ncol(output), 11L) + +# Wide-to-long melt produced one Run per intensity column (8 columns -> 8 runs) +expect_equal(length(intensity_cols), 8L) +expect_equal(data.table::uniqueN(output$Run), length(intensity_cols)) + +# Run values map to the correct Condition and BioReplicate. +# Output Run is the standardized form, so standardize the annotation Run too. +annotation_std = data.table::as.data.table(annotation) +annotation_std[, Run := MSstatsConvert:::.standardizeColnames(Run)] +run_map = merge( + unique(output[, list(Run, Condition, BioReplicate)]), + annotation_std[, list(Run, Condition, BioReplicate)], + by = "Run", suffixes = c("_out", "_annot")) +expect_equal(nrow(run_map), length(intensity_cols)) +expect_true(all(as.character(run_map$Condition_out) == + as.character(run_map$Condition_annot))) +expect_true(all(as.character(run_map$BioReplicate_out) == + as.character(run_map$BioReplicate_annot))) + +# Zero intensities became NA, not 0 +expect_false(any(output$Intensity == 0, na.rm = TRUE)) +expect_true(any(is.na(output$Intensity))) + +# IsotopeLabelType is "L" everywhere, and the annotation's own IsotopeLabelType +# column neither broke the merge nor produced a duplicate column +expect_true("IsotopeLabelType" %in% colnames(annotation)) +expect_true(all(output$IsotopeLabelType == "L")) +expect_equal(sum(colnames(output) == "IsotopeLabelType"), 1L) + +# Regression: a plain Run/Condition/BioReplicate annotation (no IsotopeLabelType, +# no Fraction) still yields a Fraction column. Fraction is supplied by +# MSstatsBalancedDesign, not by columns_to_fill, so it must appear even when the +# annotation carries none. The main fixture annotation deliberately includes +# IsotopeLabelType, so the plain three-column case is otherwise untested. +annotation_min = annotation[, list(Run, Condition, BioReplicate)] +expect_false("Fraction" %in% colnames(annotation_min)) +expect_false("IsotopeLabelType" %in% colnames(annotation_min)) +output_min = data.table::as.data.table( + SagetoMSstatsFormat(sage_raw, annotation_min, use_log_file = FALSE)) +expect_true("Fraction" %in% colnames(output_min)) +expect_true(all(output_min$Fraction == 1)) + +# FragmentIon and ProductCharge are NA +expect_true(all(is.na(output$FragmentIon))) +expect_true(all(is.na(output$ProductCharge))) + +# q-value filter is load-bearing. Isolate its effect with removeFewMeasurements +# = FALSE so the only difference between the two runs is filter_with_Qvalue. +high_q_peptides = sage_raw[q_value > 0.01, unique(peptide)] +expect_true(length(high_q_peptides) >= 1L) + +out_filter = SagetoMSstatsFormat(sage_raw, annotation, + filter_with_Qvalue = TRUE, + removeFewMeasurements = FALSE, + use_log_file = FALSE) +out_filter = data.table::as.data.table(out_filter) +out_nofilter = SagetoMSstatsFormat(sage_raw, annotation, + filter_with_Qvalue = FALSE, + removeFewMeasurements = FALSE, + use_log_file = FALSE) +out_nofilter = data.table::as.data.table(out_nofilter) + +signal_with_filter = unique(out_filter[!is.na(Intensity), PeptideSequence]) +signal_without_filter = unique(out_nofilter[!is.na(Intensity), PeptideSequence]) + +# With filter FALSE the 5 high-q peptides keep signal; with TRUE they are gone +expect_true(all(high_q_peptides %in% signal_without_filter)) +expect_false(any(high_q_peptides %in% signal_with_filter)) + +# The peptides that lose all signal when the filter is switched on are exactly +# the high-q peptides +expect_equal(sort(setdiff(signal_without_filter, signal_with_filter)), + sort(high_q_peptides)) + + +# Test SagetoMSstatsFormat on the charge-resolved fixture ------------------- +sage_cr_path = system.file("tinytest/raw_data/Sage/lfq_charge_resolved.tsv", + package = "MSstatsConvert") +annot_cr_path = system.file("tinytest/raw_data/Sage/annotation_charge_resolved.csv", + package = "MSstatsConvert") +if (nzchar(sage_cr_path) && nzchar(annot_cr_path)) { + sage_cr = data.table::fread(sage_cr_path) + annot_cr = data.table::fread(annot_cr_path) + + # Single run: keep single-measurement features so nothing is dropped + out_cr = SagetoMSstatsFormat(sage_cr, annot_cr, + removeFewMeasurements = FALSE, + use_log_file = FALSE) + out_cr = data.table::as.data.table(out_cr) + + # Real charges come through, and -1 does not appear + expect_true(all(c(2, 3) %in% unique(out_cr$PrecursorCharge))) + expect_false(-1 %in% unique(out_cr$PrecursorCharge)) + + # A peptide seen at two charges yields two distinct features + charge_counts = out_cr[, list(n_charge = data.table::uniqueN(PrecursorCharge)), + by = PeptideSequence] + expect_true(any(charge_counts$n_charge == 2L)) + two_charge_pep = charge_counts[n_charge == 2L, PeptideSequence][1] + feats = unique(out_cr[PeptideSequence == two_charge_pep, + list(PeptideSequence, PrecursorCharge)]) + expect_equal(nrow(feats), 2L) +} diff --git a/man/MSstatsClean.Rd b/man/MSstatsClean.Rd index d8cccabe..f4e11e32 100644 --- a/man/MSstatsClean.Rd +++ b/man/MSstatsClean.Rd @@ -16,6 +16,7 @@ \alias{MSstatsClean,MSstatsMetamorpheusFiles-method} \alias{MSstatsClean,MSstatsProteinProspectorFiles-method} \alias{MSstatsClean,MSstatsMZMineFiles-method} +\alias{MSstatsClean,MSstatsSageFiles-method} \title{Clean files generated by a signal processing tools.} \usage{ MSstatsClean(msstats_object, ...) @@ -83,6 +84,8 @@ MSstatsClean(msstats_object, ...) \S4method{MSstatsClean}{MSstatsProteinProspectorFiles}(msstats_object) \S4method{MSstatsClean}{MSstatsMZMineFiles}(msstats_object, mzmine_annotations, sirius_annotations = NULL) + +\S4method{MSstatsClean}{MSstatsSageFiles}(msstats_object) } \arguments{ \item{msstats_object}{object that inherits from \code{MSstatsInputFiles} class.} @@ -250,6 +253,8 @@ data.table data.table +data.table + data.table } \description{ @@ -282,6 +287,8 @@ Clean Metamorpheus files Clean Protein Prospector files Clean MZMine files + +Clean Sage files } \examples{ evidence_path = system.file("tinytest/raw_data/MaxQuant/mq_ev.csv", diff --git a/man/MSstatsInputFiles.Rd b/man/MSstatsInputFiles.Rd index eb4f8646..d5f06d2b 100644 --- a/man/MSstatsInputFiles.Rd +++ b/man/MSstatsInputFiles.Rd @@ -18,6 +18,7 @@ \alias{MSstatsMetamorpheusFiles-class} \alias{MSstatsProteinProspectorFiles-class} \alias{MSstatsMZMineFiles-class} +\alias{MSstatsSageFiles-class} \title{Class to model files that describe a single MS dataset.} \description{ Class to model files that describe a single MS dataset. @@ -51,6 +52,8 @@ MSstatsMetamorpheusFiles: class for Metamorpheus files. MSstatsProteinProspectorFiles: class for ProteinProspector files. MSstatsMZMineFiles: class for MZMine files. + +MSstatsSageFiles: class for Sage files. } \section{Slots}{ diff --git a/man/SagetoMSstatsFormat.Rd b/man/SagetoMSstatsFormat.Rd new file mode 100644 index 00000000..06c5f647 --- /dev/null +++ b/man/SagetoMSstatsFormat.Rd @@ -0,0 +1,140 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/converters_SagetoMSstatsFormat.R +\name{SagetoMSstatsFormat} +\alias{SagetoMSstatsFormat} +\title{Import Sage LFQ files} +\usage{ +SagetoMSstatsFormat( + 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, + ... +) +} +\arguments{ +\item{input}{Sage \code{lfq.tsv} report, as a \code{data.frame}/\code{data.table} or a path. +Wide format with fixed columns \code{peptide}, \code{charge}, \code{proteins}, \code{q_value}, +\code{score}, \code{spectral_angle}, followed by one intensity column per input mzML, +each headed by the run's file name.} + +\item{annotation}{\code{data.frame} with \code{Run}, \code{Condition} and \code{BioReplicate} +columns (a \code{Fraction} column may also be supplied). This argument is +\strong{required}: Sage's \code{lfq.tsv} carries no experimental design, so condition +and replicate information must be provided separately.} + +\item{useUniquePeptide}{TRUE (default) removes peptides that are assigned for more than one proteins. +We assume to use unique peptide for each protein.} + +\item{removeFewMeasurements}{TRUE (default) will remove the features that have 1 or 2 measurements across runs.} + +\item{removeProtein_with1Peptide}{TRUE will remove the proteins which have only 1 peptide and charge. FALSE is default.} + +\item{qvalue_cutoff}{Cutoff for the \code{q_value} column. Default is 0.01.} + +\item{filter_with_Qvalue}{TRUE (default) replaces intensities whose \code{q_value} +exceeds \code{qvalue_cutoff} with \code{NA} (treated as censored missing downstream); +FALSE performs no q-value filtering. See the "FDR filtering" section for +why this matters.} + +\item{use_log_file}{logical. If TRUE, information about data processing +will be saved to a file.} + +\item{append}{logical. If TRUE, information about data processing will be added +to an existing log file.} + +\item{verbose}{logical. If TRUE, information about data processing will be printed +to the console.} + +\item{log_file_path}{character. Path to a file to which information about +data processing will be saved. +If not provided, such a file will be created automatically. +If \code{append = TRUE}, has to be a valid path to a file.} + +\item{...}{additional parameters to \code{data.table::fread}.} +} +\value{ +\code{data.frame} in the MSstats required format. +} +\description{ +Converts the label-free quantification report (\code{lfq.tsv}) produced by the Sage +search engine into a \code{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 \code{dataProcess}. +} +\section{Input file}{ + +Use \code{lfq.tsv}, which Sage writes only when \code{quant.lfq: true} is set in the +search configuration. Do \strong{not} use \code{results.sage.tsv}: its \code{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}{ + +\code{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 \code{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 +\code{peptide_q_value} set to \code{0.01} will still contain rows well above \code{0.01}. +This converter applies the filter: \code{filter_with_Qvalue} defaults to \code{TRUE} +and \code{qvalue_cutoff} defaults to \code{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 \code{0.01}. Set +\code{filter_with_Qvalue = FALSE} to return unfiltered data. +} + +\section{Charge states}{ + +Sage's \code{combine_charge_states} option (default \code{true}) sums charge states and +writes \code{charge} as \code{-1} for every row, so \code{PrecursorCharge} will be \code{-1} +throughout and the feature key reduces to the peptide sequence. Setting +\code{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 (\code{.}) +while preserving hyphens and underscores. The melted \code{Run} values (the +intensity column headers) and the annotation \code{Run} values are both passed +through this same standardization before merging, so they match automatically. +For example, a run named \code{B.naive_01steady-state.mzML.gz} in the annotation +becomes \code{Bnaive_01steady-statemzMLgz}; supply the raw file name in the +annotation and the merge resolves it. Note that the \code{Run} values in the +returned table are the standardized form. +} + +\section{Shared peptides}{ + +Sage pre-joins shared proteins into a single semicolon-delimited \code{proteins} +value (e.g. \verb{sp|A|X;sp|B|Y}). Because that is one \code{ProteinName} string rather +than several, MSstats' shared-peptide removal sees a single protein per +peptide and does not treat these rows as shared. Consequently +\code{useUniquePeptide = TRUE} has no effect on peptides that Sage reports against a +shared (semicolon-joined) protein group. +} + +\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) +} + +} diff --git a/man/dot-cleanRawSage.Rd b/man/dot-cleanRawSage.Rd new file mode 100644 index 00000000..25375aba --- /dev/null +++ b/man/dot-cleanRawSage.Rd @@ -0,0 +1,26 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/clean_Sage.R +\name{.cleanRawSage} +\alias{.cleanRawSage} +\title{Clean raw Sage LFQ output} +\usage{ +.cleanRawSage(msstats_object) +} +\arguments{ +\item{msstats_object}{an object of class \code{MSstatsSageFiles}.} +} +\value{ +data.table +} +\description{ +Operates on Sage's \code{lfq.tsv} report (produced when \code{quant.lfq: true}). This +is a wide-format table with fixed columns \code{peptide}, \code{charge}, \code{proteins}, +\code{q_value}, \code{score}, \code{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 \code{NA}. +} +\keyword{internal} diff --git a/man/dot-formatLogMessage.Rd b/man/dot-formatLogMessage.Rd index 907fb77b..f3639dd1 100644 --- a/man/dot-formatLogMessage.Rd +++ b/man/dot-formatLogMessage.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/utils_logging.R \name{.formatLogMessage} \alias{.formatLogMessage} -\title{Format a log line as "LEVEL \link{timestamp} message"} +\title{Format one log line: level, timestamp, message} \usage{ .formatLogMessage(level, ...) }