分享

Seurat | 我用「TBtools」学会了「单细胞测序」数据分析

 生信药丸 2023-01-09 发布于贵州

写在前面

TBtools已经拆除了绝大部分生信分析的门槛,让众多湿实验工作者们(无代码基础)可以畅快地分析他们的数据。

不求人非常爽!

不用跟公司低效沟通分析细节非常爽!

很早之前陈师兄就跟我提过是否有兴趣整一个单细胞分析方面的插件,一方面由于我R用的非常不6,另一方面确实是懒。。。所以搁置了。。近期又提了一下,于是菜菜的我终于踏上了一段辛酸的写bug之旅。。

总的来说,第一次接触shiny开发体系,确实非常不习惯,跟之前数据库整的html+css+JavaScript+php完全不一样,特别是在debug上,菜狗如我,太难了。。。不过shiny确实在实现可交互功能的时候非常强大,能够在ui端可视化地实时控制参数进行分析或者自定义绘图。

磕磕绊绊,尽管质量有待提高,最终还是完成了这个插件,Seurat Plugin for TBtools。将Seurat的基本分析流程进行了打包与可视化,能够在ui界面上直接完成从读入标准表达矩阵到数据质控、降维,然后进行细胞分群、marker基因鉴定与可视化等基本的single cell分析流程。

Seurat Plugin for TBtools

插件安装方式

注意到这是一个 TBtools R-plugin。与其他R插件类似,如果没有安装 Rserver插件,那么需要先安装 Rserver 插件(可以直接在 TBtools Plugin Store 找到并安装)。

安装完成后,打开插件。进入方式非常简单,你只需要点一下 Start 。

1. 上传数据

输入表达量数据,目前支持两种格式的表达量文件

  1. 基于Cell Ranger输出的稀疏矩阵,包括三个文件:barcodes.tsv,genes.tsv,matrix.mtx。

  1. count矩阵。一般在GEO数据库或者其他公共数据库存储平台中(比如浙大樊龙江老师团队开发的PlantscRNAdb)下载单细胞公共数据,大部分都是count矩阵,就是这个文件可能会比较大。。。同样,我们选择count data选项,然后上传相应的文件,注意需要是csv格式,即逗号分隔符,如果是tab分割的,可以TBtools或者notepad++做一下转换。

  1. 我们也准备了一个demo data,使用的是PRJNA497883这套公共数据(拟南芥根系)。

2. 数据预处理以及创建Seurat对象

上传完数据之后,在右边的“Data Preview”中会展示count矩阵,用户可以自行浏览或者检索。由于单细胞数据量(细胞数目)一般非常大,因此我们只展示了前30个细胞的数据。

接下来对数据进行预处理,过滤掉一些不表达的基因以及基本没有表达基因的细胞,一定程度可以去除一些异常值,加速下游的计算。通常这一步我们过滤两个指标:

  1. min.cells,保留至少在几个细胞中检测到表达的基因,默认为3个细胞,可以根据自己样本测序的实际细胞数量来进行调整,目的是为了去除一些几乎没有表达的基因。

  2. min.features,保留至少检测到有多少个基因表达的细胞,默认为200个基因。如果表达细胞表达的基因过少,那么很有可能是一个异常细胞,将其过滤。


预处理完数据之后,后台自动对其构建Seurat对象。

质控

质量控制。这一步我们对数据的一些QC指标进行可视化,用户可以根据可视化图片判断自己的数据是否正常。此外,我们提供了一个接口,用户可以在QC step1中计算一些基因在所有细胞中表达的比例,例如计算线粒体基因表达比例,其是过滤异常细胞的一个常用指标。一般我们认为线粒体基因表达过高的细胞,是一个处于应激或者其他不好状态的细胞,需要将其过滤。也可以计算核糖体基因等的表达情况。

Step1:计算线粒体(或者其他基因)的表达比例

这里我们以线粒体基因为例。拟南芥基因组中(TAIR10)的线粒体基因为ATMG开头,因此我们可以用正则表达式“^ATMG”来匹配拟南芥中所有线粒体基因。如果是人类,线粒体基因则以MT-开头,小鼠以mt-开头。我们提供了两种基因检索方式:

  1. 正则匹配,匹配线粒体基因等具有特定模式的基因集,例如拟南芥中:^ATMG

  2. 用户自定义基因集,用户输入自定义的基因集进行计算,每行一个基因。

Tips:如果你研究的物种并不是拟南芥等模式物种,我们该怎么确定它的线粒体基因呢?

  • 到Ensembl数据库或者其他的物种数据库中找到该物种的注释文件,在Ensembl中通常会有一个线粒体基因注释文件,下载到本地,打开你就可以发现这个物种的线粒体基因ID是否有规律,如果有固定的开头,即可以用正则表达式去匹配,如^ATMG,如果没有,那么就TBtools提取所有的线粒体基因ID,然后使用mode2,输入线粒体基因ID集,进行计算。

Step2:QC指标的可视化

根据QC指标过滤异常细胞

在Step2中,我们对几个QC指标进行了可视化,其中包括Step1中用户指定计算的基因表达量(如线粒体)。可视化的指标如下:

  1. 每个细胞的线粒体基因比例(如果用户有在Step1中计算)

  2. 每个细胞检测到的基因数目(nFeature_RNA)

  3. 每个细胞测得的UMI总数(nCount_RNA)

接着根据线粒体基因比例每个细胞检测到的基因数目进行细胞过滤。线粒体基因表达过高的细胞一般是异常细胞,而检测到基因数目异常过高的细胞,大概率是双胞(doublets,一个文库中存在两个细胞的情况,会对后续的分群和细胞鉴定造成影响和误导,需要剔除。当然有专门的doublets预测软件,如DoubletFinder等)

注意:示例数据使用的是拟南芥根尖的公共数据,该套数据已经经过质控,因此QC指标基本正常,并不存在异常值。

查看单细胞数据的质量

我们也对特征值进行了相关性分析,并可视化。用户可以根据nFeature_RNA和nCount_RNA的相关性进行大致判断自己的数据是否正常,理论上一个细胞测得的UMI越多,其测得的基因应该也越多,因此这两个指标应该是强正相关。

数据标准化与归一化

在对数据进行降维之前,需要对数据进行处理,即三步法:Normalization,Find Variable Genes以及Scale。通常情况下,默认只是标准化2000个高变基因,运算速度更快,不影响 PCA降维和细胞分群。对数据进行scale,可以消除极高表达的基因的影响。在Seurat plugin中,用户可以直接点击start按钮进行一键化跑三步法。

PCA线性降维

对处理过后的数据进行pca降维,默认计算前50个主成分。用户点击PCA降维按钮之后,会实时可视化展示PCA结果。

细胞聚类

在进行细胞聚类之前,我们首先要确定好两个参数:1,使用的PC个数;2,Resolution(分辨率)。

  1. 对于PC数目来说,我们需要使用合适的PC数目来进行聚类,太多的PC可能会对聚类造成一定的干扰,太少的PC又不能很好的解释数据集。因此我们推荐使用Elbow Plot来确定使用的PC数目。在Elbow Plot选择处于拐点处,即趋近于水平处的PC个数,来进行聚类,这些PC能够很大程度上解释整个数据集。

  2. 对于Resolution。Resolution参数决定下游聚类所得到的分群数,对于3千左右的细胞量,seurat官方说0.4-1.2能得到较好的结果。如果数据量增大,该参数也应该适当增大。那么,到底选择多少的Resolution最为合适呢?这里我们推荐使用clustree来确定Resolution。clustree能够进行多次不同的Resolution分群,每个不同的Resolution的分群结果进行可视化,用户可以根据Resolution的变化而得到的不同细胞群数,来确定适合你数据的最佳Resolution。

个人项目经验,如果Resolution选取较小,则分群较少,会将一些相似的细胞类群合并成一个大的cluster,从而对后续的细胞类型鉴定造成困扰,丢失一些细胞类型信息。而Resolution过大,则会将cluster分的过多,加大后续细胞类型鉴定的难度。因此我们一般会将Resolution选择在中等稍微偏大一点的值,鉴定细胞类型的时候,可以将具有相同marker基因,且在umap/tSNE上距离相近的cluster归为同一个细胞类型。

确定好PC数目和Resolution之后,则可以点击按钮进行细胞聚类。

UMAP/tSNE非线性降维

细胞聚类完成之后,用户可以点击按钮进行UMAP/tSNE非线性降维,并对其进行可视化。

marker基因鉴定/差异表达分析

marker基因鉴定,本质上其实是进行差异分析,鉴定每个cluster的差异表达基因(差异上调,往往cluster的marker基因是只在该cluster中特异表达,而在其他cluster中不表达。因此鉴定marker的时候,我们通常会设置only.pos=TRUE,即只保留差异上调基因)。在该模块中,提供了三种鉴定marker基因的模式。

  1. 一键化鉴定所有cluster的marker基因(比较慢,运行时间与细胞数和cluster数成正比),即对每一个cluster都与剩余其他cluster进行差异分析;

  2. 用户指定鉴定某个cluster的基因(快速),即对选定的cluster与其他所有cluster进行差异分析;

  3. 用户指定鉴定某两个cluster之间的marker基因,实际上就是在对这两个选定cluster之间进行差异分析,选用这种模式我们就没有必要只保留差异上调基因。

marker基因可视化

鉴定好marker基因之后,用户可以在该模块中,对感兴趣的marker基因进行各种可视化,分析其是否确实在某些cluster中特异表达,而在其他cluster中不表达。

  1. 对每个cluster中的Top5 marker基因进行热图展示,用户可以选定需要展示的marker基因数目

  1. 对每个cluster中的Top5 marker基因进行Dot plot展示,用户可以选定需要展示的marker基因数目

  1. 用户输入基因或基因集,对其进行Violin plot可视化

  1. 用户输入基因或基因集,对其进行Feature Plot可视化,将marker基因的表达量映射到UMAP或者tSNE上,可以非常清晰地看到marker基因在cluster中的表达分布情况

  1. Ridge plot。

写在最后

磕磕绊绊,最终还是写完了这个插件。虽然写的确实很烂,功能也非常简单,但在写这些代码的折腾中,确确实实能感受到收获了新的东西,接触到了新的网页工具开发方式(shiny跟html和php等网站开发确实不一样。。。) 。Seurat Plugin for TBtools,确实还有许多需要改进的地方以及新增的功能,后续也会抽空慢慢将其完善,副教授说“优秀的产品都是一步一步优化的”哈哈,希望自己也能像产品一样,”一步一步优化“!

接触了shiny才发现,用其与R脚本结合,搭建实时可交互的网页/插件工具确实非常强大,后续也会慢慢摸索,尝试将单细胞分析中常用的几个包写成插件,顺便继续磨炼一下自己的能力。。我果然还是太菜了。。

博士生涯转眼已过一半,临近过年,希望各位老板舒舒服服过个好年,来年再战!

    转藏 分享 献花(0

    0条评论

    发表

    请遵守用户 评论公约

    类似文章 更多