此文章是通过学习瑞典国家生物信息学基础设施(NBIS) 所开放的单细胞分析教程加上网上所查找的资料,自身的理解所形成的,可能会有不足之处。该部分是对质控后的数据进行特征基因的选择和降维的过程。
参考来源:https://nbisweden.github.io/workshop-scRNAseq/labs/compiled/seurat/seurat_01_qc.html
感兴趣的话可以阅读原文。
降维预处理
首先,把上一步得到的质控结果读进R中,然后进行后续分析。
读取数据
suppressPackageStartupMessages({
library(Seurat)
library(cowplot)
library(ggplot2)
library(scran)
})
alldata <- readRDS("data/results/seurat_covid_qc.rds")
dim(alldata)
# [1] 18133 6374
特征基因选取
由上可以看出,经过上一步处理,还是得到了18133个基因,对所有的基因进行降维处理不仅耗时耗力,而且也没有必要,所以需要挑选出比较重要的基因进行下一步的分析。需要找到的是跨细胞高度可变的基因。
Seurat是按照基因的表达量对基因进行分类,分成不同的bins,然后选择bins中细胞间差异最大基因作为最后的重要基因,首先必须要找到那些对于细胞类群分型很重要的基因/特征,找到那些跨细胞高度可变的基因,也会为将来细胞蔟的分离做准备。
# 特征基因选择
# 首先查看基因的分布状况
alldata <- FindVariableFeatures(alldata, selection.method = "vst",
nfeatures = 2000, verbose = FALSE, assay = "RNA")
top20 <- head(VariableFeatures(alldata), 20)
LabelPoints(plot = VariableFeaturePlot(alldata), points = top20, repel = TRUE)
Z-score标准化
从上图可以看的出来,的确有部分基因,他们在细胞之间是高度可辨的,接下来就是要找到这些基因然后按照这些重要基因对细胞进行分型。
找重要基因的方法(降维)的方法有很多,常用的有PCA,t-SNE,UMAP等等,但是在进行降维处理之前,需要对基因的表达量进行Z-score标准化。因为每个基因都有不同的表达水平,而在降维处理时,如果不将其转换为相同的权重,后续将很难进行下去。
用自带的函数ScaleData
进行缩放。
alldata <- ScaleData(alldata, vars.to.regress = c("percent_mito", "nFeature_RNA"),
assay = "RNA")
PCA
可以从很多层面上对该数据进行PCA降维处理,但是在这里,只希望通过基因的表达模式来对其进行降维处理,从而找到具有生物学意义的分群。
# 计算PCA,并绘图
alldata <- RunPCA(alldata, npcs = 50, verbose = F)
plot_grid(ncol = 3, DimPlot(alldata, reduction = "pca", group.by = "orig.ident",
dims = 1:2), DimPlot(alldata, reduction = "pca", group.by = "orig.ident", dims = 3:4),
DimPlot(alldata, reduction = "pca", group.by = "orig.ident", dims = 5:6))
# 紧接着查看在每个主成分中高度可变的基因
VizDimLoadings(alldata, dims = 1:5, reduction = "pca", ncol = 5, balanced = T)
# 还可以查看每个PC对结果的解释作用的大小
ElbowPlot(alldata, reduction = "pca", ndims = 50)
可以看到前8个PC就已经包含了大量的信息,但是仍然推荐使用较多的PC,因为后面那些包含少量信息的PC可能包含稀有细胞类型的信息。
t-SNE
# 进行tSNE聚类分析并绘图
alldata <- RunTSNE(alldata, reduction = "pca", dims = 1:30, perplexity = 30, max_iter = 1000,theta = 0.5, eta = 200, num_threads = 0)
plot_grid(ncol = 3, DimPlot(alldata, reduction = "tsne", group.by = "orig.ident"))
根据t-SNE的原理,从图上可以看出,这个降维并没有把具有相同特征基因的细胞分在一起,而是不同样本之间的细胞聚在了一起,所以具有较为明显的批次效应。
UMAP
除了t-SNE与PCA之外,降维的方式还有UMAP。
alldata <- RunUMAP(alldata, reduction = "pca", dims = 1:30, n.components = 2, n.neighbors = 30,
n.epochs = 200, min.dist = 0.3, learning.rate = 1, spread = 1)
经过以上几步,基本上就完成了对可变基因的选择,接下来就是测试不同的可变基因选择方法对最后降维结果的影响。
PS:UMAP与t-SNE不一样,不会收到原始数据维数(在这里也就是基因数量)的限制,能够改变
n.components
参数从而影响最后的结果。
alldata <- RunUMAP(alldata, reduction.name = "UMAP10_on_PCA", reduction = "pca",
dims = 1:30, n.components = 10, n.neighbors = 30, n.epochs = 200, min.dist = 0.3,
learning.rate = 1, spread = 1)
降维方法比较
- 不同维度UMAP得到的降维结果比较
plot_grid(ncol = 3, DimPlot(alldata, reduction = "umap", group.by = "orig.ident") +
ggplot2::ggtitle(label = "UMAP_on_PCA"), DimPlot(alldata, reduction = "UMAP10_on_PCA",
group.by = "orig.ident", dims = 1:2) + ggplot2::ggtitle(label = "UMAP10_on_PCA"),
DimPlot(alldata, reduction = "UMAP10_on_PCA", group.by = "orig.ident", dims = 3:4) +
ggplot2::ggtitle(label = "UMAP10_on_PCA"))
也可以看出较为明显的批次效应。
- 不同降维方法结果展示
plot_grid(ncol = 3, DimPlot(alldata, reduction = "pca", group.by = "orig.ident"),
DimPlot(alldata, reduction = "tsne", group.by = "orig.ident"), DimPlot(alldata,
reduction = "umap", group.by = "orig.ident"))
缩放数据或者使用图方法降维
缩放数据进行UMAP降维
不使用全部的数据,而只使用高度可辨的基因进行UMAP降维处理。
alldata <- RunUMAP(alldata, reduction.name = "UMAP_on_ScaleData", features = alldata@assays$RNA@var.features,
assay = "RNA", n.components = 2, n.neighbors = 30, n.epochs = 200, min.dist = 0.3,
learning.rate = 1, spread = 1)
构建图方法进行UMAP降维
因为不管是t-SNE或者还是UMAP,都需要首先根据指定的距离矩阵从数据中构建graph,然后优化他们的降维方式,其实graph就是一个包含每个细胞之间距离的矩阵,所以可以先构建任意距离的graph,然后根据这个graph来进行优化。[这个点可以先去了解一些t-SNE/UMAP的基本原理,之后就好理解了]
# 构建graph(距离矩阵)
alldata <- FindNeighbors(alldata, reduction = "pca", graph.name = "SNN", assay = "RNA",
k.param = 20, features = alldata@assays$RNA@var.features)
# 在graph上运行UMAP/t-SNE
alldata <- RunUMAP(alldata, reduction.name = "UMAP_on_Graph", graph = "SNN", assay = "RNA")
比较不同降维方式的区别
p1 <- DimPlot(alldata, reduction = "umap", group.by = "orig.ident") + ggplot2::ggtitle(label = "UMAP_on_PCA")
p2 <- DimPlot(alldata, reduction = "UMAP_on_ScaleData", group.by = "orig.ident") +
ggplot2::ggtitle(label = "UMAP_on_ScaleData")
p3 <- DimPlot(alldata, reduction = "UMAP_on_Graph", group.by = "orig.ident") + ggplot2::ggtitle(label = "UMAP_on_Graph")
leg <- get_legend(p1)
gridExtra::grid.arrange(gridExtra::arrangeGrob(p1 + NoLegend() + NoAxes(), p2 + NoLegend() +
NoAxes(), p3 + NoLegend() + NoAxes(), leg, nrow = 2), ncol = 1, widths = c(1))
使用graph方法所降维的图像基本没用
展示基因
可以在降维的图中展示自己感兴趣的基因,首先要先设定自己感兴趣的基因集,然后在对每个基因进行绘图,其实就是通过特定的marker基因展示细胞类型所在位置。
Markers | Cell Type |
---|---|
CD3E | T cells |
CD3E CD4 | CD4+ T cells |
CD3E CD8A | CD8+ T cells |
GNLY, NKG7 | NK cells |
MS4A1 | B cells |
CD14, LYZ, CST3, MS4A7 | CD14+ Monocytes |
FCGR3A, LYZ, CST3, MS4A7 | FCGR3A+ Monocytes |
FCER1A, CST3 | DCs |
myfeatures <- c("CD3E", "CD4", "CD8A", "NKG7", "GNLY", "MS4A1", "CD14", "LYZ", "MS4A7",
"FCGR3A", "CST3", "FCER1A")
plot_list <- list()
for (i in myfeatures) {
plot_list[[i]] <- FeaturePlot(alldata, reduction = "umap", dims = 1:2, features = i,
ncol = 3, order = T) + NoLegend() + NoAxes() + NoGrid()
}
plot_grid(ncol = 3, plotlist = plot_list)
从上面这张图也可以很明显的看出批次效应影响了降维的结果,同一个marker基因并没有在图上有聚集的效果,证明基因对于降维结果的影响小于样本之间的影响。
保存数据
saveRDS(alldata, "data/results/covid_qc_dr.rds")