7.3 SCTransform标准化和聚类
由于我们已经通过额外的QC去除了双胞和空细胞，现在可以应用SCTransform标准化，这有利于通过提高信噪比来寻找稀有细胞群。SCTransform命令取代了NormalizeData、ScaleData和FindVariableFeatures。我们还将使用vars.to.regress变量来校正线粒体基因和细胞周期得分；我们之前的探索表明，细胞周期得分和线粒体百分比在cluster之间都没有发生很大变化，因此我们不会去除生物信号，而只会去除一些不必要的变化。
> srat <- SCTransform(srat, method = "glmGamPoi", ncells = 8824,
vars.to.regress = c("percent.mt","S.Score","G2M.Score"), verbose = F)
> srat
An object of class Seurat
56857 features across 8824 samples within 2 assays
Active assay: SCT (20256 features, 3000 variable features)
1 other assay present: RNA
2 dimensional reductions calculated: pca, umap
之后，进行标准PCA、UMAP和聚类。注意SCT现在是处在活动状态的assay。通常在使用SCTransform时使用更多PC；确切数量可以根据数据集进行调整。
> srat <- RunPCA(srat, verbose = F)
> srat <- RunUMAP(srat, dims = 1:30, verbose = F)
> srat <- FindNeighbors(srat, dims = 1:30, verbose = F)
> srat <- FindClusters(srat, verbose = F)
> table(srat[[]]$seurat_clusters)
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
1766 1554 1214 1066 927 379 328 316 301 268 249 131 107 103 98 17
> DimPlot(srat, label = T)
正确定义cluster非常重要。我们来检查一下之前提到的较小细胞群的marker——即血小板和DC。
> FeaturePlot(srat,"PPBP") &
scale_colour_gradientn(colours = rev(brewer.pal(n = 11, name = "Spectral")))
> FeaturePlot(srat,"LILRA4") &
scale_colour_gradientn(colours = rev(brewer.pal(n = 11, name = "Spectral")))
我们可以看到，在cluster 6和cluster 14之间有一簇血小板，但尚未被识别（注：这里感觉图文不对版，关注思路就好）。将FindClusters中的聚类分辨率增加到2将有助于分离血小板簇，但也会生成过多的簇。
如果有需要，我们可以手动分离一些簇。还有一些聚类方法专门用于识别稀有细胞群。让我们尝试在KNN图中使用更少的邻居，结合leiden算法（现在是scanpy中的默认算法）并稍微增加分辨率：
> srat <- FindNeighbors(srat, dims = 1:30, k.param = 15, verbose = F)
> srat <- FindClusters(srat, verbose = F, algorithm = 4, resolution = 0.9)
> table(srat[[]]$seurat_clusters)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
1738 1342 1211 1056 934 380 332 318 315 265 251 211 133 105 104 99
17 18
17 13
> DimPlot(srat, label = T)
已知cluster 16对应血小板，cluster 15对应DC（看不出来）。然后我们粗略地看一下大细胞群是什么。
> FeaturePlot(srat,"MS4A1") +
scale_colour_gradientn(colours = rev(brewer.pal(n = 11, name = "Spectral"))) +
ggtitle("MS4A1: B cells")
> FeaturePlot(srat,"LYZ") +
scale_colour_gradientn(colours = rev(brewer.pal(n = 11, name = "Spectral"))) +
ggtitle("LYZ: monocytes")
> FeaturePlot(srat,"NKG7") +
scale_colour_gradientn(colours = rev(brewer.pal(n = 11, name = "Spectral"))) +
ggtitle("NKG7: natural killers")
> FeaturePlot(srat,"CD8B") +
scale_colour_gradientn(colours = rev(brewer.pal(n = 11, name = "Spectral"))) +
ggtitle("CD8B: CD8 T cells")
> FeaturePlot(srat,"IL7R") +
scale_colour_gradientn(colours = rev(brewer.pal(n = 11, name = "Spectral"))) +
ggtitle("IL7R: CD4 T cells")
可以看到一些亚群的分离效果更好。这些将在下文进一步讨论。
7.4 差异表达和marker选择
差异表达分析使我们能够定义每个cluster特有的marker基因。它受cluster定义方式的影响，因此在筛选marker之前找到正确的聚类分辨率非常重要。如果某些cluster缺少显著的marker，则调整聚类。建议在RNA assay上进行差异表达分析，而不是SCTransform。
首先，让我们将处于活动状态的assay重新设置为“RNA”，并重新进行标准化和中心化（因为我们删除了很大一部分未通过QC的细胞）：
> DefaultAssay(srat) <- "RNA"
> srat <- NormalizeData(srat)
> srat <- FindVariableFeatures(srat, selection.method = "vst", nfeatures = 2000)
> all.genes <- rownames(srat)
> srat <- ScaleData(srat, features = all.genes)
Centering and scaling data matrix
以下函数允许通过将每个cluster与所有剩余细胞进行比较来找到每个cluster的marker，同时仅输出阳性结果。许多检验可用于定义marker，默认使用Wilcoxon秩和检验。这需要大概几分钟的时间。为了提高速度，可以增加默认的最小百分比和log2FC阈值；调整这些参数以适应你的数据。
> all.markers <- FindAllMarkers(srat, only.pos = T, min.pct = 0.5, logfc.threshold = 0.5)
快速浏览一下这些marker。
> dim(all.markers)
[1] 4898 7
> table(all.markers$cluster)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
544 131 130 74 530 145 447 183 62 261 154 84 249 152 524 117 224 887
> top3_markers <- as.data.frame(all.markers %>%
group_by(cluster) %>%
top_n(n = 3, wt = avg_log2FC))
> top3_markers
p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene
1 0.000000e+00 4.5268744 0.997 0.177 0.000000e+00 1 S100A8
2 0.000000e+00 4.2834962 0.997 0.200 0.000000e+00 1 S100A9
3 0.000000e+00 3.2865113 0.972 0.115 0.000000e+00 1 S100A12
4 0.000000e+00 1.2341390 0.911 0.372 0.000000e+00 2 TCF7
5 1.164915e-302 1.1729103 0.888 0.383 4.263707e-298 2 BCL11B
6 1.180905e-198 1.3114609 0.699 0.304 4.322230e-194 2 TRBC1
7 0.000000e+00 2.5047609 0.982 0.069 0.000000e+00 3 CD8B
8 0.000000e+00 1.7737856 0.677 0.023 0.000000e+00 3 LINC02446
9 0.000000e+00 1.6422981 0.876 0.078 0.000000e+00 3 CD8A
10 0.000000e+00 1.6336200 0.986 0.475 0.000000e+00 4 IL32
11 0.000000e+00 1.6077712 0.991 0.653 0.000000e+00 4 LTB
12 3.238201e-273 1.4391753 0.934 0.406 1.185214e-268 4 IL7R
13 0.000000e+00 2.1857867 0.900 0.193 0.000000e+00 5 APOBEC3A
14 0.000000e+00 2.1288663 0.978 0.290 0.000000e+00 5 MARCKS
15 0.000000e+00 2.1100918 0.984 0.324 0.000000e+00 5 IFITM3
16 0.000000e+00 3.3394836 0.997 0.051 0.000000e+00 6 MS4A1
17 0.000000e+00 3.2002061 0.682 0.072 0.000000e+00 6 IGKC
18 0.000000e+00 3.1635207 1.000 0.088 0.000000e+00 6 CD79A
19 0.000000e+00 3.7269130 1.000 0.139 0.000000e+00 7 FCGR3A
20 0.000000e+00 3.5624660 0.913 0.027 0.000000e+00 7 CDKN1C
21 3.857245e-241 2.8285081 1.000 0.370 1.411790e-236 7 IFITM3
22 0.000000e+00 3.8483906 0.937 0.033 0.000000e+00 8 GZMH
23 0.000000e+00 3.7596608 0.997 0.117 0.000000e+00 8 CCL5
24 0.000000e+00 3.7303640 1.000 0.170 0.000000e+00 8 NKG7
25 0.000000e+00 3.2997953 0.841 0.034 0.000000e+00 9 GZMK
26 0.000000e+00 2.6895744 0.921 0.120 0.000000e+00 9 CCL5
27 7.722065e-230 1.7610753 0.800 0.126 2.826353e-225 9 GZMA
28 0.000000e+00 5.2176026 1.000 0.073 0.000000e+00 10 GNLY
29 0.000000e+00 4.2078064 1.000 0.176 0.000000e+00 10 NKG7
30 0.000000e+00 3.8337378 0.992 0.103 0.000000e+00 10 PRF1
31 0.000000e+00 4.0114869 1.000 0.069 0.000000e+00 11 IGHM
32 0.000000e+00 3.7919296 1.000 0.072 0.000000e+00 11 IGKC
33 0.000000e+00 3.5643880 0.964 0.021 0.000000e+00 11 TCL1A
34 6.325413e-39 0.8709246 0.777 0.360 2.315164e-34 12 CCR7
35 7.346774e-38 0.8722015 0.957 0.534 2.688993e-33 12 CD3E
36 2.507370e-26 0.8702442 0.853 0.537 9.177223e-22 12 TRBC2
37 0.000000e+00 2.7732163 0.842 0.008 0.000000e+00 13 FCER1A
38 2.353340e-111 2.9986254 0.992 0.298 8.613459e-107 13 HLA-DQA1
39 3.412537e-88 2.7180890 1.000 0.495 1.249023e-83 13 HLA-DPB1
40 0.000000e+00 4.4809498 0.905 0.025 0.000000e+00 14 IGLC2
41 0.000000e+00 3.9502507 0.686 0.018 0.000000e+00 14 IGLC3
42 1.608644e-238 3.8078629 1.000 0.085 5.887798e-234 14 IGHM
43 0.000000e+00 3.5738916 0.837 0.040 0.000000e+00 15 JCHAIN
44 1.000645e-238 4.1104811 0.817 0.055 3.662459e-234 15 GZMB
45 1.605595e-139 3.5527436 0.846 0.117 5.876638e-135 15 ITM2C
46 0.000000e+00 2.8907374 0.980 0.052 0.000000e+00 16 GZMK
47 3.476722e-192 3.6726402 1.000 0.114 1.272515e-187 16 KLRB1
48 2.240730e-121 2.2655451 0.848 0.124 8.201297e-117 16 NCR3
49 1.173559e-117 8.0515428 1.000 0.030 4.295342e-113 17 IGLC1
50 3.450223e-91 9.4362396 0.941 0.036 1.262816e-86 17 IGHA1
51 1.364084e-24 7.8462060 0.765 0.097 4.992685e-20 17 IGKC
52 0.000000e+00 2.9967998 1.000 0.004 0.000000e+00 18 TYMS
53 8.837517e-20 3.8742685 1.000 0.188 3.234619e-15 18 STMN1
54 1.157218e-09 3.3464621 1.000 0.740 4.235532e-05 18 TUBA1B
有些marker的信息量不如其他marker。为了进行详细分析，最好对亚群之间的差异表达进行分析。
往期内容：
重生之我在剑桥大学学习单细胞RNA-seq分析——7. 使用Seurat进行单细胞RNA测序分析（1）
重生之我在剑桥大学学习单细胞RNA-seq分析——7. 使用Seurat进行单细胞RNA测序分析（2）
