YangLiBMBL的个人博客分享 http://blog.sciencenet.cn/u/YangLiBMBL

博文

单细胞转录组测序的完整分析流程

已有 9203 次阅读 2021-6-5 09:16 |个人分类:科研笔记|系统分类:科研笔记

    作为数学出身的生物信息学研究员,设计模型和开发算法是我们的强项。但是短板在于不得不依赖合作者提前将数据处理成我们所“认识”的格式,例如:矩阵、表格和文本,等等。众多生物数据中,转录组数据始终占据“C位”。因此,今天我来详细总结一下从拿到测序数据到完成基本数据分析(例如:整合、降维和聚类)的步骤。

    目前多数测序技术原始数据保存在多个fastq.gz文件中。每个fastq.gz文件中包含大量的barcode,可以等同为细胞。以fastq.gz文件名“S3_CKDL210002749-1a-SI_TT_A2_H2NV5DSX2_S1_L001_I1_001.fastq.gz”为例。文件名中的“S1”就是所属的Sample名称。接下来我将分步介绍数据处理的流程:

    1. 创建目录。一个Sample通常包括多个文件,因此我们的第一步就是将属于同一Sample的文件放置到一个目录下,该操作可以通过Shell、Python或者Perl等脚本语言完成。目录整理完毕,我们需要创建一个.txt文本文档来存储目录的名字,如图1所示。

image.png

图1:Sample文件目录。

    2. 获取.bam文件。我们用Cellranger软件来获取.bam文件。Cellranger可以同时对多个Sample进行处理,因此这一步建议用Linux并行任务解决。具体的命令行如下:

CellRanger count --id=$NAME --transcriptome=${MM10_Refer} --fastqs=${FastqFolder}/$NAME --sample=$NAME --localcores=16 --localmem=128

其中$NAME指的是Sample的名称,${MM10_Refer}是参考基因组.gtf文件的名称。经过这一步,每个Sample对应的目录下会生成一个.bam文件,一个.bam.bai序号文件以及一个.h5文件。其中,这个.h5文件就是每个Sample下生成的基因计数矩阵。接下来我们需要用批处理工具,例如Shell、Python或Perl把所有的.bam文件汇集到一起。

    3. 读取矩阵文件。对于下面的操作,我们转战R语言。首先,加载所需要的各种库,包括:Seurat(不用说了,大家都认识,这是单细胞领域的国际标准),cowplot是ggplot2的一个辅助库,用来绘制高质量的论文图片,dplyr(提供了大量的数据框处理工具),ggplot2(不用说了,绘图界的旗舰),patchwork(用于组合各种各样的复杂图片),qs(提供了一种比.rds更快速的变量读取格式.qsave),here(用于锁定项目空间和工作路径)和 Polychrome(提供了精美的调色盘)。加载各种库之后,我们那需要用函数here::set_here()设定当前的工作目录。Seurat对.h5文件的读取是每个Sample逐个进行的,如果你有8个Sample,那就需要用函数Seurat::Read10X_h5()分别读取8次,每个.h5文件用一个变量保存。然后,我们用函数Seurat::CreateSeuratObject()将这些读取的矩阵转化成Seurat对象。最后,用函数Seurat::merge()将这些Seurat对象合并成一个。记住,合并的同时别忘了记录每个对象来自哪个Sample。

    4. 质量控制(Quality Control,QC)。质量控制首先从统计核糖体RNA和线粒体基因的表达开始。原理很简单:就是对不同基因的名字检查它的前缀。例如:人类核糖体RNA具有前缀“RP”,有的情况下具有前缀“RPSL”;老鼠的核糖体具有前缀“Rp”或“Rpsl”。同理,人类的线粒体基因具有前缀“MT-”,而老鼠的线粒体基因具有前缀“mt-”。我们分别用函数grep()判断每个基因名称的前缀是否符合上述条件,然后返回一个bool向量。在对基因和细胞进行筛选之前,通常先用函数Seurat::VlnPlot()绘制一个小提琴图,用来展示nFeature_RNA、nCount_RNA、percent.mito和percent.ribo。其中,nFeature_RNA和nCount_RNA分别指的是每个基因下RNA分子的数量(UMI的数量)和所有基因上UMI的数量;percent.mito和percent.ribo分别表示线粒体基因上的UMI和核糖体RNA对应的UMI占所有UMI的比例。QC的目的是为了去除Outlier数据,R语言中恰好含有自动确定Outlier的函数,即:scater::isOutlier()。严格地说,QC应该分别针对不同的Sample采用不同的标准,但是为了效率,人们大都对所有的Sample采用统一的标准,这也是可行的。最后一步,用函数subset完成对Outlier的删减。为了作对照,做完QC之后,我们会再次绘制一张小提琴图。Seurat还提供了一个新的功能——去除细胞周期异质性。Seurat对每个细胞计算它处于各个细胞周期的分数,然后将细胞周期的影响去除。Seurat所用的函数是CellCycleScoring(),具体原理在此不详细介绍。

    5. 数据整合。不同Sample可能测于不同的时间和设备,因此批次效应(Batch Effect)难以避免。我们借用Seurat v3中的CCA算法将多个Sample的数据整合在一起。首先,我们用SplitObject()将这个Seurat对象根据Sample拆分成若干个对象,并保存到一个列表里。然后用Seurat::FindVariableFeatures()搜索每个Sample上的可变基因(Variable Gene),通常选取前2000个。基于这些可变基因,我们用Seurat v3中的FindIntegrationAnchors()把所有Sample“锚定”到一起。数据整合就此完成。

    6. 后续分析。下面就可以做后续分析了,例如:降维和聚类。降维常用的函数有RunPCA、RunUMAP和RunTSNE。通常,我们会先对数据用PCA降维,然后在此基础上用UMAP或TSNE降维。聚类用的函数是FindNeighbors()和FindClusters()。这两个分析大家一定很熟悉了,在此不再赘述。

    据我所知,去除细胞周期异质性的操作是Seurat后期追加的功能。因此作为生物信息领域的工作者,要牢记数据分析的流程一定要紧随潮流,不能“刻舟求剑”,否则投稿中所有不完备的操作都会被审稿人一一列举出来。

    如果想详细了解单细胞数据分析的完整细节,请参考这篇文章:https://abego.cn/2019/08/17/the-scrna-sequence-analysis-pipe-and-theory/。



https://wap.sciencenet.cn/blog-3447504-1289807.html

上一篇:肿瘤微环境基础知识汇总
下一篇:scATAC-Seq数据预处理和基础分析流程
收藏 IP: 144.121.166.*| 热度|

1 许培扬

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...
扫一扫,分享此博文

全部作者的精选博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-3-29 06:41

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部