10 RNA Processing
10.1 Normalisation and variable feature selection
Normalise the RNA assay using log-normalisation and identify the top 3000 most variable genes. Cell cycle scores (S phase and G2/M phase) are assigned to each cell using the canonical marker gene lists from Seurat. These scores are used later to assess whether cell cycle phase is a confounding source of variation in the data.
DefaultAssay(combined) <- "RNA"
combined<-NormalizeData(combined) %>%
FindVariableFeatures(nfeatures = 3000) %>%
CellCycleScoring(s.features = cc.genes.updated.2019$s.genes,
g2m.features = cc.genes.updated.2019$g2m.genes)
combined<- FindVariableFeatures(combined,nfeatures = 3000)
combined <- SCTransform(combined)
combined <- ScaleData(combined)
combined <- RunPCA(combined, npcs = 50)
combined <- RunUMAP(combined, dims = 1:20,
reduction.name = "umap.rna",
reduction.key = "UMAPRNA_")
DefaultAssay(combined)
#> [1] "SCT"10.2 PCA diagnostics
Visualise the top gene loadings for the first two principal components, the PCA plot, heatmaps for the top six PCs, and an elbow plot to guide the choice of how many PCs to retain for downstream neighbour finding.
VizDimLoadings(combined, dims = 1:2, reduction = "pca")
DimPlot(combined, reduction = "pca")
DimHeatmap(combined, dims = c(1:6), cells = 500, balanced = TRUE)
ElbowPlot(combined, ndims = 50, reduction = "pca")
10.3 Initial clustering
Find neighbours using the first 15 PCs and sweep a range of resolutions from 0.1 to 2.0. The clustree section below uses these resolution sweeps to guide the final choice of resolution.
options(future.globals.maxSize = Inf)
combined <- FindNeighbors(combined, dims = 1:15)
resolution <- 2
combined <- FindClusters(combined,
reduction.type = "umap",
resolution = seq(0.1, resolution, 0.1),
dims.use = 1:15,
save.SNN = TRUE)10.4 Naive UMAP — by sample and donor
DimPlot(combined, group.by = "sample", reduction = "umap.rna")
DimPlot(combined, group.by = "sample_donor", reduction = "umap.rna")
10.5 Clustree — RNA resolution sweep
The clustree plot traces how cluster membership changes across resolutions. Stable clusters that persist across several resolutions and avoid extensive branching are preferred. Use this to select a resolution for final annotation.
clustree(combined, prefix = 'SCT_snn_res.', show_axis = TRUE) +
theme(legend.key.size = unit(0.20, 'cm'))
10.6 Confounding Factors
Before finalising the RNA embedding, it is important to evaluate whether technical or biological covariates are driving the clustering structure rather than true cell-type differences. The three main sources of confounding assessed here are cell cycle phase, ribosomal gene expression, and immunoglobulin (IG) gene expression. Mitochondrial gene expression is also scored as a reference but is not regressed out.
10.6.1 Cell cycle phase
Visualise cell cycle phase assignment on the initial RNA UMAP. If cells separate primarily by phase rather than by cell type, phase regression should be considered.
DimPlot(combined, group.by = "Phase", reduction = "umap.rna")
10.6.2 Gene set identification
Identify mitochondrial, ribosomal, and immunoglobulin gene sets from the feature names of the RNA assay. These sets are used to compute per-cell module scores, which quantify the aggregate expression of each gene set in each cell.
- Mitochondrial genes: prefix
MT- - Ribosomal genes: prefix
RP - IG genes: prefixes
IGH,IGK,IGL
DefaultAssay(combined) <- 'RNA'
mito_genes <- rownames(combined)[grepl('^MT', rownames(combined))]
ribo_genes <- rownames(combined)[grepl('^RP', rownames(combined))]
IG_genes <- c(rownames(combined)[grepl('^IGH', rownames(combined))],
rownames(combined)[grepl('^IGK', rownames(combined))],
rownames(combined)[grepl('^IGL', rownames(combined))])10.6.3 Mitochondrial genes
mito_genes
#> [1] "MTOR" "MTOR-AS1" "MTHFR" "MTFR1L"
#> [5] "MTF1" "MTF2" "MTMR11" "MTX1"
#> [9] "MTR" "MT1HL1" "MTRNR2L11" "MTA3"
#> [13] "MTIF2" "MTHFD2" "MTLN" "MTX2"
#> [17] "MTERF4" "MTMR14" "MTRNR2L12" "MTHFD2L"
#> [21] "MTTP" "MTRNR2L13" "MTNR1A" "MTRR"
#> [25] "MTMR12" "MTREX" "MTX3" "MTCH1"
#> [29] "MTO1" "MTRES1" "MTFR2" "MTHFD1L"
#> [33] "MTRF1L" "MTURN" "MTERF1" "MTPN"
#> [37] "MTRNR2L6" "MTMR9" "MTMR7" "MTUS1"
#> [41] "MTFR1" "MTERF3" "MTDH" "MTBP"
#> [45] "MTSS1" "MTAP" "MTPAP" "MTRNR2L7"
#> [49] "MTRNR2L5" "MTG1" "MTRNR2L8" "MTCH2"
#> [53] "MTA2" "MTNR1B" "MTMR2" "MTERF2"
#> [57] "MTMR6" "MTIF3" "MTUS2" "MTUS2-AS2"
#> [61] "MTUS2-AS1" "MTRF1" "MTHFD1" "MTA1"
#> [65] "MTMR10" "MTFMT" "MTHFS" "MTRNR2L4"
#> [69] "MT4" "MT3" "MT2A" "MT1E"
#> [73] "MT1M" "MT1A" "MT1B" "MT1F"
#> [77] "MT1G" "MT1H" "MT1X" "MTSS2"
#> [81] "MTHFSD" "MTRNR2L1" "MTMR4" "MTCL1"
#> [85] "MTRNR2L3" "MTG2" "MTMR3" "MTFP1"
#> [89] "MTRNR2L10" "MTMR8" "MTM1" "MTMR1"
#> [93] "MTCP1" "MT-ND1" "MT-ND2" "MT-CO1"
#> [97] "MT-CO2" "MT-ATP8" "MT-ATP6" "MT-CO3"
#> [101] "MT-ND3" "MT-ND4L" "MT-ND4" "MT-ND5"
#> [105] "MT-ND6" "MT-CYB"10.6.4 Ribosomal genes
ribo_genes
#> [1] "RPL22" "RPL11" "RPS6KA1"
#> [4] "RPA2" "RPS8" "RPE65"
#> [7] "RPF1" "RPAP2" "RPL5"
#> [10] "RPRD2" "RPTN" "RPS27"
#> [13] "RPS6KC1" "RPS7" "RPS27A"
#> [16] "RPIA" "RPL31" "RPRM"
#> [19] "RPE" "RPL37A" "RPUSD3"
#> [22] "RPL32" "RPL15" "RPSA"
#> [25] "RPL14" "RPL29" "RPP14"
#> [28] "RPL24" "RPN1" "RPL22L1"
#> [31] "RPL39L" "RPL35A" "RPL9"
#> [34] "RPL34-AS1" "RPL34" "RPS3A"
#> [37] "RPL37" "RPS23" "RPS14"
#> [40] "RPL26L1" "RPP40" "RPP21"
#> [43] "RPS18" "RPS10-NUDT3" "RPS10"
#> [46] "RPL10A" "RPL7L1" "RPF2"
#> [49] "RPS12" "RPS6KA2" "RPS6KA2-IT1"
#> [52] "RPS6KA2-AS1" "RPA3" "RP9"
#> [55] "RP1L1" "RP1" "RPS20"
#> [58] "RPL7" "RPL30" "RPL8"
#> [61] "RPS6" "RPP25L" "RPL35"
#> [64] "RPL12" "RPL7A" "RPP38-DT"
#> [67] "RPP38" "RPS24" "RPP30"
#> [70] "RPARP-AS1" "RPEL1" "RPLP2"
#> [73] "RPL27A" "RPS13" "RPS6KA4"
#> [76] "RPS6KB2" "RPS6KB2-AS1" "RPS3"
#> [79] "RPS25" "RPUSD4" "RPAP3"
#> [82] "RPS26" "RPL41" "RPL6"
#> [85] "RPH3A" "RPLP0" "RPL21"
#> [88] "RPGRIP1" "RPL10L" "RPS29"
#> [91] "RPL36AL" "RPS6KL1" "RPS6KA5"
#> [94] "RPUSD2" "RPAP1" "RPS27L"
#> [97] "RPL4" "RPLP1" "RPP25"
#> [100] "RPS17" "RPUSD1" "RPL3L"
#> [103] "RPS2" "RPS15A" "RPGRIP1L"
#> [106] "RPL13" "RPH3AL" "RPA1"
#> [109] "RPAIN" "RPL26" "RPL23A"
#> [112] "RPL23" "RPL19" "RPL27"
#> [115] "RPRML" "RPS6KB1" "RPL38"
#> [118] "RPTOR" "RPRD1A" "RPL17"
#> [121] "RPS15" "RPL36" "RPS28"
#> [124] "RPL18A" "RPS16" "RPS19"
#> [127] "RPL18" "RPL13A" "RPS11"
#> [130] "RPS9" "RPL28" "RPS5"
#> [133] "RPN2" "RPRD1B" "RPS21"
#> [136] "RPL3" "RPS19BP1" "RPS6KA3"
#> [139] "RPGR" "RP2" "RPS4X"
#> [142] "RPS6KA6" "RPA4" "RPL36A"
#> [145] "RPL39" "RPL10" "RPS4Y1"
#> [148] "RPS4Y2"10.6.5 IG genes
IG_genes
#> [1] "IGHEP2" "IGHMBP2" "IGHA2"
#> [4] "IGHE" "IGHG4" "IGHG2"
#> [7] "IGHGP" "IGHA1" "IGHEP1"
#> [10] "IGHG1" "IGHG3" "IGHD"
#> [13] "IGHM" "IGHJ6" "IGHJ3P"
#> [16] "IGHJ5" "IGHJ4" "IGHJ3"
#> [19] "IGHJ2P" "IGHJ2" "IGHJ1"
#> [22] "IGHD7-27" "IGHJ1P" "IGHD1-26"
#> [25] "IGHD6-25" "IGHD5-24" "IGHD4-23"
#> [28] "IGHD3-22" "IGHD2-21" "IGHD1-20"
#> [31] "IGHD6-19" "IGHD5-18" "IGHD4-17"
#> [34] "IGHD3-16" "IGHD2-15" "IGHD1-14"
#> [37] "IGHD6-13" "IGHD5-12" "IGHD4-11"
#> [40] "IGHD3-10" "IGHD3-9" "IGHD2-8"
#> [43] "IGHD1-7" "IGHD6-6" "IGHD5-5"
#> [46] "IGHD4-4" "IGHD3-3" "IGHD2-2"
#> [49] "IGHD1-1" "IGHV6-1" "IGHVII-1-1"
#> [52] "IGHV1-2" "IGHVIII-2-1" "IGHV1-3"
#> [55] "IGHV4-4" "IGHV7-4-1" "IGHV2-5"
#> [58] "IGHVIII-5-1" "IGHVIII-5-2" "IGHV3-6"
#> [61] "IGHV3-7" "IGHV3-64D" "IGHV5-10-1"
#> [64] "IGHV3-11" "IGHVIII-11-1" "IGHV1-12"
#> [67] "IGHV3-13" "IGHVIII-13-1" "IGHV1-14"
#> [70] "IGHV3-15" "IGHVII-15-1" "IGHV3-16"
#> [73] "IGHVIII-16-1" "IGHV1-17" "IGHV1-18"
#> [76] "IGHV3-19" "IGHV3-20" "IGHV3-21"
#> [79] "IGHV3-22" "IGHVII-22-1" "IGHVIII-22-2"
#> [82] "IGHV3-23" "IGHV1-24" "IGHV3-25"
#> [85] "IGHVIII-25-1" "IGHV2-26" "IGHVIII-26-1"
#> [88] "IGHVII-26-2" "IGHV7-27" "IGHV4-28"
#> [91] "IGHVII-28-1" "IGHV3-32" "IGHV3-30"
#> [94] "IGHVII-30-1" "IGHV3-30-2" "IGHV4-31"
#> [97] "IGHVII-30-21" "IGHV3-29" "IGHV3-33"
#> [100] "IGHVII-33-1" "IGHV3-33-2" "IGHV4-34"
#> [103] "IGHV7-34-1" "IGHV3-35" "IGHV3-36"
#> [106] "IGHV3-37" "IGHV3-38" "IGHVIII-38-1"
#> [109] "IGHV4-39" "IGHV7-40" "IGHVII-40-1"
#> [112] "IGHV3-41" "IGHV3-42" "IGHV3-43"
#> [115] "IGHVII-43-1" "IGHVIII-44" "IGHVIV-44-1"
#> [118] "IGHVII-44-2" "IGHV1-45" "IGHV1-46"
#> [121] "IGHVII-46-1" "IGHV3-47" "IGHVIII-47-1"
#> [124] "IGHV3-48" "IGHV3-49" "IGHVII-49-1"
#> [127] "IGHV3-50" "IGHV5-51" "IGHVIII-51-1"
#> [130] "IGHVII-51-2" "IGHV3-52" "IGHV3-53"
#> [133] "IGHVII-53-1" "IGHV3-54" "IGHV4-55"
#> [136] "IGHV7-56" "IGHV3-57" "IGHV1-58"
#> [139] "IGHV4-59" "IGHV3-60" "IGHVII-60-1"
#> [142] "IGHV4-61" "IGHV3-62" "IGHVII-62-1"
#> [145] "IGHV3-63" "IGHV3-64" "IGHV3-65"
#> [148] "IGHVII-65-1" "IGHV3-66" "IGHV1-67"
#> [151] "IGHVII-67-1" "IGHVIII-67-2" "IGHVIII-67-3"
#> [154] "IGHVIII-67-4" "IGHV1-68" "IGHV1-69"
#> [157] "IGHV2-70D" "IGHV3-69-1" "IGHV1-69-2"
#> [160] "IGHV1-69D" "IGHV2-70" "IGHV3-71"
#> [163] "IGHV3-72" "IGHV3-73" "IGHV3-74"
#> [166] "IGHVII-74-1" "IGHV3-75" "IGHV3-76"
#> [169] "IGHVIII-76-1" "IGHV5-78" "IGHVII-78-1"
#> [172] "IGHV3-79" "IGHV4-80" "IGHV7-81"
#> [175] "IGHVIII-82" "IGHV1OR15-9" "IGHV1OR15-2"
#> [178] "IGHV3OR15-7" "IGHD5OR15-5A" "IGHD4OR15-4A"
#> [181] "IGHD3OR15-3A" "IGHD2OR15-2A" "IGHD1OR15-1A"
#> [184] "IGHV1OR15-6" "IGHD5OR15-5B" "IGHD4OR15-4B"
#> [187] "IGHD3OR15-3B" "IGHD2OR15-2B" "IGHD1OR15-1B"
#> [190] "IGHV1OR15-1" "IGHV1OR15-3" "IGHV4OR15-8"
#> [193] "IGHV1OR15-4" "IGHV1OR16-1" "IGHV1OR16-3"
#> [196] "IGHV3OR16-9" "IGHV2OR16-5" "IGHV3OR16-15"
#> [199] "IGHV3OR16-6" "IGHV1OR16-2" "IGHV3OR16-10"
#> [202] "IGHV1OR16-4" "IGHV3OR16-8" "IGHV3OR16-12"
#> [205] "IGHV3OR16-13" "IGHV3OR16-11" "IGHV3OR16-7"
#> [208] "IGHV1OR21-1" "IGKV1OR1-1" "IGKV3OR2-268"
#> [211] "IGKC" "IGKJ5" "IGKJ4"
#> [214] "IGKJ3" "IGKJ2" "IGKJ1"
#> [217] "IGKV4-1" "IGKV5-2" "IGKV7-3"
#> [220] "IGKV2-4" "IGKV1-5" "IGKV1-6"
#> [223] "IGKV3-7" "IGKV1-8" "IGKV1-9"
#> [226] "IGKV2-10" "IGKV3-11" "IGKV1-12"
#> [229] "IGKV1-13" "IGKV2-14" "IGKV3-15"
#> [232] "IGKV1-16" "IGKV1-17" "IGKV2-18"
#> [235] "IGKV2-19" "IGKV3-20" "IGKV6-21"
#> [238] "IGKV1-22" "IGKV2-23" "IGKV2-24"
#> [241] "IGKV3-25" "IGKV2-26" "IGKV1-27"
#> [244] "IGKV2-28" "IGKV2-29" "IGKV2-30"
#> [247] "IGKV3-31" "IGKV1-32" "IGKV1-33"
#> [250] "IGKV3-34" "IGKV1-35" "IGKV2-36"
#> [253] "IGKV1-37" "IGKV2-38" "IGKV1-39"
#> [256] "IGKV2-40" "IGKV2D-40" "IGKV1D-39"
#> [259] "IGKV2D-38" "IGKV1D-37" "IGKV2D-36"
#> [262] "IGKV1D-35" "IGKV3D-34" "IGKV1D-33"
#> [265] "IGKV1D-32" "IGKV3D-31" "IGKV2D-30"
#> [268] "IGKV2D-29" "IGKV2D-28" "IGKV1D-27"
#> [271] "IGKV2D-26" "IGKV3D-25" "IGKV2D-24"
#> [274] "IGKV2D-23" "IGKV1D-22" "IGKV6D-21"
#> [277] "IGKV3D-20" "IGKV2D-19" "IGKV2D-18"
#> [280] "IGKV6D-41" "IGKV1D-17" "IGKV1D-16"
#> [283] "IGKV3D-15" "IGKV2D-14" "IGKV1D-13"
#> [286] "IGKV1D-12" "IGKV3D-11" "IGKV2D-10"
#> [289] "IGKV1D-42" "IGKV1D-43" "IGKV1D-8"
#> [292] "IGKV3D-7" "IGKV1OR2-118" "IGKV1OR2-1"
#> [295] "IGKV2OR2-1" "IGKV2OR2-2" "IGKV1OR2-3"
#> [298] "IGKV1OR2-9" "IGKV2OR2-10" "IGKV2OR2-7D"
#> [301] "IGKV3OR2-5" "IGKV1OR2-6" "IGKV2OR2-7"
#> [304] "IGKV2OR2-8" "IGKV1OR2-11" "IGKV1OR2-108"
#> [307] "IGKV1OR9-2" "IGKV1OR-2" "IGKV1OR9-1"
#> [310] "IGKV1OR-3" "IGKV1OR10-1" "IGKV1OR22-5"
#> [313] "IGKV2OR22-4" "IGKV2OR22-3" "IGKV3OR22-2"
#> [316] "IGKV1OR22-1" "IGLV8OR8-1" "IGLJCOR18"
#> [319] "IGLON5" "IGLVI-70" "IGLV4-69"
#> [322] "IGLVI-68" "IGLV10-54" "IGLV10-67"
#> [325] "IGLVIV-66-1" "IGLVV-66" "IGLVIV-65"
#> [328] "IGLVIV-64" "IGLVI-63" "IGLV1-62"
#> [331] "IGLV8-61" "IGLV4-60" "IGLVIV-59"
#> [334] "IGLVV-58" "IGLV6-57" "IGLVI-56"
#> [337] "IGLV11-55" "IGLVIV-53" "IGLV5-52"
#> [340] "IGLV1-51" "IGLV1-50" "IGLV9-49"
#> [343] "IGLV5-48" "IGLV1-47" "IGLV7-46"
#> [346] "IGLV5-45" "IGLV1-44" "IGLV7-43"
#> [349] "IGLVI-42" "IGLVVII-41-1" "IGLV1-41"
#> [352] "IGLV1-40" "IGLVI-38" "IGLV5-37"
#> [355] "IGLV1-36" "IGLV7-35" "IGLV2-34"
#> [358] "IGLV2-33" "IGLV3-32" "IGLV3-31"
#> [361] "IGLV3-30" "IGLV3-29" "IGLV2-28"
#> [364] "IGLV3-27" "IGLV3-26" "IGLVVI-25-1"
#> [367] "IGLV3-25" "IGLV3-24" "IGLV2-23"
#> [370] "IGLVVI-22-1" "IGLV3-22" "IGLV3-21"
#> [373] "IGLVI-20" "IGLV3-19" "IGLV2-18"
#> [376] "IGLV3-17" "IGLV3-16" "IGLV3-15"
#> [379] "IGLV2-14" "IGLV3-13" "IGLV3-12"
#> [382] "IGLV2-11" "IGLV3-10" "IGLV3-9"
#> [385] "IGLV2-8" "IGLV3-7" "IGLV3-6"
#> [388] "IGLV2-5" "IGLV3-4" "IGLV4-3"
#> [391] "IGLV3-2" "IGLV3-1" "IGLJ1"
#> [394] "IGLC1" "IGLJ2" "IGLC2"
#> [397] "IGLJ3" "IGLC3" "IGLJ4"
#> [400] "IGLC4" "IGLJ5" "IGLC5"
#> [403] "IGLJ6" "IGLC6" "IGLJ7"
#> [406] "IGLC7" "IGLL1" "IGLVIVOR22-1"
#> [409] "IGLCOR22-1" "IGLCOR22-2" "IGLVIVOR22-2"10.7 Module scores
Compute per-cell module scores for each gene set. These scores are stored in
the metadata as mito_genes1, ribo_genes2, and IG_genes3 and can be
overlaid on the UMAP to evaluate spatial distribution of each signal.
genes_regress <- list(mito_genes = mito_genes,
ribo_genes = ribo_genes,
IG_genes = IG_genes)
combined <- AddModuleScore(object = combined,
features = genes_regress,
ctrl = 5,
name = c("mito_genes", "ribo_genes", "IG_genes"),
search = TRUE)10.7.1 Mitochondrial gene module score
FeaturePlot(combined, features = "mito_genes1", label = TRUE, repel = TRUE,
reduction = "umap.rna") +
scale_colour_gradientn(colours = c("lightblue", "beige", "red"))
10.7.2 Ribosomal gene module score
FeaturePlot(combined, features = "ribo_genes2", label = TRUE, repel = TRUE,
reduction = "umap.rna") +
scale_colour_gradientn(colours = c("lightblue", "beige", "red"))
10.7.3 IG gene module score
FeaturePlot(combined, features = "IG_genes3", label = TRUE, repel = TRUE,
reduction = "umap.rna") +
scale_colour_gradientn(colours = c("lightblue", "beige", "red"))
10.8 Regression of ribosomal genes
Ribosomal gene expression can dominate the variance in immune cell datasets,
particularly in plasma cells and activated B cells, and cause biologically
unrelated cell types to cluster together. The RNA data are rescaled with
ribosomal module score regressed out, a new PCA is computed, and a new UMAP
embedding (umap.rna_regressedRibo) is generated for comparison.
10.9 Post-regression evaluation
Overlay cell cycle phase and each module score on the ribo-regressed UMAP to confirm that ribosomal signal has been reduced without distorting the biologically meaningful structure.
10.9.2 Mitochondrial gene module score
FeaturePlot(combined, features = "mito_genes1", label = TRUE, repel = TRUE,
reduction = "umap.rna_regressedRibo") +
scale_colour_gradientn(colours = c("lightblue", "beige", "red"))
10.9.3 Ribosomal gene module score
FeaturePlot(combined, features = "ribo_genes2", label = TRUE, repel = TRUE,
reduction = "umap.rna_regressedRibo") +
scale_colour_gradientn(colours = c("lightblue", "beige", "red"))
10.9.4 IG gene module score
FeaturePlot(combined, features = "IG_genes3", label = TRUE, repel = TRUE,
reduction = "umap.rna_regressedRibo") +
scale_colour_gradientn(colours = c("lightblue", "beige", "red"))