Single-cell RNA in R
This page displays the Quarto source; it does not run the analysis.
templates/single_cell/guide_r.qmd · Download .qmd · GitHub source
---
title: "Single-cell RNA-seq in R with Seurat"
format:
html:
toc: true
code-fold: false
execute:
warning: false
message: false
---
> Status: draft. Check species-specific gene labels, QC limits, package versions, and study design before use.
Input: a gene-by-cell count matrix and one metadata row per cell. The main path
creates and filters a Seurat object, clusters cells, and finds cluster markers
with MAST. A sample-level pseudobulk path supports condition comparisons; the
remaining analyses are optional.
## 1. Load one count matrix and its cell metadata
The matrix has genes in rows and cells in columns. `metadata.tsv` has one row
per cell and contains `cell_id` and `sample_id`.
```{r}
library(Seurat)
library(qs2)
library(readr)
counts <- qs2::qs_read("input/counts.qs2")
metadata <- readr::read_tsv("input/metadata.tsv", show_col_types = FALSE)
stopifnot("cell_id" %in% names(metadata), "sample_id" %in% names(metadata))
stopifnot(!anyDuplicated(colnames(counts)), !anyDuplicated(metadata$cell_id))
stopifnot(setequal(colnames(counts), metadata$cell_id))
rownames(metadata) <- metadata$cell_id
metadata$cell_id <- NULL
metadata <- metadata[colnames(counts), , drop = FALSE]
object <- Seurat::CreateSeuratObject(
counts = counts,
meta.data = metadata,
min.cells = 3,
min.features = 200
)
object[["percent.mt"]] <- Seurat::PercentageFeatureSet(object, pattern = "^MT-")
print(object)
```
`min.cells` and `min.features` are object-construction filters. They are not
final QC thresholds. Change the mitochondrial pattern for the organism and
gene annotation.
## 2. Inspect and apply cell QC
These example cutoffs are starting values. Inspect each sample before choosing
them; a pooled cutoff can hide sample-specific quality problems.
```{r}
min_features <- 200
max_features <- Inf
min_counts <- 500
max_counts <- Inf
max_percent_mt <- 20
metadata <- object[[]]
metadata$qc_pass <- with(
metadata,
!is.na(nFeature_RNA) &
!is.na(nCount_RNA) &
!is.na(percent.mt) &
nFeature_RNA >= min_features &
nFeature_RNA <= max_features &
nCount_RNA >= min_counts &
nCount_RNA <= max_counts &
percent.mt <= max_percent_mt
)
print(table(metadata$sample_id, metadata$qc_pass, useNA = "ifany"))
VlnPlot(object, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"),
group.by = "sample_id", ncol = 3, pt.size = 0)
object <- subset(object, cells = rownames(metadata)[metadata$qc_pass])
```
## 3. Call doublets by capture library (optional)
Enable this cell after confirming that `sample_id` identifies independent
capture libraries. Review calls by library before removing cells.
```{r}
#| eval: false
library(SingleCellExperiment)
library(scDblFinder)
sce <- Seurat::as.SingleCellExperiment(object, assay = "RNA")
sce <- scDblFinder::scDblFinder(sce, samples = sce$sample_id)
object$doublet_class <- SummarizedExperiment::colData(sce)$scDblFinder.class
object$doublet_score <- SummarizedExperiment::colData(sce)$scDblFinder.score
table(object$doublet_class, object$sample_id, useNA = "ifany")
# Filter only after reviewing per-library calls:
# object <- subset(object, subset = doublet_class == "singlet")
```
## 4. Normalize, reduce dimensions, and cluster
This guide uses log normalization as the default route. SCTransform is a
separate choice; do not run both normalization routes in sequence.
```{r}
object <- Seurat::NormalizeData(object, normalization.method = "LogNormalize",
scale.factor = 10000)
object <- Seurat::FindVariableFeatures(object, selection.method = "vst",
nfeatures = 3000)
object <- Seurat::ScaleData(object, features = Seurat::VariableFeatures(object))
object <- Seurat::RunPCA(object, features = Seurat::VariableFeatures(object),
npcs = 30)
object <- Seurat::FindNeighbors(object, reduction = "pca", dims = 1:30)
object <- Seurat::FindClusters(object, resolution = 0.5, random.seed = 1234)
object <- Seurat::RunUMAP(object, reduction = "pca", dims = 1:30, seed.use = 1234)
DimPlot(object, group.by = "seurat_clusters", label = TRUE)
```
To use SCTransform instead, enable this cell and skip the log-normalization cell above:
```{r}
#| eval: false
object <- Seurat::SCTransform(object, assay = "RNA", vst.flavor = "v2")
object <- Seurat::RunPCA(object, assay = "SCT", npcs = 30)
```
Harmony is optional batch correction. Add it only when the study design and
diagnostics support integration, then use the `harmony` reduction for
neighbors, clustering, and UMAP.
```{r}
#| eval: false
object <- harmony::RunHarmony(object, group.by.vars = "sample_id",
reduction = "pca", dims.use = 1:30)
object <- Seurat::FindNeighbors(object, reduction = "harmony", dims = 1:30)
object <- Seurat::FindClusters(object, resolution = 0.5, random.seed = 1234)
object <- Seurat::RunUMAP(object, reduction = "harmony", dims = 1:30,
reduction.name = "umap_harmony", seed.use = 1234)
```
## 5. Find cluster markers with MAST
This is cell-level marker discovery. MAST adjusts for the listed cell-level
covariates; add other prespecified technical covariates such as `batch` when
available. These covariates do not make cells independent biological
replicates. Use the pseudobulk comparison below for condition-level inference.
```{r}
library(MAST)
Seurat::Idents(object) <- "seurat_clusters"
marker_covariates <- c("nCount_RNA", "percent.mt")
stopifnot(all(marker_covariates %in% names(object[[]])))
markers <- Seurat::FindAllMarkers(
object,
assay = "RNA",
test.use = "MAST",
latent.vars = marker_covariates,
logfc.threshold = 0.25,
only.pos = TRUE
)
head(markers[order(markers$p_val_adj), ], 20)
```
## 6. Test condition differences by sample (optional)
Use one row per biological sample in `sample_metadata.tsv`, with `sample_id`
and `condition` plus any prespecified sample-level covariates. Each sample must
belong to one condition. The example expects `control` and `case`; edit the
factor levels and contrast for other labels. The cell-type comparison below
aggregates raw counts within sample and cluster before testing.
```{r}
#| eval: false
library(edgeR)
library(limma)
sample_metadata <- readr::read_tsv("input/sample_metadata.tsv", show_col_types = FALSE)
stopifnot(all(c("sample_id", "condition") %in% names(sample_metadata)))
stopifnot(!anyDuplicated(sample_metadata$sample_id))
stopifnot(!anyNA(sample_metadata$sample_id), !anyNA(sample_metadata$condition))
stopifnot(all(unique(object$sample_id) %in% sample_metadata$sample_id))
sample_covariates <- setdiff(names(sample_metadata), c("sample_id", "condition"))
sample_metadata <- sample_metadata[match(object$sample_id, sample_metadata$sample_id), , drop = FALSE]
rownames(sample_metadata) <- colnames(object)
object <- Seurat::AddMetaData(
object,
sample_metadata[, setdiff(names(sample_metadata), "sample_id"), drop = FALSE]
)
object$cell_type <- as.character(object$seurat_clusters)
print(table(object$sample_id, object$cell_type))
pseudobulk <- Seurat::AggregateExpression(
object,
assays = "RNA",
return.seurat = TRUE,
group.by = c("sample_id", "condition", "cell_type", sample_covariates)
)
cell_type_to_test <- "0"
pseudobulk <- subset(pseudobulk, subset = cell_type == cell_type_to_test)
pb_counts <- as.matrix(Seurat::GetAssayData(pseudobulk, assay = "RNA", layer = "counts"))
pb_metadata <- pseudobulk[[]][colnames(pb_counts), , drop = FALSE]
stopifnot(all(as.character(pb_metadata$condition) %in% c("control", "case")))
pb_metadata$condition <- factor(pb_metadata$condition, levels = c("control", "case"))
print(table(pb_metadata$condition))
stopifnot(!anyDuplicated(pb_metadata$sample_id))
stopifnot(all(table(pb_metadata$condition) >= 2))
model_covariates <- intersect("batch", sample_covariates)
# For paired samples, add "donor_id" only when the design is paired:
# model_covariates <- c(model_covariates, "donor_id")
design <- stats::model.matrix(
stats::reformulate(c("condition", model_covariates), intercept = FALSE),
data = pb_metadata
)
stopifnot(qr(design)$rank == ncol(design))
dge <- edgeR::DGEList(counts = pb_counts)
keep <- edgeR::filterByExpr(dge, design = design)
stopifnot(any(keep))
dge <- dge[keep, , keep.lib.sizes = FALSE]
dge <- edgeR::normLibSizes(dge)
fit <- edgeR::glmQLFit(dge, design)
contrast <- limma::makeContrasts(conditioncase - conditioncontrol, levels = design)
test <- edgeR::glmQLFTest(fit, contrast = contrast)
pseudobulk_results <- edgeR::topTags(test, n = nrow(dge))$table
pseudobulk_results$gene_id <- rownames(pseudobulk_results)
head(pseudobulk_results[order(pseudobulk_results$FDR), ], 20)
dir.create("output", recursive = TRUE, showWarnings = FALSE)
readr::write_tsv(pseudobulk_results, "output/pseudobulk_rna_de.tsv")
```
Cells are not biological replicates. Two samples per condition is only a
minimal code check, not evidence of adequate replication. Check sample counts
and the design before interpreting this result; a full-rank design cannot
resolve a condition completely confounded with batch. Add a paired donor term
only for a paired design. MAST markers above remain useful for describing
clusters.
## 7. Score one directional gene signature (optional)
The TSV has `gene` and `direction` columns; direction is `UP` or `DOWN`.
```{r}
#| eval: false
library(UCell)
signature <- readr::read_tsv("input/signature.tsv", show_col_types = FALSE)
genes_up <- intersect(signature$gene[signature$direction == "UP"],
rownames(object[["RNA"]]))
genes_down <- intersect(signature$gene[signature$direction == "DOWN"],
rownames(object[["RNA"]]))
stopifnot(length(genes_up) > 0, length(genes_down) > 0)
object <- UCell::AddModuleScore_UCell(
object,
features = list(signature_UP = genes_up, signature_DOWN = genes_down),
name = NULL
)
object$signature_score <- object$signature_UP_UCell - object$signature_DOWN_UCell
FeaturePlot(object, features = "signature_score")
```
## 8. Estimate cell potency from raw counts (optional)
CytoTRACE2 uses raw RNA counts. Set `species` to match the data. Its score is a
potency estimate, not a lineage trace or a replacement for trajectory inference.
```{r}
#| eval: false
library(CytoTRACE2)
counts <- Seurat::GetAssayData(object, assay = "RNA", layer = "counts")
cytotrace <- CytoTRACE2::cytotrace2(counts, is_seurat = FALSE,
species = "human", ncores = 1)
object$cytotrace2_score <- cytotrace$CytoTRACE2_Score[colnames(object)]
object$cytotrace2_potency <- cytotrace$CytoTRACE2_Potency[colnames(object)]
FeaturePlot(object, features = "cytotrace2_score", reduction = "umap")
```
## 9. Infer pseudotime, then test expression along it (optional)
Choose the starting cluster from the study design. Slingshot infers curves
from the PCA coordinates and cluster labels. tradeSeq then tests count changes
along the selected curve; neither step proves cell lineage or causal direction.
```{r}
#| eval: false
library(SingleCellExperiment)
library(slingshot)
sce <- Seurat::as.SingleCellExperiment(object, assay = "RNA")
reducedDim(sce, "PCA") <- Seurat::Embeddings(object, "pca")
sce <- slingshot::slingshot(sce, clusterLabels = "seurat_clusters",
reducedDim = "PCA", start.clus = "0")
pseudotime <- slingshot::slingPseudotime(sce)[, 1]
weights <- slingshot::slingCurveWeights(sce)[, 1]
object$pseudotime_1 <- pseudotime[colnames(object)]
object$lineage_weight_1 <- weights[colnames(object)]
```
```{r}
#| eval: false
library(tradeSeq)
counts <- Seurat::GetAssayData(object, assay = "RNA", layer = "counts")
fit <- tradeSeq::fitGAM(
counts = counts,
pseudotime = matrix(object$pseudotime_1, ncol = 1,
dimnames = list(colnames(object), "lineage_1")),
cellWeights = matrix(object$lineage_weight_1, ncol = 1,
dimnames = list(colnames(object), "lineage_1")),
nknots = 5,
verbose = FALSE
)
trajectory_results <- tradeSeq::associationTest(fit)
head(trajectory_results[order(trajectory_results$pvalue), ], 20)
```
## 10. Infer co-expression modules in one population (optional)
hdWGCNA is a network analysis, not a score for a supplied gene set. Choose
the target population for the network; pooling distinct cell types can make
modules reflect cell composition. See the [hdWGCNA tutorial](https://smorabit.github.io/hdWGCNA/articles/basic_tutorial.html).
```{r}
#| eval: false
library(hdWGCNA)
object <- hdWGCNA::SetupForWGCNA(object, features = Seurat::VariableFeatures(object),
wgcna_name = "wgcna")
object <- hdWGCNA::MetacellsByGroups(object, group.by = c("sample_id", "seurat_clusters"),
reduction = "pca", k = 25, max_shared = 10)
object <- hdWGCNA::NormalizeMetacells(object)
network_group <- "0"
stopifnot(network_group %in% as.character(object$seurat_clusters))
object <- hdWGCNA::SetDatExpr(object,
group_name = network_group,
group.by = "seurat_clusters", assay = "RNA",
layer = "data")
object <- hdWGCNA::TestSoftPowers(object, networkType = "signed hybrid")
hdWGCNA::PlotSoftPowers(object)
object <- hdWGCNA::ConstructNetwork(object, soft_power = NULL,
networkType = "signed hybrid",
randomSeed = 1234)
object <- hdWGCNA::ModuleEigengenes(object, assay = "RNA")
object <- hdWGCNA::ModuleConnectivity(object, group.by = "seurat_clusters",
group_name = network_group)
```
## 11. Save the object
The object carries the retained cells and any optional results added above.
```{r}
dir.create("output", recursive = TRUE, showWarnings = FALSE)
qs2::qs_save(object, "output/seurat_analysis.qs2")
```