diff --git a/DESCRIPTION b/DESCRIPTION index 83fdc00..8a6629b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -15,8 +15,9 @@ Description: MSstats package provide tools for preprocessing, summarization and License: Artistic-2.0 Encoding: UTF-8 RoxygenNote: 7.3.3 -Imports: +Imports: arrow, + data.table, DBI, dplyr, MSstats, diff --git a/NAMESPACE b/NAMESPACE index 9ec4823..8312349 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -10,5 +10,9 @@ importFrom(MSstats,groupComparison) importFrom(MSstatsConvert,MSstatsClean) importFrom(MSstatsConvert,MSstatsImport) importFrom(MSstatsConvert,MSstatsMakeAnnotation) +importFrom(data.table,":=") +importFrom(data.table,.SD) +importFrom(data.table,setDT) +importFrom(data.table,setnames) importFrom(utils,head) importFrom(utils,sessionInfo) diff --git a/R/clean_spectronaut.R b/R/clean_spectronaut.R index 1eeb89c..f57bbce 100644 --- a/R/clean_spectronaut.R +++ b/R/clean_spectronaut.R @@ -1,3 +1,6 @@ +#' @importFrom data.table := .SD setDT setnames +NULL + #' @keywords internal reduceBigSpectronaut <- function(input_file, output_path, intensity="F.NormalizedPeakArea", @@ -5,29 +8,108 @@ reduceBigSpectronaut <- function(input_file, output_path, filter_by_identified = FALSE, filter_by_qvalue = TRUE, qvalue_cutoff = 0.01, - calculateAnomalyScores=FALSE, - anomalyModelFeatures=c()) { + calculateAnomalyScores=FALSE, + anomalyModelFeatures=c(), + block_size = 16L * 1024L * 1024L) { + block_size <- as.integer(block_size) + stopifnot(length(block_size) == 1L, !is.na(block_size), block_size > 0L) + if (grepl("csv", input_file)) { - delim = "," + delim <- "," } else if (grepl("tsv|xls", input_file)) { - delim = "\t" + delim <- "\t" } else { delim <- ";" } - spec_chunk <- function(x, pos) cleanSpectronautChunk(x, - output_path, - intensity, - filter_by_excluded, - filter_by_identified, - filter_by_qvalue, - qvalue_cutoff, - pos, - calculateAnomalyScores, - anomalyModelFeatures) - readr::read_delim_chunked(input_file, - readr::DataFrameCallback$new(spec_chunk), - delim = delim, - chunk_size = 1e6) + + # Columns cleanSpectronautChunk actually consumes; Arrow's + # convert_options$include_columns drops everything else at parse time so + # we never materialize the ~35 unused columns Spectronaut exports. + needed_cols <- c("R.FileName", "R.Condition", "R.Replicate", + "PG.ProteinAccessions", "EG.ModifiedSequence", + "FG.LabeledSequence", "FG.Charge", + "F.FrgIon", "F.Charge", + "EG.Identified", "F.ExcludedFromQuantification", + "F.FrgLossType", "PG.Qvalue", "EG.Qvalue", + intensity) + if (calculateAnomalyScores) { + needed_cols <- c(needed_cols, anomalyModelFeatures) + } + + # Arrow's CSV reader replaces readr::read_delim_chunked. Arrow releases + # per-batch state as soon as a batch is consumed, so peak memory is + # bounded by one record batch instead of growing with the dataset (readr + # keeps a string-interning pool that accumulates across chunks). The + # `delim` switch above already covers comma / tab / semicolon variants; + # Arrow's CSV reader handles all three the same way through + # CsvParseOptions$delimiter. + parse_opts <- arrow::CsvParseOptions$create(delimiter = delim) + convert_opts <- arrow::CsvConvertOptions$create() + read_opts <- arrow::CsvReadOptions$create(block_size = block_size) + + ds <- arrow::open_dataset( + input_file, + format = "csv", + parse_options = parse_opts, + convert_options = convert_opts, + read_options = read_opts + ) + + # Project to only the columns cleanSpectronautChunk consumes. Scanner + # projection (not CsvConvertOptions$include_columns, which collides with the + # open_dataset schema layer) is the mechanism that works with the dataset + # API. intersect() keeps only columns actually present so a partial export + # still runs instead of erroring, matching cleanSpectronautChunk's handling. + present_cols <- intersect(needed_cols, names(ds)) + # Scanner applies the projection; ToRecordBatchReader() then yields one + # projected batch at a time on demand, keeping peak memory to a single batch. + scanner <- arrow::Scanner$create(ds, projection = present_cols) + reader <- scanner$ToRecordBatchReader() + + t_start <- Sys.time() + pos <- 1L + batch_idx <- 0L + repeat { + batch <- reader$read_next_batch() + if (is.null(batch)) break + chunk_df <- as.data.frame(batch) + cleanSpectronautChunk(chunk_df, + output_path, + intensity, + filter_by_excluded, + filter_by_identified, + filter_by_qvalue, + qvalue_cutoff, + pos, + calculateAnomalyScores, + anomalyModelFeatures) + pos <- pos + nrow(chunk_df) + batch_idx <- batch_idx + 1L + + if (batch_idx %% 1000L == 0L) { + elapsed <- as.numeric(Sys.time() - t_start, units = "secs") + rate <- (pos - 1L) / elapsed + message(sprintf( + "[reduceBigSpectronaut] %d batches | %s rows | %.1fk rows/s | %.0fs elapsed", + batch_idx, + format(pos - 1L, big.mark = ","), + rate / 1000, + elapsed)) + } + + rm(batch, chunk_df) + } + + if (batch_idx %% 1000L != 0L) { + elapsed <- as.numeric(Sys.time() - t_start, units = "secs") + rate <- (pos - 1L) / elapsed + message(sprintf( + "[reduceBigSpectronaut] done: %d batches | %s rows | %.1fk rows/s | %.0fs elapsed", + batch_idx, + format(pos - 1L, big.mark = ","), + rate / 1000, + elapsed)) + } } #' @keywords internal @@ -38,85 +120,108 @@ cleanSpectronautChunk = function(input, output_path, filter_by_qvalue = TRUE, qvalue_cutoff = 0.01, pos = NULL, - calculateAnomalyScores=FALSE, + calculateAnomalyScores=FALSE, anomalyModelFeatures=c()) { + data.table::setDT(input) + all_cols <- c("R.FileName", "R.Condition", "R.Replicate", "PG.ProteinAccessions", "EG.ModifiedSequence", "FG.LabeledSequence", "FG.Charge", "F.FrgIon", "F.Charge", "EG.Identified", "F.ExcludedFromQuantification", "F.FrgLossType", "PG.Qvalue", "EG.Qvalue", intensity) - - if (calculateAnomalyScores){ - all_cols <- c(all_cols, anomalyModelFeatures) - } - - cols <- intersect(all_cols, colnames(input)) - input <- dplyr::select(input, all_of(cols)) - input <- dplyr::rename_with(input, .fn = MSstatsConvert:::.standardizeColnames) - new_names <- c("Run", "Condition", "BioReplicate", "ProteinName", "PeptideSequence", "LabeledSequence", "PrecursorCharge", "FragmentIon", "ProductCharge", "Identified", "Excluded", "FFrgLossType", "PGQvalue", "EGQvalue", "Intensity") - if (calculateAnomalyScores){ + if (calculateAnomalyScores) { + all_cols <- c(all_cols, anomalyModelFeatures) new_names <- c(new_names, MSstatsConvert:::.standardizeColnames(anomalyModelFeatures)) } - - # non_standardized = - old_names <- MSstatsConvert:::.standardizeColnames(all_cols) - names(old_names) <- new_names - old_names <- old_names[old_names %in% colnames(input)] - - input <- dplyr::rename(input, !!old_names) - input <- dplyr::mutate(input, Intensity = as.numeric(Intensity)) - - if (is.character(dplyr::pull(dplyr::collect(head(dplyr::select(input, Excluded))), Excluded))) { - input <- dplyr::mutate(input, Excluded = Excluded == "True") + + present_orig <- intersect(all_cols, colnames(input)) + if (length(present_orig) == 0L) { + stop(sprintf( + paste0("cleanSpectronautChunk: none of the expected Spectronaut ", + "columns were found in the input batch. ", + "Expected any of: %s. Found: %s. ", + "Check that the file is comma-delimited and that the ", + "Spectronaut export uses the standard column names."), + paste(all_cols, collapse = ", "), + paste(colnames(input), collapse = ", "))) + } + input <- input[, present_orig, with = FALSE] + + # Two-step rename matching the MSstatsConvert family pattern: standardize + # all column names, then map standardized -> MSstats final names. + data.table::setnames(input, MSstatsConvert:::.standardizeColnames(colnames(input))) + std_to_msstats <- stats::setNames(new_names, + MSstatsConvert:::.standardizeColnames(all_cols)) + data.table::setnames(input, + old = names(std_to_msstats), + new = unname(std_to_msstats), + skip_absent = TRUE) + + input[, Intensity := as.numeric(Intensity)] + + if (is.character(input[["Excluded"]])) { + input[, Excluded := Excluded == "True"] + } + if ("Identified" %in% colnames(input) && is.character(input[["Identified"]])) { + input[, Identified := Identified == "True"] } - if (is.element("Identified", colnames(input))) { - if (is.character(dplyr::pull(dplyr::collect(head(dplyr::select(input, Identified))), Identified))) { - input <- dplyr::mutate(input, Identified = Identified == "True") + + # Fail fast if a filter's source column is absent, rather than letting + # data.table raise a cryptic "object not found" mid-run. The message reports + # the original Spectronaut column name (what the user controls in the export). + msstats_to_spectronaut <- stats::setNames(all_cols, new_names) + require_filter_cols <- function(cols, filter_name) { + missing <- setdiff(cols, colnames(input)) + if (length(missing) > 0L) { + stop(sprintf( + paste0("cleanSpectronautChunk: %s needs Spectronaut column(s) %s, ", + "which were not found in the input export. Found: %s."), + filter_name, + paste(msstats_to_spectronaut[missing], collapse = ", "), + paste(colnames(input), collapse = ", "))) } } - + if (filter_by_excluded) { - input <- dplyr::mutate( - input, Intensity = dplyr::if_else(Excluded, NA_real_, Intensity)) - + require_filter_cols("Excluded", "filter_by_excluded") + input[Excluded == TRUE, Intensity := NA_real_] } - if (filter_by_identified) { - input <- dplyr::mutate( - input, Intensity = dplyr::if_else(Identified, Intensity, NA_real_)) + require_filter_cols("Identified", "filter_by_identified") + input[Identified == FALSE, Intensity := NA_real_] } - if (filter_by_qvalue) { - input <- dplyr::mutate( - input, - Intensity = dplyr::if_else(EGQvalue < qvalue_cutoff, Intensity, NA_real_)) - input <- dplyr::mutate( - input, - Intensity = dplyr::if_else(PGQvalue < qvalue_cutoff, Intensity, NA_real_)) + require_filter_cols(c("EGQvalue", "PGQvalue"), "filter_by_qvalue") + # Preserve dplyr::if_else semantics: rows with NA q-values become NA. + input[is.na(EGQvalue) | EGQvalue >= qvalue_cutoff, Intensity := NA_real_] + input[is.na(PGQvalue) | PGQvalue >= qvalue_cutoff, Intensity := NA_real_] } - - input <- dplyr::filter(input, FFrgLossType == "noloss") - if (is.element("LabeledSequence", colnames(input))) { - input <- dplyr::mutate(input, IsLabeled = grepl("Lys8", LabeledSequence) | grepl("Arg10", LabeledSequence)) - input <- dplyr::mutate(input, IsotopeLabelType := dplyr::if_else(IsLabeled, "H", "L")) + + require_filter_cols("FFrgLossType", "the noloss fragment-loss filter") + input <- input[FFrgLossType == "noloss"] + + if ("LabeledSequence" %in% colnames(input)) { + input[, IsotopeLabelType := ifelse( + grepl("Lys8", LabeledSequence) | grepl("Arg10", LabeledSequence), + "H", "L")] } else { - input <- dplyr::mutate(input, IsotopeLabelType = "L") + input[, IsotopeLabelType := "L"] } - - select_cols = c("ProteinName", "PeptideSequence", "PrecursorCharge", "FragmentIon", - "ProductCharge", "IsotopeLabelType", "Run", "BioReplicate", "Condition", - "Intensity") - if (calculateAnomalyScores){ - select_cols = c(select_cols, - MSstatsConvert:::.standardizeColnames(anomalyModelFeatures)) + + select_cols <- c("ProteinName", "PeptideSequence", "PrecursorCharge", "FragmentIon", + "ProductCharge", "IsotopeLabelType", "Run", "BioReplicate", "Condition", + "Intensity") + if (calculateAnomalyScores) { + select_cols <- c(select_cols, + MSstatsConvert:::.standardizeColnames(anomalyModelFeatures)) } - - input <- dplyr::select(input, select_cols) + + input <- input[, select_cols, with = FALSE] .writeChunkToFile(input, output_path, pos) NULL } diff --git a/R/converters.R b/R/converters.R index 13b4383..668edb4 100644 --- a/R/converters.R +++ b/R/converters.R @@ -117,6 +117,50 @@ bigFragPipetoMSstatsFormat <- function(input_file, output_file_name, #' @param filter_by_identified if TRUE, will filter by the `EG.Identified` column. #' @param filter_by_qvalue if TRUE, will filter by EG.Qvalue and PG.Qvalue columns. #' @param qvalue_cutoff cutoff which will be used for q-value filtering. +#' @param block_size Arrow CSV reader block size in bytes: the amount of raw +#' CSV that is read and parsed at a time (roughly one record batch per block). +#' Defaults to 16 MiB (`16L * 1024L * 1024L`). Two things matter when tuning it: +#' \itemize{ +#' \item It must be larger than your single widest row. If you see +#' `Invalid: straddling object straddles two block boundaries` on +#' extra-wide exports, raise it (e.g. `64L * 1024L * 1024L`). +#' \item Larger blocks reduce parsing overhead, but peak memory grows +#' roughly linearly with `block_size`, so it is also the main lever on +#' out-of-memory risk. Gains taper off quickly because per-batch time is +#' dominated by row cleaning, not by block size. +#' } +#' To speed up processing without risking a crash, double the value +#' (32 -> 64 -> 128 MiB) while watching memory use, and stop once throughput +#' stops improving. Keep it comfortably below free RAM -- and below free RAM +#' divided by the number of jobs if you run several conversions in parallel. +#' +#' @details +#' `bigSpectronauttoMSstatsFormat` consumes the following Spectronaut export +#' columns. +#' +#' Always required: +#' \itemize{ +#' \item `R.FileName`, `R.Condition`, `R.Replicate` -> Run, Condition, BioReplicate +#' \item `PG.ProteinAccessions` -> ProteinName +#' \item `EG.ModifiedSequence` -> PeptideSequence +#' \item `FG.Charge` -> PrecursorCharge; `F.FrgIon` -> FragmentIon; `F.Charge` -> ProductCharge +#' \item the intensity column (`F.NormalizedPeakArea` by default; see `intensity`) -> Intensity +#' \item `F.FrgLossType` -> used by the always-on `noloss` fragment filter +#' } +#' +#' Optional: +#' \itemize{ +#' \item `FG.LabeledSequence` -> if present, sets IsotopeLabelType ("H" for +#' Lys8/Arg10 labels, otherwise "L"); if absent, IsotopeLabelType defaults to "L". +#' } +#' +#' Required only when the matching filter is enabled: +#' \itemize{ +#' \item `EG.Identified` -> `filter_by_identified` +#' \item `F.ExcludedFromQuantification` -> `filter_by_excluded` +#' \item `EG.Qvalue` and `PG.Qvalue` -> `filter_by_qvalue` +#' } +#' A missing required column raises an informative error naming the column. #' #' @export #' @@ -143,14 +187,16 @@ bigSpectronauttoMSstatsFormat <- function(input_file, output_file_name, aggregate_psms = FALSE, filter_few_obs = FALSE, remove_annotation = FALSE, - calculateAnomalyScores=FALSE, + calculateAnomalyScores=FALSE, anomalyModelFeatures=c(), - connection = NULL) { + connection = NULL, + block_size = 16L * 1024L * 1024L) { reduced_file <- .prefixedPath("reduce_output_", output_file_name) reduceBigSpectronaut(input_file, reduced_file, intensity, filter_by_excluded, filter_by_identified, filter_by_qvalue, qvalue_cutoff, - calculateAnomalyScores, anomalyModelFeatures) + calculateAnomalyScores, anomalyModelFeatures, + block_size = block_size) msstats_data <- MSstatsPreprocessBig( input_file = reduced_file, output_file_name = output_file_name, diff --git a/man/bigSpectronauttoMSstatsFormat.Rd b/man/bigSpectronauttoMSstatsFormat.Rd index 01706ef..22d4c2d 100644 --- a/man/bigSpectronauttoMSstatsFormat.Rd +++ b/man/bigSpectronauttoMSstatsFormat.Rd @@ -20,7 +20,8 @@ bigSpectronauttoMSstatsFormat( remove_annotation = FALSE, calculateAnomalyScores = FALSE, anomalyModelFeatures = c(), - connection = NULL + connection = NULL, + block_size = 16L * 1024L * 1024L ) } \arguments{ @@ -63,6 +64,23 @@ using dataProcess function. Only applicable to sparklyr backend.} \item{connection}{Connection to a spark instance created with the `spark_connect` function from `sparklyr` package.} + +\item{block_size}{Arrow CSV reader block size in bytes: the amount of raw +CSV that is read and parsed at a time (roughly one record batch per block). +Defaults to 16 MiB (`16L * 1024L * 1024L`). Two things matter when tuning it: +\itemize{ + \item It must be larger than your single widest row. If you see + `Invalid: straddling object straddles two block boundaries` on + extra-wide exports, raise it (e.g. `64L * 1024L * 1024L`). + \item Larger blocks reduce parsing overhead, but peak memory grows + roughly linearly with `block_size`, so it is also the main lever on + out-of-memory risk. Gains taper off quickly because per-batch time is + dominated by row cleaning, not by block size. +} +To speed up processing without risking a crash, double the value +(32 -> 64 -> 128 MiB) while watching memory use, and stop once throughput +stops improving. Keep it comfortably below free RAM -- and below free RAM +divided by the number of jobs if you run several conversions in parallel.} } \value{ either arrow object or sparklyr table that can be optionally collected @@ -71,6 +89,34 @@ into memory by using dplyr::collect function. \description{ Convert out-of-memory Spectronaut files to MSstats format. } +\details{ +`bigSpectronauttoMSstatsFormat` consumes the following Spectronaut export +columns. + +Always required: +\itemize{ + \item `R.FileName`, `R.Condition`, `R.Replicate` -> Run, Condition, BioReplicate + \item `PG.ProteinAccessions` -> ProteinName + \item `EG.ModifiedSequence` -> PeptideSequence + \item `FG.Charge` -> PrecursorCharge; `F.FrgIon` -> FragmentIon; `F.Charge` -> ProductCharge + \item the intensity column (`F.NormalizedPeakArea` by default; see `intensity`) -> Intensity + \item `F.FrgLossType` -> used by the always-on `noloss` fragment filter +} + +Optional: +\itemize{ + \item `FG.LabeledSequence` -> if present, sets IsotopeLabelType ("H" for + Lys8/Arg10 labels, otherwise "L"); if absent, IsotopeLabelType defaults to "L". +} + +Required only when the matching filter is enabled: +\itemize{ + \item `EG.Identified` -> `filter_by_identified` + \item `F.ExcludedFromQuantification` -> `filter_by_excluded` + \item `EG.Qvalue` and `PG.Qvalue` -> `filter_by_qvalue` +} +A missing required column raises an informative error naming the column. +} \examples{ converted_data <- bigSpectronauttoMSstatsFormat( system.file("extdata", "spectronaut_input.csv", package = "MSstatsBig"), diff --git a/man/dot-prefixedPath.Rd b/man/dot-prefixedPath.Rd new file mode 100644 index 0000000..be036ac --- /dev/null +++ b/man/dot-prefixedPath.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/utils.R +\name{.prefixedPath} +\alias{.prefixedPath} +\title{Build an intermediate output path by prefixing only the basename.} +\usage{ +.prefixedPath(prefix, path) +} +\arguments{ +\item{prefix}{Character scalar prepended to the basename.} + +\item{path}{Output file path supplied by the caller.} +} +\value{ +Character scalar. +} +\description{ +Naive `paste0(prefix, output_file_name)` corrupts paths that contain a +directory (`subdir/out.csv` → `topN_subdir/out.csv`, +`/tmp/out.csv` → `topN_/tmp/out.csv`). Splitting via dirname/basename keeps +the directory component intact so intermediate files land beside the final +output. +} +\keyword{internal} diff --git a/tests/testthat/test-converters.R b/tests/testthat/test-converters.R index 78f6da3..89b3862 100644 --- a/tests/testthat/test-converters.R +++ b/tests/testthat/test-converters.R @@ -94,6 +94,258 @@ test_that("bigSpectronauttoMSstatsFormat works correctly", { unlink(paste0("reduce_output_", output_file), recursive = TRUE, force = TRUE) }) +make_spectronaut_input <- function(n = 1L, ...) { + base <- data.frame( + R.FileName = rep("run1", n), + R.Condition = "A", + R.Replicate = 1L, + PG.ProteinAccessions = "P1", + EG.ModifiedSequence = paste0("PEP", seq_len(n)), + FG.LabeledSequence = paste0("PEP", seq_len(n)), + FG.Charge = 2L, + F.FrgIon = "y1", + F.Charge = 1L, + EG.Identified = "True", + F.ExcludedFromQuantification = "False", + F.FrgLossType = "noloss", + PG.Qvalue = 0.001, + EG.Qvalue = 0.001, + F.NormalizedPeakArea = seq_len(n) * 100, + stringsAsFactors = FALSE + ) + overrides <- list(...) + for (col in names(overrides)) base[[col]] <- overrides[[col]] + base +} + +test_that("cleanSpectronautChunk produces the expected MSstats schema", { + input <- make_spectronaut_input(n = 1L) + output_file <- tempfile(fileext = ".csv") + on.exit(unlink(output_file, force = TRUE), add = TRUE) + + MSstatsBig:::cleanSpectronautChunk(input, output_file, pos = 1L) + result <- readr::read_csv(output_file, show_col_types = FALSE) + + expected_cols <- c("ProteinName", "PeptideSequence", "PrecursorCharge", "FragmentIon", + "ProductCharge", "IsotopeLabelType", "Run", "BioReplicate", "Condition", + "Intensity") + expect_setequal(colnames(result), expected_cols) + expect_equal(nrow(result), 1L) + expect_equal(result$Intensity, 100) + expect_equal(result$IsotopeLabelType, "L") +}) + +test_that("cleanSpectronautChunk filter_by_excluded sets Intensity to NA on excluded rows", { + input <- make_spectronaut_input( + n = 2L, + F.ExcludedFromQuantification = c("True", "False"), + F.NormalizedPeakArea = c(100, 200) + ) + output_file <- tempfile(fileext = ".csv") + on.exit(unlink(output_file, force = TRUE), add = TRUE) + + MSstatsBig:::cleanSpectronautChunk(input, output_file, pos = 1L, + filter_by_excluded = TRUE, + filter_by_qvalue = FALSE) + result <- readr::read_csv(output_file, show_col_types = FALSE) + result <- result[order(result$PeptideSequence), ] + + expect_true(is.na(result$Intensity[1])) + expect_equal(result$Intensity[2], 200) +}) + +test_that("cleanSpectronautChunk filter_by_identified sets Intensity to NA on unidentified rows", { + input <- make_spectronaut_input( + n = 2L, + EG.Identified = c("True", "False"), + F.NormalizedPeakArea = c(100, 200) + ) + output_file <- tempfile(fileext = ".csv") + on.exit(unlink(output_file, force = TRUE), add = TRUE) + + MSstatsBig:::cleanSpectronautChunk(input, output_file, pos = 1L, + filter_by_identified = TRUE, + filter_by_qvalue = FALSE) + result <- readr::read_csv(output_file, show_col_types = FALSE) + result <- result[order(result$PeptideSequence), ] + + expect_equal(result$Intensity[1], 100) + expect_true(is.na(result$Intensity[2])) +}) + +test_that("cleanSpectronautChunk filter_by_qvalue NA-aware semantics match dplyr::if_else", { + input <- make_spectronaut_input( + n = 3L, + EG.Qvalue = c(0.001, 0.5, NA_real_), + F.NormalizedPeakArea = c(100, 200, 300) + ) + output_file <- tempfile(fileext = ".csv") + on.exit(unlink(output_file, force = TRUE), add = TRUE) + + MSstatsBig:::cleanSpectronautChunk(input, output_file, pos = 1L, + filter_by_qvalue = TRUE, + qvalue_cutoff = 0.01) + result <- readr::read_csv(output_file, show_col_types = FALSE) + result <- result[order(result$PeptideSequence), ] + + # PEP1: EGQvalue below cutoff -> kept + expect_equal(result$Intensity[1], 100) + # PEP2: EGQvalue above cutoff -> NA + expect_true(is.na(result$Intensity[2])) + # PEP3: EGQvalue NA -> NA (preserves dplyr::if_else semantics) + expect_true(is.na(result$Intensity[3])) +}) + +test_that("cleanSpectronautChunk drops rows where FFrgLossType != noloss", { + input <- make_spectronaut_input( + n = 3L, + F.FrgLossType = c("noloss", "H2O", "noloss") + ) + output_file <- tempfile(fileext = ".csv") + on.exit(unlink(output_file, force = TRUE), add = TRUE) + + MSstatsBig:::cleanSpectronautChunk(input, output_file, pos = 1L, + filter_by_qvalue = FALSE) + result <- readr::read_csv(output_file, show_col_types = FALSE) + + expect_equal(nrow(result), 2L) + expect_setequal(result$PeptideSequence, c("PEP1", "PEP3")) +}) + +test_that("cleanSpectronautChunk fails fast with a clear error when a filter column is absent", { + output_file <- tempfile(fileext = ".csv") + on.exit(unlink(output_file, force = TRUE), add = TRUE) + + # FFrgLossType filter always runs: dropping the source column must error, + # naming the original Spectronaut column, not raise a cryptic data.table error. + input <- make_spectronaut_input(n = 1L) + input$F.FrgLossType <- NULL + expect_error( + MSstatsBig:::cleanSpectronautChunk(input, output_file, pos = 1L, + filter_by_qvalue = FALSE), + regexp = "F.FrgLossType") + + # q-value filter only errors when it is actually requested. + input_noq <- make_spectronaut_input(n = 1L) + input_noq$EG.Qvalue <- NULL + expect_error( + MSstatsBig:::cleanSpectronautChunk(input_noq, output_file, pos = 1L, + filter_by_qvalue = TRUE), + regexp = "EG.Qvalue") + + # ...and is tolerated when the q-value filter is switched off. + expect_error( + MSstatsBig:::cleanSpectronautChunk(make_spectronaut_input(n = 1L), + output_file, pos = 1L, + filter_by_qvalue = FALSE), + NA) +}) + +test_that("reduceBigSpectronaut projects out columns cleanSpectronautChunk does not use", { + # Real Spectronaut columns plus two junk columns that are not in needed_cols. + input <- make_spectronaut_input(n = 2L) + input$JunkA <- "zzz" + input$JunkB <- 999L + input_file <- tempfile(fileext = ".csv") + output_file <- tempfile() + on.exit({ + unlink(input_file, force = TRUE) + unlink(output_file, recursive = TRUE, force = TRUE) + }, add = TRUE) + readr::write_csv(input, input_file) + + # Capture the columns of each batch actually handed to the cleaner. The + # Scanner projection should have dropped JunkA/JunkB before this point. + captured <- new.env(parent = emptyenv()) + captured$cols <- NULL + stub(reduceBigSpectronaut, "cleanSpectronautChunk", + function(input, ...) { + captured$cols <- colnames(input) + NULL + }) + + reduceBigSpectronaut(input_file, output_file) + + expect_false("JunkA" %in% captured$cols) + expect_false("JunkB" %in% captured$cols) + # Needed columns still make it through. + expect_true("F.NormalizedPeakArea" %in% captured$cols) + expect_true("PG.ProteinAccessions" %in% captured$cols) +}) + +test_that("reduceBigSpectronaut rejects invalid block_size values", { + input_file <- tempfile(fileext = ".csv") + writeLines("a,b\n1,2", input_file) + output_file <- tempfile() + on.exit({ + unlink(input_file, force = TRUE) + unlink(output_file, recursive = TRUE, force = TRUE) + }, add = TRUE) + + expect_error(reduceBigSpectronaut(input_file, output_file, block_size = -1L), + regexp = "block_size") + expect_error(reduceBigSpectronaut(input_file, output_file, block_size = 0L), + regexp = "block_size") + expect_error(reduceBigSpectronaut(input_file, output_file, block_size = NA_integer_), + regexp = "block_size") + expect_error(reduceBigSpectronaut(input_file, output_file, block_size = c(1L, 2L)), + regexp = "block_size") + expect_error(suppressWarnings( + reduceBigSpectronaut(input_file, output_file, block_size = "16MB") + ), regexp = "block_size") +}) + +test_that("bigSpectronauttoMSstatsFormat plumbs block_size through to reduceBigSpectronaut", { + captured <- new.env(parent = emptyenv()) + captured$block_size <- NULL + + spy_reduce <- function(input_file, output_path, intensity, filter_by_excluded, + filter_by_identified, filter_by_qvalue, qvalue_cutoff, + calculateAnomalyScores, anomalyModelFeatures, + block_size = 16L * 1024L * 1024L) { + captured$block_size <- block_size + msstats_data <- data.frame( + ProteinName = "P1", PeptideSequence = "PEPTIDE", PrecursorCharge = 2, + FragmentIon = "frag1", ProductCharge = 1, + IsotopeLabelType = "L", Condition = "A", BioReplicate = 1, + Run = "run1", Intensity = 100 + ) + readr::write_csv(msstats_data, output_path) + } + + input_file <- "dummy_spectro_input.csv" + + # Default forwards 16 MiB. + stub(bigSpectronauttoMSstatsFormat, "reduceBigSpectronaut", spy_reduce) + output_file_default <- tempfile(fileext = ".csv") + on.exit({ + unlink(output_file_default, recursive = TRUE, force = TRUE) + unlink(paste0("reduce_output_", basename(output_file_default)), + recursive = TRUE, force = TRUE) + }, add = TRUE) + bigSpectronauttoMSstatsFormat( + input_file = input_file, + output_file_name = output_file_default, + backend = "arrow" + ) + expect_identical(captured$block_size, 16L * 1024L * 1024L) + + # Override forwards the user's value. + output_file_override <- tempfile(fileext = ".csv") + on.exit({ + unlink(output_file_override, recursive = TRUE, force = TRUE) + unlink(paste0("reduce_output_", basename(output_file_override)), + recursive = TRUE, force = TRUE) + }, add = TRUE) + bigSpectronauttoMSstatsFormat( + input_file = input_file, + output_file_name = output_file_override, + backend = "arrow", + block_size = 8L * 1024L * 1024L + ) + expect_identical(captured$block_size, 8L * 1024L * 1024L) +}) + # test_that("bigDIANNtoMSstatsFormat works with real MSstatsConvert tinytest data", { # input_file <- "/Users/rudhikshah/NorthEasternContractWork/MSstatsConvert/inst/tinytest/raw_data/DIANN/diann_input.tsv" # annotation_file <- "/Users/rudhikshah/NorthEasternContractWork/MSstatsConvert/inst/tinytest/raw_data/DIANN/annotation.csv"