Bioinformatics Coding Guide

# OWN-DATA TEMPLATE: use instead of list(preT, postT) integration.
# Input: one saved, independently QC'd and SCT-normalized object per library.
# Change these example paths/labels to your actual design.
library(Seurat)
samples <- read.csv("sample-sheet.csv", stringsAsFactors = FALSE)
required <- c("library_id", "sample_id", "donor_id", "condition", "batch", "counts_rds")
stopifnot(all(required %in% names(samples)), nrow(samples) >= 2)
stopifnot(!anyDuplicated(samples$library_id), all(file.exists(samples$counts_rds)))
stopifnot(all(vapply(samples[required], function(x) all(!is.na(x) & nzchar(trimws(x))), logical(1))))
objects <- lapply(seq_len(nrow(samples)), function(i) {
  obj <- readRDS(samples$counts_rds[i])
  stopifnot(inherits(obj, "Seurat"), all(c("RNA", "SCT") %in% Assays(obj)))
  obj <- RenameCells(obj, add.cell.id = samples$library_id[i])
  obj$library_id <- samples$library_id[i]
  obj$sample_id <- samples$sample_id[i]
  obj$donor_id <- samples$donor_id[i]
  obj$condition <- samples$condition[i]
  obj$batch <- samples$batch[i]
  obj
})
names(objects) <- samples$library_id
features <- SelectIntegrationFeatures(object.list = objects, nfeatures = 3000)
prepared <- PrepSCTIntegration(object.list = objects, anchor.features = features)
anchors <- FindIntegrationAnchors(object.list = prepared,
                                 normalization.method = "SCT",
                                 anchor.features = features)
combined <- IntegrateData(anchorset = anchors, normalization.method = "SCT")
# This replaces the two-object integration block; continue PCA/clustering on combined.
# Preserve RNA counts and all metadata. Evaluate batch correction biologically.
print(table(combined$sample_id, combined$condition))
saveRDS(combined, "integrated_with_sample_metadata.rds")

# OWN-DATA TEMPLATE: run after cell annotation, NOT on the workshop's
# one-sample-per-condition demo. Requires independent biological replication.
# Add installation separately if absent: BiocManager::install("DESeq2")
library(Seurat)
library(DESeq2)
TARGET_CELLTYPE <- "REPLACE_WITH_YOUR_CELL_LABEL"
PAIRED <- TRUE  # TRUE only when each donor contributes both conditions
stopifnot(TARGET_CELLTYPE != "REPLACE_WITH_YOUR_CELL_LABEL")
# Assign your validated annotations to obj40$celltype before this template.
needed <- c("sample_id", "donor_id", "condition", "celltype")
stopifnot(all(needed %in% colnames(obj40[[]])))
selected <- colnames(obj40)[obj40$celltype == TARGET_CELLTYPE & !is.na(obj40$celltype)]
stopifnot(length(selected) > 0)
ct <- subset(obj40, cells = selected)
# With Seurat v5 split RNA layers, join them before accessing counts.
if (inherits(ct[["RNA"]], "Assay5")) ct <- JoinLayers(ct, assay = "RNA")
raw_counts <- GetAssayData(ct, assay = "RNA", layer = "counts")
# Count-based DE needs valid count input. Fractional ambient-corrected counts
# require an explicit reviewed choice of count source/correction method.
values <- if (inherits(raw_counts, "sparseMatrix")) raw_counts@x else as.vector(raw_counts)
stopifnot(all(is.finite(values)), all(values >= 0))
if (any(abs(values - round(values)) > 1e-8)) {
  stop("RNA counts are fractional: review count source; do not round silently.")
}
meta <- unique(ct[[]][, c("sample_id", "donor_id", "condition"), drop = FALSE])
stopifnot(!anyDuplicated(meta$sample_id))
# Cells are aggregated within each specimen, not across all controls/treatments.
pb <- AggregateExpression(ct, assays = "RNA", group.by = "sample_id",
                          return.seurat = TRUE)
counts <- GetAssayData(pb, assay = "RNA", layer = "counts")
stopifnot(!is.null(pb$sample_id), ncol(counts) == length(pb$sample_id))
idx <- match(pb$sample_id, meta$sample_id)
stopifnot(!anyNA(idx))
coldata <- meta[idx, , drop = FALSE]
rownames(coldata) <- colnames(counts)
coldata$condition <- factor(coldata$condition, levels = c("control", "treated"))
stopifnot(!anyNA(coldata$condition))  # replace factor levels for your real groups
coldata$donor_id <- factor(coldata$donor_id)
if (any(tapply(as.character(coldata$condition), coldata$donor_id,
               function(x) anyDuplicated(x) > 0))) {
  stop("Multiple specimens per donor-condition: review repeated-measures design.")
}
if (PAIRED) {
  donor_table <- table(coldata$donor_id, coldata$condition)
  if (any(donor_table != 1)) stop("Paired template requires one specimen per donor per condition.")
  if (nrow(donor_table) < 2) stop("No replicated donor pairs available.")
  design_formula <- ~ donor_id + condition
} else {
  if (any(table(coldata$donor_id) != 1)) stop("Unpaired template requires independent donors.")
  if (any(table(coldata$condition) < 2)) stop("No biological replication in one or more groups.")
  design_formula <- ~ condition
}
# These checks prevent obvious invalid designs; they are not a power calculation.
# If batch and condition are confounded, a model cannot separate their effects.
mm <- model.matrix(design_formula, data = coldata)
stopifnot(qr(mm)$rank == ncol(mm), nrow(mm) > ncol(mm))
keep <- rowSums(counts >= 10) >= 2  # example low-count filter: review for your design
counts <- counts[keep, , drop = FALSE]
stopifnot(nrow(counts) > 0)
dds <- DESeqDataSetFromMatrix(countData = as.matrix(counts),
                              colData = coldata, design = design_formula)
dds <- DESeq(dds)
result <- results(dds, contrast = c("condition", "treated", "control"))
write.csv(as.data.frame(result), "pseudobulk_treated_vs_control.csv")
# Review sample-level QC, cell numbers per specimen, effect sizes and adjusted p-values.
# Missing cell types are not zero-expression pseudobulk replicates.