重生之我在剑桥大学学习单细胞RNA-seq分析——7. 使用Seurat进行单细胞RNA测序分析（3）

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）

重生之我在剑桥大学学习单细胞RNA-seq分析——7. 使用Seurat进行单细胞RNA测序分析（3）

推荐阅读更多精彩内容