pbmc5k_donor4_cluster_de

# Build Tessera's PBMC cluster differential-expression toy directly from a
# 10x Genomics filtered feature-barcode matrix directory:
#   barcodes.tsv.gz · features.tsv.gz · matrix.mtx.gz
#
# Usage:
#   Rscript --vanilla pbmc5k_donor4_cluster_de.R <10x-directory> [output.csv]
#
# Clusters 0, 1 and 2 are compared independently with the same reference,
# cluster 3. The upstream QC, normalization, PCA and Louvain clustering are
# intentionally identical to pbmc5k_donor4_embeddings.R.

args <- commandArgs(trailingOnly = TRUE)
if (length(args) < 1L) {
  stop("Provide the 10x filtered feature-barcode matrix directory", call. = FALSE)
}

input <- args[[1L]]
output <- if (length(args) >= 2L) args[[2L]] else "pbmc5k_donor4_cluster_de.csv"
required_files <- c("barcodes.tsv.gz", "features.tsv.gz", "matrix.mtx.gz")
missing_files <- required_files[!file.exists(file.path(input, required_files))]
if (length(missing_files)) {
  stop("Missing 10x files: ", paste(missing_files, collapse = ", "), call. = FALSE)
}

set.seed(1)
counts <- Seurat::Read10X(
  data.dir = input,
  gene.column = 2,
  cell.column = 1,
  unique.features = TRUE,
  strip.suffix = FALSE
)
pbmc <- Seurat::CreateSeuratObject(
  counts = counts,
  project = "pbmc5k_donor4",
  min.cells = 3,
  min.features = 40
)

pbmc[["percent_mito"]] <- Seurat::PercentageFeatureSet(pbmc, pattern = "^MT-")
pbmc <- subset(
  pbmc,
  subset = nCount_RNA >= 500 & nCount_RNA <= 15000 &
    nFeature_RNA >= 200 & nFeature_RNA <= 5000 & percent_mito < 10
)
pbmc <- Seurat::NormalizeData(
  pbmc,
  normalization.method = "LogNormalize",
  scale.factor = 10000,
  verbose = FALSE
)
pbmc <- Seurat::FindVariableFeatures(
  pbmc,
  selection.method = "vst",
  nfeatures = 2000,
  verbose = FALSE
)
pbmc <- Seurat::ScaleData(
  pbmc,
  features = Seurat::VariableFeatures(pbmc),
  verbose = FALSE
)
pbmc <- Seurat::RunPCA(
  pbmc,
  features = Seurat::VariableFeatures(pbmc),
  npcs = 30,
  seed.use = 1,
  verbose = FALSE
)
pbmc <- Seurat::FindNeighbors(
  pbmc,
  reduction = "pca",
  dims = 1:20,
  k.param = 20,
  verbose = FALSE
)
pbmc <- Seurat::FindClusters(
  pbmc,
  resolution = 0.8,
  algorithm = 1,
  random.seed = 1,
  verbose = FALSE
)
SeuratObject::Idents(pbmc) <- "seurat_clusters"

target_clusters <- c("0", "1", "2")
reference_cluster <- "3"
cluster_counts <- table(SeuratObject::Idents(pbmc))
required_clusters <- c(target_clusters, reference_cluster)
if (!all(required_clusters %in% names(cluster_counts))) {
  stop("Expected clusters 0, 1, 2 and 3 were not reproduced", call. = FALSE)
}

results <- lapply(target_clusters, function(target_cluster) {
  markers <- Seurat::FindMarkers(
    pbmc,
    ident.1 = target_cluster,
    ident.2 = reference_cluster,
    assay = "RNA",
    slot = "data",
    test.use = "wilcox",
    logfc.threshold = 0,
    min.pct = 0.01,
    only.pos = FALSE,
    verbose = FALSE
  )
  markers$gene <- rownames(markers)
  comparison <- paste0("cluster_", target_cluster, "_vs_", reference_cluster)
  data.frame(
    result_id = paste(comparison, markers$gene, sep = "::"),
    comparison = comparison,
    target_cluster = as.integer(target_cluster),
    reference_cluster = as.integer(reference_cluster),
    n_target = unname(cluster_counts[[target_cluster]]),
    n_reference = unname(cluster_counts[[reference_cluster]]),
    gene = markers$gene,
    avg_log2fc = markers$avg_log2FC,
    pct_target = markers$pct.1,
    pct_reference = markers$pct.2,
    p_value = markers$p_val,
    p_adj = markers$p_val_adj,
    row.names = NULL,
    check.names = FALSE
  )
})

frozen <- do.call(rbind, results)
frozen <- frozen[order(frozen$target_cluster, frozen$p_adj, -abs(frozen$avg_log2fc), frozen$gene), ]
write.csv(frozen, output, row.names = FALSE)