Merge pull request #170 from Alamar-Biosciences/release/v1.4.2

Release v1.4.2

Author: jcbeer

.gitignore

@@ -44,10 +44,14 @@ vignettes/*.pdf
 
 # Example HTML file generated by skeleton.Rmd
 inst/rmarkdown/templates/nulisaseq/skeleton/skeleton.html
+inst/rmarkdown/templates/nulisaseq/skeleton/outputFiles
 inst/doc
 /doc/
 /Meta/
 
+# outputFiles generated by tests
+tests/testthat/fixtures/outputFiles
+
 # Book build artifacts (don't commit these)
 vignettes/book/_book/
 vignettes/book/_bookdown_files/
@@ -61,3 +65,5 @@ vignettes/book/_main.Rmd
 CLAUDE.md
 *.bak
 docs/
+.claude/
+

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: NULISAseqR
 Title: Analysis of NULISAseq Data
-Version: 1.4.1
+Version: 1.4.2
 Authors@R: 
     person("Alamar Biosciences Bioinformatics Team", , , "bioinfo@alamarbio.com", role = c("aut", "cre", "cph"),
            comment = "Contributors: Dwight Kuo, Joanne C. Beer, Eliza Chai, Sumedh Sankhe, Kasun Buddika")
@@ -62,7 +62,8 @@ Suggests:
     renv,
     rmarkdown,
     shiny,
-    testthat (>= 3.0.0)
+    testthat (>= 3.0.0),
+    withr
 VignetteBuilder: 
     knitr
 Config/testthat/edition: 3

NEWS.md

@@ -1,3 +1,37 @@
+# NULISAseqR 1.4.2 (2026-02-15)
+
+## Changes
+
+### Enhancements
+* **render_QC_report()** - Improved function parameter ordering and defaults:
+  - `xml_files` parameter moved to first position for more intuitive usage
+  - Added default values for `output_filename` ("NULISAseq_QC_Report.html") and `output_dir` (current working directory)
+  - Added default value for `dataDir` (current working directory)
+  - Simplified `Rmd_input_file` path construction using `system.file()`
+* **lod()** - Enhanced documentation and parameter handling:
+  - Improved parameter ordering (moved `data_matrix` before `blanks`)
+  - Enhanced roxygen documentation with clearer return value descriptions
+  - Added filtering to ensure `targetNoOutlierDetection` only includes targets present in `data_matrix`
+
+### Bug Fixes
+* **quantifiability()** - Fixed sample subsetting issue that could cause errors when sample lists don't match between AQ data and sample information:
+  - Now uses `intersect()` to find common samples between `Data_AQ_aM` and `SampleNames`
+  - Correctly calculates sample counts for overall and subgroup quantifiability 
+  - Prevents errors when processing data with mismatched sample lists
+* **loadNULISAseq()** - Added calculation of `LOD_pgmL` (limit of detection in pg/mL units) from XML data for AQ assays
+* **targetBoxplot()** - Fixed parameter naming in `lod()` function call to use `data_matrix=` explicitly
+
+### Testing
+* **New comprehensive test suites** added to ensure code quality and reliability:
+  - `test-importNULISAseq.R` - Tests for `importNULISAseq()` function with and without NULISAseqAQ package, including fallback mode validation and AQ data consistency checks
+  - `test-reverse-curve.R` - Tests for reverse curve target handling, including correlation validation, data transformation verification, and NPQ value consistency between `loadNULISAseq()` and `importNULISAseq()`
+  - `test-writeNULISAseq.R` - Tests for Excel output generation with both RQ-only and AQ data, including validation of sheet structure, column names, and specific data values
+* **Test infrastructure improvements**:
+  - Moved test fixtures from `inst/rmarkdown/templates/nulisaseq/skeleton/` to `tests/testthat/fixtures/` for better organization
+  - Removed unnecessary `.gitignore` file from skeleton template directory
+
+---
+
 # NULISAseqR 1.4.1 (2026-01-16)
 
 ## Changes

R/lod.R

@@ -22,13 +22,14 @@ outliers_index_mad <- function(col, min_blanks = 4, threshold = 2.5) {
 #'
 #' Calculates limit of detection (LoD) for each target based on the
 #' negative controls (blanks). LoD = mean(blanks) + 3*SD(blanks).
+#' LoD is typically calculated on IC-IPC normalized reads (unlogged).
 #' Designates data as either above or below LoD. 
 #' Option to specify minimum count threshold for detectability.
 #'
-#' @param blanks Column indices or column names of the blanks in the
-#' data_matrix.
 #' @param data_matrix The Data matrix output from readNULISAseq.R
 #' or normalized data from normalization functions.
+#' @param blanks Column indices or column names of the blanks in the
+#' data_matrix.
 #' @param min_count Optional count threshold to apply in addition
 #' to the LoD. Default is 0.
 #' @param min_blank_no Optional numeric parameter defining the minimum number of
@@ -45,14 +46,21 @@ outliers_index_mad <- function(col, min_blanks = 4, threshold = 2.5) {
 #' Lists samples/targets that should not be reported 
 #'
 #'
-#' @return A list.
-#' @param LOD Vector of limits of detection.
-#' @param aboveLOD Logical matrix indicating whether counts are 
-#' above or below LoD for that target.
+#' @return A list containing:
+#' \item{LOD}{Vector of limits of detection.}
+#' \item{aboveLOD}{Logical matrix indicating whether counts are 
+#' above or below LoD for that target.}
 #'
 #' @export
 #' 
-lod <- function(data_matrix, blanks, min_count = 0, min_blank_no = 4, mad_threshold = 2.5, ignore_target_blank = NULL, targetNoOutlierDetection = NULL, match_matrix = NULL) {
+lod <- function(data_matrix, 
+                blanks, 
+                min_count = 0, 
+                min_blank_no = 4, 
+                mad_threshold = 2.5, 
+                ignore_target_blank = NULL, 
+                targetNoOutlierDetection = NULL, 
+                match_matrix = NULL) {
   # Determine blank names if blank indices are provided
   if (is.numeric(blanks)) {
     blank_names <- colnames(data_matrix)[blanks]
@@ -120,8 +128,10 @@ lod <- function(data_matrix, blanks, min_count = 0, min_blank_no = 4, mad_thresh
 
   # Calculate the blank_mean and blank_sd assuming no outlier detection (blank_meanCount, blank_sdCount)
   # replace blank_mean and blank_sd with blank_meanCount / blank_sdCount for targets where no outlier 
-  # deteection is desired
+  # detection is desired
   if(!is.null(targetNoOutlierDetection)){
+    # filter targetNoOutlierDetection to include only targets in the data_matrix
+    targetNoOutlierDetection <- targetNoOutlierDetection[targetNoOutlierDetection %in% rownames(data_matrix)]
     blank_meanCount <- rowMeans(blank_data, na.rm = TRUE )
     blank_sdCount <- apply(blank_data, 1, sd, na.rm = TRUE)
     blank_mean[targetNoOutlierDetection] <- blank_meanCount[targetNoOutlierDetection]

R/mergeNULISAseq.R

@@ -108,7 +108,7 @@ process_loadNULISAseq <- function(data) {
         processed_data$aqParams <- processed_data$AQ$targetAQ_param %>%
           tibble::as_tibble() %>%
           select(!any_of(c("Encrypted", "Decrypted")))
-        
+
         logger::log_info("aqParams found, creating pg/mL matrix")
         processed_data$Data_AQ_pgmL <- processed_data[["Data_AQ"]] %>%
           tibble::as_tibble(rownames = "targetName") %>%
@@ -122,7 +122,7 @@ process_loadNULISAseq <- function(data) {
             names_to = "sampleName") %>%
           mutate(
             raw = dplyr::if_else(is.na(raw),
-                          raw, 
+                          raw,
                           NULISAseqAQ::unit_convert_am_conc(raw, MW_kDa))) %>%
           select(-MW_kDa) %>%
           tidyr::pivot_wider(
@@ -131,13 +131,27 @@ process_loadNULISAseq <- function(data) {
             values_from = "raw") %>%
           column_to_rownames("targetName") %>%
           as.matrix()
-        
+
         logger::log_info("Data_AQlog2_pg/mL matrix created")
-        
+
         processed_data$Data_AQlog2_pgmL <- log2(processed_data$Data_AQ_pgmL)
-        processed_data$AQ_unit <- "pg/mL" 
+        processed_data$AQ_unit <- "pg/mL"
       } else{
-        logger::log_info("Molecular weights in kDa not available, only aM units available")
+        logger::log_info("Molecular weights in kDa not available, checking for pre-existing pg/mL data")
+        # Still create aqParams from targetAQ_param even without MW_kDa
+        processed_data$aqParams <- processed_data$AQ$targetAQ_param %>%
+        tibble::as_tibble()
+
+        # Check if pre-existing Data_AQ (pg/mL) is available from fallback XML
+        if(!is.null(processed_data$AQ[["Data_AQ"]])){
+          logger::log_info("Using pre-existing pg/mL matrix from input data")
+          processed_data$Data_AQ_pgmL <- processed_data$AQ$Data_AQ %>%
+            rename_cols(., names_df = names_df)
+          processed_data$Data_AQlog2_pgmL <- log2(processed_data$Data_AQ_pgmL)
+          processed_data$AQ_unit <- "pg/mL"
+        } else {
+          logger::log_info("No pre-existing pg/mL matrix found; only aM units will be available")
+        }
       }
     }
   }
@@ -304,23 +318,57 @@ process_loadNULISAseq <- function(data) {
   processed_data["normed_untransformedReverse"] <- NULL # now called Data_Reverse
   processed_data["AQ"] <- NULL
 
-  matrices <-  c("Data_IC","Data_IClog2","Data_raw","Data_rawlog2","aboveLOD", "Data_AQ", "Data_AQlog2", "withinDR", "Data_Reverse", "Data_Reverselog2","Data_AQ_pgmL","Data_AQlog2_pgmL")
+  matrices <-  c("Data_IC","Data_IClog2","Data_raw","Data_rawlog2","aboveLOD", "Data_AQ", "Data_AQlog2", "Data_Reverse", "Data_Reverselog2","Data_AQ_pgmL","Data_AQlog2_pgmL")
   matrices <- matrices[matrices %in% names(processed_data)]
-  
+
+  # Track dropped samples for reporting
+  dropped_samples_info <- list()
+
   for (i in matrices) {
-    logger::log_info("Removing samples with all NaN or NA values -- ",i)
-    missing <- apply(processed_data[[i]], 2, function(x) all(is.nan(x))  || all(is.na(x)))
-    
-    if(i %in% c("Data_AQ", "Data_AQlog2", "Data_AQ_pgmL","Data_AQlog2_pgmL", "withinDR")){
+    logger::log_info("Removing samples with all NaN or NA values -- ", i)
+    missing <- apply(processed_data[[i]], 2, function(x) all(is.nan(x)) || all(is.na(x)))
+
+    if (sum(missing) > 0) {
+      dropped_sample_names <- colnames(processed_data[[i]])[missing]
+      dropped_samples_info[[i]] <- dropped_sample_names
+    }
+
+    if(i %in% c("Data_AQ", "Data_AQlog2", "Data_AQ_pgmL", "Data_AQlog2_pgmL")){
       targets <- processed_data[["aqParams"]] %>%
         pull(targetName)
     } else{
       targets <- processed_data[["targets"]] %>%
-        pull(targetName) 
+        pull(targetName)
     }
-    
+
     processed_data[[i]] <- processed_data[[i]][targets, !missing]
   }
+
+  # Ensure all AQ matrices have consistent rows and columns, then align withinDR
+  aq_matrices <- c("Data_AQ", "Data_AQlog2", "Data_AQ_pgmL", "Data_AQlog2_pgmL")
+  aq_matrices <- aq_matrices[aq_matrices %in% names(processed_data)]
+
+  if (length(aq_matrices) >= 1) {
+    common_rows <- Reduce(union, lapply(aq_matrices, function(m) rownames(processed_data[[m]])))
+    common_cols <- Reduce(union, lapply(aq_matrices, function(m) colnames(processed_data[[m]])))
+  
+    # Align all AQ matrices to common rows and columns
+    common_rows <- sort(common_rows)
+    for (m in aq_matrices) {
+      processed_data[[m]] <- processed_data[[m]][common_rows, common_cols, drop = FALSE]
+    }
+
+    logger::log_info("AQ matrices aligned to ", length(common_rows), " targets and ", length(common_cols), " samples")
+
+    # Align withinDR to match AQ matrices
+    if ("withinDR" %in% names(processed_data)) {
+      processed_data$withinDR <- processed_data$withinDR[common_rows, common_cols, drop = FALSE]
+      logger::log_info("withinDR aligned to AQ matrices")
+    }
+  }
+
+  # Store dropped samples info for upstream reporting
+  processed_data$droppedSamples <- dropped_samples_info
 
   processed_data[["samples"]] <- processed_data[["samples"]] %>%
     filter(sampleName %in% colnames(processed_data[["Data_IC"]])) %>%
@@ -686,6 +734,61 @@ mergeNULISAseq <- function(dataList, fileNameList, sample_group_covar = "SAMPLE_
     }
   }
 
+  # Collect and report dropped samples from all plates
+  all_dropped_samples <- list()
+  for (plate in names(dataList)) {
+    if (!is.null(dataList[[plate]]$droppedSamples) && length(dataList[[plate]]$droppedSamples) > 0) {
+      for (matrix_name in names(dataList[[plate]]$droppedSamples)) {
+        dropped <- dataList[[plate]]$droppedSamples[[matrix_name]]
+        if (length(dropped) > 0) {
+          all_dropped_samples[[length(all_dropped_samples) + 1]] <- data.frame(
+            plateID = plate,
+            matrix = matrix_name,
+            sampleName = dropped,
+            stringsAsFactors = FALSE
+          )
+        }
+      }
+    }
+  }
+
+  # Emit warning if any samples were dropped
+  if (length(all_dropped_samples) > 0) {
+    dropped_samples_df <- do.call(rbind, all_dropped_samples)
+
+    # Create summary warning message
+    warning_summary <- dropped_samples_df %>%
+      dplyr::group_by(plateID, matrix) %>%
+      dplyr::summarize(
+        n = dplyr::n(),
+        samples = paste(sampleName, collapse = ", "),
+        .groups = "drop"
+      )
+
+    warning_messages <- sapply(seq_len(nrow(warning_summary)), function(i) {
+      row <- warning_summary[i, ]
+      sprintf("Plate: '%s', Data matrix: '%s': %d sample(s) dropped (%s)",
+              row$plateID, row$matrix, row$n, row$samples)
+    })
+
+    # Identify unique samples dropped from Data_AQ matrices (affects quantifiability)
+    aq_dropped <- dropped_samples_df %>%
+      dplyr::filter(grepl("^Data_AQ", matrix)) %>%
+      dplyr::pull(sampleName) %>%
+      unique()
+
+    quant_message <- ""
+    if (length(aq_dropped) > 0) {
+      quant_message <- sprintf("\nSamples not included in quantifiability calculation: %s",
+                               paste(aq_dropped, collapse = ", "))
+    }
+
+    warning(paste0("Samples with all NaN or NA values were removed:\n",
+                   paste("  -", warning_messages, collapse = "\n"),
+                   quant_message),
+            call. = FALSE)
+  }
+
   # Extract control sample names and SampleNames from merged samples dataframe
   IPC_samples <- samples$sampleName[samples$sampleType == 'IPC']
   SC_samples <- samples$sampleName[samples$sampleType == 'SC']
@@ -728,7 +831,7 @@ mergeNULISAseq <- function(dataList, fileNameList, sample_group_covar = "SAMPLE_
   }
   # Add dataMatrix and unit
   return_list <- c(return_list, dataMatrix, unit = unit)
-  
+
   # return the output
   return(return_list)
 }
@@ -843,9 +946,7 @@ mergeNULISAseq <- function(dataList, fileNameList, sample_group_covar = "SAMPLE_
 #'           \item{\code{inconsistent_targets}: Placeholder for targets inconsistent across plates/runs (NULL if none)}
 #'           \item{\code{detectability}: Data frame of detectability by target and sample matrix}
 #'           \item{\code{quantifiability}: Data frame of quantifiability by target and sample matrix}
-#'           \item{\code{Data_raw}, \code{Data_rawlog2}: Raw counts and log2-transformed counts (matrix, targets × samples)}
-#'           \item{\code{Data_IC}, \code{Data_IClog2}: Internal control–normalized data (linear and log2)}
-#'           \item{\code{Data_Reverse}, \code{Data_Reverselog2}: Reverse-transformed IC-IPC normalized data (linear and log2)}
+#'           \item{\code{Data_raw}: Raw counts (matrix, targets × samples)}
 #'           \item{\code{Data_AQ_aM}, \code{Data_AQlog2_aM}: Absolute quantitation in attomolar units (linear and log2), 
 #'           if AQ data available}
 #'           \item{\code{Data_AQ_pgmL}, \code{Data_AQlog2_pgmL}: Absolute quantitation in pg/mL units (linear and log2), 
@@ -2176,7 +2277,7 @@ filter_target_qc_by_mode <- function(data, AQ = FALSE) {
 #'     \item{\code{targets} - Target metadata data frame}
 #'     \item{\code{samples} - Sample metadata data frame}
 #'     \item{\code{ExecutionDetails} - Execution details list}
-#'     \item{\code{Data_Reverselog2} or \code{Data_NPQ} - NPQ data matrix for RQ}
+#'     \item{\code{Data_IClog2} or \code{Data_NPQ} - NPQ data matrix for RQ}
 #'     \item{\code{Data_raw} - Raw count matrix for RQ}
 #'     \item{\code{Data_AQ} or \code{Data_AQ_aM} - Absolute quantification matrix in aM (optional)}
 #'     \item{\code{Data_AQ_pgmL} - Absolute quantification matrix in pg/mL (optional)}
@@ -2258,12 +2359,12 @@ format_wide_to_long <- function(merged, AQ = FALSE, exclude_sample_cols = "plate
 
   } else {
     # Check which NPQ data column is available
-    NPQ_data_used <- if("Data_Reverselog2" %in% names(merged)) {
-      "Data_Reverselog2"
+    NPQ_data_used <- if("Data_IClog2" %in% names(merged)) {
+      "Data_IClog2"
     } else if("Data_NPQ" %in% names(merged)) {
-      "Data_NPQ" 
+      "Data_NPQ"
     } else {
-      stop("Neither 'Data_Reverselog2' nor 'Data_NPQ' found in the merged data")
+      stop("Neither 'Data_IClog2' nor 'Data_NPQ' found in the merged data")
     }
 
     data_long <- convert_to_long(data_matrix = merged[[NPQ_data_used]], data_col ="NPQ")

R/quantifiability.R

@@ -143,7 +143,8 @@ quantifiability <- function(runs,
         DR$LLOQ[DR$LOD_aM > DR$LLOQ & !is.na(DR$LOD_aM) & !is.na(DR$LLOQ)] <- DR$LOD_aM[DR$LOD_aM > DR$LLOQ & !is.na(DR$LOD_aM) & !is.na(DR$LLOQ)]
 
         # calculate overall quantifiability
-        AQ_quant <- x$AQ$Data_AQ_aM[,x$SampleNames, drop=FALSE]
+        intersect_samples <- intersect(colnames(x$AQ$Data_AQ_aM), x$SampleNames)
+        AQ_quant <- x$AQ$Data_AQ_aM[,intersect_samples, drop=FALSE] 
         AQ_quant <- merge(AQ_quant, DR, by.x='row.names', by.y='targetName')
         rownames(AQ_quant) <- AQ_quant[,1]
         AQ_quant <- AQ_quant[,2:ncol(AQ_quant)]
@@ -156,15 +157,17 @@ quantifiability <- function(runs,
         })
 
         AQ_quant_output_columns <- c('overall')
-        n_samples <- c(overall=length(x$SampleNames))
+        n_samples <- c(overall=length(intersect_samples)) 
 
         # calculate subgroup quantifiability
         if(!is.null(sampleGroupCovar)){
           if(sampleGroupCovar %in% colnames(x$samples)){
             subgroup_names <- unique(x$samples[x$samples$sampleType=='Sample', sampleGroupCovar])
 
             for(i in 1:length(subgroup_names)){
-              subgroup_samples <- x$samples$sampleName[x$samples[,sampleGroupCovar]==subgroup_names[i]]
+              # Filter samples to those in intersect_samples AND belonging to this subgroup
+              samples_in_subgroup <- x$samples$sampleName[x$samples[,sampleGroupCovar] == subgroup_names[i]]
+              subgroup_samples <- intersect(samples_in_subgroup, intersect_samples)
               subgroup_sample_data <- AQ_quant[,subgroup_samples, drop=FALSE]
               AQ_quant[,as.character(subgroup_names[i])] <- NA
               AQ_quant_output_columns <- c(AQ_quant_output_columns, as.character(subgroup_names[i]))