PBMCs from 8 donors, control vs interferon-beta (Kang et al. 2018)
This is one complete single-cell analysis on a laptop, from a count matrix to pathways. Every step says what we do and why we do it. Most Seurat tutorials cover one step each, so you have to stitch them together yourself. Here everything is in one place and in order.
The data. Blood cells (PBMCs) from 8 people. Each person’s cells were split in two: one half left alone (ctrl) and one half treated with interferon-beta (stim). That gives us 8 control samples and 8 treated samples, and every treated sample has a matching control from the same person. It was run on 10x Genomics Chromium.
The question. Which genes and pathways change in each cell type when cells are treated with interferon-beta?
The steps.
Want to run it yourself? The setup steps cover installing R and RStudio, the packages, and running this file. No experience needed.
library(Seurat)
library(SingleCellExperiment)
library(muscData) # the example data
library(celda) # decontX, for ambient RNA
library(scDblFinder) # doublets
library(scuttle) # QC helper
library(harmony) # integration
library(DESeq2) # differential expression
library(fgsea) # pathway analysis
library(msigdbr) # gene sets
library(ggplot2)
library(dplyr)
set.seed(42) # so you get the same results every time
dir.create("results", showWarnings = FALSE)What: download the count matrix. Rows are genes, columns are cells, and each value is how many RNA molecules of that gene were counted in that cell.
Why: everything else starts from these raw counts. The first download takes a minute, and after that the file is cached.
sce <- Kang18_8vs8()
sceclass: SingleCellExperiment
dim: 35635 29065
metadata(0):
assays(1): counts
rownames(35635): MIR1302-10 FAM138A ... MT-ND6 MT-CYB
rowData names(2): ENSEMBL SYMBOL
colnames(29065): AAACATACAATGCC-1 AAACATACATTTCC-1 ... TTTGCATGGTTTGG-1
TTTGCATGTCTTAC-1
colData names(5): ind stim cluster cell multiplets
reducedDimNames(1): TSNE
mainExpName: NULL
altExpNames(0):A few genes appear more than once, for example CD79A_ENSG00000105369 and CD79A_ENSG00000272985. The extra copies come from alternative versions of some chromosome regions and have almost no counts. We keep the copy with the most counts for each gene and give it its plain name, so that names like CD79A work in plots and match the pathway gene sets later.
total_counts <- rowSums(counts(sce))
sce <- sce[order(total_counts, decreasing = TRUE), ]
sce <- sce[!duplicated(rowData(sce)$SYMBOL), ]
rownames(sce) <- rowData(sce)$SYMBOL
nrow(sce)[1] 32938The table that comes with the cells tells us who each cell came from (ind), whether it was treated (stim), and what the original authors called it (cell). It also has a multiplets column, which we come back to in step 3.
# how many cells per person and condition
table(sce$ind, sce$stim)
ctrl stim
101 1027 1363
107 645 641
1015 3138 2703
1016 2180 1998
1039 552 797
1244 2250 1755
1256 2451 2275
1488 2376 2914We make one sample name per person and condition, like donor101-ctrl. This is the unit we compare later.
sce$sample <- paste0("donor", sce$ind, "-", sce$stim)In this experiment all control cells were run in one 10x lane and all treated cells in another. A lane is called a capture. Ambient RNA and doublets happen inside one capture, so steps 2 to 4 are done separately for each capture. Here the capture is the same as stim.
What: when cells are loaded, some float-around RNA from broken cells ends up in every droplet. This is called ambient RNA or soup. We estimate how much each cell has and subtract it.
Why: without this, genes from very common cells show up in every cell type. For example, monocyte genes can appear in T cells. This makes clusters look mixed and can create false differences later.
Which tool: if you have the raw_feature_bc_matrix from Cell Ranger (all droplets, including empty ones), SoupX is the usual choice. This dataset only has the cells, with no empty droplets, so we use decontX. It works out the contamination from the cells alone.
sce <- decontX(sce, batch = sce$stim)
# contamination per cell, from 0 (none) to 1 (all soup)
summary(sce$decontX_contamination) Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000771 0.0094743 0.0276354 0.0709858 0.0789133 0.8500350 ggplot(as.data.frame(colData(sce)), aes(decontX_contamination)) +
geom_histogram(bins = 50) +
facet_wrap(~stim) +
labs(x = "estimated contamination", y = "cells") +
theme_bw()
Most cells have a small amount of contamination. We keep the corrected counts, rounded to whole numbers because later steps expect counts.
corrected <- round(decontXcounts(sce))What: sometimes two cells land in the same droplet. That droplet looks like one cell that expresses the genes of both, for example a T cell and a B cell at once. These are doublets.
Why: doublets can form fake clusters that look like a new cell type. We want to remove them before clustering.
How: scDblFinder makes fake doublets by adding random pairs of cells together, then learns what doublets look like and scores every real cell. It runs on each capture separately because doublets can only form inside one capture.
sce_clean <- SingleCellExperiment(list(counts = corrected), colData = colData(sce))
sce_clean <- scDblFinder(sce_clean, samples = "stim")
|
| | 0%
|
|=================================== | 50%
|
|======================================================================| 100%table(sce_clean$scDblFinder.class)
singlet doublet
25759 3306 This dataset has a bonus. The authors mixed cells from 8 people in each lane and used each person’s DNA differences (with a tool called demuxlet) to find droplets containing two people. That is in the multiplets column. Let’s compare the two.
table(scDblFinder = sce_clean$scDblFinder.class, demuxlet = sce_clean$multiplets) demuxlet
scDblFinder ambs doublet singlet
singlet 1150 1507 23102
doublet 67 1662 1577demuxlet can only see doublets made of two different people, and scDblFinder can only see doublets made of two different cell types, so they catch different droplets. We remove a cell if either tool calls it a doublet. demuxlet’s ambs (ambiguous) cells are also removed, because we can’t tell which person they came from.
What: remove cells that look broken or empty.
Why: a cell with very few genes is probably a broken cell or an empty droplet. A cell with a lot of mitochondrial RNA is often a dying cell that has lost most of its other RNA. These cells cluster together and look like a separate cell type when they are not.
First we put the data into Seurat, which we use for the rest of the analysis. We also keep the uncorrected counts as a second assay, so you can compare later.
obj <- CreateSeuratObject(counts = corrected, meta.data = as.data.frame(colData(sce_clean)))
obj[["uncorrected"]] <- CreateAssay5Object(counts = counts(sce))
obj[["percent.mt"]] <- PercentageFeatureSet(obj, pattern = "^MT-")
objAn object of class Seurat
65876 features across 29065 samples within 2 assays
Active assay: RNA (32938 features, 0 variable features)
1 layer present: counts
1 other assay present: uncorrectedThe three numbers we look at:
nCount_RNA: total RNA molecules in the cellnFeature_RNA: number of different genes found in the cellpercent.mt: percent of RNA from mitochondrial genesVlnPlot(obj, features = c("nCount_RNA", "nFeature_RNA", "percent.mt"),
group.by = "stim", pt.size = 0, ncol = 3)
How to pick cutoffs. Don’t copy numbers from another dataset. Instead, flag cells that are far from the middle of their own capture. “Far” here means more than 3 MADs (median absolute deviations, a robust version of the standard deviation). We use log values because counts are skewed.
For mitochondrial percent we only flag the high side. We also never set that cutoff below 10%, because in healthy blood cells mito is low and a pure MAD rule would throw out good cells.
low_counts <- isOutlier(obj$nCount_RNA, nmads = 3, log = TRUE, type = "lower", batch = obj$stim)
high_counts <- isOutlier(obj$nCount_RNA, nmads = 3, log = TRUE, type = "higher", batch = obj$stim)
low_genes <- isOutlier(obj$nFeature_RNA, nmads = 3, log = TRUE, type = "lower", batch = obj$stim)
high_mito <- isOutlier(obj$percent.mt, nmads = 3, type = "higher", batch = obj$stim)
high_mito <- high_mito & obj$percent.mt > 10
obj$qc_fail <- low_counts | high_counts | low_genes | high_mito
obj$doublet <- obj$scDblFinder.class == "doublet" | obj$multiplets != "singlet"How many cells each rule removes, per capture:
qc_table <- data.frame(
capture = c("ctrl", "stim"),
cells = as.vector(table(obj$stim)),
low_counts = as.vector(tapply(low_counts, obj$stim, sum)),
high_counts = as.vector(tapply(high_counts, obj$stim, sum)),
low_genes = as.vector(tapply(low_genes, obj$stim, sum)),
high_mito = as.vector(tapply(high_mito, obj$stim, sum)),
doublets = as.vector(tapply(obj$doublet, obj$stim, sum)),
kept = as.vector(tapply(!obj$qc_fail & !obj$doublet, obj$stim, sum))
)
qc_table capture cells low_counts high_counts low_genes high_mito doublets kept
1 ctrl 14619 6 48 155 0 3103 11417
2 stim 14446 10 32 145 0 2860 11497write.csv(qc_table, "results/qc_table.csv", row.names = FALSE)obj$removed <- obj$qc_fail | obj$doublet
ggplot(obj@meta.data, aes(nCount_RNA, nFeature_RNA, colour = removed)) +
geom_point(size = 0.2, alpha = 0.4) +
scale_x_log10() + scale_y_log10() +
facet_wrap(~stim) +
theme_bw()
Now remove them.
obj <- subset(obj, subset = removed == FALSE)
objAn object of class Seurat
65876 features across 22914 samples within 2 assays
Active assay: RNA (32938 features, 0 variable features)
1 layer present: counts
1 other assay present: uncorrectedWhat: make cells comparable. Some cells were sequenced more deeply than others, so their counts are all higher. Normalizing divides each cell by its total, multiplies by 10,000, and takes the log.
Why: without this, the biggest difference between cells would just be how deeply each one was sequenced.
We then pick the 2,000 most variable genes. These carry most of the information about cell types, and using only them makes the next steps faster and less noisy. Scaling puts every gene on the same range so that very highly expressed genes don’t dominate. PCA then squeezes those 2,000 genes into 30 summary dimensions.
obj <- NormalizeData(obj)
obj <- FindVariableFeatures(obj, nfeatures = 2000)
obj <- ScaleData(obj)
obj <- RunPCA(obj, npcs = 50)The elbow plot shows how much each PCA dimension explains. Past about 20 to 30 the line flattens, so 30 is a safe choice.
ElbowPlot(obj, ndims = 50)
What: line up the two captures so the same cell type from each sits together.
Why: interferon changes the expression of hundreds of genes in every cell. On top of that, the two conditions were run in different lanes. So before integration, a treated T cell and a control T cell end up far apart, and we can’t name cell types across both.
First look without integration. Always do this first. If cells already mix, you don’t need integration.
obj <- RunUMAP(obj, dims = 1:30, reduction.name = "umap_unintegrated")
DimPlot(obj, reduction = "umap_unintegrated", group.by = "stim") + ggtitle("Before integration")
Each cell type splits into a control island and a treated island. We use Harmony, which nudges the PCA dimensions until the two groups mix within each cell type.
obj <- RunHarmony(obj, group.by.vars = "stim")
obj <- RunUMAP(obj, reduction = "harmony", dims = 1:30)
DimPlot(obj, reduction = "umap", group.by = "stim") + ggtitle("After Harmony")
Important: Harmony only changes the embedding (the coordinates used for clustering and UMAP). The counts stay untouched. We never use integrated values to test genes.
What: group cells that have similar expression.
Why: clusters are the starting point for naming cell types.
FindNeighbors links each cell to the cells most like it. FindClusters cuts that network into groups. The resolution controls how many groups you get: higher means more, smaller clusters. 0.5 is a common start for PBMCs.
obj <- FindNeighbors(obj, reduction = "harmony", dims = 1:30)
obj <- FindClusters(obj, resolution = 0.5)Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 22914
Number of edges: 797751
Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9038
Number of communities: 13
Elapsed time: 2 secondsDimPlot(obj, reduction = "umap", label = TRUE) + NoLegend()
Check that every cluster has cells from both conditions and from all donors. A cluster from only one person deserves a closer look.
table(obj$seurat_clusters, obj$stim)
ctrl stim
0 2776 2738
1 1977 1938
2 1651 1730
3 992 1036
4 1023 1003
5 826 806
6 578 608
7 503 556
8 426 403
9 221 252
10 248 202
11 190 190
12 6 35What: for each cluster, find genes that are higher there than in all other clusters.
Why: marker genes tell you what each cluster is. For example, CD3E marks T cells and MS4A1 marks B cells.
We use the normalized RNA values from step 5 and the Wilcoxon test. The p-values here only describe the clusters. Each cell is not an independent sample, so don’t use them as real statistics.
obj <- JoinLayers(obj)
markers <- FindAllMarkers(obj, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.5)
write.csv(markers, "results/cluster_markers.csv", row.names = FALSE)
# top 5 per cluster
top5 <- markers %>% group_by(cluster) %>% slice_max(avg_log2FC, n = 5)
top5# A tibble: 65 × 7
# Groups: cluster [13]
p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene
<dbl> <dbl> <dbl> <dbl> <dbl> <fct> <chr>
1 0 9.83 0.336 0.003 0 0 CCL7
2 0 7.12 0.717 0.023 0 0 CCL2
3 0 7.11 0.451 0.013 0 0 CCL8
4 0 6.83 0.663 0.012 0 0 S100A8
5 0 6.73 0.689 0.016 0 0 S100A9
6 0 2.38 0.273 0.055 0 1 LEF1
7 0 2.13 0.699 0.201 0 1 SELL
8 1.56e-281 1.98 0.322 0.094 5.15e-277 1 AES
9 0 1.92 0.725 0.229 0 1 LTB
10 2.47e-228 1.89 0.284 0.083 8.14e-224 1 GIMAP1
# ℹ 55 more rowsWell-known PBMC marker genes, as a dot plot. Bigger dots mean more cells express the gene, and darker means higher expression.
known_markers <- c("CD3E", "IL7R", "CCR7", "CD8A", "GNLY", "NKG7",
"MS4A1", "CD79A", "CD14", "LYZ", "FCGR3A",
"FCER1A", "CST3", "PPBP")
DotPlot(obj, features = known_markers) + RotatedAxis()
What: give each cluster a cell type name.
Why: the rest of the analysis is done per cell type.
How: first make an automatic guess from known marker genes, then check it by eye.
We write down a short list of well-known marker genes for each blood cell type. AddModuleScore gives every cell a score for each list: how much higher those genes are in that cell than random genes with similar expression. Each cell gets the name of its highest-scoring list.
cell_type_markers <- list(
"CD14+ Monocytes" = c("CD14", "LYZ", "S100A8", "S100A9"),
"FCGR3A+ Monocytes" = c("FCGR3A", "MS4A7"),
"CD4+ T cells" = c("CD3D", "IL7R"),
"CD8+ T cells" = c("CD8A", "CD8B"),
"NK cells" = c("GNLY", "NKG7", "KLRF1"),
"B cells" = c("MS4A1", "CD79A"),
"Dendritic cells" = c("FCER1A", "CD1C"),
"pDC" = c("LILRA4", "CLEC4C"),
"Megakaryocytes" = c("PPBP", "PF4")
)
obj <- AddModuleScore(obj, features = cell_type_markers, name = "score")
# the scores are saved as score1, score2, ... in the same order as the list
scores <- obj@meta.data[, paste0("score", 1:length(cell_type_markers))]
colnames(scores) <- names(cell_type_markers)
# each cell gets the name with the highest score
obj$marker_guess <- colnames(scores)[apply(scores, 1, which.max)]Now take the most common guess in each cluster.
label_table <- table(obj$seurat_clusters, obj$marker_guess)
proposed <- data.frame(
cluster = rownames(label_table),
n_cells = rowSums(label_table),
marker_guess = colnames(label_table)[apply(label_table, 1, which.max)],
agreement = round(apply(label_table, 1, max) / rowSums(label_table), 2)
)
proposed$celltype <- proposed$marker_guess
proposed cluster n_cells marker_guess agreement celltype
0 0 5514 CD14+ Monocytes 0.70 CD14+ Monocytes
1 1 3915 CD4+ T cells 0.74 CD4+ T cells
2 2 3381 CD4+ T cells 0.66 CD4+ T cells
3 3 2028 NK cells 0.46 NK cells
4 4 2026 B cells 0.84 B cells
5 5 1632 NK cells 0.98 NK cells
6 6 1186 CD4+ T cells 0.56 CD4+ T cells
7 7 1059 FCGR3A+ Monocytes 0.98 FCGR3A+ Monocytes
8 8 829 CD4+ T cells 0.49 CD4+ T cells
9 9 473 Megakaryocytes 0.63 Megakaryocytes
10 10 450 CD14+ Monocytes 0.52 CD14+ Monocytes
11 11 380 B cells 0.58 B cells
12 12 41 pDC 0.71 pDCwrite.csv(proposed, "results/annotation_proposed.csv", row.names = FALSE)agreement is the share of cells in the cluster with that guess. Low agreement means the cluster is mixed or the marker lists don’t separate it well, so check the dot plot and the top markers from step 8.
An automatic guess is only as good as its marker lists. A cluster with markers that aren’t on any list, for example stressed cells with high heat-shock genes, still gets the closest name. That is why a person checks every label.
There are other ways to make the automatic guess. Reference-based tools such as SingleR or Azimuth compare your cells to a labelled dataset. They need extra packages and a reference download, and they make mistakes too.
Your turn. Open results/annotation_proposed.csv, check each label against the dot plot and markers, change the celltype column if needed, and save it as annotation_approved.csv next to this file. If that file exists, the notebook uses it.
if (file.exists("annotation_approved.csv")) {
annotation <- read.csv("annotation_approved.csv")
} else {
annotation <- proposed
}
obj$celltype <- annotation$celltype[match(obj$seurat_clusters, annotation$cluster)]As a check, compare with the labels from the original paper.
table(ours = obj$celltype, paper = obj$cell) paper
ours B cells CD14+ Monocytes CD4 T cells CD8 T cells
B cells 2343 1 3 0
CD14+ Monocytes 0 5049 1 0
CD4+ T cells 37 11 7875 151
CD8+ T cells 0 0 1352 1594
Dendritic cells 0 128 0 0
FCGR3A+ Monocytes 0 28 0 0
Megakaryocytes 32 0 235 39
NK cells 0 0 3 64
pDC 2 0 0 0
paper
ours Dendritic cells FCGR3A+ Monocytes Megakaryocytes NK cells
B cells 35 20 3 1
CD14+ Monocytes 9 422 27 6
CD4+ T cells 7 15 16 13
CD8+ T cells 0 0 1 267
Dendritic cells 279 43 0 0
FCGR3A+ Monocytes 0 1028 2 1
Megakaryocytes 1 2 136 28
NK cells 0 1 0 1564
pDC 38 1 0 0DimPlot(obj, reduction = "umap", group.by = "celltype", label = TRUE, repel = TRUE) + NoLegend() +
DimPlot(obj, reduction = "umap", group.by = "stim")
What: for each cell type, find genes that go up or down with interferon.
Why pseudobulk: the replicates in this experiment are the 8 people, not the thousands of cells. Cells from the same person are more alike than cells from different people. If you test cell by cell, you pretend you have thousands of replicates and get far too many “significant” genes.
How: for each cell type, add up the raw counts of all cells from each sample. That gives one column per sample, like a small bulk RNA-seq experiment. Then run DESeq2, the standard bulk tool.
Because each person has both a control and a treated sample, the design is ~ ind + stim. This compares treated to control within each person, which removes differences between people.
celltypes <- unique(obj$celltype)
de_results <- list()
for (ct in celltypes) {
cells <- subset(obj, subset = celltype == ct)
# only use samples with at least 10 cells of this type
cell_counts <- table(cells$sample)
keep_samples <- names(cell_counts)[cell_counts >= 10]
# need at least 2 samples in each condition, or there is nothing to compare
n_ctrl <- sum(grepl("-ctrl$", keep_samples))
n_stim <- sum(grepl("-stim$", keep_samples))
if (n_ctrl < 2 | n_stim < 2) {
cat(paste0("Skipping ", ct, ": only ", n_ctrl, " control and ", n_stim,
" treated samples have 10 or more cells.\n"))
next
}
cells <- subset(cells, subset = sample %in% keep_samples)
# add up counts per sample
pb <- AggregateExpression(cells, assays = "RNA", group.by = "sample")$RNA
pb <- round(as.matrix(pb))
# sample table: which person and condition each column is
info <- data.frame(sample = colnames(pb))
info$ind <- factor(cells$ind[match(info$sample, cells$sample)])
info$stim <- factor(cells$stim[match(info$sample, cells$sample)], levels = c("ctrl", "stim"))
# drop genes with almost no counts
pb <- pb[rowSums(pb >= 10) >= 2, ]
dds <- DESeqDataSetFromMatrix(pb, colData = info, design = ~ ind + stim)
dds <- DESeq(dds, quiet = TRUE)
res <- as.data.frame(results(dds, contrast = c("stim", "stim", "ctrl")))
res$gene <- rownames(res)
res$celltype <- ct
de_results[[ct]] <- res
}Skipping pDC: only 0 control and 0 treated samples have 10 or more cells.de_all <- bind_rows(de_results)
write.csv(de_all, "results/pseudobulk_DE.csv", row.names = FALSE)How many genes change in each cell type (adjusted p-value under 0.05):
de_summary <- de_all %>%
group_by(celltype) %>%
summarise(genes_tested = n(), significant = sum(padj < 0.05, na.rm = TRUE)) %>%
arrange(desc(significant))
de_summary# A tibble: 8 × 3
celltype genes_tested significant
<chr> <int> <int>
1 CD14+ Monocytes 5395 3140
2 FCGR3A+ Monocytes 2622 1325
3 CD4+ T cells 6406 925
4 Dendritic cells 1666 925
5 B cells 3857 688
6 CD8+ T cells 4488 481
7 NK cells 2588 469
8 Megakaryocytes 875 113Interferon is a strong signal, so expect many significant genes. Known interferon genes like ISG15, IFI6 and MX1 should be near the top.
de_all %>% filter(celltype == de_summary$celltype[1]) %>% arrange(padj) %>% head(10) baseMean log2FoldChange lfcSE stat pvalue
RPL7 658.5996 -2.159763 0.08081148 -26.72594 2.351751e-157
CCL8 7025.4913 9.964247 0.37702120 26.42888 6.383319e-154
RPL6 387.3201 -2.633680 0.10085513 -26.11349 2.562012e-150
ISG20 2884.5713 6.009875 0.23034865 26.09034 4.692700e-150
IL1RN 877.8711 6.722971 0.26839444 25.04885 1.797067e-138
SSB 347.6087 3.624881 0.14989222 24.18325 3.338337e-129
TMSB10 2429.2409 2.299250 0.09606988 23.93310 1.385970e-126
IFITM2 471.5635 4.951345 0.21324804 23.21871 2.946911e-119
DYNLT1 901.9296 2.753866 0.11880313 23.18008 7.232393e-119
NT5C3A 208.5972 6.258726 0.27367520 22.86918 9.419615e-116
padj gene celltype
RPL7 1.268769e-153 RPL7 CD14+ Monocytes
CCL8 1.721900e-150 CCL8 CD14+ Monocytes
RPL6 4.607351e-147 RPL6 CD14+ Monocytes
ISG20 6.329279e-147 ISG20 CD14+ Monocytes
IL1RN 1.939035e-135 IL1RN CD14+ Monocytes
SSB 3.001721e-126 SSB CD14+ Monocytes
TMSB10 1.068187e-123 TMSB10 CD14+ Monocytes
IFITM2 1.987323e-116 IFITM2 CD14+ Monocytes
DYNLT1 4.335418e-116 DYNLT1 CD14+ Monocytes
NT5C3A 5.081882e-113 NT5C3A CD14+ MonocytesWhat: find groups of genes (pathways) that move together.
Why: a list of hundreds of genes is hard to read. Pathways turn it into a few ideas, like “interferon response is up”.
How: GSEA (with the fgsea package). Rank every tested gene by the DESeq2 statistic, from most up to most down, and ask whether the genes of a pathway pile up at the top or bottom. It uses all genes, not just the significant ones, so it also works when few genes pass the cutoff.
Gene sets come from MSigDB: Hallmark (50 broad, well-curated sets) and Reactome (detailed pathways).
hallmark <- msigdbr(species = "human", collection = "H")
reactome <- msigdbr(species = "human", collection = "C2", subcollection = "CP:REACTOME")
genesets <- rbind(hallmark, reactome)
genesets <- split(genesets$gene_symbol, genesets$gs_name)
length(genesets)[1] 1889gsea_results <- list()
for (ct in names(de_results)) {
res <- de_results[[ct]]
res <- res[!is.na(res$stat), ]
ranks <- res$stat
names(ranks) <- res$gene
ranks <- sort(ranks, decreasing = TRUE)
gs <- fgsea(genesets, ranks, minSize = 15, maxSize = 500)
gs$celltype <- ct
gsea_results[[ct]] <- gs
}
gsea_all <- bind_rows(gsea_results)
gsea_all$leadingEdge <- sapply(gsea_all$leadingEdge, paste, collapse = ";")
write.csv(gsea_all, "results/gsea.csv", row.names = FALSE)Top pathways for the cell type with the most significant genes. Positive NES means the pathway goes up with interferon.
top_ct <- de_summary$celltype[1]
top_paths <- gsea_all %>% filter(celltype == top_ct) %>% arrange(padj) %>% head(20)
ggplot(top_paths, aes(NES, reorder(pathway, NES), fill = padj < 0.05)) +
geom_col() +
labs(y = NULL, title = paste("Top pathways, stim vs ctrl:", top_ct)) +
theme_bw() +
theme(axis.text.y = element_text(size = 7))
HALLMARK_INTERFERON_ALPHA_RESPONSE and HALLMARK_INTERFERON_GAMMA_RESPONSE should be at the top. If they are, the whole pipeline worked.
saveRDS(obj, "results/pbmc_ifnb_annotated.rds")
sessionInfo()R version 4.4.2 (2024-10-31)
Platform: aarch64-apple-darwin20
Running under: macOS 26.6.2
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
locale:
[1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
time zone: America/New_York
tzcode source: internal
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] future_1.70.0 dplyr_1.2.0
[3] ggplot2_4.0.2 msigdbr_26.1.1
[5] fgsea_1.32.4 DESeq2_1.46.0
[7] harmony_1.2.4 Rcpp_1.1.1
[9] scuttle_1.16.0 scDblFinder_1.20.2
[11] celda_1.22.1 Matrix_1.7-5
[13] muscData_1.20.0 ExperimentHub_2.14.0
[15] AnnotationHub_3.14.0 BiocFileCache_2.14.0
[17] dbplyr_2.5.1 SingleCellExperiment_1.28.1
[19] SummarizedExperiment_1.36.0 Biobase_2.66.0
[21] GenomicRanges_1.58.0 GenomeInfoDb_1.42.3
[23] IRanges_2.40.1 S4Vectors_0.44.0
[25] BiocGenerics_0.52.0 MatrixGenerics_1.18.1
[27] matrixStats_1.5.0 Seurat_5.4.0
[29] SeuratObject_5.3.0 sp_2.2-1
loaded via a namespace (and not attached):
[1] spatstat.sparse_3.1-0 bitops_1.0-9 httr_1.4.8
[4] RColorBrewer_1.1-3 doParallel_1.0.17 tools_4.4.2
[7] sctransform_0.4.3 utf8_1.2.6 R6_2.6.1
[10] lazyeval_0.2.2 uwot_0.2.4 withr_3.0.2
[13] gridExtra_2.3 progressr_0.18.0 cli_3.6.5
[16] spatstat.explore_3.8-0 fastDummies_1.7.5 enrichR_3.4
[19] labeling_0.4.3 S7_0.2.1 spatstat.data_3.1-9
[22] ggridges_0.5.7 pbapply_1.7-4 Rsamtools_2.22.0
[25] dbscan_1.2.7 dichromat_2.0-0.1 scater_1.34.1
[28] parallelly_1.46.1 WriteXLS_6.8.0 limma_3.62.2
[31] RSQLite_2.4.6 generics_0.1.4 BiocIO_1.16.0
[34] combinat_0.0-9 ica_1.0-3 spatstat.random_3.4-5
[37] dendextend_1.19.1 ggbeeswarm_0.7.3 abind_1.4-8
[40] lifecycle_1.0.5 yaml_2.3.12 edgeR_4.4.2
[43] SparseArray_1.6.2 Rtsne_0.17 grid_4.4.2
[46] blob_1.3.0 promises_1.5.0 dqrng_0.4.1
[49] crayon_1.5.3 miniUI_0.1.2 lattice_0.22-9
[52] beachmat_2.22.0 cowplot_1.2.0 KEGGREST_1.46.0
[55] metapod_1.14.0 pillar_1.11.1 knitr_1.51
[58] rjson_0.2.23 xgboost_3.2.1.1 future.apply_1.20.2
[61] codetools_0.2-20 fastmatch_1.1-8 glue_1.8.0
[64] spatstat.univar_3.1-7 data.table_1.18.2.1 vctrs_0.7.2
[67] png_0.1-9 spam_2.11-3 gtable_0.3.6
[70] assertthat_0.2.1 cachem_1.1.0 xfun_0.57
[73] S4Arrays_1.6.0 mime_0.13 RcppEigen_0.3.4.0.2
[76] survival_3.8-6 iterators_1.0.14 statmod_1.5.1
[79] bluster_1.16.0 fitdistrplus_1.2-6 ROCR_1.0-12
[82] nlme_3.1-168 bit64_4.6.0-1 filelock_1.0.3
[85] RcppAnnoy_0.0.23 irlba_2.3.7 vipor_0.4.7
[88] KernSmooth_2.23-26 otel_0.2.0 colorspace_2.1-2
[91] DBI_1.3.0 ggrastr_1.0.2 tidyselect_1.2.1
[94] bit_4.6.0 compiler_4.4.2 curl_7.0.0
[97] BiocNeighbors_2.0.1 ggdendro_0.2.0 DelayedArray_0.32.0
[100] plotly_4.12.0 rtracklayer_1.66.0 scales_1.4.0
[103] lmtest_0.9-40 rappdirs_0.3.4 stringr_1.6.0
[106] digest_0.6.39 goftest_1.2-3 presto_1.0.0
[109] spatstat.utils_3.2-2 rmarkdown_2.30 RhpcBLASctl_0.23-42
[112] XVector_0.46.0 htmltools_0.5.9 pkgconfig_2.0.3
[115] fastmap_1.2.0 rlang_1.1.7 htmlwidgets_1.6.4
[118] UCSC.utils_1.2.0 shiny_1.13.0 farver_2.1.2
[121] zoo_1.8-15 jsonlite_2.0.0 BiocParallel_1.40.2
[124] BiocSingular_1.22.0 RCurl_1.98-1.18 magrittr_2.0.4
[127] GenomeInfoDbData_1.2.13 dotCall64_1.2 patchwork_1.3.2
[130] viridis_0.6.5 reticulate_1.45.0 stringi_1.8.7
[133] pROC_1.19.1 MCMCprecision_0.4.2 zlibbioc_1.52.0
[136] MASS_7.3-65 plyr_1.8.9 parallel_4.4.2
[139] listenv_0.10.1 ggrepel_0.9.8 deldir_2.0-4
[142] Biostrings_2.74.1 splines_4.4.2 tensor_1.5.1
[145] locfit_1.5-9.12 igraph_2.2.2 spatstat.geom_3.7-2
[148] RcppHNSW_0.6.0 reshape2_1.4.5 ScaledMatrix_1.14.0
[151] BiocVersion_3.20.0 XML_3.99-0.23 evaluate_1.0.5
[154] scran_1.34.0 BiocManager_1.30.27 foreach_1.5.2
[157] httpuv_1.6.17 RANN_2.6.2 tidyr_1.3.2
[160] purrr_1.2.1 polyclip_1.10-7 scattermore_1.2
[163] rsvd_1.0.5 xtable_1.8-8 restfulr_0.0.16
[166] RSpectra_0.16-2 later_1.4.8 viridisLite_0.4.3
[169] tibble_3.3.1 memoise_2.0.1 beeswarm_0.4.0
[172] AnnotationDbi_1.68.0 GenomicAlignments_1.42.0 cluster_2.1.8.2
[175] globals_0.19.1 The advanced tutorial starts from raw sequencing files (FASTQ) on a computing cluster, runs Cell Ranger, and analyses single-nucleus data from human kidney. It uses SoupX with the raw droplet matrix and SCTransform normalization.