重生之我在劍橋大學(xué)學(xué)習(xí)單細(xì)胞RNA-seq分析——7. 使用Seurat進(jìn)行單細(xì)胞RNA測序分析（3）

7.3 SCTransform標(biāo)準(zhǔn)化和聚類
由于我們已經(jīng)通過額外的QC去除了雙胞和空細(xì)胞，現(xiàn)在可以應(yīng)用SCTransform標(biāo)準(zhǔn)化，這有利于通過提高信噪比來尋找稀有細(xì)胞群。SCTransform命令取代了NormalizeData、ScaleData和FindVariableFeatures。我們還將使用vars.to.regress變量來校正線粒體基因和細(xì)胞周期得分；我們之前的探索表明，細(xì)胞周期得分和線粒體百分比在cluster之間都沒有發(fā)生很大變化，因此我們不會(huì)去除生物信號(hào)，而只會(huì)去除一些不必要的變化。

> 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

之后，進(jìn)行標(biāo)準(zhǔn)PCA、UMAP和聚類。注意SCT現(xiàn)在是處在活動(dòng)狀態(tài)的assay。通常在使用SCTransform時(shí)使用更多PC；確切數(shù)量可以根據(jù)數(shù)據(jù)集進(jìn)行調(diào)整。

> 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非常重要。我們來檢查一下之前提到的較小細(xì)胞群的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之間有一簇血小板，但尚未被識(shí)別（注：這里感覺圖文不對版，關(guān)注思路就好）。將FindClusters中的聚類分辨率增加到2將有助于分離血小板簇，但也會(huì)生成過多的簇。
如果有需要，我們可以手動(dòng)分離一些簇。還有一些聚類方法專門用于識(shí)別稀有細(xì)胞群。讓我們嘗試在KNN圖中使用更少的鄰居，結(jié)合leiden算法（現(xiàn)在是scanpy中的默認(rèn)算法）并稍微增加分辨率：

> 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對應(yīng)血小板，cluster 15對應(yīng)DC（看不出來）。然后我們粗略地看一下大細(xì)胞群是什么。

> 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")

可以看到一些亞群的分離效果更好。這些將在下文進(jìn)一步討論。
7.4 差異表達(dá)和marker選擇
差異表達(dá)分析使我們能夠定義每個(gè)cluster特有的marker基因。它受cluster定義方式的影響，因此在篩選marker之前找到正確的聚類分辨率非常重要。如果某些cluster缺少顯著的marker，則調(diào)整聚類。建議在RNA assay上進(jìn)行差異表達(dá)分析，而不是SCTransform。
首先，讓我們將處于活動(dòng)狀態(tài)的assay重新設(shè)置為“RNA”，并重新進(jìn)行標(biāo)準(zhǔn)化和中心化（因?yàn)槲覀儎h除了很大一部分未通過QC的細(xì)胞）：

> 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

以下函數(shù)允許通過將每個(gè)cluster與所有剩余細(xì)胞進(jìn)行比較來找到每個(gè)cluster的marker，同時(shí)僅輸出陽性結(jié)果。許多檢驗(yàn)可用于定義marker，默認(rèn)使用Wilcoxon秩和檢驗(yàn)。這需要大概幾分鐘的時(shí)間。為了提高速度，可以增加默認(rèn)的最小百分比和log2FC閾值；調(diào)整這些參數(shù)以適應(yīng)你的數(shù)據(jù)。

> 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。為了進(jìn)行詳細(xì)分析，最好對亞群之間的差異表達(dá)進(jìn)行分析。

往期內(nèi)容：
重生之我在劍橋大學(xué)學(xué)習(xí)單細(xì)胞RNA-seq分析——7. 使用Seurat進(jìn)行單細(xì)胞RNA測序分析（1）
重生之我在劍橋大學(xué)學(xué)習(xí)單細(xì)胞RNA-seq分析——7. 使用Seurat進(jìn)行單細(xì)胞RNA測序分析（2）

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡書系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容