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

写在前面

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脚本结合,搭建实时可交互的网页/插件工具确实非常强大,后续也会慢慢摸索,尝试将单细胞分析中常用的几个包写成插件,顺便继续磨炼一下自己的能力。。 我果然还是太菜了。。

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

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 194,491评论 5 459
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 81,856评论 2 371
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 141,745评论 0 319
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 52,196评论 1 263
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 61,073评论 4 355
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 46,112评论 1 272
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 36,531评论 3 381
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 35,215评论 0 253
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 39,485评论 1 290
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 34,578评论 2 309
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 36,356评论 1 326
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 32,215评论 3 312
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 37,583评论 3 299
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 28,898评论 0 17
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,174评论 1 250
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 41,497评论 2 341
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 40,697评论 2 335

推荐阅读更多精彩内容