Data Splitting for Count Data

Introduction

Simulate a raw count matrix with no cluster structure.

set.seed(1)
raw_counts <- matrix(rbinom(n = 100000, size = 10, prob = 0.5), ncol = 100)
rownames(raw_counts) <- paste0("gene", 1:nrow(raw_counts))
colnames(raw_counts) <- paste0("sample", 1:ncol(raw_counts))

Count Splitting (Binomial Thinning)

Split the counts evenly (epsilon = 0.5) by drawing from a Binomial distribution where the size parameter is the original raw count.

epsilon <- 0.5
train_counts <- matrix(
    rbinom(n = length(raw_counts), size = raw_counts, prob = epsilon),
    nrow = nrow(raw_counts),
    ncol = ncol(raw_counts)
)

# Assign gene and sample names to the training set
rownames(train_counts) <- rownames(raw_counts)
colnames(train_counts) <- colnames(raw_counts)

# Define the test matrix as the remainder of the counts
test_counts <- raw_counts - train_counts

1. Discover Clusters on Training Data

Use the train_counts to perform unsupervised learning without “using up” our data for downstream inference. Here, we transpose train_counts to cluster the samples (columns) into 2 groups using K-means.

set.seed(2)
kmeans_res <- kmeans(t(train_counts), centers = 2)
clusters <- kmeans_res$cluster

2. Differential Abundance on Test Data

Perform differential abundance analysis on test_counts using clusters discovered from training_counts without “double dipping.”

p_values <- apply(test_counts, 1, function(x) {
    t.test(x ~ clusters)$p.value
})

# Visualize the distribution of p-values
hist(p_values, main = "Histogram of p-values", xlab = "p-value", col = "lightblue")