NBIS单细胞教程：降维（二）

此文章是通过学习瑞典国家生物信息学基础设施(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")