M sstats big/work/20260513 arrow chunking followups

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,

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)

R/clean_spectronaut.R

@@ -1,33 +1,115 @@
+#' @importFrom data.table := .SD setDT setnames
+NULL
+
 #' @keywords internal
 reduceBigSpectronaut <- function(input_file, output_path,
                                  intensity="F.NormalizedPeakArea",
                                  filter_by_excluded = FALSE,
                                  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
 }

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,