diff --git a/DESCRIPTION b/DESCRIPTION
index 8e031ec..4d2a743 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: splikit
 Title: Analysing RNA Splicing in Single-Cell RNA Sequencing Data
-Version: 2.3.2
+Version: 2.3.3
 Authors@R:
     person("Arsham", "Mikaeili Namini", , "arsham.mikaeilinamini@mail.mcgill.ca", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0002-9453-6951"))
diff --git a/NAMESPACE b/NAMESPACE
index f42032c..cfac55d 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,6 +1,7 @@
 # Generated by roxygen2: do not edit by hand
 
 S3method(print,splikit_junction_plot)
+S3method(print,splikit_pseudo_correlation_result)
 export(SplikitObject)
 export(find_variable_events)
 export(find_variable_genes)
diff --git a/NEWS.md b/NEWS.md
index 1ed7220..d112d4e 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,27 @@
+# splikit 2.3.3
+
+## Gene-expression input selection
+
+* `make_gene_count()` now separates count-matrix selection from barcode
+  filtering. `matrix_source` selects the `raw` or `filtered` directory, while
+  `matrix_file` supports either the standard `matrix.mtx` or an alternative
+  Matrix Market file. By default, the function reads from `raw/`, prefers
+  STARsolo's `UniqueAndMult-EM.mtx` when available, and falls back to
+  `raw/matrix.mtx`. Explicit arguments can still force any supported source or
+  filename.
+
+## Pseudo-correlation result export
+
+* `get_pseudo_correlation()` and `SplikitObject$getPseudoCorrelation()` now
+  return a structured `splikit_pseudo_correlation_result`. Its `statistics`
+  component retains the existing per-event summaries and empirical inference,
+  while `null_distribution` contains one exportable row for every retained
+  event and permutation. The long table preserves the event and permutation
+  identifiers and replaces the previous `null_draws` matrix attribute.
+* The pooled null values are provided for descriptive analysis and export;
+  event-level empirical p-values continue to use each event's own permutation
+  draws.
+
 # splikit 2.3.2
 
 ## New Features
@@ -199,4 +223,4 @@ library(splikit)
 * Further performance optimizations
 
 ---
-*For questions or issues, please visit: https://github.com/Arshammik/splikit/issues*
\ No newline at end of file
+*For questions or issues, please visit: https://github.com/Arshammik/splikit/issues*
diff --git a/R/SplikitObject.R b/R/SplikitObject.R
index e8d7b97..af2de12 100644
--- a/R/SplikitObject.R
+++ b/R/SplikitObject.R
@@ -227,9 +227,8 @@ SplikitObject <- R6::R6Class("SplikitObject",
     #'   permutations for reproducibility without disturbing the global RNG
     #'   (default: NULL).
     #'
-    #' @return A data.table with event names, pseudo_correlation, null_distribution,
-    #'   null_sd, n_perm_valid, emp_pvalue, and emp_padj, plus the full null draws
-    #'   in attr(result, "null_draws").
+    #' @return A `splikit_pseudo_correlation_result` with per-event statistics,
+    #'   a long event-by-permutation null distribution, and computation metadata.
     getPseudoCorrelation = function(ZDB_matrix, metric = "CoxSnell",
                                      suppress_warnings = TRUE,
                                      permutation_count = 100L,
@@ -266,7 +265,7 @@ SplikitObject <- R6::R6Class("SplikitObject",
 
       # Warn about NA removal
       n_before <- nrow(self$m1)
-      n_after <- nrow(result)
+      n_after <- nrow(result$statistics)
       if (n_before > n_after) {
         n_removed <- n_before - n_after
         message("Removed ", n_removed, " event(s) with NA values (",
diff --git a/R/general_tools.R b/R/general_tools.R
index b1bca5c..3482ccf 100644
--- a/R/general_tools.R
+++ b/R/general_tools.R
@@ -33,7 +33,18 @@ utils::globalVariables(c("unique_mapped", "i", "j", "x_1", "x_tot", "x_2", "grou
 #'   caller's global RNG stream is saved and restored, so passing a seed does not
 #'   perturb downstream randomness. Defaults to `NULL` (use the ambient RNG state).
 #'
-#' @return A `data.table` with the following columns:
+#' @return An S3 object of class `splikit_pseudo_correlation_result` with three
+#' components:
+#' \describe{
+#'   \item{statistics}{A per-event `data.table` with the columns described below.}
+#'   \item{null_distribution}{A long `data.table` containing one row per retained
+#'   event and permutation, with columns `event`, `permutation`, and
+#'   `null_pseudo_correlation`. This table is suitable for direct export or pooled
+#'   descriptive analyses; event-level p-values remain based on each event's own
+#'   permutation draws.}
+#'   \item{metadata}{A list recording the permutation settings and null-draw counts.}
+#' }
+#' The `statistics` component contains:
 #' \describe{
 #'   \item{event}{The event names from `ZDB_matrix` rownames.}
 #'   \item{pseudo_correlation}{The computed pseudo R² correlation values using the specified metric.}
@@ -43,8 +54,6 @@ utils::globalVariables(c("unique_mapped", "i", "j", "x_1", "x_tot", "x_2", "grou
 #'   \item{emp_pvalue}{Two-sided empirical p-value, `(b + 1) / (n_perm_valid + 1)`, where `b` is the number of valid null draws whose absolute value is greater than or equal to the absolute observed correlation.}
 #'   \item{emp_padj}{Benjamini-Hochberg adjusted `emp_pvalue` across the retained events.}
 #' }
-#' The full matrix of per-event null draws (events x `permutation_count`) is attached
-#' as `attr(result, "null_draws")` for custom downstream analyses.
 #'
 #' @examples
 #' \donttest{
@@ -83,9 +92,10 @@ utils::globalVariables(c("unique_mapped", "i", "j", "x_1", "x_tot", "x_2", "grou
 #'                                                  permutation_count = 20,
 #'                                                  permutation_seed = 1)
 #'
-#' # Per-event empirical p-values and the full null draws
-#' head(pseudo_r_square_sparse[, c("event", "pseudo_correlation", "emp_pvalue")])
-#' dim(attr(pseudo_r_square_sparse, "null_draws"))
+#' # Per-event empirical p-values and the exportable long null distribution
+#' head(pseudo_r_square_sparse$statistics[,
+#'   c("event", "pseudo_correlation", "emp_pvalue")])
+#' dim(pseudo_r_square_sparse$null_distribution)
 #'
 #' # Example using Nagelkerke R-squared instead of Cox-Snell
 #' pseudo_r_square_nagel <- get_pseudo_correlation(ZDB_matrix, m1_inclusion, m2_exclusion,
@@ -291,15 +301,65 @@ get_pseudo_correlation <- function(ZDB_matrix, m1_inclusion = NULL, m2_exclusion
   # Benjamini-Hochberg adjusted empirical p-values across the retained events.
   results$emp_padj <- stats::p.adjust(results$emp_pvalue, method = "BH")
 
-  # Attach the full per-event null draws for custom downstream analyses.
-  if (nrow(null_mat) > 0L) {
-    rownames(null_mat) <- results$event
-    colnames(null_mat) <- paste0("perm_", seq_len(permutation_count))
-  }
-  data.table::setattr(results, "null_draws", null_mat)
+  # Expose every retained event-by-permutation draw in a long, exportable table.
+  # R matrices are column-major, so this ordering groups all retained events for
+  # permutation 1, followed by all retained events for permutation 2, and so on.
+  null_distribution <- data.table::data.table(
+    event = rep(results$event, times = permutation_count),
+    permutation = rep(seq_len(permutation_count), each = nrow(results)),
+    null_pseudo_correlation = as.vector(null_mat)
+  )
+
+  metadata <- list(
+    metric = metric,
+    permutation_count = permutation_count,
+    permutation_seed = permutation_seed,
+    n_events_input = nrow(ZDB_matrix),
+    n_events_retained = nrow(results),
+    n_null_draws = nrow(null_distribution),
+    n_null_valid = sum(!is.na(null_distribution$null_pseudo_correlation)),
+    pooled_null_usage = "descriptive"
+  )
+
+  result <- structure(
+    list(
+      statistics = results,
+      null_distribution = null_distribution,
+      metadata = metadata
+    ),
+    class = c("splikit_pseudo_correlation_result", "list")
+  )
 
   if (verbose) message("Computation completed successfully.")
-  return(results)
+  return(result)
+}
+
+#' Print a pseudo-correlation result
+#'
+#' Prints a short preview of the per-event statistics and reports the size of
+#' the long null distribution without printing every permutation draw.
+#'
+#' @param x A `splikit_pseudo_correlation_result` object.
+#' @param n Number of statistic rows to preview. Defaults to `6L`.
+#' @param ... Additional arguments passed to the `data.table` print method.
+#'
+#' @return `x`, invisibly.
+#' @method print splikit_pseudo_correlation_result
+#' @export
+print.splikit_pseudo_correlation_result <- function(x, n = 6L, ...) {
+  if (!is.numeric(n) || length(n) != 1L || is.na(n) || n < 0) {
+    stop("n must be a single non-negative integer.", call. = FALSE)
+  }
+  n <- as.integer(n)
+
+  cat("Pseudo-correlation result\n")
+  cat("Statistics:", nrow(x$statistics), "retained event(s)\n")
+  print(utils::head(x$statistics, n), ...)
+  cat(
+    "Null distribution:", nrow(x$null_distribution), "row(s);",
+    x$metadata$n_null_valid, "valid draw(s)\n"
+  )
+  invisible(x)
 }
 
 #' Calculate Row-wise Variance for Dense or Sparse Matrices
diff --git a/R/star_solo_processing.R b/R/star_solo_processing.R
index 25c1845..4d03965 100644
--- a/R/star_solo_processing.R
+++ b/R/star_solo_processing.R
@@ -1122,8 +1122,10 @@ make_eventdata_plus <- function(eventdata, GTF_file_direction) {
 #' Constructs sparse gene expression matrices from one or more directories containing 10X Genomics-style output.
 #' The function supports barcode filtering using either an external whitelist or the internally provided filtered barcode file.
 #'
-#' @param expression_dirs A character vector or list of strings. Each element must be a path to a directory containing the
-#'   gene expression matrix files: \code{matrix.mtx}, \code{barcodes.tsv}, and \code{features.tsv} (or \code{genes.tsv}).
+#' @param expression_dirs A character vector or list of strings. Each element
+#'   must point to a gene-expression output directory containing \code{raw}
+#'   and/or \code{filtered} subdirectories with the matrix, barcodes, and
+#'   features files.
 #'
 #' @param sample_ids A character vector or list of unique sample identifiers, one for each element in \code{expression_dirs}.
 #'   These are used to name outputs in the returned list when multiple samples are provided.
@@ -1136,6 +1138,17 @@ make_eventdata_plus <- function(eventdata, GTF_file_direction) {
 #'   the function will attempt to use the default filtered barcode list from the input directory.
 #'   If \code{FALSE}, no internal filtration will be applied unless a whitelist is explicitly provided.
 #' @param verbose Logical. If \code{TRUE}, prints progress and informational messages. Default is \code{FALSE}.
+#' @param matrix_source Character. Directory below each \code{expression_dirs}
+#'   entry from which the count matrix, barcodes, and features are read. One of
+#'   \code{"raw"}, \code{"filtered"}, or \code{"auto"}. The default is
+#'   \code{"raw"}. The compatibility option \code{"auto"} preserves the
+#'   historical directory behavior: it uses
+#'   \code{"filtered"} when \code{use_internal_whitelist = TRUE} and
+#'   \code{"raw"} otherwise.
+#' @param matrix_file Character. Name of the Matrix Market count file inside the
+#'   selected source directory. The default, \code{"auto"}, uses STARsolo's
+#'   \code{"UniqueAndMult-EM.mtx"} when available and otherwise falls back to
+#'   \code{"matrix.mtx"}. Supply a filename explicitly to force that file.
 #'
 #' @return
 #' If a single sample is provided, returns a sparse matrix of class \code{"dgCMatrix"} with genes as rows and barcodes as columns.
@@ -1150,11 +1163,28 @@ make_eventdata_plus <- function(eventdata, GTF_file_direction) {
 #' If neither an external whitelist nor an internal filtered barcode file is available, all barcodes from the
 #' raw matrix will be retained.
 #'
+#' Matrix selection and barcode filtering are independent. For example,
+#' \code{matrix_source = "raw"} with \code{use_internal_whitelist = TRUE}
+#' reads the raw count matrix and then retains barcodes listed in
+#' \code{filtered/barcodes.tsv}. In contrast,
+#' \code{matrix_source = "filtered"} reads the already-filtered matrix directly.
+#' With the defaults, \code{raw/UniqueAndMult-EM.mtx} is preferred and
+#' \code{raw/matrix.mtx} is used as a fallback.
+#'
 #' @section Dependencies:
 #' Requires the \pkg{Matrix} package for sparse matrix handling and potentially \pkg{data.table} for efficient I/O.
 #'
 #' @export
-make_gene_count <- function(expression_dirs, sample_ids, whitelist_barcodes = NULL, use_internal_whitelist = TRUE, verbose = FALSE) {
+make_gene_count <- function(expression_dirs, sample_ids, whitelist_barcodes = NULL,
+                            use_internal_whitelist = TRUE, verbose = FALSE,
+                            matrix_source = c("raw", "filtered", "auto"),
+                            matrix_file = "auto") {
+
+  matrix_source <- match.arg(matrix_source)
+  if (!is.character(matrix_file) || length(matrix_file) != 1L ||
+      is.na(matrix_file) || !nzchar(matrix_file)) {
+    stop("matrix_file must be a single non-empty filename.", call. = FALSE)
+  }
 
   # Handle single sample input by converting to lists
   if (!is.list(expression_dirs)) expression_dirs <- list(expression_dirs)
@@ -1172,22 +1202,50 @@ make_gene_count <- function(expression_dirs, sample_ids, whitelist_barcodes = NU
 
   # Helper function to process one sample
   process_ex_sample <- function(expression_dir, sample_id, whitelist_barcode) {
-    # Determine directory (filtered or raw)
-    data_dir <- if (use_internal_whitelist) paste0(expression_dir, "/filtered") else paste0(expression_dir, "/raw")
+    # "auto" retains the pre-existing coupling for backward compatibility.
+    selected_source <- if (identical(matrix_source, "auto")) {
+      if (use_internal_whitelist) "filtered" else "raw"
+    } else {
+      matrix_source
+    }
+    data_dir <- file.path(expression_dir, selected_source)
+
+    selected_matrix_file <- matrix_file
+    if (identical(matrix_file, "auto")) {
+      candidate_files <- c("UniqueAndMult-EM.mtx", "matrix.mtx")
+      available <- candidate_files[
+        file.exists(file.path(data_dir, candidate_files))
+      ]
+      if (length(available) == 0L) {
+        stop(
+          "No expression matrix file found for sample ", sample_id,
+          " in ", data_dir, ". Tried: ",
+          paste(candidate_files, collapse = ", "),
+          call. = FALSE
+        )
+      }
+      selected_matrix_file <- available[[1L]]
+    }
 
     # Define paths
-    expression_matrix_dir <- paste0(data_dir, "/matrix.mtx")
-    expression_barcodes_dir <- paste0(data_dir, "/barcodes.tsv")
-    expression_features_dir <- paste0(data_dir, "/features.tsv")
-    filtered_barcodes_dir <- paste0(expression_dir, "/filtered/barcodes.tsv")  # For barcode filtration
+    expression_matrix_dir <- file.path(data_dir, selected_matrix_file)
+    expression_barcodes_dir <- file.path(data_dir, "barcodes.tsv")
+    expression_features_dir <- file.path(data_dir, "features.tsv")
+    filtered_barcodes_dir <- file.path(expression_dir, "filtered", "barcodes.tsv")
 
     # Check for required files
-    if (!file.exists(expression_matrix_dir)) stop("No expression matrix file found for sample: ", sample_id, call. = FALSE)
+    if (!file.exists(expression_matrix_dir)) {
+      stop("No expression matrix file found for sample ", sample_id,
+           ": ", expression_matrix_dir, call. = FALSE)
+    }
     if (!file.exists(expression_barcodes_dir)) stop("No barcodes file found for sample: ", sample_id, call. = FALSE)
     if (!file.exists(expression_features_dir)) stop("No features file found for sample: ", sample_id, call. = FALSE)
 
     # Read gene expression data
-    if (verbose) message("|-- Processing gene expression data for sample: ", sample_id)
+    if (verbose) {
+      message("|-- Processing gene expression data for sample: ", sample_id,
+              " (", selected_source, "/", selected_matrix_file, ")")
+    }
     g_mtx <- Matrix::readMM(expression_matrix_dir)
     g_brc <- data.table::fread(expression_barcodes_dir, header = FALSE, showProgress = FALSE)$V1
     g_feature <- data.table::fread(expression_features_dir, header = FALSE, showProgress = FALSE)
diff --git a/README.md b/README.md
index 3a50d40..a5457f2 100644
--- a/README.md
+++ b/README.md
@@ -6,7 +6,9 @@
 [![CRAN status](https://www.r-pkg.org/badges/version/splikit)](https://CRAN.R-project.org/package=splikit)
 [![R-CMD-check](https://github.com/csglab/splikit/actions/workflows/R-CMD-check.yml/badge.svg)](https://github.com/csglab/splikit/actions/workflows/R-CMD-check.yml)
 [![Documentation](https://img.shields.io/badge/Docs-Learn%20More-blue.svg)](https://csglab.github.io/splikit/)
-[![CRAN downloads](https://cranlogs.r-pkg.org/badges/grand-total/splikit)](https://cran.r-project.org/package=splikit)
+[![CRAN status](https://www.r-pkg.org/badges/version/splikit)](https://CRAN.R-project.org/package=splikit) [![GitHub version](https://img.shields.io/github/v/release/csglab/splikit?label=GitHub)](https://github.com/csglab/splikit/releases/latest)
+
+
 
 ## **Requirements**
 
@@ -100,4 +102,3 @@ This project is licensed under the MIT License – see the `LICENSE.md` file for
 
 The author is aware of the package’s limitations and potential breakpoints. It was developed under limited knowledge, time, and resources, and is provided with the hope that it will be useful. Feedback and contributions from the community are warmly welcomed. If you encounter any issues or have suggestions, please [open an issue](https://github.com/csglab/splikit/issues/new).
 
-
diff --git a/RELEASENOTES.md b/RELEASENOTES.md
index d4beb2e..5dbe3e7 100644
--- a/RELEASENOTES.md
+++ b/RELEASENOTES.md
@@ -6,6 +6,24 @@ canonical, CRAN-facing changelog see `NEWS.md`.
 
 ---
 
+## 2.3.3 — 2026-07-27
+
+**Summary — Exportable pseudo-correlation null results.**
+
+- Replaced the pseudo-correlation result's `null_draws` matrix attribute with a
+  structured S3 result containing `statistics`, a long `null_distribution`, and
+  computation `metadata`.
+- The long table has one row per retained event and permutation and preserves
+  both identifiers for reliable export and reconstruction.
+- Per-event empirical p-values and adjusted p-values are unchanged; the pooled
+  null column is descriptive and is not used as a global event-level null.
+- Added independent `matrix_source` and `matrix_file` controls to
+  `make_gene_count()`, including direct support for raw or filtered matrices
+  and STARsolo files such as `UniqueAndMult-EM.mtx`. The defaults now read
+  `raw/UniqueAndMult-EM.mtx` when available and fall back to `raw/matrix.mtx`.
+
+---
+
 ## 2.3.2 — 2026-06-09
 
 **Author:** Arsham Mikaeili Namini
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 91533a6..5524cf7 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -82,6 +82,7 @@ reference:
   desc: Additional analysis utilities
   contents:
   - get_pseudo_correlation
+  - print.splikit_pseudo_correlation_result
   - get_silhouette_mean
 
 - title: Data
diff --git a/man/SplikitObject.Rd b/man/SplikitObject.Rd
index f3993b2..ab3d3c5 100644
--- a/man/SplikitObject.Rd
+++ b/man/SplikitObject.Rd
@@ -228,9 +228,8 @@ permutations for reproducibility without disturbing the global RNG
 \if{html}{\out{</div>}}
 }
 \subsection{Returns}{
-A data.table with event names, pseudo_correlation, null_distribution,
-null_sd, n_perm_valid, emp_pvalue, and emp_padj, plus the full null draws
-in attr(result, "null_draws").
+A \code{splikit_pseudo_correlation_result} with per-event statistics,
+a long event-by-permutation null distribution, and computation metadata.
 }
 }
 \if{html}{\out{<hr>}}
diff --git a/man/get_pseudo_correlation.Rd b/man/get_pseudo_correlation.Rd
index 407a935..b34f18c 100644
--- a/man/get_pseudo_correlation.Rd
+++ b/man/get_pseudo_correlation.Rd
@@ -43,7 +43,18 @@ caller's global RNG stream is saved and restored, so passing a seed does not
 perturb downstream randomness. Defaults to \code{NULL} (use the ambient RNG state).}
 }
 \value{
-A \code{data.table} with the following columns:
+An S3 object of class \code{splikit_pseudo_correlation_result} with three
+components:
+\describe{
+\item{statistics}{A per-event \code{data.table} with the columns described below.}
+\item{null_distribution}{A long \code{data.table} containing one row per retained
+event and permutation, with columns \code{event}, \code{permutation}, and
+\code{null_pseudo_correlation}. This table is suitable for direct export or pooled
+descriptive analyses; event-level p-values remain based on each event's own
+permutation draws.}
+\item{metadata}{A list recording the permutation settings and null-draw counts.}
+}
+The \code{statistics} component contains:
 \describe{
 \item{event}{The event names from \code{ZDB_matrix} rownames.}
 \item{pseudo_correlation}{The computed pseudo R² correlation values using the specified metric.}
@@ -53,8 +64,6 @@ A \code{data.table} with the following columns:
 \item{emp_pvalue}{Two-sided empirical p-value, \code{(b + 1) / (n_perm_valid + 1)}, where \code{b} is the number of valid null draws whose absolute value is greater than or equal to the absolute observed correlation.}
 \item{emp_padj}{Benjamini-Hochberg adjusted \code{emp_pvalue} across the retained events.}
 }
-The full matrix of per-event null draws (events x \code{permutation_count}) is attached
-as \code{attr(result, "null_draws")} for custom downstream analyses.
 }
 \description{
 This function calculates a pseudo R²-like correlation metric using a beta-binomial model
@@ -99,9 +108,10 @@ pseudo_r_square_sparse <- get_pseudo_correlation(ZDB_matrix, m1_inclusion, m2_ex
                                                  permutation_count = 20,
                                                  permutation_seed = 1)
 
-# Per-event empirical p-values and the full null draws
-head(pseudo_r_square_sparse[, c("event", "pseudo_correlation", "emp_pvalue")])
-dim(attr(pseudo_r_square_sparse, "null_draws"))
+# Per-event empirical p-values and the exportable long null distribution
+head(pseudo_r_square_sparse$statistics[,
+  c("event", "pseudo_correlation", "emp_pvalue")])
+dim(pseudo_r_square_sparse$null_distribution)
 
 # Example using Nagelkerke R-squared instead of Cox-Snell
 pseudo_r_square_nagel <- get_pseudo_correlation(ZDB_matrix, m1_inclusion, m2_exclusion,
diff --git a/man/make_gene_count.Rd b/man/make_gene_count.Rd
index 476c1fb..ed85a08 100644
--- a/man/make_gene_count.Rd
+++ b/man/make_gene_count.Rd
@@ -9,12 +9,16 @@ make_gene_count(
   sample_ids,
   whitelist_barcodes = NULL,
   use_internal_whitelist = TRUE,
-  verbose = FALSE
+  verbose = FALSE,
+  matrix_source = c("raw", "filtered", "auto"),
+  matrix_file = "auto"
 )
 }
 \arguments{
-\item{expression_dirs}{A character vector or list of strings. Each element must be a path to a directory containing the
-gene expression matrix files: \code{matrix.mtx}, \code{barcodes.tsv}, and \code{features.tsv} (or \code{genes.tsv}).}
+\item{expression_dirs}{A character vector or list of strings. Each element
+must point to a gene-expression output directory containing \code{raw}
+and/or \code{filtered} subdirectories with the matrix, barcodes, and
+features files.}
 
 \item{sample_ids}{A character vector or list of unique sample identifiers, one for each element in \code{expression_dirs}.
 These are used to name outputs in the returned list when multiple samples are provided.}
@@ -28,6 +32,19 @@ the function will attempt to use the default filtered barcode list from the inpu
 If \code{FALSE}, no internal filtration will be applied unless a whitelist is explicitly provided.}
 
 \item{verbose}{Logical. If \code{TRUE}, prints progress and informational messages. Default is \code{FALSE}.}
+
+\item{matrix_source}{Character. Directory below each \code{expression_dirs}
+entry from which the count matrix, barcodes, and features are read. One of
+\code{"raw"}, \code{"filtered"}, or \code{"auto"}. The default is
+\code{"raw"}. The compatibility option \code{"auto"} preserves the
+historical directory behavior: it uses
+\code{"filtered"} when \code{use_internal_whitelist = TRUE} and
+\code{"raw"} otherwise.}
+
+\item{matrix_file}{Character. Name of the Matrix Market count file inside the
+selected source directory. The default, \code{"auto"}, uses STARsolo's
+\code{"UniqueAndMult-EM.mtx"} when available and otherwise falls back to
+\code{"matrix.mtx"}. Supply a filename explicitly to force that file.}
 }
 \value{
 If a single sample is provided, returns a sparse matrix of class \code{"dgCMatrix"} with genes as rows and barcodes as columns.
@@ -45,6 +62,14 @@ output by tools like CellRanger.
 
 If neither an external whitelist nor an internal filtered barcode file is available, all barcodes from the
 raw matrix will be retained.
+
+Matrix selection and barcode filtering are independent. For example,
+\code{matrix_source = "raw"} with \code{use_internal_whitelist = TRUE}
+reads the raw count matrix and then retains barcodes listed in
+\code{filtered/barcodes.tsv}. In contrast,
+\code{matrix_source = "filtered"} reads the already-filtered matrix directly.
+With the defaults, \code{raw/UniqueAndMult-EM.mtx} is preferred and
+\code{raw/matrix.mtx} is used as a fallback.
 }
 \section{Dependencies}{
 
diff --git a/man/print.splikit_pseudo_correlation_result.Rd b/man/print.splikit_pseudo_correlation_result.Rd
new file mode 100644
index 0000000..62b33bb
--- /dev/null
+++ b/man/print.splikit_pseudo_correlation_result.Rd
@@ -0,0 +1,22 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/general_tools.R
+\name{print.splikit_pseudo_correlation_result}
+\alias{print.splikit_pseudo_correlation_result}
+\title{Print a pseudo-correlation result}
+\usage{
+\method{print}{splikit_pseudo_correlation_result}(x, n = 6L, ...)
+}
+\arguments{
+\item{x}{A \code{splikit_pseudo_correlation_result} object.}
+
+\item{n}{Number of statistic rows to preview. Defaults to \code{6L}.}
+
+\item{...}{Additional arguments passed to the \code{data.table} print method.}
+}
+\value{
+\code{x}, invisibly.
+}
+\description{
+Prints a short preview of the per-event statistics and reports the size of
+the long null distribution without printing every permutation draw.
+}
diff --git a/tests/testthat/test-SplikitObject.R b/tests/testthat/test-SplikitObject.R
index d3a3adf..f9e8b61 100644
--- a/tests/testthat/test-SplikitObject.R
+++ b/tests/testthat/test-SplikitObject.R
@@ -662,10 +662,14 @@ test_that("SplikitObject$getPseudoCorrelation executes correctly", {
   # Keep permutation_count small here: this runs on the full 2000 x 2000 toy data.
   res <- obj$getPseudoCorrelation(ZDB_matrix = ZDB_matrix, suppress_warnings = TRUE,
                                   permutation_count = 3, permutation_seed = 1)
-  expect_true(data.table::is.data.table(res))
+  expect_s3_class(res, "splikit_pseudo_correlation_result")
+  expect_true(data.table::is.data.table(res$statistics))
+  expect_true(data.table::is.data.table(res$null_distribution))
   # With pure random data, some might be NA and removed, check columns
   expect_true(all(c("event", "pseudo_correlation", "null_distribution",
-                    "null_sd", "n_perm_valid", "emp_pvalue", "emp_padj") %in% names(res)))
-  expect_true(all(res$emp_pvalue > 0 & res$emp_pvalue <= 1))
+                    "null_sd", "n_perm_valid", "emp_pvalue", "emp_padj") %in%
+                  names(res$statistics)))
+  expect_true(all(res$statistics$emp_pvalue > 0 &
+                  res$statistics$emp_pvalue <= 1))
+  expect_equal(nrow(res$null_distribution), nrow(res$statistics) * 3L)
 })
-
diff --git a/tests/testthat/test-general_tools.R b/tests/testthat/test-general_tools.R
index 30222b1..0f95cd0 100644
--- a/tests/testthat/test-general_tools.R
+++ b/tests/testthat/test-general_tools.R
@@ -1,5 +1,5 @@
 test_that("get_pseudo_correlation works on valid data", {
-  n_events <- 5
+  n_events <- 10
   n_cells <- 100
   set.seed(42)
   
@@ -12,27 +12,36 @@ test_that("get_pseudo_correlation works on valid data", {
   
   # Test with default arguments
   res <- get_pseudo_correlation(ZDB, m1, m2)
-  expect_true(data.table::is.data.table(res))
-  expect_equal(nrow(res), n_events)
-  expect_true(all(c("event", "pseudo_correlation", "null_distribution") %in% names(res)))
-  expect_type(res$pseudo_correlation, "double")
+  expect_s3_class(res, "splikit_pseudo_correlation_result")
+  expect_named(res, c("statistics", "null_distribution", "metadata"))
+  expect_true(data.table::is.data.table(res$statistics))
+  expect_true(data.table::is.data.table(res$null_distribution))
+  expect_equal(nrow(res$statistics), n_events)
+  expect_true(all(c("event", "pseudo_correlation", "null_distribution") %in%
+                  names(res$statistics)))
+  expect_named(res$null_distribution,
+               c("event", "permutation", "null_pseudo_correlation"))
+  expect_type(res$statistics$pseudo_correlation, "double")
+  expect_output(print(res), "Null distribution: 1000 row\\(s\\)")
   
   # Test with Nagelkerke metric
   res_nagel <- get_pseudo_correlation(ZDB, m1, m2, metric = "Nagelkerke")
-  expect_true(data.table::is.data.table(res_nagel))
-  expect_equal(nrow(res_nagel), n_events)
+  expect_s3_class(res_nagel, "splikit_pseudo_correlation_result")
+  expect_equal(nrow(res_nagel$statistics), n_events)
   
   # Test with sparse matrices
   m1_sp <- Matrix::Matrix(m1, sparse = TRUE)
   m2_sp <- Matrix::Matrix(m2, sparse = TRUE)
   res_sp <- get_pseudo_correlation(ZDB, m1_sp, m2_sp)
-  expect_equal(res_sp$pseudo_correlation, res$pseudo_correlation)
+  expect_equal(res_sp$statistics$pseudo_correlation,
+               res$statistics$pseudo_correlation)
   
   # Test with SplikitObject wrapper
   m1_obj <- list(m1 = m1_sp, m2 = m2_sp)
   class(m1_obj) <- "SplikitObject"
   res_obj <- get_pseudo_correlation(m1_obj, m1_inclusion = ZDB)
-  expect_equal(res_obj$pseudo_correlation, res$pseudo_correlation)
+  expect_equal(res_obj$statistics$pseudo_correlation,
+               res$statistics$pseudo_correlation)
 })
 
 test_that("get_pseudo_correlation permutation null is valid and reproducible", {
@@ -47,33 +56,61 @@ test_that("get_pseudo_correlation permutation null is valid and reproducible", {
 
   res <- get_pseudo_correlation(ZDB, m1, m2,
                                 permutation_count = 50, permutation_seed = 123)
+  stats <- res$statistics
+  draws <- res$null_distribution
 
-  # New columns are present and well formed
-  expect_true(all(c("null_sd", "n_perm_valid", "emp_pvalue", "emp_padj") %in% names(res)))
-  expect_true(all(res$emp_pvalue > 0 & res$emp_pvalue <= 1))
-  expect_true(all(res$emp_padj >= 0 & res$emp_padj <= 1))
-  expect_true(all(res$n_perm_valid <= 50L))
+  # Statistic columns are present and well formed
+  expect_true(all(c("null_sd", "n_perm_valid", "emp_pvalue", "emp_padj") %in%
+                  names(stats)))
+  expect_true(all(stats$emp_pvalue > 0 & stats$emp_pvalue <= 1))
+  expect_true(all(stats$emp_padj >= 0 & stats$emp_padj <= 1))
+  expect_true(all(stats$n_perm_valid <= 50L))
   # Smallest achievable empirical p-value is 1 / (n_perm_valid + 1)
-  expect_true(all(res$emp_pvalue >= 1 / (res$n_perm_valid + 1) - 1e-9))
+  expect_true(all(stats$emp_pvalue >= 1 / (stats$n_perm_valid + 1) - 1e-9))
 
-  # Full null draws attached as an attribute, correctly shaped
-  draws <- attr(res, "null_draws")
-  expect_true(is.matrix(draws))
-  expect_equal(ncol(draws), 50L)
-  expect_equal(nrow(draws), nrow(res))
+  # The long null table has one row per retained event and permutation and
+  # preserves the matrix's permutation-major ordering.
+  expect_equal(nrow(draws), nrow(stats) * 50L)
+  expect_equal(draws$event, rep(stats$event, times = 50L))
+  expect_equal(draws$permutation, rep(seq_len(50L), each = nrow(stats)))
+
+  null_mat <- matrix(draws$null_pseudo_correlation,
+                     nrow = nrow(stats), ncol = 50L)
+  expect_equal(rowSums(!is.na(null_mat)), stats$n_perm_valid)
+  expect_equal(rowMeans(null_mat, na.rm = TRUE), stats$null_distribution)
+  expect_equal(apply(null_mat, 1L, stats::sd, na.rm = TRUE), stats$null_sd)
+
+  reconstructed_p <- vapply(seq_len(nrow(stats)), function(i) {
+    valid <- !is.na(null_mat[i, ])
+    (sum(abs(null_mat[i, valid]) >= abs(stats$pseudo_correlation[i])) + 1) /
+      (sum(valid) + 1)
+  }, numeric(1))
+  expect_equal(reconstructed_p, stats$emp_pvalue)
+
+  expect_equal(res$metadata$n_events_input, n_events)
+  expect_equal(res$metadata$n_events_retained, nrow(stats))
+  expect_equal(res$metadata$n_null_draws, nrow(draws))
+  expect_equal(res$metadata$n_null_valid,
+               sum(!is.na(draws$null_pseudo_correlation)))
+  expect_identical(res$metadata$pooled_null_usage, "descriptive")
 
   # Same seed -> identical null; observed correlation is seed-independent
   res2 <- get_pseudo_correlation(ZDB, m1, m2,
                                  permutation_count = 50, permutation_seed = 123)
+  expect_equal(res$statistics$null_distribution,
+               res2$statistics$null_distribution)
+  expect_equal(res$statistics$emp_pvalue, res2$statistics$emp_pvalue)
+  expect_equal(res$statistics$pseudo_correlation,
+               res2$statistics$pseudo_correlation)
   expect_equal(res$null_distribution, res2$null_distribution)
-  expect_equal(res$emp_pvalue, res2$emp_pvalue)
-  expect_equal(res$pseudo_correlation, res2$pseudo_correlation)
 
   # Different seed -> different draws, but the same observed correlation
   res3 <- get_pseudo_correlation(ZDB, m1, m2,
                                  permutation_count = 50, permutation_seed = 999)
-  expect_false(isTRUE(all.equal(attr(res, "null_draws"), attr(res3, "null_draws"))))
-  expect_equal(res$pseudo_correlation, res3$pseudo_correlation)
+  expect_false(isTRUE(all.equal(res$null_distribution$null_pseudo_correlation,
+                                res3$null_distribution$null_pseudo_correlation)))
+  expect_equal(res$statistics$pseudo_correlation,
+               res3$statistics$pseudo_correlation)
 
   # permutation_seed must not perturb the caller's global RNG stream
   set.seed(7); a <- runif(3)
@@ -86,9 +123,19 @@ test_that("get_pseudo_correlation permutation null is valid and reproducible", {
   # equals the single column of draws
   res1 <- get_pseudo_correlation(ZDB, m1, m2,
                                  permutation_count = 1, permutation_seed = 5)
-  expect_equal(ncol(attr(res1, "null_draws")), 1L)
-  expect_equal(unname(res1$null_distribution),
-               unname(attr(res1, "null_draws")[, 1]))
+  expect_equal(nrow(res1$null_distribution), nrow(res1$statistics))
+  expect_equal(unname(res1$statistics$null_distribution),
+               unname(res1$null_distribution$null_pseudo_correlation))
+
+  # Both component tables can be exported and imported without losing rows or
+  # values. Metadata intentionally remains an R list rather than CSV columns.
+  statistics_file <- tempfile(fileext = ".csv")
+  null_file <- tempfile(fileext = ".csv")
+  on.exit(unlink(c(statistics_file, null_file)), add = TRUE)
+  data.table::fwrite(stats, statistics_file)
+  data.table::fwrite(draws, null_file)
+  expect_equal(data.table::fread(statistics_file), stats)
+  expect_equal(data.table::fread(null_file), draws)
 
   # Input validation
   expect_error(get_pseudo_correlation(ZDB, m1, m2, permutation_count = 0),
diff --git a/tests/testthat/test-make_gene_count.R b/tests/testthat/test-make_gene_count.R
new file mode 100644
index 0000000..17116e1
--- /dev/null
+++ b/tests/testthat/test-make_gene_count.R
@@ -0,0 +1,166 @@
+make_gene_count_fixture <- function(include_em = TRUE) {
+  root <- tempfile("splikit-gene-")
+  raw_dir <- file.path(root, "raw")
+  filtered_dir <- file.path(root, "filtered")
+  dir.create(raw_dir, recursive = TRUE)
+  dir.create(filtered_dir, recursive = TRUE)
+
+  raw_counts <- Matrix::Matrix(
+    matrix(c(1, 2, 3, 4, 5, 6), nrow = 2), sparse = TRUE
+  )
+  em_counts <- Matrix::Matrix(
+    matrix(c(11, 12, 13, 14, 15, 16), nrow = 2), sparse = TRUE
+  )
+  filtered_counts <- raw_counts[, c(1, 3), drop = FALSE]
+
+  Matrix::writeMM(raw_counts, file.path(raw_dir, "matrix.mtx"))
+  if (include_em) {
+    Matrix::writeMM(em_counts, file.path(raw_dir, "UniqueAndMult-EM.mtx"))
+  }
+  Matrix::writeMM(filtered_counts, file.path(filtered_dir, "matrix.mtx"))
+
+  data.table::fwrite(
+    data.table::data.table(barcode = c("AA", "BB", "CC")),
+    file.path(raw_dir, "barcodes.tsv"), col.names = FALSE
+  )
+  data.table::fwrite(
+    data.table::data.table(barcode = c("AA", "CC")),
+    file.path(filtered_dir, "barcodes.tsv"), col.names = FALSE
+  )
+
+  features <- data.table::data.table(
+    id = c("gene1", "gene2"),
+    name = c("Gene 1", "Gene 2"),
+    type = "Gene Expression"
+  )
+  data.table::fwrite(features, file.path(raw_dir, "features.tsv"),
+                     sep = "\t", col.names = FALSE)
+  data.table::fwrite(features, file.path(filtered_dir, "features.tsv"),
+                     sep = "\t", col.names = FALSE)
+
+  list(
+    root = root,
+    raw_counts = raw_counts,
+    em_counts = em_counts,
+    filtered_counts = filtered_counts
+  )
+}
+
+test_that("make_gene_count selects raw, filtered, and custom count matrices", {
+  fixture <- make_gene_count_fixture()
+  on.exit(unlink(fixture$root, recursive = TRUE), add = TRUE)
+
+  raw_result <- make_gene_count(
+    fixture$root, "sample1",
+    use_internal_whitelist = FALSE,
+    matrix_source = "raw",
+    matrix_file = "matrix.mtx"
+  )
+  expect_equal(unname(as.matrix(raw_result)),
+               unname(as.matrix(fixture$raw_counts)))
+  expect_identical(colnames(raw_result), c("AA-sample1", "BB-sample1", "CC-sample1"))
+
+  filtered_result <- make_gene_count(
+    fixture$root, "sample1",
+    use_internal_whitelist = FALSE,
+    matrix_source = "filtered"
+  )
+  expect_equal(unname(as.matrix(filtered_result)),
+               unname(as.matrix(fixture$filtered_counts)))
+  expect_identical(colnames(filtered_result), c("AA-sample1", "CC-sample1"))
+
+  em_result <- make_gene_count(
+    fixture$root, "sample1",
+    use_internal_whitelist = FALSE,
+    matrix_source = "raw",
+    matrix_file = "UniqueAndMult-EM.mtx"
+  )
+  expect_equal(unname(as.matrix(em_result)),
+               unname(as.matrix(fixture$em_counts)))
+})
+
+test_that("matrix source is independent from barcode filtering", {
+  fixture <- make_gene_count_fixture()
+  on.exit(unlink(fixture$root, recursive = TRUE), add = TRUE)
+
+  result <- make_gene_count(
+    fixture$root, "sample1",
+    use_internal_whitelist = TRUE,
+    matrix_source = "raw",
+    matrix_file = "matrix.mtx"
+  )
+
+  expect_equal(unname(as.matrix(result)),
+               unname(as.matrix(fixture$raw_counts[, c(1, 3)])))
+  expect_identical(colnames(result), c("AA-sample1", "CC-sample1"))
+
+  em_result <- make_gene_count(
+    fixture$root, "sample1",
+    use_internal_whitelist = TRUE,
+    matrix_source = "raw",
+    matrix_file = "UniqueAndMult-EM.mtx"
+  )
+  expect_equal(unname(as.matrix(em_result)),
+               unname(as.matrix(fixture$em_counts[, c(1, 3)])))
+  expect_identical(colnames(em_result), c("AA-sample1", "CC-sample1"))
+})
+
+test_that("defaults prefer raw EM counts and fall back to raw matrix.mtx", {
+  fixture <- make_gene_count_fixture()
+  on.exit(unlink(fixture$root, recursive = TRUE), add = TRUE)
+
+  default_result <- make_gene_count(fixture$root, "sample1")
+  expect_equal(unname(as.matrix(default_result)),
+               unname(as.matrix(fixture$em_counts[, c(1, 3)])))
+  expect_identical(colnames(default_result), c("AA-sample1", "CC-sample1"))
+
+  fallback <- make_gene_count_fixture(include_em = FALSE)
+  on.exit(unlink(fallback$root, recursive = TRUE), add = TRUE)
+  fallback_result <- make_gene_count(fallback$root, "sample1")
+  expect_equal(unname(as.matrix(fallback_result)),
+               unname(as.matrix(fallback$raw_counts[, c(1, 3)])))
+})
+
+test_that("auto matrix source preserves historical directory selection", {
+  fixture <- make_gene_count_fixture()
+  on.exit(unlink(fixture$root, recursive = TRUE), add = TRUE)
+
+  auto_filtered <- make_gene_count(
+    fixture$root, "sample1",
+    use_internal_whitelist = TRUE,
+    matrix_source = "auto",
+    matrix_file = "matrix.mtx"
+  )
+  auto_raw <- make_gene_count(
+    fixture$root, "sample1",
+    use_internal_whitelist = FALSE,
+    matrix_source = "auto",
+    matrix_file = "matrix.mtx"
+  )
+
+  expect_equal(unname(as.matrix(auto_filtered)),
+               unname(as.matrix(fixture$filtered_counts)))
+  expect_equal(unname(as.matrix(auto_raw)),
+               unname(as.matrix(fixture$raw_counts)))
+})
+
+test_that("make_gene_count validates matrix selection", {
+  fixture <- make_gene_count_fixture()
+  on.exit(unlink(fixture$root, recursive = TRUE), add = TRUE)
+
+  expect_error(
+    make_gene_count(fixture$root, "sample1", matrix_source = "unknown"),
+    "arg.*should be one of"
+  )
+  expect_error(
+    make_gene_count(fixture$root, "sample1", matrix_file = ""),
+    "single non-empty filename"
+  )
+  expect_error(
+    make_gene_count(
+      fixture$root, "sample1",
+      matrix_source = "raw", matrix_file = "missing.mtx"
+    ),
+    "missing\\.mtx"
+  )
+})
diff --git a/vignettes/methods.Rmd b/vignettes/methods.Rmd
index f9e47f7..bb1de14 100644
--- a/vignettes/methods.Rmd
+++ b/vignettes/methods.Rmd
@@ -25,7 +25,7 @@ Splice junctions are grouped into *local junction variants* — junctions sharin
 
 ## Event-covariate association
 
-`get_pseudo_correlation()` fits a per-event binomial logistic GLM of the inclusion ratio on a target covariate by iteratively reweighted least squares, and reports a Cox-Snell / Nagelkerke pseudo-R-squared computed from the residual deviance. This quantifies how strongly each event tracks the covariate (e.g. a cluster label or a gene's expression). A per-event empirical null is built by permuting the cells (columns) of the covariate `permutation_count` times (default 100) and recomputing the statistic; this yields a two-sided empirical p-value (`emp_pvalue`) and a Benjamini-Hochberg adjustment (`emp_padj`). The permutations can be made reproducible with `permutation_seed` without disturbing the global RNG stream.
+`get_pseudo_correlation()` fits a per-event binomial logistic GLM of the inclusion ratio on a target covariate by iteratively reweighted least squares, and reports a Cox-Snell / Nagelkerke pseudo-R-squared computed from the residual deviance. This quantifies how strongly each event tracks the covariate (e.g. a cluster label or a gene's expression). A per-event empirical null is built by permuting the cells (columns) of the covariate `permutation_count` times (default 100) and recomputing the statistic; this yields a two-sided empirical p-value (`emp_pvalue`) and a Benjamini-Hochberg adjustment (`emp_padj`). The permutations can be made reproducible with `permutation_seed` without disturbing the global RNG stream. Results are returned as a structured object with a per-event `statistics` table and a long, exportable `null_distribution` table containing every retained event-by-permutation draw; the combined null values are descriptive and do not replace the event-specific null used for inference.
 
 ## Implementation
 
diff --git a/vignettes/splikit_manual.Rmd b/vignettes/splikit_manual.Rmd
index fdb1f6e..aadc4ee 100644
--- a/vignettes/splikit_manual.Rmd
+++ b/vignettes/splikit_manual.Rmd
@@ -545,12 +545,31 @@ summary(cor_values)  # Often negative correlation
 
 ```r
 gene_expression <- make_gene_count(
-  expression_dirs = "/path/to/10X/output",
+  expression_dirs = "/path/to/STARsolo/Gene",
   sample_ids = "Sample1",
   use_internal_whitelist = TRUE
 )
 ```
 
+`matrix_source` controls whether counts, features, and barcodes are read from
+the `raw/` or `filtered/` directory independently of barcode filtering. By
+default, `make_gene_count()` reads from `raw/`, uses
+`UniqueAndMult-EM.mtx` when that file is available, and otherwise falls back to
+`raw/matrix.mtx`. With `use_internal_whitelist = TRUE`, either matrix is then
+restricted to the barcodes in `filtered/barcodes.tsv`.
+
+To force the standard filtered matrix instead:
+
+```r
+gene_expression_filtered <- make_gene_count(
+  expression_dirs = "/path/to/STARsolo/Gene",
+  sample_ids = "Sample1",
+  matrix_source = "filtered",
+  matrix_file = "matrix.mtx",
+  use_internal_whitelist = FALSE
+)
+```
+
 #### `make_velo_count()`
 
 **Purpose**: Processes spliced/unspliced matrices for RNA velocity analysis.
@@ -694,7 +713,13 @@ get_pseudo_correlation(ZDB_matrix, m1_inclusion, m2_exclusion,
 - `permutation_count`: Number of cell-permutation null draws per event (default `100`); use `1` for the original single-permutation behaviour, and `1000`+ for reliable empirical FDR
 - `permutation_seed`: Optional seed for reproducible permutations (does not disturb the global RNG)
 
-**Output columns**: `event`, `pseudo_correlation`, `null_distribution` (mean of the per-event null draws), `null_sd`, `n_perm_valid`, `emp_pvalue` (two-sided empirical p-value), and `emp_padj` (BH-adjusted). The full per-event null draws are attached as `attr(result, "null_draws")`.
+**Output**: A structured result with `statistics`, `null_distribution`, and
+`metadata`. `statistics` contains `event`, `pseudo_correlation`,
+`null_distribution` (the per-event null mean), `null_sd`, `n_perm_valid`,
+`emp_pvalue` (two-sided event-wise empirical p-value), and `emp_padj`
+(BH-adjusted). The long `null_distribution` table contains `event`,
+`permutation`, and `null_pseudo_correlation`, with one row per retained event
+and permutation for direct export.
 
 **Use Cases**:
 - Correlating splicing patterns with gene expression
@@ -721,13 +746,20 @@ rownames(external_data) <- rownames(m1_matrix)
 correlations <- get_pseudo_correlation(external_data, m1_matrix, m2_matrix)
 
 # Examine results
-head(correlations[order(-abs(pseudo_correlation))])
+head(correlations$statistics[order(-abs(pseudo_correlation))])
 
-# Compare to null distribution
-plot(correlations$pseudo_correlation, correlations$null_distribution,
+# Compare observed correlations to their per-event null means
+plot(correlations$statistics$pseudo_correlation,
+     correlations$statistics$null_distribution,
      xlab = "Observed Correlation", ylab = "Null Distribution",
      main = "Pseudo-correlation Analysis")
 abline(0, 1, col = "red")
+
+# Export the two rectangular result tables independently
+statistics_file <- tempfile(fileext = ".csv")
+null_file <- tempfile(fileext = ".csv")
+data.table::fwrite(correlations$statistics, statistics_file)
+data.table::fwrite(correlations$null_distribution, null_file)
 ```
 
 ### `get_rowVar()`
@@ -885,7 +917,7 @@ if ("gene_expression" %in% names(toy_m1m2_data)) {
   )
   
   print("Top correlated events:")
-  print(head(correlations[order(-abs(pseudo_correlation))], 5))
+  print(head(correlations$statistics[order(-abs(pseudo_correlation))], 5))
 }
 
 cat("\n=== Analysis Complete ===\n")