四 确定维度
如何确定使用多少PCA呢？以下2种方法科研作为参考
4.1 JackStrawPlot
使用JackStrawPlot()函数对PCA结果进行可视化，"重要"的PC 将会显示在虚线上方且P值较低，本示例中PC10-PC12后显著性下降较明显。较耗时。
# NOTE: This process can take a long time for big datasets, comment out for expediency. More # approximate techniques such as those implemented in ElbowPlot() can be used to reduce # computation time pbmc <- JackStraw(pbmc, num.replicate = 100) pbmc <- ScoreJackStraw(pbmc, dims = 1:20) JackStrawPlot(pbmc, dims = 1:15)
4.2 Elbow图
展示每个主成分对数据方差的解释情况（百分比表示），并进行排序。发现第9个主成分是一个拐点，后续的主成分(PC)变化不大。
ElbowPlot(pbmc)
A：以上结果“供参考”；
B：Seurat官网鼓励用户使用不同数量的PC（10、15，甚至50）重复下游分析，虽然结果通常没有显著差异。
C：建议设置此参数时偏高一些，较少维度进行下游分析可能会对结果产生一些负面影响。
4.3 显著相关基因
这个也可以作为后面分析选择基因的一个参考。
#Returns a set of genes, based on the JackStraw analysis, that have statistically significant associations with a set of PCs. ?PCASigGenes head(PCASigGenes(pbmc,pcs.use=2,pval.cut = 0.7)) [1] "PPBP" "LYZ" "S100A9" "IGLL5" "GNLY" "FTL"
五 聚类
基于PCA结果进行聚类；
大约3K细胞的单细胞数据集，将resolution参数设置在0.4-1.2之间，数据集增加resolution 值对应增加。这个参数的设置目前没有找到标准，有清楚的小伙伴欢迎后台告知，谢谢。
pbmc <- FindNeighbors(pbmc, dims = 1:10) pbmc <- FindClusters(pbmc, resolution = 0.5) # Look at cluster IDs of the first 5 cells head(Idents(pbmc), 5) head(pbmc@active.ident,5) #同上 AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1 AAACCGTGCTTCCG-1 AAACCGTGTATGCG-1 0 3 2 1 6 Levels: 0 1 2 3 4 5 6 7 8
5.1 查看指定cluster的cell
head(subset(as.data.frame(pbmc@active.ident),pbmc@active.ident=="2")) pbmc@active.identAAACATTGATCAGC-1 2 AAACGCACTGGTAC-1 2 AAAGCCTGTATGCG-1 2 AAATCAACTCGCAA-1 2 AAATGTTGCCACAA-1 2 AACACGTGCAGAGG-1 2
5.2 提取某一cluster细胞
subpbmc<-subset(x = pbmc,idents="2") subpbmc head(subpbmc@active.ident,5) AAACATTGATCAGC-1 AAACGCACTGGTAC-1 AAAGCCTGTATGCG-1 AAATCAACTCGCAA-1 AAATGTTGCCACAA-1 2 2 2 2 2 Levels: 2
六 非线性降维聚类
建议使用与聚类分析相同的pc维度进行非线性降维度分析，如tSNE和UMAP，并可视化展示。
6.1 UMAP
# If you haven't installed UMAP, you can do so via reticulate::py_install(packages = # 'umap-learn') pbmc <- RunUMAP(pbmc, dims = 1:10) # note that you can set `label = TRUE` or use the LabelClusters function to help label # individual clusters DimPlot(pbmc, reduction = "umap")
6.2 tSNE
pbmc <- RunTSNE(pbmc, dims = 1:10) head(pbmc@reductions$tsne@cell.embeddings) DimPlot(pbmc, reduction = "tsne")
6.3 比较
# note that you can set `label = TRUE` or use the LabelClusters function to help label # individual clusters plot1<-DimPlot(pbmc, reduction = "umap",label = TRUE) plot2<-DimPlot(pbmc, reduction = "tsne",label = TRUE) plot1 + plot2
可以看出，两者的降维可视化的结构是一致的，UMAP方法相对更加紧凑。
七 差异表达基因
Seurat可以通过FindMarkers函数 和 FindAllMarkers函数寻找不同cluster的差异表达基因。
min.pct参数:设定在两个细胞群中任何一个被检测到的百分比，通过此设定不检测很少表达基因来缩短程序运行时间。默认0.1
thresh.test参数：设定在两个细胞群中基因差异表达量。可以设置为0 ，程序运行时间会更长。
max.cells.per.ident参数：每个类群细胞抽样设置；也可以缩短程序运行时间。
7.1 特定cluster
是one-others的差异分析方法，由ident.1来制定cluster，本例就是cluster2与其余的cluster来做比较。
# find all markers of cluster 2 cluster2.markers <- FindMarkers(pbmc, ident.1 = 2, min.pct = 0.25) head(cluster2.markers, n = 5) ## p_val avg_log2FC pct.1 pct.2 p_val_adj ## IL32 2.593535e-91 1.2154360 0.949 0.466 3.556774e-87 ## LTB 7.994465e-87 1.2828597 0.981 0.644 1.096361e-82 ## CD3D 3.922451e-70 0.9359210 0.922 0.433 5.379250e-66 ## IL7R 1.130870e-66 1.1776027 0.748 0.327 1.550876e-62 ## LDHB 4.082189e-65 0.8837324 0.953 0.614 5.598314e-61
7.2 指定cluster
本例为cluster5 和 cluster0_cluster3的差异
# find all markers distinguishing cluster 5 from clusters 0 and 3 cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25) head(cluster5.markers, n = 5) p_val avg_log2FC pct.1 pct.2 p_val_adjFCGR3A 8.331882e-208 4.261784 0.975 0.040 1.142634e-203 CFD 1.932644e-198 3.423863 0.938 0.036 2.650429e-194 IFITM3 2.710023e-198 3.876058 0.975 0.049 3.716525e-194 CD68 1.069778e-193 3.013656 0.926 0.035 1.467094e-189 RP11-290F20.3 4.218926e-190 2.722303 0.840 0.016 5.785835e-186
7.3 FindAllMarkers
每个cluster分别与其他所有cluster进行比较，展示各个cluster的前2个基因
# find markers for every cluster compared to all remaining cells, report only the positive ones pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25) pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC) # A tibble: 18 x 7# Groups: cluster [9] p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene <dbl> <dbl> <dbl> <dbl> <dbl> <fct> <chr> 1 1.88e-117 1.08 0.913 0.588 2.57e-113 0 LDHB 2 5.01e- 85 1.34 0.437 0.108 6.88e- 81 0 CCR7 3 0 5.57 0.996 0.215 0 1 S100A9 4 0 5.48 0.975 0.121 0 1 S100A8 5 1.93e- 80 1.27 0.981 0.65 2.65e- 76 2 LTB 6 2.91e- 58 1.27 0.667 0.251 3.98e- 54 2 CD2 7 0 4.31 0.939 0.042 0 3 CD79A 8 1.06e-269 3.59 0.623 0.022 1.45e-265 3 TCL1A 9 3.60e-221 3.21 0.984 0.226 4.93e-217 4 CCL5 10 4.27e-176 3.01 0.573 0.05 5.85e-172 4 GZMK 11 3.51e-184 3.31 0.975 0.134 4.82e-180 5 FCGR3A 12 2.03e-125 3.09 1 0.315 2.78e-121 5 LST1 13 3.17e-267 4.83 0.961 0.068 4.35e-263 6 GZMB 14 1.04e-189 5.28 0.961 0.132 1.43e-185 6 GNLY 15 1.48e-220 3.87 0.812 0.011 2.03e-216 7 FCER1A 16 1.67e- 21 2.87 1 0.513 2.28e- 17 7 HLA-DPB1 17 7.73e-200 7.24 1 0.01 1.06e-195 8 PF4 18 3.68e-110 8.58 1 0.024 5.05e-106 8 PPBP
?FindAllMarkers查看更多参数。
cluster0.markers <- FindMarkers(pbmc, ident.1 = 0, logfc.threshold = 0.25,
test.use = "roc", # 检验的方式
only.pos = TRUE) #只输出pos的基因
其中test.use 可选参数有：wilcox（默认），bimod，roc，t，negbinom，poisson，LR，MAST，DESeq2。
7.4 可视化
可以通过seurat的内置函数查看重点基因的基本情况：
1）VlnPlot: 基因表达概率分布
VlnPlot(pbmc, features = c("MS4A1", "CD79A")) # you can plot raw counts as well VlnPlot(pbmc, features = c("NKG7", "PF4"), slot = "counts", log = TRUE)
2）FeaturePlot：在tSNE 或 PCA图中展示重点基因的表达
FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))
3）DoHeatmap()查看重点基因细胞和cluster的表达热图
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC) DoHeatmap(pbmc, features = top10$gene) + NoLegend()
此外，还建议使用RidgePlot()、CellScatter()和DotPlot()查看数据集的情况。
这也可以作为后续手动注释的一些参考。
