Skip to content

Commit e4b78fe

Browse files
tonywu1999tonywu1999
authored
feat(protein turnover): Add functionality to add heavy isotope label type to Spectronaut converter (#122)
Useful for protein turnover
1 parent 2b96a3d commit e4b78fe

12 files changed

Lines changed: 570 additions & 71 deletions

‎R/clean_Spectronaut.R‎

Lines changed: 178 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -4,60 +4,211 @@
44
#' @return `data.table`
55
#' @keywords internal
66
.cleanRawSpectronaut = function(msstats_object, intensity,
7-
calculateAnomalyScores,
8-
anomalyModelFeatures) {
7+
calculateAnomalyScores,
8+
anomalyModelFeatures,
9+
peptideSequenceColumn = "EG.ModifiedSequence",
10+
heavyLabels = NULL) {
911
FFrgLossType = FExcludedFromQuantification = NULL
10-
12+
1113
spec_input = getInputFile(msstats_object, "input")
12-
.validateSpectronautInput(spec_input)
14+
.validateSpectronautInput(spec_input, peptideSequenceColumn)
15+
spec_input = .addSpectronautColumnsIfMissing(spec_input)
1316
spec_input = spec_input[FFrgLossType == "noloss", ]
14-
15-
f_charge_col = .findAvailable(c("FCharge", "FFrgZ"), colnames(spec_input))
17+
18+
f_charge_col = .findAvailable(c("FCharge", "FFrgZ"), colnames(spec_input), fall_back = "FCharge")
1619
pg_qval_col = .findAvailable(c("PGQvalue"), colnames(spec_input))
17-
interference_col = .findAvailable(c("FPossibleInterference"),
20+
interference_col = .findAvailable(c("FPossibleInterference"),
1821
colnames(spec_input))
19-
exclude_col = .findAvailable(c("FExcludedFromQuantification"),
22+
exclude_col = .findAvailable(c("FExcludedFromQuantification"),
2023
colnames(spec_input))
21-
intensity_column_mapping = c(
22-
"PeakArea" = "FPeakArea",
23-
"NormalizedPeakArea" = "FNormalizedPeakArea",
24-
"MS1Quantity" = "FGMS1Quantity"
25-
)
26-
intensity = match.arg(intensity, names(intensity_column_mapping))
27-
intensity_column = intensity_column_mapping[[intensity]]
28-
cols = c("PGProteinGroups", "EGModifiedSequence", "FGCharge", "FFrgIon",
29-
f_charge_col, "RFileName", "RCondition", "RReplicate",
24+
intensity_col = .resolveSpectronautIntensityColumn(intensity, colnames(spec_input))
25+
peptide_col = .findAvailable(.standardizeColnames(peptideSequenceColumn),
26+
colnames(spec_input))
27+
protein_col = .findAvailable(c("PGProteinGroups", "PGProteinAccessions"),
28+
colnames(spec_input), fall_back = "PGProteinGroups")
29+
30+
cols = c(protein_col, peptide_col, "FGCharge", "FFrgIon",
31+
f_charge_col, "RFileName", "RCondition", "RReplicate",
3032
"EGQvalue", pg_qval_col, interference_col, exclude_col,
31-
intensity_column)
33+
intensity_col)
3234
if (calculateAnomalyScores){
3335
cols = c(cols, anomalyModelFeatures)
3436
}
3537
cols = intersect(cols, colnames(spec_input))
3638
spec_input = spec_input[, cols, with = FALSE]
39+
3740
data.table::setnames(
38-
spec_input,
39-
c("PGProteinGroups", "EGModifiedSequence", "FGCharge", "FFrgIon",
40-
f_charge_col, "RFileName", intensity_column,
41+
spec_input,
42+
c(protein_col, peptide_col, "FGCharge", "FFrgIon",
43+
f_charge_col, "RFileName", intensity_col,
4144
"RCondition", "RReplicate"),
4245
c("ProteinName", "PeptideSequence", "PrecursorCharge", "FragmentIon",
43-
"ProductCharge", "Run", "Intensity", "Condition", "BioReplicate"),
46+
"ProductCharge", "Run", "Intensity", "Condition", "BioReplicate"),
4447
skip_absent = TRUE)
48+
49+
spec_input = .assignSpectronautIsotopeLabelType(
50+
spec_input, heavyLabels)
51+
4552
.logSuccess("Spectronaut", "clean")
4653
spec_input
4754
}
4855

4956
#' Helper method to validate input has necessary columns
5057
#' @param spec_input dataframe input
58+
#' @param peptideSequenceColumn character, name of the column containing peptide
59+
#' sequences, passed from user
5160
#' @noRd
52-
.validateSpectronautInput = function(spec_input) {
53-
required_columns = c(
54-
"FFrgLossType", "FExcludedFromQuantification", "PGProteinGroups",
55-
"FGCharge")
61+
.validateSpectronautInput = function(spec_input, peptideSequenceColumn) {
62+
# Only FGCharge is truly required
63+
required_columns = c("FGCharge")
5664
missing_columns = setdiff(required_columns, colnames(spec_input))
5765
if (length(missing_columns) > 0) {
58-
msg = paste("The following columns are missing from the input data:",
66+
msg = paste("The following columns are missing from the input data:",
5967
paste(missing_columns, sep = ", ", collapse = ", "))
6068
getOption("MSstatsLog")("ERROR", msg)
6169
stop(msg)
6270
}
71+
# Ensure at least one protein name column is present
72+
if (!any(c("PGProteinGroups", "PGProteinAccessions") %in% colnames(spec_input))) {
73+
msg = paste("The following columns are missing from the input data:",
74+
"PGProteinGroups")
75+
getOption("MSstatsLog")("ERROR", msg)
76+
stop(msg)
77+
}
78+
# Ensure at least one protein name column is present
79+
if (!(.standardizeColnames(peptideSequenceColumn) %in% colnames(spec_input))) {
80+
msg = paste("The following column are missing from the input data:",
81+
peptideSequenceColumn)
82+
getOption("MSstatsLog")("ERROR", msg)
83+
stop(msg)
84+
}
85+
}
86+
87+
#' Add synthetic columns that are absent in protein turnover Spectronaut reports.
88+
#'
89+
#' Spectronaut's protein turnover export omits several columns that the
90+
#' standard MSstats pipeline expects. Rather than requiring callers to patch
91+
#' the data frame before passing it in (the "hacks" documented in the design
92+
#' document), we synthesize sensible defaults here so the rest of the pipeline
93+
#' sees a consistent schema.
94+
#'
95+
#' Columns synthesized when absent:
96+
#' \describe{
97+
#' \item{FFrgLossType}{"noloss" — the filter \code{spec_input[FFrgLossType == "noloss"]}
98+
#' keeps all rows, which is the correct behavior when the column is missing.}
99+
#' \item{FExcludedFromQuantification}{FALSE — no rows are excluded.}
100+
#' \item{FFrgIon}{NA — fragment ion identity is not available at the precursor
101+
#' level used by MS1-based protein turnover workflows.}
102+
#' \item{FCharge}{NA — product charge is not applicable when fragment-level
103+
#' columns are absent.}
104+
#' }
105+
#'
106+
#' @param spec_input `data.table` with standardized column names.
107+
#' @return `data.table` with missing columns added.
108+
#' @keywords internal
109+
#' @noRd
110+
.addSpectronautColumnsIfMissing = function(spec_input) {
111+
if (!("FFrgLossType" %in% colnames(spec_input))) {
112+
spec_input[, FFrgLossType := "noloss"]
113+
}
114+
if (!("FExcludedFromQuantification" %in% colnames(spec_input))) {
115+
spec_input[, FExcludedFromQuantification := FALSE]
116+
}
117+
if (!("FFrgIon" %in% colnames(spec_input))) {
118+
spec_input[, FFrgIon := NA_character_]
119+
}
120+
if (!("FCharge" %in% colnames(spec_input))) {
121+
spec_input[, FCharge := NA_integer_]
122+
}
123+
spec_input
124+
}
125+
126+
127+
#' Resolve the Spectronaut intensity column from user input.
128+
#'
129+
#' Accepts either a legacy enum alias or a raw (standardized) column name.
130+
#' The legacy aliases map to their canonical standardized column names so that
131+
#' old code continues to work unchanged. When the user passes a raw column
132+
#' name (e.g. \code{"FGMS1Quantity"} or \code{"FGMS2Quantity"}), it is used
133+
#' directly after verifying that the column exists.
134+
#'
135+
#' @param intensity Character scalar: enum alias or raw column name.
136+
#' @param available_cols Character vector of available standardized column names.
137+
#' @return The resolved standardized column name.
138+
#' @keywords internal
139+
#' @noRd
140+
.resolveSpectronautIntensityColumn = function(intensity, available_cols) {
141+
legacy_mapping = c(
142+
"PeakArea" = "FPeakArea",
143+
"NormalizedPeakArea" = "FNormalizedPeakArea"
144+
)
145+
146+
if (intensity %in% names(legacy_mapping)) {
147+
resolved = legacy_mapping[[intensity]]
148+
} else {
149+
resolved = .standardizeColnames(intensity)
150+
}
151+
152+
if (!(resolved %in% available_cols)) {
153+
stop(paste0(
154+
"Intensity column '", intensity, "' not found in input data. ", collapse = ", ")
155+
)
156+
}
157+
resolved
158+
}
159+
160+
161+
#' Assign IsotopeLabelType based on heavy label detection.
162+
#'
163+
#' When \code{heavyLabel} is provided, each row is classified as heavy
164+
#' (\code{"H"}), light (\code{"L"}), or unlabeled (\code{NA}) by inspecting
165+
#' the labeled sequence column for the presence of the label tag.
166+
#'
167+
#' In Spectronaut protein turnover reports, heavy peptides appear in
168+
#' \code{FG.LabeledSequence} with a bracketed modification, e.g.
169+
#' \code{_PEPTIDEK[Lys6]_}. Any sequence that contains
170+
#' \code{[<heavyLabel>]} is classified as heavy; all others are light.
171+
#' Sequences that do not have amino acids that can carry the label
172+
#' are classified as \code{NA}. For example, if \code{heavyLabels} is
173+
#' \code{"Lys6"}, then \code{PEPTIDEZ} is classified as NA since it
174+
#' has no lysine residues that could be labeled.
175+
#' When \code{heavyLabel} is \code{NULL} the column is left untouched so
176+
#' that the downstream \code{columns_to_fill} default of \code{"L"} applies,
177+
#' preserving backwards compatibility.
178+
#'
179+
#' @param spec_input `data.table` after column renaming.
180+
#' @param heavyLabels Character scalar heavy label name (e.g. \code{"Lys6"}),
181+
#' or \code{NULL}.
182+
#' @return `data.table` with \code{IsotopeLabelType} column added or updated.
183+
#' @keywords internal
184+
#' @noRd
185+
.assignSpectronautIsotopeLabelType = function(spec_input, heavyLabels) {
186+
IsotopeLabelType = PeptideSequence = StrippedSequence = NULL
187+
if (is.null(heavyLabels)) {
188+
return(spec_input)
189+
}
190+
191+
bare_amino_acids = sub("\\[.*\\]", "", heavyLabels)
192+
bare_amino_acids_pattern = paste0(bare_amino_acids, collapse = "|")
193+
heavy_pattern = paste0(heavyLabels, collapse = "|")
194+
heavy_brackets_escaped_pattern = paste(
195+
gsub("([\\[\\]])", "\\\\\\1", heavy_pattern, perl = TRUE),
196+
collapse = "|"
197+
)
198+
199+
spec_input[, StrippedSequence := gsub("\\[.*?\\]", "", PeptideSequence)]
200+
201+
spec_input[, IsotopeLabelType := data.table::fcase(
202+
grepl(heavy_brackets_escaped_pattern, PeptideSequence, perl = TRUE), "H",
203+
grepl(bare_amino_acids_pattern, StrippedSequence, perl = TRUE), "L",
204+
default = NA_character_
205+
)]
206+
207+
spec_input[, StrippedSequence := NULL]
208+
209+
for (i in seq_along(heavyLabels)) {
210+
escaped = gsub("([\\[\\]])", "\\\\\\1", heavyLabels[i], perl = TRUE)
211+
spec_input[, PeptideSequence := gsub(escaped, bare_amino_acids[i], PeptideSequence, perl = TRUE)]
212+
}
213+
spec_input
63214
}

‎R/converters_SpectronauttoMSstatsFormat.R‎

Lines changed: 43 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,28 @@
11
#' Import Spectronaut files
2-
#'
2+
#'
33
#' @param input name of Spectronaut output, which is long-format. ProteinName, PeptideSequence, PrecursorCharge, FragmentIon, ProductCharge, IsotopeLabelType, Condition, BioReplicate, Run, Intensity, F.ExcludedFromQuantification are required. Rows with F.ExcludedFromQuantification=True will be removed.
44
#' @param annotation name of 'annotation.txt' data which includes Condition, BioReplicate, Run. If annotation is already complete in Spectronaut, use annotation=NULL (default). It will use the annotation information from input.
5-
#' @param intensity 'PeakArea'(default) uses not normalized MS2 peak area. 'NormalizedPeakArea' uses MS2 peak area normalized by Spectronaut.
6-
#' 'MS1Quantity' uses MS1 level quantification, which should be used if MS2 is unreliable.
5+
#' @param intensity Intensity column to use. Accepts legacy enum values
6+
#' \code{'PeakArea'} (default, uses F.PeakArea), \code{'NormalizedPeakArea'}
7+
#' (uses F.NormalizedPeakArea). Can also be any raw
8+
#' Spectronaut column name passed as a string (e.g.
9+
#' \code{"FG.MS1Quantity"}); the column name is standardized internally.
10+
#' For protein turnover workflows the recommended default is
11+
#' \code{"FG.MS1Quantity"}.
12+
#' @param peptideSequenceColumn Name of the Spectronaut column that contains the
13+
#' peptide sequence. Defaults to \code{"EG.ModifiedSequence"}. The value is
14+
#' standardized internally (dots and spaces removed) before column lookup.
15+
#' @param heavyLabels Character list identifying the heavy isotope labels as it
16+
#' appears inside square brackets in the peptide sequence column, e.g.
17+
#' \code{c("Lys6")} matches peptides containing \code{[Lys6]}.
18+
#' \code{c("Lys6", "Arg10")} matches peptides containing either \code{[Lys6]} or \code{[Arg10]}.
19+
#' Supports any novel label name reported by Spectronaut (e.g. \code{"Leu6"},
20+
#' \code{"Phe10"}). When provided, peptides are
21+
#' classified as heavy (\code{IsotopeLabelType = "H"}), light
22+
#' (\code{IsotopeLabelType = "L"}), or unlabeled
23+
#' (\code{IsotopeLabelType = NA}) based on its labeled sequence. When
24+
#' \code{NULL} (default) all peptides receive \code{IsotopeLabelType = "L"}.
25+
#' Useful for protein turnover experiments.
726
#' @param excludedFromQuantificationFilter Remove rows with F.ExcludedFromQuantification=TRUE Default is TRUE.
827
#' @param filter_with_Qvalue FALSE(default) will not perform any filtering. TRUE will filter out the intensities that have greater than qvalue_cutoff in EG.Qvalue column. Those intensities will be replaced with zero and will be considered as censored missing values for imputation purpose.
928
#' @param qvalue_cutoff Cutoff for EG.Qvalue. default is 0.01.
@@ -33,8 +52,10 @@
3352
#' head(spectronaut_imported)
3453
#'
3554
SpectronauttoMSstatsFormat = function(
36-
input, annotation = NULL,
37-
intensity = c('PeakArea', 'NormalizedPeakArea', 'MS1Quantity'),
55+
input, annotation = NULL,
56+
intensity = 'PeakArea',
57+
peptideSequenceColumn = "EG.ModifiedSequence",
58+
heavyLabels = NULL,
3859
excludedFromQuantificationFilter = TRUE,
3960
filter_with_Qvalue = FALSE, qvalue_cutoff = 0.01,
4061
useUniquePeptide = TRUE, removeFewMeasurements=TRUE,
@@ -47,10 +68,11 @@ SpectronauttoMSstatsFormat = function(
4768
use_log_file = TRUE, append = FALSE, verbose = TRUE,
4869
log_file_path = NULL, ...
4970
) {
71+
5072
validation_config = list(
51-
input = input,
52-
annotation = annotation,
53-
intensity = intensity,
73+
input = input,
74+
annotation = annotation,
75+
intensity = intensity,
5476
excludedFromQuantificationFilter = excludedFromQuantificationFilter,
5577
filter_with_Qvalue = filter_with_Qvalue,
5678
qvalue_cutoff = qvalue_cutoff,
@@ -84,8 +106,10 @@ SpectronauttoMSstatsFormat = function(
84106
"MSstats", "Spectronaut", ...)
85107

86108
input = MSstatsConvert::MSstatsClean(input, intensity = intensity,
87-
calculateAnomalyScores,
88-
anomalyModelFeatures)
109+
calculateAnomalyScores,
110+
anomalyModelFeatures,
111+
peptideSequenceColumn = peptideSequenceColumn,
112+
heavyLabels = heavyLabels)
89113
annotation = MSstatsConvert::MSstatsMakeAnnotation(input, annotation)
90114

91115
pq_filter = list(score_column = "PGQvalue",
@@ -113,20 +137,24 @@ SpectronauttoMSstatsFormat = function(
113137
drop_column = TRUE
114138
)
115139

116-
feature_columns = c("PeptideSequence", "PrecursorCharge",
140+
feature_columns = c("PeptideSequence", "PrecursorCharge",
117141
"FragmentIon", "ProductCharge")
142+
143+
fill_isotope_label_type = if (is.null(heavyLabels))
144+
list("IsotopeLabelType" = "L") else list()
145+
118146
input = MSstatsConvert::MSstatsPreprocess(
119-
input,
120-
annotation,
147+
input,
148+
annotation,
121149
feature_columns,
122150
remove_shared_peptides = useUniquePeptide,
123151
remove_single_feature_proteins = removeProtein_with1Feature,
124152
feature_cleaning = list(remove_features_with_few_measurements = removeFewMeasurements,
125153
summarize_multiple_psms = summaryforMultipleRows),
126-
score_filtering = list(pgq = pq_filter,
154+
score_filtering = list(pgq = pq_filter,
127155
psm_q = qval_filter),
128156
exact_filtering = list(excluded_quant = excluded_quant_filter),
129-
columns_to_fill = list("IsotopeLabelType" = "L"),
157+
columns_to_fill = fill_isotope_label_type,
130158
anomaly_metrics = anomalyModelFeatures)
131159
input[, Intensity := ifelse(Intensity == 0, NA, Intensity)]
132160

‎R/utils_balanced_design.R‎

Lines changed: 23 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -94,19 +94,35 @@
9494
measurement_col, group_col)
9595
result[, 2:4, with = FALSE]
9696
} else {
97-
labels = unique(input[["IsotopeLabelType"]])
97+
labels = na.omit(unique(input[["IsotopeLabelType"]]))
9898
groups = unique(input[[group_col]])
9999
by_group = vector("list", length(groups))
100100
measurements = unique(input[[measurement_col]])
101101
for (group_id in seq_along(groups)) {
102102
group = groups[group_id]
103103
group_filter = input[[group_col]] == group
104-
by_group[[group_id]] = data.table::as.data.table(
105-
expand.grid(
106-
labels = labels,
107-
features = unique(input[[feature_col]][group_filter]),
108-
measurements = unique(input[[measurement_col]][group_filter])
109-
))
104+
non_na_filter = group_filter & !is.na(input[["IsotopeLabelType"]])
105+
na_filter = group_filter & is.na(input[["IsotopeLabelType"]])
106+
non_na_features = unique(input[[feature_col]][non_na_filter])
107+
na_label_features = unique(input[[feature_col]][na_filter])
108+
parts = list()
109+
if (length(non_na_features) > 0) {
110+
parts[[length(parts) + 1]] = data.table::as.data.table(
111+
expand.grid(
112+
labels = labels,
113+
features = non_na_features,
114+
measurements = unique(input[[measurement_col]][group_filter])
115+
))
116+
}
117+
if (length(na_label_features) > 0) {
118+
parts[[length(parts) + 1]] = data.table::as.data.table(
119+
expand.grid(
120+
labels = NA,
121+
features = na_label_features,
122+
measurements = unique(input[[measurement_col]][group_filter])
123+
))
124+
}
125+
by_group[[group_id]] = data.table::rbindlist(parts)
110126
by_group[[group_id]]$group = group
111127
}
112128
result = data.table::rbindlist(by_group)

0 commit comments

Comments
 (0)