# Splikit Manual for Single-Cell Splicing Analysis

## Table of Contents

1. Introduction
2. Installation
3. Core Concepts
4. Getting Started
5. R6 Object-Oriented Interface
6. Main Workflow Functions
7. Feature Selection Functions
8. Utility Functions
9. Complete Workflow Example
10. Advanced Usage
11. Troubleshooting
12. Function Reference

## Introduction

Splikit is a comprehensive R package designed for analyzing alternative splicing events in single-cell RNA sequencing data. Developed by the Computational and Statistical Genomics Laboratory at McGill University, splikit provides a complete toolkit for processing STARsolo splicing junction output, identifying highly variable splicing events, and performing downstream analyses.

### What Does Splikit Do?

Splikit transforms raw splicing junction data into meaningful biological insights through several key capabilities:

- Junction Data Processing: Converts STARsolo splicing junction output into structured matrices suitable for analysis
- Inclusion/Exclusion Matrix Generation: Creates M1 (inclusion) and M2 (exclusion) matrices that capture splicing behavior across cells
- Variable Event Detection: Identifies highly variable splicing events using advanced statistical methods
- Gene Expression Integration: Processes standard gene expression data alongside splicing information
- Statistical Analysis: Provides tools for correlation analysis, variance computation, and clustering evaluation ### How splikit get these data? We grouped splice junctions according to the coordinates of their intronic regions into “local junction variants” (LJVs). An LJV is defined as a set of junctions that share either the first or the last coordinate.

For each junction within an LJV, we extracted the junction abundance count for every cell, which we refer to as the M1 count. The M2 count for a given junction in each cell is defined as the sum of the counts of its alternative junctions—that is, all other junctions within the same LJV. In other words, for each junction and cell, the number of reads supporting that junction is contrasted with the total number of reads supporting its alternative junctions. For example, if \(( M1_{ij} )\) represents the count for junction \(j\) in cell \(i\), then the M2 count can be written as:

    M2_{ij} = \sum_{\substack{k \in J \\ k \neq j}} M1_{ik}

where \(J\) denotes the set of all junctions in the LJV. The grouping method is determined by whether the first or the last coordinate is used, with an appended `-E` or `-S` added to the event IDs accordingly. This approach may result in a single junction receiving two different measurements in the M2 counts (in the exclusion matrix) while having identical measurements in the M1 matrix (inclusion matrix).

We also address sample-specific junctions. If a junction is present in only a subset of samples, a vector of zeros is assigned to the M1 matrix for the samples where the junction is absent. The M2 measurements are then computed accordingly. This figure illustrates two exemplary LJVs.

Traditional single-cell RNA-seq analysis focuses primarily on gene expression levels, but alternative splicing represents an additional layer of cellular complexity. Splikit enables researchers to:

1. Capture Splicing Diversity: Identify cell-type-specific splicing patterns that might be missed by gene expression alone
2. Integrate Multiple Modalities: Combine splicing and expression data for comprehensive cellular characterization
3. Scale Efficiently: Handle large single-cell datasets with optimized C++ implementations
4. Ensure Reproducibility: Follow standardized workflows with well-documented functions

## Installation

### Prerequisites

Splikit requires several dependencies that are automatically installed:

    # Core dependencies (installed automatically)
    required_packages <- c("Rcpp", "RcppEigen", "RcppArmadillo", "Matrix", "data.table")

### Installation from Source

    # Install from GitHub (replace with actual repository)
    devtools::install_github("csglab/splikit")

    # Load the package
    library(splikit)

### Verify Installation

    # Load example data to test functionality
    toy_data <- load_toy_SJ_object()
    print(names(toy_data))

## Core Concepts

Understanding splikit requires familiarity with several key concepts in splicing analysis:

### Splicing Junction Events

A splicing junction represents the connection between two exons after intron removal. In single-cell data, we observe:

- Inclusion Events: When a particular exon or junction is included in the mature transcript
- Exclusion Events: When a particular exon or junction is skipped during splicing

### Matrix Representations

Splikit uses two primary matrix types:

M1 Matrix (Inclusion Matrix): - Rows represent splicing events - Columns represent individual cells - Values indicate the number of reads supporting inclusion of each event

M2 Matrix (Exclusion Matrix): - Same dimensions as M1 - Values indicate the number of reads supporting exclusion of each event - Calculated based on the total splicing activity at shared coordinate groups

### Coordinate Groups

Splikit groups splicing junctions that share start or end coordinates, enabling the calculation of exclusion events. This grouping is essential because:

- Multiple junctions can share the same donor (5’) or acceptor (3’) site
- The total splicing activity at these sites helps quantify exclusion events
- Groups with only one member cannot generate exclusion events

### Variable Event Detection

The package identifies highly variable splicing events using statistical methods:

- Deviance-based Approach: Models splicing ratios using beta-binomial or similar distributions
- Variance Stabilization: Applies transformations to identify events with genuine biological variability
- Library-wise Calculation: Accounts for technical variation between samples or libraries

## Getting Started

Let’s begin with a simple example using the provided toy dataset:

    library(splikit)

    # Load the toy splicing junction data
    toy_sj_data <- load_toy_SJ_object()

    # Examine the structure
    str(toy_sj_data, max.level = 2)

    # The toy data contains junction abundance information for multiple samples
    # Each sample has 'eventdata' (metadata) and 'junction_ab' (sparse matrix)

### Quick Start Workflow

Here’s the minimal workflow to get from raw data to variable events:

    # Step 1: Create M1 matrix
    m1_result <- make_m1(toy_sj_data)
    m1_matrix <- m1_result$m1_inclusion_matrix
    event_data <- m1_result$event_data

    # Step 2: Create M2 matrix
    m2_matrix <- make_m2(m1_matrix, event_data)

    # Step 3: Find highly variable events
    variable_events <- find_variable_events(m1_matrix, m2_matrix)

    # View top variable events
    head(variable_events[order(-sum_deviance)])

## R6 Object-Oriented Interface

Starting with version 2.0.0, splikit introduces a modern R6 object-oriented interface through the `SplikitObject` class. This provides a more intuitive and streamlined workflow with method chaining support.

### Creating a SplikitObject

    library(splikit)

    # Load toy data for demonstration
    toy_sj_data <- load_toy_SJ_object()

    # Create SplikitObject from junction abundance data
    spl <- SplikitObject$new(junction_ab_object = toy_sj_data)

    # View object summary
    print(spl)

### Complete Workflow with Method Chaining

The R6 interface supports method chaining for a clean, pipeline-style workflow:

    # Complete analysis pipeline with method chaining
    spl <- SplikitObject$new(junction_ab_object = toy_sj_data)$
      makeM1(min_counts = 1)$
      makeM2(n_threads = 4, verbose = TRUE)$
      findVariableEvents(min_row_sum = 30, n_threads = 4)

    # Access results
    m1_matrix <- spl$m1
    m2_matrix <- spl$m2
    variable_events <- spl$variableEvents
    event_data <- spl$eventData

### Step-by-Step Workflow

For more control, you can call methods individually:

    # Create object
    spl <- SplikitObject$new(junction_ab_object = toy_sj_data)

    # Step 1: Create M1 matrix
    spl$makeM1(min_counts = 1, verbose = TRUE)
    print(dim(spl$m1))  # Check dimensions

    # Step 2: Create M2 matrix with multi-threading
    spl$makeM2(
      n_threads = 4,           # Number of OpenMP threads
      use_cpp = TRUE,          # Use optimized C++ implementation
      verbose = TRUE
    )
    print(dim(spl$m2))  # Same dimensions as M1

    # Step 3: Find highly variable events
    spl$findVariableEvents(
      min_row_sum = 30,
      n_threads = 4
    )

    # View top variable events
    top_events <- spl$variableEvents[order(-sum_deviance)][1:10]
    print(top_events)

### Subsetting SplikitObjects

SplikitObjects support flexible subsetting by events, cells, or both:

    # Subset by events (rows)
    top_100_events <- spl$variableEvents$events[1:100]
    spl_subset <- spl$subset(events = top_100_events)

    # Subset by cells (columns)
    selected_cells <- colnames(spl$m1)[1:500]
    spl_cells <- spl$subset(cells = selected_cells)

    # Subset both events and cells
    spl_combined <- spl$subset(
      events = top_100_events,
      cells = selected_cells
    )

    # Check dimensions
    print(dim(spl_subset$m1))
    print(dim(spl_cells$m1))
    print(dim(spl_combined$m1))

### Working with Existing Matrices

You can also create a SplikitObject from pre-computed matrices:

    # Load pre-computed M1/M2 data
    toy_m1m2 <- load_toy_M1_M2_object()

    # Create SplikitObject from existing matrices
    spl <- SplikitObject$new(
      m1_matrix = toy_m1m2$m1,
      m2_matrix = toy_m1m2$m2
    )

    # Directly find variable events (M1 and M2 already computed)
    spl$findVariableEvents(min_row_sum = 50, n_threads = 4)

### Accessing SplikitObject Fields

The SplikitObject provides access to all data through fields:

    # Core matrices
    spl$m1              # Inclusion matrix (events × cells)
    spl$m2              # Exclusion matrix (events × cells)

    # Metadata
    spl$eventData       # Event metadata data.table
    spl$variableEvents  # Variable events results

    # Original data
    spl$junctionAbObject  # Raw junction abundance data

### R6 vs Function-Based Interface Comparison

Both interfaces produce identical results. Choose based on your workflow preference:

R6 Interface (recommended for new users):

    spl <- SplikitObject$new(junction_ab_object = toy_sj_data)$
      makeM1(min_counts = 1)$
      makeM2(n_threads = 4)$
      findVariableEvents(min_row_sum = 30)

    m1 <- spl$m1
    m2 <- spl$m2
    hve <- spl$variableEvents

Function-Based Interface (traditional):

    m1_result <- make_m1(toy_sj_data, min_counts = 1)
    m2 <- make_m2(m1_result$m1_inclusion_matrix, m1_result$event_data, n_threads = 4)
    hve <- find_variable_events(m1_result$m1_inclusion_matrix, m2, min_row_sum = 30)

### Performance Features

The SplikitObject methods support the same performance optimizations as the underlying functions:

- Multi-threading: Set `n_threads` parameter for OpenMP parallelization
- Memory-efficient M2: Version 2.0.0 includes a completely rewritten `make_m2` implementation with ~30x memory reduction
- C++ optimization: Set `use_cpp = TRUE` (default) for maximum performance

    # Check if OpenMP is available
    check_openmp_enabled()

    # Use all available cores for large datasets
    n_cores <- parallel::detectCores() - 1

    spl$makeM2(n_threads = n_cores, use_cpp = TRUE, verbose = TRUE)
    spl$findVariableEvents(n_threads = n_cores, min_row_sum = 50)

## Main Workflow Functions

### `make_junction_ab()`

Purpose: Processes STARsolo splicing junction output into structured R objects.

Function Signature:

    make_junction_ab(STARsolo_SJ_dirs, white_barcode_lists = NULL,
                     sample_ids, use_internal_whitelist = TRUE,
                     verbose = FALSE, keep_multi_mapped_junctions = FALSE, ...)

Parameters:

- `STARsolo_SJ_dirs`: Character vector of paths to STARsolo SJ directories
- `white_barcode_lists`: Optional list of barcode whitelists for filtering
- `sample_ids`: Unique identifiers for each sample
- `use_internal_whitelist`: Whether to use STARsolo’s internal filtered barcodes (default TRUE)
- `verbose`: Logical for detailed progress reporting (default FALSE)
- `keep_multi_mapped_junctions`: Whether to keep multi-mapped junctions (default FALSE)
- `...`: Additional parameters for future extensions

Expected Directory Structure:

    STARsolo_output/
    ├── SJ/
    │   ├── raw/
    │   │   ├── matrix.mtx      # Junction count matrix
    │   │   ├── barcodes.tsv    # Cell barcodes
    │   │   └── features.tsv    # Junction features
    │   └── filtered/
    │       └── barcodes.tsv    # Filtered barcodes
    └── SJ.out.tab              # STAR junction output

Returns: A named list where each element corresponds to a sample and contains: - `eventdata`: Data.table with junction metadata (chromosome, coordinates, strand) - `junction_ab`: Sparse matrix of junction abundances (junctions × cells)

Example:

    # Single sample processing
    sj_dirs <- "/path/to/STARsolo/output/SJ"
    sample_names <- "Sample1"

    junction_data <- make_junction_ab(
      STARsolo_SJ_dirs = sj_dirs,
      sample_ids = sample_names,
      use_internal_whitelist = TRUE
    )

    # Multiple samples with custom barcodes
    sj_dirs <- c("/path/to/sample1/SJ", "/path/to/sample2/SJ")
    sample_names <- c("Control", "Treatment")
    custom_barcodes <- list(
      c("AAACCTGAGCATCGCA", "AAACCTGAGCGAACGC"),  # Control barcodes
      c("AAACCTGAGGAATGGT", "AAACCTGAGGGCTCTC")   # Treatment barcodes
    )

    junction_data <- make_junction_ab(
      STARsolo_SJ_dirs = sj_dirs,
      white_barcode_lists = custom_barcodes,
      sample_ids = sample_names,
      use_internal_whitelist = FALSE
    )

Understanding the Output:

The `eventdata` component contains crucial information for each junction: - `chr`, `start`, `end`, `strand`: Genomic coordinates - `start_cor_id`, `end_cor_id`: Coordinate identifiers for grouping - `row_names_mtx`: Unique junction identifiers - `sample_id`: Sample origin

### `make_m1()`

Purpose: Transforms junction abundance data into the inclusion matrix (M1) with coordinate-based grouping.

Function Signature:

    make_m1(junction_ab_object, min_counts = 1, verbose = FALSE)

Parameters: - `junction_ab_object`: Output from `make_junction_ab()` - `min_counts`: Minimum count threshold for filtering events (default 1) - `verbose`: Logical for detailed progress reporting (default FALSE)

The Coordinate Grouping Process:

This function performs sophisticated grouping of junctions that share coordinates:

1. Start Coordinate Grouping: Junctions sharing the same 5’ splice site are grouped
2. End Coordinate Grouping: Junctions sharing the same 3’ splice site are grouped
3. Matrix Expansion: Each group generates new rows in the M1 matrix
4. Event Naming: Groups are labeled with “_S” (start) or “_E” (end) suffixes

Returns: - `m1_inclusion_matrix`: Sparse matrix with expanded events based on coordinate groups - `event_data`: Enhanced metadata with group information

Example:

    # Process junction data into M1 matrix
    m1_result <- make_m1(junction_data)

    # Examine the structure
    print(dim(m1_result$m1_inclusion_matrix))  # Events × Cells
    print(nrow(m1_result$event_data))          # Number of events

    # Check the grouping information
    table(m1_result$event_data$group_count)    # Distribution of group sizes

Understanding M1 Structure:

The M1 matrix represents inclusion evidence: - Each row corresponds to a specific splicing event (possibly grouped) - Each column corresponds to a single cell - Values represent the number of reads supporting inclusion of that event - Events with “_S” suffix represent start coordinate groups - Events with “_E” suffix represent end coordinate groups

### `make_m2()`

Purpose: Creates the exclusion matrix (M2) from the M1 matrix and event data with intelligent memory management.

Function Signature:

    make_m2(m1_inclusion_matrix, eventdata, batch_size = 5000,
            memory_threshold = 2e9, force_fast = FALSE,
            multi_thread = FALSE, n_threads = 1, use_cpp = TRUE,
            verbose = FALSE)

Parameters: - `m1_inclusion_matrix`: Output from `make_m1()` - `eventdata`: Event metadata from `make_m1()` - `batch_size`: Number of groups to process per batch when using batched mode (default 5000) - `memory_threshold`: Maximum rows before switching to batched processing (default 2e9) - `force_fast`: Force fast processing regardless of memory estimates (default FALSE) - `multi_thread`: Enable parallel processing for batched operations (default FALSE) - `n_threads`: Number of OpenMP threads for C++ implementation (default 1) - `use_cpp`: Use optimized C++ implementation (default TRUE, recommended) - `verbose`: Logical for detailed progress reporting (default FALSE)

Performance Notes (v2.0.0): The C++ implementation has been completely rewritten for version 2.0.0 with significant improvements: - ~8x faster execution through two-pass CSC format construction - ~30x lower memory usage by eliminating intermediate dense matrices - Memory usage now scales with output size only: O(nnz_output) + O(n_groups) - Previous versions used O(n_groups × n_cells) for intermediate storage

The Exclusion Calculation Logic:

For each coordinate group, M2 values are calculated as:

    M2[event, cell] = Total_Group_Activity[group, cell] - M1[event, cell]

This means: 1. For each coordinate group, sum all inclusion reads across group members 2. For each individual event, subtract its inclusion reads from the group total 3. The remainder represents “exclusion” evidence for that specific event

Example:

    # Create M2 matrix
    m2_matrix <- make_m2(m1_result$m1_inclusion_matrix, m1_result$event_data)

    # Compare dimensions
    print(dim(m1_result$m1_inclusion_matrix))
    print(dim(m2_matrix))  # Should be identical

    # Examine the relationship between M1 and M2
    cor_values <- diag(cor(as.matrix(m1_result$m1_inclusion_matrix[1:100, 1:50]),
                          as.matrix(m2_matrix[1:100, 1:50])))
    summary(cor_values)  # Often negative correlation

Interpreting M2 Values:

- High M2 values indicate strong evidence for exclusion of that specific event
- M2 values of zero occur when an event represents all activity in its coordinate group
- The M1/M2 ratio provides information about splicing ratios

### Additional Data Processing Functions

#### `make_gene_count()`

Purpose: Processes standard gene expression matrices alongside splicing data.

    gene_expression <- make_gene_count(
      expression_dirs = "/path/to/10X/output",
      sample_ids = "Sample1",
      use_internal_whitelist = TRUE
    )

#### `make_velo_count()`

Purpose: Processes spliced/unspliced matrices for RNA velocity analysis.

    velocity_data <- make_velo_count(
      velocyto_dirs = "/path/to/velocyto/output",
      sample_ids = "Sample1",
      merge_counts = TRUE  # Combine spliced/unspliced across samples
    )

## Feature Selection Functions

### `find_variable_events()`

Purpose: Identifies highly variable splicing events using deviance-based statistics.

Function Signature:

    find_variable_events(m1_matrix, m2_matrix, min_row_sum = 50,
                         n_threads = 1, verbose = TRUE, ...)

Parameters: - `m1_matrix`, `m2_matrix`: Inclusion and exclusion matrices - `min_row_sum`: Minimum total reads per event for inclusion in analysis - `n_threads`: Number of threads for parallel processing (requires OpenMP) - `verbose`: Whether to print progress messages

The Statistical Method:

The function calculates deviance statistics for each event:

1. Library-wise Processing: Analyzes each sample/library separately
2. Ratio Modeling: Models the M1/(M1+M2) ratio for each event
3. Deviance Calculation: Computes how much each event deviates from expected behavior
4. Summation: Combines deviance scores across all libraries

Returns: A data.table with columns: - `events`: Event identifiers - `sum_deviance`: Combined deviance score across all libraries

Example:

    # Load toy data
    toy_data <- load_toy_M1_M2_object()
    m1_matrix <- toy_data$m1
    m2_matrix <- toy_data$m2

    # Find variable events with different parameters
    hve_default <- find_variable_events(m1_matrix, m2_matrix)
    hve_strict <- find_variable_events(m1_matrix, m2_matrix, min_row_sum = 100)

    # Compare results
    print(nrow(hve_default))  # Number of events analyzed
    print(nrow(hve_strict))   # Fewer events with stricter filter

    # Top variable events
    top_events <- hve_default[order(-sum_deviance)][1:20]
    print(top_events)

Interpreting Results:

- Higher `sum_deviance` values indicate more variable splicing behavior
- Events with high variability might represent:
  - Cell-type-specific splicing patterns
  - Response to environmental conditions
  - Technical artifacts (require validation)

### `find_variable_genes()`

Purpose: Identifies highly variable genes using either variance stabilization or deviance methods.

Function Signature:

    find_variable_genes(gene_expression_matrix, method = "vst",
                       n_threads = 1, verbose = TRUE, ...)

Parameters: - `gene_expression_matrix`: Sparse gene expression matrix - `method`: Either “vst” (variance stabilizing transformation) or “sum_deviance” - `n_threads`: Number of threads for parallel processing (only for sum_deviance method) - `verbose`: Logical for progress reporting (default TRUE)

Method Comparison:

VST Method (Default): - Implements Seurat-style variance stabilization - Faster computation through C++ implementation - Good for identifying highly expressed variable genes

Sum Deviance Method: - Similar to the splicing event approach - Models gene expression using negative binomial deviances - Better for lowly expressed genes

Example:

    # Using toy gene expression data
    toy_data <- load_toy_M1_M2_object()
    gene_expr <- toy_data$gene_expression

    # VST method (faster)
    hvg_vst <- find_variable_genes(gene_expr, method = "vst")
    print(head(hvg_vst[order(-standardize_variance)]))

    # Deviance method (more comprehensive)
    hvg_dev <- find_variable_genes(gene_expr, method = "sum_deviance")
    print(head(hvg_dev[order(-sum_deviance)]))

    # Compare methods
    length(intersect(hvg_vst[1:500]$events, hvg_dev[1:500]$events))

## Utility Functions

### `get_pseudo_correlation()`

Purpose: Computes pseudo-R² correlation between external data and splicing behavior using beta-binomial modeling.

Function Signature:

    get_pseudo_correlation(ZDB_matrix, m1_inclusion, m2_exclusion,
                          metric = "CoxSnell", suppress_warnings = TRUE)

Parameters: - `ZDB_matrix`: External data matrix (e.g., gene expression, protein levels) - must be dense - `m1_inclusion`, `m2_exclusion`: Can be either sparse or dense matrices - `metric`: R² metric to use - “CoxSnell” (default) or “Nagelkerke” - `suppress_warnings`: Whether to suppress computational warnings

Use Cases: - Correlating splicing patterns with gene expression - Identifying splicing events that respond to external stimuli - Validating splicing predictions against protein data

Example:

    # Create example data
    toy_sj <- load_toy_SJ_object()
    m1_obj <- make_m1(toy_sj)
    m1_matrix <- as.matrix(m1_obj$m1_inclusion_matrix)  # Convert to dense
    m2_matrix <- as.matrix(make_m2(m1_obj$m1_inclusion_matrix, m1_obj$event_data))

    # Create dummy external data (e.g., gene expression)
    external_data <- matrix(
      rnorm(nrow(m1_matrix) * ncol(m1_matrix)),
      nrow = nrow(m1_matrix),
      ncol = ncol(m1_matrix)
    )
    rownames(external_data) <- rownames(m1_matrix)

    # Calculate pseudo-correlations
    correlations <- get_pseudo_correlation(external_data, m1_matrix, m2_matrix)

    # Examine results
    head(correlations[order(-abs(pseudo_correlation))])

    # Compare to null distribution
    plot(correlations$pseudo_correlation, correlations$null_distribution,
         xlab = "Observed Correlation", ylab = "Null Distribution",
         main = "Pseudo-correlation Analysis")
    abline(0, 1, col = "red")

### `get_rowVar()`

Purpose: Efficiently computes row-wise variance for dense or sparse matrices.

Function Signature:

    get_rowVar(M, verbose = FALSE)

Example:

    # Dense matrix
    dense_mat <- matrix(rnorm(1000), nrow = 100)
    var_dense <- get_rowVar(dense_mat)

    # Sparse matrix
    library(Matrix)
    sparse_mat <- rsparsematrix(100, 10, density = 0.1)
    var_sparse <- get_rowVar(sparse_mat, verbose = TRUE)

    # Compare with base R (for dense matrices)
    var_base <- apply(dense_mat, 1, var)
    all.equal(var_dense, var_base)

### `get_silhouette_mean()`

Purpose: Computes average silhouette width for clustering evaluation.

Function Signature:

    get_silhouette_mean(X, cluster_assignments, n_threads = 1)

Example:

    # Create example clustering data
    set.seed(42)
    pc_matrix <- matrix(rnorm(1000 * 15), nrow = 1000, ncol = 15)
    clusters <- sample(1:5, 1000, replace = TRUE)

    # Calculate silhouette score
    n_threads <- parallel::detectCores()
    silhouette_score <- get_silhouette_mean(pc_matrix, clusters, n_threads)

    print(paste("Average silhouette score:", round(silhouette_score, 3)))

## Complete Workflow Example

Here’s a comprehensive example that demonstrates the entire splikit workflow:

### Scenario: Analyzing Cell Type-Specific Splicing

    library(splikit)
    library(Matrix)
    library(data.table)

    # Step 1: Load and examine data
    cat("=== Loading Data ===\n")
    toy_sj_data <- load_toy_SJ_object()
    toy_m1m2_data <- load_toy_M1_M2_object()

    print("Available samples:")
    print(names(toy_sj_data))

    # Step 2: Process junction data into M1/M2 matrices
    cat("\n=== Creating Inclusion/Exclusion Matrices ===\n")
    m1_result <- make_m1(toy_sj_data)
    m1_matrix <- m1_result$m1_inclusion_matrix
    event_data <- m1_result$event_data

    m2_matrix <- make_m2(m1_matrix, event_data)

    print(paste("M1 matrix dimensions:", paste(dim(m1_matrix), collapse = " x ")))
    print(paste("M2 matrix dimensions:", paste(dim(m2_matrix), collapse = " x ")))

    # Step 3: Identify highly variable events
    cat("\n=== Finding Variable Events ===\n")
    variable_events <- find_variable_events(
      m1_matrix = m1_matrix,
      m2_matrix = m2_matrix,
      min_row_sum = 30,  # Lower threshold for toy data
      verbose = TRUE
    )

    print(paste("Found", nrow(variable_events), "variable events"))

    # Step 4: Examine top variable events
    top_events <- variable_events[order(-sum_deviance)][1:10]
    print("Top 10 most variable events:")
    print(top_events)

    # Step 5: Calculate splicing ratios for top events
    cat("\n=== Analyzing Splicing Ratios ===\n")
    top_event_names <- top_events$events[1:5]
    top_m1 <- m1_matrix[top_event_names, ]
    top_m2 <- m2_matrix[top_event_names, ]

    # Calculate PSI (Percent Spliced In) values
    psi_values <- top_m1 / (top_m1 + top_m2)
    psi_values[is.nan(psi_values)] <- 0

    # Examine PSI distribution
    for (i in 1:nrow(psi_values)) {
      event_name <- rownames(psi_values)[i]
      psi_vector <- as.numeric(psi_values[i, ])
      psi_vector <- psi_vector[psi_vector > 0]  # Remove zeros

      if (length(psi_vector) > 10) {
        cat(sprintf("Event %s: PSI mean=%.3f, sd=%.3f, n=%d\n",
                    event_name, mean(psi_vector), sd(psi_vector), length(psi_vector)))
      }
    }

    # Step 6: Gene expression analysis (if available)
    cat("\n=== Gene Expression Analysis ===\n")
    if ("gene_expression" %in% names(toy_m1m2_data)) {
      gene_expr <- toy_m1m2_data$gene_expression

      hvg_results <- find_variable_genes(
        gene_expression_matrix = gene_expr,
        method = "sum_deviance"
      )

      print(paste("Found", nrow(hvg_results), "genes for analysis"))
      print("Top 5 most variable genes:")
      print(head(hvg_results[order(-sum_deviance)], 5))
    }

    # Step 7: Advanced analysis - pseudo-correlations
    cat("\n=== Pseudo-correlation Analysis ===\n")
    if ("gene_expression" %in% names(toy_m1m2_data)) {
      # Select subset for demonstration
      n_events <- min(100, nrow(m1_matrix))
      n_cells <- min(50, ncol(m1_matrix))

      m1_subset <- as.matrix(m1_matrix[1:n_events, 1:n_cells])
      m2_subset <- as.matrix(m2_matrix[1:n_events, 1:n_cells])

      # Create matching gene expression subset
      gene_subset <- as.matrix(gene_expr[1:n_events, 1:n_cells])
      rownames(gene_subset) <- rownames(m1_subset)

      # Calculate correlations
      correlations <- get_pseudo_correlation(
        ZDB_matrix = gene_subset,
        m1_inclusion = m1_subset,
        m2_exclusion = m2_subset
      )

      print("Top correlated events:")
      print(head(correlations[order(-abs(pseudo_correlation))], 5))
    }

    cat("\n=== Analysis Complete ===\n")

### Advanced Workflow: Multi-sample Comparison

    # Simulate multi-sample analysis
    analyze_multi_sample <- function() {
      cat("=== Multi-Sample Splicing Analysis ===\n")

      # Load data
      toy_data <- load_toy_M1_M2_object()
      m1_matrix <- toy_data$m1
      m2_matrix <- toy_data$m2

      # Extract sample information from cell barcodes
      cell_barcodes <- colnames(m1_matrix)
      sample_ids <- sub("^.{16}-(.*$)", "\\1", cell_barcodes)
      unique_samples <- unique(sample_ids)

      print(paste("Found", length(unique_samples), "samples"))
      print(table(sample_ids))

      # Analyze each sample separately
      sample_results <- list()

      for (sample in unique_samples) {
        cat(sprintf("\nAnalyzing sample: %s\n", sample))

        # Subset matrices for this sample
        sample_cells <- which(sample_ids == sample)
        m1_sample <- m1_matrix[, sample_cells, drop = FALSE]
        m2_sample <- m2_matrix[, sample_cells, drop = FALSE]

        # Find variable events for this sample
        if (ncol(m1_sample) > 5) {  # Minimum cells for analysis
          var_events <- find_variable_events(
            m1_sample, m2_sample,
            min_row_sum = 10,  # Lower threshold for subsets
            verbose = FALSE
          )

          sample_results[[sample]] <- var_events[order(-sum_deviance)][1:20]
          print(paste("Top variable event:", var_events[1]$events))
        }
      }

      # Compare top events across samples
      if (length(sample_results) > 1) {
        cat("\n=== Cross-Sample Comparison ===\n")

        all_top_events <- unique(unlist(lapply(sample_results, function(x) x$events)))
        print(paste("Total unique top events:", length(all_top_events)))

        # Find common highly variable events
        event_counts <- table(unlist(lapply(sample_results, function(x) x$events[1:10])))
        common_events <- names(event_counts)[event_counts >= 2]

        if (length(common_events) > 0) {
          print("Events variable in multiple samples:")
          print(common_events)
        }
      }

      return(sample_results)
    }

    # Run multi-sample analysis
    multi_sample_results <- analyze_multi_sample()

## Advanced Usage

### Performance Optimization

#### Parallel Processing

Splikit supports parallel processing through OpenMP for computationally intensive functions:

    # Check OpenMP availability
    openmp_available <- check_openmp_enabled()
    print(paste("OpenMP enabled:", openmp_available))

    if (openmp_available) {
      # Use multiple threads for large datasets
      n_threads <- parallel::detectCores() - 1  # Leave one core free

      variable_events <- find_variable_events(
        m1_matrix, m2_matrix,
        n_threads = n_threads
      )
    }

#### Memory Management

For large datasets, consider these strategies:

    # Process data in chunks for very large matrices
    process_large_dataset <- function(m1_matrix, m2_matrix, chunk_size = 1000) {
      n_events <- nrow(m1_matrix)
      results <- list()

      for (i in seq(1, n_events, chunk_size)) {
        end_idx <- min(i + chunk_size - 1, n_events)
        cat(sprintf("Processing events %d to %d\n", i, end_idx))

        m1_chunk <- m1_matrix[i:end_idx, , drop = FALSE]
        m2_chunk <- m2_matrix[i:end_idx, , drop = FALSE]

        chunk_result <- find_variable_events(m1_chunk, m2_chunk, verbose = FALSE)
        results[[length(results) + 1]] <- chunk_result
      }

      # Combine results
      do.call(rbind, results)
    }

### Integration with Other Packages

#### Seurat Integration

    library(Seurat)

    # Create Seurat object with splicing data
    create_splicing_seurat <- function(m1_matrix, m2_matrix, metadata = NULL) {
      # Calculate PSI values
      psi_matrix <- m1_matrix / (m1_matrix + m2_matrix)
      psi_matrix[is.nan(psi_matrix)] <- 0

      # Create Seurat object
      seurat_obj <- CreateSeuratObject(
        counts = psi_matrix,
        project = "SplicingAnalysis",
        meta.data = metadata
      )

      # Add M1 and M2 as additional assays
      seurat_obj[["M1"]] <- CreateAssayObject(counts = m1_matrix)
      seurat_obj[["M2"]] <- CreateAssayObject(counts = m2_matrix)

      return(seurat_obj)
    }

#### SingleCellExperiment Integration

    library(SingleCellExperiment)

    # Create SCE object with multiple modalities
    create_splicing_sce <- function(gene_expr, m1_matrix, m2_matrix) {
      # Ensure consistent cell names
      common_cells <- intersect(colnames(gene_expr), colnames(m1_matrix))

      gene_expr <- gene_expr[, common_cells]
      m1_matrix <- m1_matrix[, common_cells]
      m2_matrix <- m2_matrix[, common_cells]

      # Create main SCE object with gene expression
      sce <- SingleCellExperiment(
        assays = list(counts = gene_expr),
        colData = DataFrame(cell_id = common_cells)
      )

      # Add splicing information as alternative experiments
      splicing_sce <- SingleCellExperiment(
        assays = list(inclusion = m1_matrix, exclusion = m2_matrix)
      )

      altExp(sce, "splicing") <- splicing_sce

      return(sce)
    }

### Custom Analysis Functions

#### Differential Splicing Analysis

    # Function to identify differentially spliced events between conditions
    find_differential_splicing <- function(m1_matrix, m2_matrix, cell_groups,
                                         group1, group2, min_cells = 10) {
      # Subset cells for comparison
      cells_group1 <- which(cell_groups == group1)
      cells_group2 <- which(cell_groups == group2)

      if (length(cells_group1) < min_cells || length(cells_group2) < min_cells) {
        stop("Insufficient cells in one or both groups")
      }

      # Calculate PSI values for each group
      psi1 <- m1_matrix[, cells_group1] / (m1_matrix[, cells_group1] + m2_matrix[, cells_group1])
      psi2 <- m1_matrix[, cells_group2] / (m1_matrix[, cells_group2] + m2_matrix[, cells_group2])

      # Handle NaN values
      psi1[is.nan(psi1)] <- 0
      psi2[is.nan(psi2)] <- 0

      # Calculate statistics for each event
      results <- data.table(
        event = rownames(m1_matrix),
        mean_psi_group1 = rowMeans(psi1),
        mean_psi_group2 = rowMeans(psi2),
        delta_psi = rowMeans(psi2) - rowMeans(psi1)
      )

      # Add statistical tests (simplified - use appropriate tests for your data)
      results[, p_value := apply(cbind(psi1, psi2), 1, function(row) {
        group1_vals <- row[1:length(cells_group1)]
        group2_vals <- row[(length(cells_group1) + 1):length(row)]

        # Remove zeros for testing
        group1_vals <- group1_vals[group1_vals > 0]
        group2_vals <- group2_vals[group2_vals > 0]

        if (length(group1_vals) < 3 || length(group2_vals) < 3) {
          return(1)  # Not enough data
        }

        tryCatch({
          wilcox.test(group1_vals, group2_vals)$p.value
        }, error = function(e) 1)
      })]

      # Adjust p-values
      results[, p_adjusted := p.adjust(p_value, method = "BH")]

      return(results[order(abs(delta_psi), decreasing = TRUE)])
    }

## Troubleshooting

### Common Issues and Solutions

#### Issue 1: Memory Errors with Large Datasets

Symptom: R crashes or reports “cannot allocate vector of size X”

Solutions:

    # 1. Increase memory limit (Windows)
    memory.limit(size = 8000)  # 8GB

    # 2. Process data in chunks
    # 3. Use more selective filtering
    variable_events <- find_variable_events(
      m1_matrix, m2_matrix,
      min_row_sum = 100  # Higher threshold reduces events
    )

    # 4. Subsample cells for initial exploration
    n_cells_subset <- min(1000, ncol(m1_matrix))
    cells_to_keep <- sample(ncol(m1_matrix), n_cells_subset)
    m1_subset <- m1_matrix[, cells_to_keep]
    m2_subset <- m2_matrix[, cells_to_keep]

#### Issue 2: STARsolo Directory Structure Problems

Symptom: “No abundance matrix found” or similar file not found errors

Solution:

    # Verify directory structure
    check_starsolo_structure <- function(starsolo_dir) {
      required_files <- c(
        "raw/matrix.mtx",
        "raw/barcodes.tsv",
        "../../SJ.out.tab"
      )

      for (file in required_files) {
        full_path <- file.path(starsolo_dir, file)
        if (!file.exists(full_path)) {
          cat("Missing file:", full_path, "\n")
          return(FALSE)
        }
      }

      cat("Directory structure is correct\n")
      return(TRUE)
    }

    # Check before processing
    check_starsolo_structure("/path/to/STARsolo/SJ")

#### Issue 3: Barcode Mismatch Issues

Symptom: “All barcodes were removed after trimming”

Solutions:

    # 1. Examine barcode formats
    raw_barcodes <- data.table::fread("raw/barcodes.tsv", header = FALSE)$V1
    filter_barcodes <- data.table::fread("filtered/barcodes.tsv", header = FALSE)$V1

    head(raw_barcodes)
    head(filter_barcodes)

    # 2. Check for format differences
    # Sometimes barcodes have different suffixes (-1, -2, etc.)

    # 3. Use more permissive matching
    flexible_barcode_match <- function(raw_barcodes, filter_barcodes) {
      # Strip suffixes for matching
      raw_stripped <- sub("-\\d+$", "", raw_barcodes)
      filter_stripped <- sub("-\\d+$", "", filter_barcodes)

      matches <- raw_barcodes[raw_stripped %in% filter_stripped]
      return(matches)
    }

#### Issue 4: Low Event Detection

Symptom: Very few variable events identified

Solutions:

    # 1. Lower minimum row sum threshold
    variable_events <- find_variable_events(
      m1_matrix, m2_matrix,
      min_row_sum = 10  # Much lower threshold
    )

    # 2. Check data quality
    print("M1 matrix summary:")
    print(summary(rowSums(m1_matrix)))

    print("M2 matrix summary:")
    print(summary(rowSums(m2_matrix)))

    # 3. Examine coordinate groups
    event_data_summary <- event_data[, .N, by = group_count]
    print("Group size distribution:")
    print(event_data_summary)

### Performance Benchmarking

    # Function to benchmark splikit performance
    benchmark_splikit <- function(m1_matrix, m2_matrix) {
      library(microbenchmark)

      cat("=== Splikit Performance Benchmark ===\n")

      # Test different thread counts
      if (check_openmp_enabled()) {
        thread_counts <- c(1, 2, 4, 8)
        thread_counts <- thread_counts[thread_counts <= parallel::detectCores()]

        for (n_threads in thread_counts) {
          cat(sprintf("Testing with %d threads:\n", n_threads))

          time_taken <- system.time({
            result <- find_variable_events(
              m1_matrix, m2_matrix,
              n_threads = n_threads,
              verbose = FALSE
            )
          })

          cat(sprintf("  Time: %.2f seconds\n", time_taken[3]))
          cat(sprintf("  Events found: %d\n", nrow(result)))
        }
      }

      # Memory usage
      cat("\nMemory usage:\n")
      print(pryr::object_size(m1_matrix))
      print(pryr::object_size(m2_matrix))
    }

## Function Reference

### Core Workflow Functions

Function | Purpose | Input | Output
--- | --- | --- | ---
`make_junction_ab()` | Process STARsolo output | Directories, sample IDs | Junction abundance objects
`make_m1()` | Create inclusion matrix | Junction abundance | M1 matrix + metadata
`make_m2()` | Create exclusion matrix | M1 matrix + metadata | M2 matrix
`make_gene_count()` | Process gene expression | 10X directories | Gene expression matrices
`make_velo_count()` | Process velocity data | Velocyto directories | Spliced/unspliced matrices

### Feature Selection Functions

Function | Purpose | Method | Output
--- | --- | --- | ---
`find_variable_events()` | Variable splicing events | Deviance-based | Events + deviance scores
`find_variable_genes()` | Variable genes | VST or deviance | Genes + variability metrics

### Utility Functions

Function | Purpose | Input | Output
--- | --- | --- | ---
`get_pseudo_correlation()` | Pseudo-R² correlation | External data + M1/M2 | Correlation scores
`get_rowVar()` | Row variance | Dense/sparse matrix | Variance vector
`get_silhouette_mean()` | Clustering validation | Data + clusters | Silhouette score

### Data Loading Functions

Function | Purpose | Output
--- | --- | ---
`load_toy_SJ_object()` | Load example junction data | Junction abundance object
`load_toy_M1_M2_object()` | Load example M1/M2 data | M1/M2 matrices + gene expression

### Technical Functions

Function | Purpose | Output
--- | --- | ---
`check_openmp_enabled()` | Check parallel processing | Boolean
C++ functions | Low-level computations | Various

## Conclusion

Splikit provides a comprehensive framework for single-cell splicing analysis, from raw STARsolo output to biological insights. The package’s strength lies in its integrated workflow that handles the complexities of splicing data while providing flexibility for advanced users.

Key advantages of splikit include:

- Complete workflow coverage from data processing to statistical analysis
- Efficient implementation with C++ acceleration for performance-critical functions
- Flexible design that accommodates various experimental setups and data types
- Integration capabilities with popular single-cell analysis frameworks

For additional support, examples, and updates, visit the package documentation and GitHub repository. The splikit development team welcomes feedback and contributions from the community.

This manual was generated for splikit version 2.0.0. For the most up-to-date information, please refer to the package documentation and vignettes.