查看: 5194|回复: 5

NBIS系列单细胞转录组数据分析实战(二):数据降维可视化

[复制链接]

版主

Rank: 10Rank: 10Rank: 10

主题
34
注册时间
2016.5.11
在线时间
658 小时

活跃会员突出贡献


发表于 2021.4.4 11:08:26 | 显示全部楼层 |阅读模式
本帖最后由 bioinfomics 于 2021.4.4 11:07 编辑
本文首发于微信公众号“bioinfomics”:https://mp.weixin.qq.com/s/SIn1SlybuTM-Zrd9d9d1Hw

第二节:数据降维可视化
加载所需的R包和质控过滤后的数据
首先,我们加载分析所需的R包和上一节中质控过滤后的数据集。
[Bash shell] 纯文本查看 复制代码
# 加载所需的R包
suppressPackageStartupMessages({
    library(Seurat)
    library(cowplot)
    library(ggplot2)
    library(scran)
})

# 加载质控过滤后的数据
alldata <- readRDS("data/results/seurat_covid_qc.rds")

# 查看数据信息
alldata
#An object of class Seurat 
#18121 features across 5532 samples within 1 assay 
#Active assay: RNA (18121 features, 2000 variable features)
# 2 dimensional reductions calculated: pca, umap

特征选择筛选高变基因
接下来,我们需要筛选出数据集中细胞之间高度可变的特征/基因,用于后续的降维分析,以区分不同的细胞类型。
[AppleScript] 纯文本查看 复制代码
alldata <- FindVariableFeatures(alldata, selection.method = "vst", assay = "RNA"
                                nfeatures = 2000, verbose = FALSE)
top20 <- head(VariableFeatures(alldata), 20)
head(top20)
# [1] "IGLC3"  "CXCL10" "IFI27"  "PTGDS"  "LYPD2"  "STMN1"

# 高变基因可视化
LabelPoints(plot = VariableFeaturePlot(alldata), points = top20, repel = TRUE)

Z-score数据归一化
接下来,我们对数据进行z-score归一化处理。由于每个基因具有不同的表达水平,这意味着具有较高表达值的基因将会具有较高的变异,这些变异也更容易被PCA捕获。因此,在进行PCA降维之前,我们需要以某种方式赋予每个基因类似的权重(请参见下文)。常见的做法是在进行PCA之前对每个基因进行中心化(centering)和缩放(scaling)。这种精确的缩放比例称为Z-score归一化,对于PCA,聚类和绘制热图都是非常有用的。

此外,我们还可以使用回归的方法来从数据集中删除任何不想要的变异源,例如cell cycle,sequencing depth,percent mitocondria。我们可以通过将这些参数作为模型中的协变量进行广义线性回归来实现,然后将模型的残差作为“回归数据”。
[AppleScript] 纯文本查看 复制代码
alldata <- ScaleData(alldata, vars.to.regress = c("percent_mito", "nFeature_RNA"), assay = "RNA")

PCA线性降维
[AppleScript] 纯文本查看 复制代码
# 使用RunPCA函数进行PCA线性降维
alldata <- RunPCA(alldata, npcs = 50, verbose = F)


然后,我们绘制前几个主成分。
[AppleScript] 纯文本查看 复制代码
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))

为了确定哪些基因或元数据对每个PC的贡献最大,我们可以使用VizDimLoadings函数查看相应基因对每个PC的贡献信息。
[AppleScript] 纯文本查看 复制代码
VizDimLoadings(alldata, dims = 1:5, reduction = "pca", ncol = 5, balanced = T)

我们还可以使用ElbowPlot函数绘制每个PC解释的方差量,以确定后续所要使用的PC数。
[AppleScript] 纯文本查看 复制代码
ElbowPlot(alldata, reduction = "pca", ndims = 50)

如图所示,我们可以看到排名前10位的PC保留了大量变异信息,而其他的PC所包含的信息则较少。但是,仍然建议使用更多的PC,因为它们可能包含有关稀有细胞类型(例如此数据集中的血小板和DC)的信息。

t-SNE降维可视化
接下来,我们进行t-SNE降维可视化处理。
[AppleScript] 纯文本查看 复制代码
# 使用RunTSNE函数进行t-SNE降维可视化
alldata <- RunTSNE(alldata, reduction = "pca", dims = 1:30, 
                   perplexity = 30, max_iter = 1000, 
                   theta = 0.5, eta = 200, num_threads = 0)
# see ?Rtsne and ?RunTSNE for more info

现在,我们可以绘制t-SNE降维后的数据,可以清楚地看到数据集中存在一定的批次效应。
[AppleScript] 纯文本查看 复制代码
DimPlot(alldata, reduction = "tsne", group.by = "orig.ident")


UMAP降维可视化

[AppleScript] 纯文本查看 复制代码
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)
# see ?RunUMAP for more info

DimPlot(alldata, reduction = "umap", group.by = "orig.ident")



Task:测试设置不同的参数
UMAP降维的另一个有用之处是,它不受数据维数的限制(与tSNE不同),我们可以简单地改变n.components参数的大小。
[AppleScript] 纯文本查看 复制代码
# we can add in additional reductions, by defulat they are named 'pca', 'umap',
# 'tsne' etc. But we can specify alternative names with reduction.name
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)
# see ?RunUMAP for more info


现在,我们可以绘制UMAP降维后的结果。尽管不像tSNE那样明显,但我们仍然可以看到数据中不同批次的影响很大。
[AppleScript] 纯文本查看 复制代码
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"))

现在,我们可以并排绘制PCA,UMAP和tSNE降维的结果进行比较。在这里,我们可以得出结论,我们的数据集中存在一定的批处理效应,需要在进行后续聚类和差异基因表达分析之前对其进行校正。
[AppleScript] 纯文本查看 复制代码
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"))

使用ScaledData和graphs进行数据降维
尽管在PCA线性降维后进行第二次降维(即tSNE或UMAP)将是一种标准方法(因为它允许更高的计算效率),但实际上这些选择是可变的。下面,我们将展示其他两个常见的做法,例如直接在缩放数据(用于PCA)上运行,或在根据缩放数据构建的图聚类上运行。从现在开始,我们将仅在UMAP上可视化降维的结果,但对于t-SNE同样如此。

使用ScaledData进行UMAP可视化
要对scaling后的数据进行t-SNE或UMAP降维可视化,首先需要选择要使用的变量数目。这是因为,那些有助于细胞类型分离的维数将最终掩盖这些差异。另一个原因是因为用所有的基因/特征一起运行也将花费更长的时间,或者可能在计算上不可行。因此,我们将使用高度可变的基因的缩放数据。
[AppleScript] 纯文本查看 复制代码
# 指定features参数设置所要使用的变异基因
alldata <- RunUMAP(alldata, reduction.name = "UMAP_on_ScaleData", 
                   features = alldata@assays$[url=mailto:RNA@var.features]RNA@var.features[/url], assay = "RNA", 
                   n.components = 2, n.neighbors = 30, n.epochs = 200, 
                   min.dist = 0.3, learning.rate = 1, spread = 1)

使用图聚类进行UMAP可视化
要在图聚类上进行t-SNE或UMAP降维可视化,我们首先需要根据数据构建图形。实际上,t-SNE和UMAP都首先使用指定的距离矩阵从数据中构建图形,然后优化嵌入。由于图形只是一个矩阵,其中包含不同细胞之间的距离,因此,我们可以使用所需的任何其他距离度量方法来运行UMAP或tSNE。Euclidean和Correlation距离通常是最常用的距离度量。
[AppleScript] 纯文本查看 复制代码
# Build Graph
alldata <- FindNeighbors(alldata, reduction = "pca", graph.name = "SNN", 
                         assay = "RNA", k.param = 20, 
                         features = alldata@assays$[url=mailto:RNA@var.features]RNA@var.features[/url])
# Run UMAP on a graph
alldata <- RunUMAP(alldata, reduction.name = "UMAP_on_Graph", graph = "SNN", assay = "RNA")

现在,我们可以在PCA,ScaledData和Graph上进行比较,绘制UMAP可视化图。
[AppleScript] 纯文本查看 复制代码
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))

可视化感兴趣基因的表达
接下来,我们对一些不同细胞类型的marker基因进行可视化展示。
MarkersCell Type
CD3ET cells
CD3E CD4CD4+ T cells
CD3E CD8ACD8+ T cells
GNLY, NKG7NK cells
MS4A1B cells
CD14, LYZ, CST3, MS4A7CD14+ Monocytes
FCGR3A, LYZ, CST3, MS4A7FCGR3A+ Monocytes
FCER1A, CST3DCs
[AppleScript] 纯文本查看 复制代码
myfeatures <- c("CD3E", "CD4", "CD8A", "NKG7", "GNLY", "MS4A1", 
                "CD14", "LYZ", "MS4A7", "FCGR3A", "CST3", "FCER1A")
# 使用FeaturePlot函数进行可视化展示
FeaturePlot(alldata, reduction = "umap", dims = 1:2, 
           features = myfeatures, ncol = 3, order = T)



保存降维可视化后的数据
[AppleScript] 纯文本查看 复制代码
saveRDS(alldata, "data/results/covid_qc_dr.rds")


本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
加油啦。
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2020.11.21
在线时间
62 小时

发表于 2021.4.4 12:39:08 | 显示全部楼层
坚持就是胜利!
回复

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
3
注册时间
2017.9.8
在线时间
65 小时

发表于 2021.4.5 10:59:53 | 显示全部楼层
加油,加油!
回复

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
注册时间
2021.3.29
在线时间
7 小时

发表于 2021.4.5 18:44:59 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
3
注册时间
2017.9.8
在线时间
65 小时

发表于 2021.4.12 11:05:51 | 显示全部楼层
加油,加油!
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
0
注册时间
2020.2.12
在线时间
69 小时

发表于 2021.4.17 13:12:33 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

快速回复 返回顶部 返回列表