单细胞转录组测序分析--初探Seurat

单细胞转录组测序--闲言碎语

时代发展的步伐总是毫不留情的将你甩在身后，连车尾灯都看不见。当你还在沉迷于普通转录组数据挖掘时，已经有人悄悄的搞上单细胞了。单细胞转录组测序，顾名思义就是在单个细胞的分辨率基础上去研究细胞内的基因表达等，其主要目的是为了研究不同细胞类型的基因表达异质性，从而解决相关生物学问题。谈到单细胞就不得不提一下当下火爆的10x Genomics服务商了，具体参见10x Genomics。本篇文章暂时不介绍10x，主要介绍单细胞转录组数据分析软件Seurat。
Seurat软件是一个R包，可以说是单细胞转录组测序分析的明星软件，很多单细胞测序文章都会引用该软件，引用次数也是杠杠的，而且也有详细的在线教程。本文也主要是根据其教程介绍一下使用Seurat软件分析一个样本的单细胞转录组数据的步骤及注意事项，供大家讨论。

360截图20190220211931988.jpg

Seurat软件介绍

导入分析需要的包

library(Seurat)
library(dplyr)
library(magrittr)

Seurat软件提供了很友好的函数可以直接读取10x Genomics的输出结果

matrix <- Read10X_h5("filtered_gene_bc_matrices_h5.h5")#该函数直接读取10x的h5文件，
#并返回一个稀疏矩阵
#也可以使用下面这个函数导入数据
matrix <- Read10X("data/matrix_dir")#data/matrix_dir目录下包含文件
#matrix.mtx, genes.tsv, and barcodes.tsv

导入文件后便可以创建Seurat对象

filter_10x_object <-CreateSeuratObject(raw.data = matrix,min.cells = 3, 
min.genes = 200, project = "d615")#该函数还具有数据过滤的功能，
#基因至少在min.cells个细胞中表达，每个细胞中至少表达min.genes个基因，
#可以根据具体情况决定这些阈值

创建完Seurat对象后，Seurat将数据保存在不同的slot中，如filter_10x_object@raw.data, filter_10x_objectt@data, filter_10x_object@meta.data, filter_10x_object@ident，其中raw.data存放的是每个细胞中每个gene的原始UMI数据，data存放的是gene的表达量，meta.data存放的是每个细胞的统计数据如UMI数目，gene数目等，ident此时存放的是project信息。

由于技术原因，一个GEM中可能会包含2个或多个细胞，也可能不包含细胞，这时候可以通过观察每个barcode中的基因数目或UMI数目来判断。

GenePlot(object = filter_10x_object, gene1 = "nUMI", gene2 = "nGene")

1.jpg

上图展示的是每个barcode中的基因数目和UMI数目的关系，一般二者都成正相关关系，有个别barcode的基因数目和UMI数目过高，有可能就是包含2个细胞的GEM，可以考虑在后续分析中将其过滤掉。
我们不仅仅可以观察每个barcode的基因数目，还可以计算每个barcode中的线粒体基因含量等，从而更加仔细的观察数据的质量。

mito.genes <- grep(pattern = "^MT-", x = rownames(x = filter_10x_object@data), 
value = TRUE)#统计线粒体基因
percent.mito <- Matrix::colSums(filter_10x_object@raw.data[mito.genes, ])/
Matrix::colSums(filter_10x_object@raw.data)#计算线粒体基因含量
filter_10x_object <- AddMetaData(object = filter_10x_object, metadata = percent.mito,
col.name = "percent.mito")#将线粒体基因含量添加到meta.data中
VlnPlot(object = filter_10x_object, features.plot = c("nGene", "nUMI", 
"percent.mito"), nCol = 3)#绘制小提琴图

2.jpg

这张图片展示了每个barcode中基因数目、UMI数目以及线粒体基因含量的分布情况，根据上述2张图片就可以大致确定是否需要过滤哪些数据进行后续分析。
Seurat提供了一个很好用的数据过滤函数：

filter_10x_object_filter <- FilterCells(object = filter_10x_object, subset.names =
 c("nGene", "percent.mito"), low.thresholds = c(200, -Inf), 
high.thresholds = c(2500, 0.05))#过滤掉小于low.thresholds或大于high.thresholds的数据，
#这个参数值对应subset.names的参数

以上就是数据的预处理过程了，接下来就进入正式的分析阶段，包括数据的标准化、归一化、数据降维以及聚类分析等。

#数据标准化，即表达量的计算
filter_10x_object_norm <- NormalizeData(object = filter_10x_object_filter, 
normalization.method = "LogNormalize", scale.factor = 10000)#在每个细胞内，
#对每个基因进行标准化：该基因的UMI除以该细胞内转录物的总UMI并乘以10000，
#然后在使用log1p对这些值进行自然对数转换（标准化的数值存放在data slot中）

#计算高度变化的基因
#（并不是所有基因都有有效的信息，寻找高度变化的基因既可以节省计算资源，又能概括样本信息）
filter_10x_object_norm_var <- FindVariableGenes(object = filter_10x_object_norm,
 mean.function = ExpMean, dispersion.function = LogVMR, 
x.low.cutoff = 0.0125, x.high.cutoff = 3, y.cutoff = 0.5)#具体 cutoff值可以根据下图调整

3.jpg

FindVariableGenes算法：首先计算基因的平均表达量，然后计算基因的离散度；接下来根据平均表达值将基因分成20块并计算每块的离散度的Z值。
如上图：横坐标代表基因的平均表达量，纵坐标代表基因的离散度的Z值，标有基因名的点就是由函数中的cutoff值决定的，改变cutoff值，这些标记也会随之改变。
数据的线性回归、中心化和比例化：对数据进行线性回归分析，去除不想要的变异源。
中心化：首先计算基因A在所有细胞中的平均表达量，然后分别将每个细胞中基因A的表达值减去平均值。
比例化：在中心化的基础上，首先计算基因A在所有细胞中的中心化值后的标准差，然后分别将每个细胞中基因A的中心化值除以标准差。这些步骤都在一个函数中完成。

filter_10x_object_norm_var_scal <- ScaleData(object = filter_10x_object_norm_var, 
vars.to.regress = c("nUMI", "percent.mito"))
#vars.to.regress中的参数就是需要去除的变异来源
#该函数会先进行先行回归分析，然后再进行数据的中心化和比例化，并将值存在scale.data中

单细胞转录组测序产生的数据是数万个基因在数万个细胞中的表达情况，属于典型的高维数据。如果把1个基因视为1个坐标轴的话，那么一个细胞的空间位置就是在数万个坐标轴中的定位，这样的话相同细胞类型的细胞就应该挨在一起，我们就可以根据细胞的空间位置判断细胞亚群了。可是我最多也就认识三维坐标啊，咋办，能不能把这些高维数据投影到二维坐标呢，那就交给PCA和t-SNE吧。PCA和tSNE都是数据降维分析方法，PCA属于线性降维，tSNE属于非线性降维。我们先执行PCA分析，使高维数据的信息最大程度保留在低维数据中，PCA分析利用的是保存在scale.data的值。

filter_10x_object_norm_var_scal_pca<- RunPCA(object= filter_10x_object_norm_var_scal,
 pc.genes = filter_10x_object_norm_var_scal@var.genes, pcs.compute=20)#利用高可变基因
#计算20个主成分，并将数据保存在@dr$pca中

执行完PCA分析后，就要根据PCA得分来进行聚类分析了，但是在进行聚类分析之前，需要选择使用对少个主成分进行计算。每个主成分实际上代表的是相关基因集的信息，因此确定多少个主成分是一个重要的步骤，我们可以根据PCElbowPlot函数来判断。

PCElbowPlot(object = filter_10x_object_norm_var_scal_pca)

4.jpg

从上图可以看到，拐点出现在10-15之间，我们可以选择15来进行聚类分析。Seurat采用的是基于图形的聚类方法，即利用PCA空间中的欧几里德距离构造一个KNN图（数学好的可以留下来帮忙讲讲）。

filter_10x_object_norm_var_scal_pca_cluster <- FindClusters(object =
 filter_10x_object_norm_var_scal_pca, reduction.type = "pca", 
dims.use = 1:15, resolution = 1.0, print.output = 0, save.SNN = TRUE)
#resolution会影响聚类的个数

好了，到此我们就知道了我们的数据中有多少种细胞亚群了，怎么可以少得了图片展示呢。超棒的可视化方法tSNE要上场了。tSNE的目标是将在高维空间中具有相似局部邻域的细胞，在低维空间中放在一起。

filter_10x_object_norm_var_scal_pca_cluster_tsne <- RunTSNE(object = 
filter_10x_object_norm_var_scal_pca_cluster, dims.use = 1:15, do.fast = TRUE)
TSNEPlot(object = filter_10x_object_norm_var_scal_pca_cluster_tsne,do.label = TRUE, 
pt.size = 0.2)#绘图

5.jpg

既然我们知道了有多少种细胞亚群，那么是不是就要分析一下这些亚群间的差异性呢，交给FindAllMarkers吧。FindAllMarkers能够同时计算所有亚群的差异性（分别计算每个亚群与剩下的所有细胞的差异性）。

pbmc.markers <- FindAllMarkers(object = 
filter_10x_object_norm_var_scal_pca_cluster_tsne, 
only.pos = TRUE, min.pct = 0.25, thresh.use = 0.25)#根据data slot的值进行计算

得到差异表达基因后，当然要进行展示了。

6.jpg
7.jpg

好了，剩下的就是进行生物学知识挖掘了，例如根据这些差异基因推断细胞类型啊之类的。
关于单个样本的单细胞转录组数据分析就介绍到这儿了，那多个样本的分析会有什么不同呢，我们下次再说吧。

单细胞转录组数据分析的小伙伴一起来讨论吧