一.单细胞转录组简介¶
单细胞 RNA 测序(Single cell RNA sequencing,scRNA-seq)是一种在单细胞水平上利用 RNA 测序对特定细胞群体进行基因表达谱定量的高通量实验技术。待测组织经过单细胞分离、RNA 提取、逆转录、文库构建和测序,便可利用数据分析获得多个细胞的基因表达谱。
1.单细胞测序与普通转录组测序的区别¶
普通转录组使用细胞混合物组成的样品进行测序,因此只能估计基因在细胞群中的平均表达水平,没有考虑样本中各个细胞的基因表达的异质性。无法分析早期发育组织或复杂组织的异质系统(如大脑组织等)。为克服这一限制,开发了单细胞水平的转录组测序技术(scRNA-seq)。
2.单细胞测序的原理¶
单细胞测序技术以单个细胞作为对象,通过对单个细胞遗传物质均匀扩增,标记建库后进行测序,最后对单个细胞基因组或转录组展开数据分析,其技术原理主要包括单细胞分离、扩增测序和数据分析3方面。 单细胞分离:主要包括荧光激活细胞分选法( flow- activated cell sorting ,FACS) 以及微流控分选法( microfluidics) 等。市场上,较成熟的商业单细胞测序公司主要有 10X Genomis 公司 的Chromium( 液滴法) 及 BD 公司的Rhapsody( 微孔法)。
(1) 10×genomics技术¶
首先在凝胶微珠上种上特定的DNA片段,DNA片段由三部分组成:Barcode、UMI、PolyT组成。Barcode是16个碱基的长度。一共有400万种Barcode,一个微珠是对应于一种Barcode,通过这400万种Barcode,可以把凝胶微珠给区分开。UMI是一段随机序列,也就是说每一个DNA分子,都有自己的UMI序列。10个碱基长的UMI,有100万种序列的变化(4^10 = 1,048,576),UMI的作用是为了区分哪些哪些reads是来自于一个原始cDNA分子,区分基因片段重复还是duplication及区分是真实的SNP位点还是PCR产生的突变。
3'端文库的构建¶
通过10×genomics仪器将单个细胞与单个凝胶微珠通过油相混在一起,形成油包水的小微滴,接下来把细胞膜破掉,让细胞当中的mRNA游离出来。游离出来的mRNA与小液滴中的水相混合,也就是和逆转录酶、结合在凝胶微珠上的核酸引物、以及dNTP底物相接触。接着,发生逆转录反应。mRNA与凝胶微珠上带标签的DNA分子相结合,在逆转录酶的作用下,逆转录出cDNA来。 把这个乳浊液当中所有的水相抽出来,也就是把所有带了标签的cDNA分子都抽出来,再把这些cDNA分子都加上接头,经过PCR扩增,做成illumina的测序文库,放到Illumina的测序仪上进行测序。
5'端文库的构建¶
与3’端文库的构建过程相似 小提示:barcode和UMI是单细胞中常见的名词,barcode就像是每个凝胶微珠的身份证号码;而UMI则像是每个DNA标签分子的身份证号码。
使用scRNA-seq,可以估计基因在细胞群的表达水平分布。 scRNA-seq 非常适合研究细胞群的异质性。例如,识别组织的细胞类型,定义不同细胞类型的"转录指纹",研究细胞分化,探索疾病或环境因素导致的细胞组成变化等。 单细胞测序的经典流程:分离单细胞(核),将RNA转换为cDNA,准备测序文库(illumina)后测序。 SMART-seq2是一种低通量单细胞测序法,提供全长转录组的定量,适合研究小细胞组(如差异isoform使用,低表达转录本的特征)。 10x Chromium是一种高通量方法,使用UMIs进行定量,适合研究高度异质组织和大量的细胞样本。
二. 如何开始单细胞RNA测序数据分析¶
在过去的 5 年中,单细胞方法能够在一次实验中监测数千个单个细胞的基因和蛋白质表达、遗传和表观遗传变化。随着测量的改进以及反应和测序成本的降低,这些数据集的大小正在迅速增加。关键瓶颈仍然是对单细胞实验产生的大量信息的分析。在这篇评论中,我们对分析管道进行了简化概述,因为它们通常在当今的现场使用。我们的目标是让研究人员从单细胞分析起步,以了解挑战和最常用的分析工具。此外,我们希望使其他人能够了解单细胞数据集的典型读数如何在已发表的文献中呈现。
1.背景¶
基于下一代测序的单细胞基因表达分析在1992年首次被描述。2015年,开发了基于encapsulation和barcode的分析。在过去的5年里,单细胞分析变得大众化,大多数学术机构都有专门的核心设施来进行单细胞表达、表观基因组或其他分析。新的统计分析方法得到了迅速发展。一些分析平台也已经开发出来,如Seurat,它是用R 语言写的(https://satijalab.org/seurat/get_started. html) ,用Python写的Scanpy (https://scanpy.readthedocs.io/en/ stable/tutorials.html)。在这里,我们回顾了基本的分析工具和概念。我们关注于10× Genomics数据,因为它们更常用(图1)。本文主要基于两个案例研究教程(https://www. github.com/theislab/single-cell-tutorial 和http://scrnaseqcourse.cog.sanger.ac.uk/website/index.html)。
2.数据矩阵生成和质量控制¶
单细胞分析的一个关键技术进步是barcode的发展,它允许大规模并行化,同时保持成本最低。barcode被添加到 在逆转录过程中的RNA分子,允许识别单个细胞和独特的分子。分析的第一步是生成一个数据矩阵,由转录本数据库从原始测序文件生成一个barcode(细胞)。对于10×基因组学数据,CellRanger(表1总结了文中提到的工具,方法和数据库)是最常用的管道,包括测序读数与基因组的注释和比对,将比对的读数注释到基因,和量化基因。替代方案包括,例如唯一分子标识(unique molecular identifier , UMI)工具、zUMIs、kallisto、STAR和STARsolo (https://github.com/alexdobin/STAR/blob/master/docs/STARsolo.md)。 每个barcode可以代表一个单细胞。一个doublets,或者一个不含细胞,但含有周围RNA的“空”液滴。需要提到的一个重要问题是,标准的测序管道将测序数据与转录组进行比对,例如处理后的成熟mRNA。然而,单核RNA数据或表观基因组数据(通过测序测定转座酶可获得的染色质,ATAC-seq)应与全基因组比对,因为细胞核主要含有premRNA,其包括内含子区域。Raw read计数通常也会过滤掉在极少数细胞中检测到的基因,从而有效地减少数据矩阵的大小。分析管道的下一步是质量控制(QC),例如识别每个barcode的counts数,每个barcode的基因数,以及每个barcode的线粒体基因counts 的比例(图2)。 低基因数量和高比例的线粒体reads通常表明细胞质量差。然而,一些细胞,包括肾近端和远端曲小管细胞,含有丰富的线粒体。异常高的read和基因count可能代表doublets。现在有几种doublet检测工具可用,包括DoubletDecon,Scrublet和DoubletFinder(图2)。然而,doublet检测中的关键问题是含有来自例如上皮和间充质来源的标记基因的过渡细胞可能被标记为doublet,有时导致假阳性检测。此外,这些工具只能识别低同型的双链体,即转录相似的细胞在基因表达空间中聚集在复合细胞类型单链体之间形成的双链体。控制环境RNA污染也很重要。环境RNA是存在于单细胞溶液中的RNA,在包裹过程中被整合到油滴中。我们通常使用SoupX,它可以从空液滴中估计周围的RNA污染(图2)。另一个包是CellBender,它可以消除来自周围环境的RNA分子和随机barcode交换的count(原始)基于UMI的单细胞RNA测序(scRNA-seq)的count 矩阵。
3.标准化¶
单细胞数据需要不同类型和级别的标准化 (图1)。例如,总测序read count会改变原始count,因此基因count应该按总count深度调整。一种常用的方法假设每个细胞具有相同的转录本初始数量,简单地将数据规范化为每百万分数。Scran使用基于pooling的大小因子估计和线性回归来对数据进行归一化,除了Seurat使用的简单对数归一化之外,它也是最流行的方法之一。其他方法已经开发,如SCtransform, SCnorm和BayNorm。归一化后,数据进行log(x + 1)转换。通常从数据中回归出与细胞周期相关的变化,它包含在Seurat或Scanpy的标准分析平台中。该平台还允许其他技术或生物学变异的回归。
4.批处理效果校正和数据集成¶
大多数情况下,会生成多个数据集,需要额外的批处理和数据集成方法(图3)。包含多个不同实验和不同方法的大型数据集通常使用非线性方法进行集成。在Seurat中,有一个基于参考的集成选项,它使用典型相关分析或互反主成分分析(PCA)。Scanorama是Scanpy中使用的另一种流行且性能良好的方法。最近,Harmony越来越受欢迎,并迅速成为单细胞数据集最常用的集成方法。第一步,利用PCA衍生的嵌入矩阵和批处理元数据进行缩放,使每个单元都有一个长度参数。然后,利用正则k-means聚类对尺度数据初始化聚类质心。最后,通过迭代地将特定于批次的质心拉到聚类质心,直到收敛,消除批效应。基因组实验关系的LinkedInference 通过整合非负矩阵分解识别共享的和数据集特定的因素。通过UMI的数量归一化后,基因表达被缩放但不居中。不同的集成方法可能会显示不同的结果。一般情况下,我们期望不同实验的相同细胞类型能够整合在一起,特别是控制细胞应该与多个实验保持一致。感兴趣的读者可以参考Tran等人和Chen等人最近发表的两篇优秀论文,它们提供了一些支持基因组和谐、关联推理的最佳证据实验关系和Seurat关于批效应校正。
5.可视化和聚类¶
可视化的第一步是在保留信息基因(1000–5000)并过滤掉其他基因时进行特征选择,这在Seurat和Scanpy中均已实现(图4和5)。可视化是在低维空间中总结数据集以观察模式的一种尝试。通常,降维是通过线性和非线性方法来实现的。PCA是聚类和轨迹推理的基础,是一种线性变换,在完整的PCA中保持细胞之间的欧氏距离。在常用的Seurat管道中,PCA用于预处理阶段。可以将PC投影到技术和生物协变量中以了解其性能。使用基于置换测试的jackstraw方法,对top PC和通过“elbow”启发式方法选择的PC数量汇总了PCA(图4)。 单细胞数据可视化主要使用其他非线性降维方法,如t分布随机邻域嵌入。该方法专注于以全局结构为代价捕获局部相似性。均匀近似和投影(UMAP)方法由于其速度而受到欢迎。UMAP似乎能更好地捕捉底层数据结构,并能在两个维度以上总结数据;因此,它现在最常用于单细胞数据可视化。UMAP和t分布随机邻域嵌入的一个关键限制是它们对用户自定义参数有很强的依赖性,并且结果对这些参数有很高的敏感性。最需要注意的是,这两种可视化都不能保持细胞-细胞之间的距离,因此结果嵌入不应该被下游分析直接使用(图4)。 根据基因表达的相似性形成的细胞群是分析的第一个直接结果。细胞聚类允许根据基因表达的相似性分组细胞来推断细胞类型。聚类是一种基于距离矩阵的无监督机器学习过程。community中默认的聚类方法是使用单细胞K近邻方法进行Louvain community检测。单元在图中表示为节点。每个细胞与它的K个最相似的细胞相连接,这些细胞通常是使用pc简化表达式空间上的欧氏距离获得的。单元格在图中表示为节点。每个单元格连接到其K个最相似的单元格,这些单元格通常使用PC缩减表达式空间上的欧几里得距离获得。一个关键问题是,用户决定Louvain聚类的分辨率,而分辨率决定数据集中识别的簇或单元类型的数量。我们建议执行子聚类,例如从初始数据集设置某些聚类,然后在不使用其他单元类型的情况下重新聚类。这允许在细胞类型内出现更精细,更细粒度的数据结构。包含在Leidenbase软件包中的Leiden community检测算法是Louvain算法的替代方案,在Monocle轨迹分析中用作默认值(见下文)。新的聚类方法使用神经网络和人工智能,例如,单细胞聚类的深度嵌入使用深度神经网络,网络权值和初始聚类从一个自动编码器获得。Clusters(簇)并不一定意味着细胞类型。这一点非常重要,因为用户定义的集群解析参数决定了观察到的集群的数量。注释和聚类是强链接的。聚类和注释以迭代方式进行,这是耗时的。目前,关于最佳聚类参数尚未达成共识。因此,可以接受多个版本的聚类和相同数据的解释。Wilcoxon 秩和检验用于通过组间表达差异对基因进行排序。经典的单元类型注释使用外部数据集,该数据集被认为是基本事实。越来越多的肾脏细胞类型注释外部数据集包括Susztak Lab,Humphreys Lab,Tabula Muris, Human Cell Atlas, Renal Epithelial Cell Ontology webpage, 和 ImmGen Consortium。最近,自动化细胞注释工具已经开发出来,如Garnett、SingleR、CHETAH和MOANA,他们提供了一种更全面和概率的细胞身份标注方法。相同细胞类型的标记基因在不同的数据集可能不同。
6.细胞水平分析:细胞组分变化,分解和轨迹分析¶
细胞组分的变化(数据集中每种细胞类型的比例)与疾病状态有很强的相关性,这是单细胞分析最简单的结果之一。这些数字可以提供条件之间的相对估计,但由于单细胞库制备过程中细胞捕获的偏差,从单细胞数据推断的细胞分数可能不准确。此外,来自肾皮质的样本中,近端小管细胞的比例比来自髓质的样本高。为了推断bulk RNA-seq数据的细胞类型组成,MuSiC是最近开发的一种以单细胞表达数据为参考的批量组织细胞类型反卷积方法。MuSiC使用加权非负最小二乘回归估计细胞类型比例。可供选择的方法包括CIBERSORT、BSEQ-sc和BisqueRNA。对样本之间的细胞识别簇比例变化的统计检验是相互依赖的,而且,因为当一个细胞识别簇的比例变化时,所有其他细胞识别簇的比例也会变化。或者,基于排列的统计测试方法可以用于差异比例分析,其中cluster比例与总细胞的随机比例进行比较。 细胞多样性不能充分地描述一个离散的分类系统,如聚类。轨迹分析捕捉转变期间细胞的显著特征,例如沿着几个时间点的器官发育期间,或疾病状态、细胞历史或拓扑信息之间的变化。驱动观察到的异质性的生物过程是连续的。因此,捕捉细胞特性、分支分化过程或生物学功能的渐进、不同步变化之间的转变需要基因表达的动态模型。Monocle是一种机器学习方法,用于重建每个细胞从一种状态转换到另一种状态时必须执行的基因表达变化序列。它基于逆图嵌入,一种高度可扩展的非线性流形学习技术。在该方法学习转换( transition)路径或轨迹之后,它将每个细胞沿着它放置在正确的位置,这被称为伪时间,即细胞在生物过程中移动多远的度量。一种分析细胞历史的新方法在最近开发的RNA速度分析中得到了应用,如包velocyto。RNA速度是基因表达状态的时间衍生物,可以通过在常见的scRNAseq协议中区分未剪接和剪接的mRNA来直接估计。RNA速度是一种高维向量,以小时为时间尺度预测单个细胞的未来状态。TradeSeq基于称为Slingshot的先前方法,在简单轨迹分析方面优于其他方法。另一个有用的软件包是PHATE,它是一种可视化方法,使用数据点之间的信息几何距离捕捉局部和全局非线性结构。推断出的轨迹不一定代表一个生物过程,应该收集进一步的证据来源来解释从这些方法得到的轨迹(图6)。
7. 基因水平分析:差异表达,基因调控网络,驱动途径和细胞-细胞相互作用¶
差异表达(DE)分析是通过包括技术和生物协变量的未校正数据进行的。Seurat使用不同的模型进行DE分析(图6)。MAST使用障碍模型来解释辍学。为了将scRNA-seq数据集信息与其他表型变量相关联,基于回归的模型可以组合几个样本及其相关的表型特征,以将某些细胞类型(例如近端肾小管细胞)中的基因表达变化与相应的定量测量表型(例如GFR,Albuminuria)。尽管DE测试工具通常允许用户灵活地合并混杂器,但用户必须警惕哪些变量被添加到模型中。例如,在大多数单细胞实验设置中,样本和条件协变量是混淆的,因为在多种条件下获得单个样本几乎不可能。基因水平分析还可以与基因集富集分析方法相结合,例如基因集富集分析或加权相关网络分析。 为了解释DE结果,我们通常根据参与共同生物过程的基因进行分组。生物过程标签存储在数据库中,例如MSigDB,基因本体论或《京都议定书》的基因和基因组百科全书以及Reactome数据库。尽管我们需要记住,某些通路成员的基因表达富集可能并不一定与通路活性相关,但基因列表上注释的富集可以使用大量工具进行测试,这些工具已经在其他地方进行了审查和比较。单细胞分析领域的最新发展是使用配对基因标签来进行配体-受体分析。在这里,细胞群之间的相互作用是从受体及其同源配体的表达推断。配体-受体对标签可以从最近的数据库获得,例如CellPhoneDB或Connectome,并用于使用统计模型解释跨簇的高表达基因。
8.单细胞分辨率下的基因调控¶
通过测序获得转座酶可及染色质的单核分析(snATAC-seq)可以通过染色质可及性图谱分析单细胞的表观基因组景观(图7)。已经为snATAC-seq分析开发了多种工具。最着名的是Ren实验室开发的SnapATAC,Satija实验室开发的Signac和Greenleaf实验室开发的ArchR。我们更喜欢SnapATAC,这是一种非线性降维方法。在CellRanger中通过barcode-by-cell 矩阵生成barcode之后,我们通过使用SnapATAC将片段二值化为均匀大小的逐个cell-by-bin来预处理矩阵。QC步骤包括过滤读深过低或过高的低质量细胞或doublet,并移除基因组黑名单区中的读码。重要的QC标准是转录起始位点的富集、峰中reads的比例和启动子区reads的比例。SnapATAC利用扩散图识别异质组织中的细胞类型。从扩散图中得到的低维嵌入作为输入,以消除批处理效应。然后用Louvain算法进行聚类,选取k个最近邻算法中的k个值作为输入。为了进行聚类注释,从选定的肾细胞类型特异性标记基因中生成细胞基因活性评分矩阵。预定义的启动子区域,(例如,来自Ensembl调控构建)或基因体12kb区域用于整合所有与基因转录本重叠的片段。为了从每种细胞类型调用峰,将从相同细胞类型获得的所有片段聚合以构建伪批量ATAC数据集和MACS2,分别针对每种细胞类型进行。ArchR实现了一种改进的方法,在独立样本上调用峰值,然后保留可重现的峰值。Fisher在edgeR中的精确测试在细胞簇之间进行测试,以揭示每种细胞类型的不同可访问区域。为了鉴定不同细胞类型中富集的基序,HOMER或chromVAR可用于snATAC-seq数据的转录因子分析,尽管遗传背景将严重影响哪些转录因子基序富集。 为了研究开放染色质变化如何与细胞分化和细胞命运决定相关联,使用Monocle3进行轨迹分析,通过使用潜在语义索引和UMAP可视化来减少维度。为了了解开放染色质和目的基因的表达变化,通过分析在Cicero中实现的两个峰的可达性,进行了峰-峰相关性研究。该策略聚合相似的单元以获得一组“元单元”,并解决snATAC-seq数据中的稀疏性问题。或者,可以将峰推算到GREAT中以识别最近的基因。
9.Webtools和数据集¶
在过去的几年中,已经生成了大量的人类和小鼠肾脏数据集。原始数据集通常可从GEO下载。作为人类细胞图谱项目和人类生物分子图谱计划的一部分,将提供大量综合参考人类肾脏注释。肾脏精准医学项目旨在为各种人类肾脏疾病状况生成数据集。重建肾脏联盟将分析开发的人类肾脏样本和体外分化的肾脏类器官。此外,一些研究人员已经为小型单个实验数据集生成了可视化工具。Humphreys实验室的工具包站点允许快速可视化他们广泛的数据(http://humphreyslab.com/SingleCell/).)。McMahon和Kim实验室使用VisCello将发育和成年小鼠的数据可视化,同时比较雄性和雌性小鼠的差异。使用VisCello平台,我们的实验室可视化开发,成人,健康和疾病小鼠模型数据(http://susztaklab.com/VisCello/)35)和开放的染色质表观基因组数据,这也是可用于相同的时间点(http://susztaklab.com/developing_adult_kidney/igv/).)。这些站点不允许进行全面分析,并且聚类参数(有些主观)是固定的,但它们对于查找和比较非常有用。数据分析自动化的另一个重要发展是Satija实验室的独立分析应用程序(http://azimuth.satijalab.org/app/azimuth)),允许将单个数据集投影到参考数据集,其中已在单个细胞中同时测量了RNA和表面蛋白表达。研究者可以上传本地生成的数据集,该软件包自动执行上述所有步骤,并与人体血液参考数据聚类。
三. 单细胞转录组的上游分析¶
1.从BCL到FASTQ¶
(1) 软件安装和介绍¶
Cellranger mkfastq 管道可用于将 BCL 文件解码为单个库的 FASTQ 文件。如果测序提供程序已经完成了这一步,则可以直接使用每个库的 FASTQ 文件进行数据分析。cellranger mkfastq的本质是调用bcl2fastq生成bcl2fastq,并生成额外的10x样本信息。Bcl2fastq是 illumina开发的bcl到fastq的转换程序。cellranger下载安装地址如下: https://support.10xgenomics.com/single-cell-gene-expression/software/downloads/latest? 安装cellranger 参考:https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/tutorial_in
tar -xzf cellranger-6.1.1.tar.gz
echo 'export PATH=~/softwares/cellranger-6.1.2:$PATH' >> ~/.bashrc
source ~/.bashrc
cellranger testrun --id=tiny
#https://anaconda.org/dranew/bcl2fastq
conda install -c dranew bcl2fastq
(2) 案例¶
公司提供的bcl原始测序文件,本次教学使用使用cellranger的测试数据
curl -O https://cf.10xgenomics.com/supp/cell-exp/cellranger-tiny-bcl-1.2.0.tar.gz
curl -O https://cf.10xgenomics.com/supp/cell-exp/cellranger-tiny-bcl-simple-1.2.0.csv
tar -zxvf cellranger-tiny-bcl-1.2.0.tar.gz
tree -L 2 cellranger-tiny-bcl-1.2.0/
cat cellranger-tiny-bcl-simple-1.2.0.csv
查看帮助文档
cellranger mkfastq --help
cellranger mkfastq --id=tiny-bcl --run=cellranger-tiny-bcl-1.2.0 --csv=cellranger-tiny-bcl-simple-1.2.0.csv
tree -L 2 tiny-bcl/
运行时间根据系统资源的不同而不同,但是它不应该占用多于几分钟的时间。现在转到 fastq _ path 目录。
cd tiny-bcl/outs/fastq_path/
cellranger mkfastq
管道的输出作为 cellranger count
的输入。
2. 从SRA到FASTQ¶
(1) 软件安装和介绍¶
SRAtoolkit是NCBI提供的SRA文件处理工具集, SRA文件是NCBI的SRA数据库数据的储存格式,许多公开的scRNA-seq数据都会上传到该数据库。SRAtoolkit将NCBI的SRA数据库中SRA文件转换为FastQ文件。 在conda的环境中安装SRAtoolkit,地址:https://anaconda.org/daler/sratoolkit
conda install -c daler sratoolkit
curl -O https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/3.0.0/sratoolkit.3.0.0-centos_linux64.tar.gz
tar zxf sratoolkit.3.0.0-centos_linux64.tar.gz
##注意解压位置
echo 'export PATH=~/sratoolkit.3.0.0-centos_linux64/bin:$PATH' >> ~/.bashrc
source ~/.bashrc
vdb-config -i
(2) 案例¶
数据地址:https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA483959&o=acc_s%3Aa 点击Total的“Accession List”下载整个研究的SRA文件编号; 点击Total的“Metadata”下载网页中SRA文件描述。
cd /mnt/f/Linux/scRNASeqTest/
mkdir sra_data
cd sra_data/
curl -O https://sra-pub-run-odp.s3.amazonaws.com/sra/SRR7722937/SRR7722937
curl -O https://sra-pub-run-odp.s3.amazonaws.com/sra/SRR7722938/SRR7722938
curl -O https://sra-pub-run-odp.s3.amazonaws.com/sra/SRR7722939/SRR7722939
curl -O https://sra-pub-run-odp.s3.amazonaws.com/sra/SRR7722940/SRR7722940
curl -O https://sra-pub-run-odp.s3.amazonaws.com/sra/SRR7722941/SRR7722941
curl -O https://sra-pub-run-odp.s3.amazonaws.com/sra/SRR7722942/SRR7722942
fastq-dump
命令可将sra数据转换成fastq。
fastq-dump --gzip --split-files -X 25000 -O fastqs_path data/SRR7722937
fastq-dump --gzip --split-files -O fastqs_path data/SRR7722937
fastq-dump --gzip --split-files -O fastqs_path data/SRR7722938
fastq-dump --gzip --split-files -O fastqs_path data/SRR7722939
fastq-dump --gzip --split-files -O fastqs_path data/SRR7722940
fastq-dump --gzip --split-files -O fastqs_path data/SRR7722941
fastq-dump --gzip --split-files -O fastqs_path data/SRR7722942
for i in `cat SRR_Acc_List.txt`
do
mv fastqs_path/${i}_1.fastq.gz fastqs_path/${i}_S1_L001_I1_001.fastq.gz
mv fastqs_path/${i}_2.fastq.gz fastqs_path/${i}_S1_L001_R1_001.fastq.gz
mv fastqs_path/${i}_3.fastq.gz fastqs_path/${i}_S1_L001_R2_001.fastq.gz
done
3.从FASTQ到count矩阵¶
cellranger计数管道将FASTQ文件中的测序结果与参考转录组进行比对,并生成一个.cloupe
文件,用于在Loupe Browser中进行可视化和分析,同时还生成了一些与其他公开工具兼容的输出,用于进一步分析。
接下来,需要一个参考基因组数据。从FASTQ文件的下载页面可以看到这些是人类细胞。在10x Genomics支持网站上有几个预构建的人类参考转录组包。下载最新的包并解压缩它。
#人
curl -O https://cf.10xgenomics.com/supp/cell-exp/refdata-cellranger-GRCh38-3.0.0.tar.gz
#小鼠
curl -O https://cf.10xgenomics.com/supp/cell-exp/refdata-gex-mm10-2020-A.tar.gz
tar zxvf refdata-gex-GRCh38-2020-A.tar.gz
curl -O https://cf.10xgenomics.com/supp/cell-exp/refdata-gex-mm10-2020-A.tar.gz
tar -xvf pbmc_1k_v3_fastqs.tar
pbmc_1k_v3_fastqs/
pbmc_1k_v3_fastqs/pbmc_1k_v3_S1_L001_R2_001.fastq.gz
pbmc_1k_v3_fastqs/pbmc_1k_v3_S1_L002_I1_001.fastq.gz
pbmc_1k_v3_fastqs/pbmc_1k_v3_S1_L001_R1_001.fastq.gz
pbmc_1k_v3_fastqs/pbmc_1k_v3_S1_L002_R1_001.fastq.gz
pbmc_1k_v3_fastqs/pbmc_1k_v3_S1_L002_R2_001.fastq.gz
pbmc_1k_v3_fastqs/pbmc_1k_v3_S1_L001_I1_001.fastq.gz
cellranger count
了。
cellranger count --id=run_count_1kpbmcs \
--fastqs= /mnt/f/Linux/run_cellranger_count/pbmc_1k_v3_fastqs \
--sample=pbmc_1k_v3 \
--transcriptome= /mnt/f/Linux/genomeAnno/hsa/refdata-gex-GRCh38-2020-A \
--expect-cells=1000 \
--localcores=6 \
--localmem=8 \
--nosecondary
cellranger count --id=cellranger_count \
--transcriptome=/mnt/f/Linux/genomeAnno/hsa/refdata-gex-GRCh38-2020-A \
--fastqs=fastq_path/H35KCBCXY/test_sample \
--sample=test_sample \
--expect-cells=1000 \
--localcores=6 \
--localmem=8 \
--nosecondary
cellranger count --id=cellranger_count \
--transcriptome=/mnt/f/Linux/genomeAnno/hsa/refdata-gex-GRCh38-2020-A \
--fastqs=/mnt/f/Linux/sradata/sra_data/fastqs_path \
--sample=SRR7722937 \
--expect-cells=1000 \
--localcores=6 \
--localmem=8 \
--nosecondary
4.STAR/STARsolo:单细胞FastQ文件比对¶
STAR是一款用于高通量转录组比对的软件,是Encode计划和TCGA计划使用的官方比对软件,基于STAR分析转录组的流程,通常称为STAR Method, STARsolo是STAR专门用于单细胞转录组定量的模块,原理类似于cellranger count, STAR官方宣称速度比cellranger快10倍。 安装:https://anaconda.org/bioconda/star 在conda环境中安装
#要使用 conda 安装此包,请运行以下命令之一
conda install -c bioconda star
conda install -c bioconda/label/cf201901 star
STAR --runMode genomeGenerate \
--genomeDir star_hg38_index \
--genomeFastaFiles refdata-gexGRCh38-2020-A/fasta/genome.fa \
--runThreadN 8
--sjdbGTFfile refdata-gex-GRCh38-2020-A/genes/genes.gtf
基因组fasta下载地址: https://asia.ensembl.org/Homo_sapiens/Info/Index 基因注释gtf下载地址: https://asia.ensembl.org/Homo_sapiens/Info/Index 以SRR7722939为例 • --genomeDir: STAR索引所在目录 • --readFilesCommand zcat:使用zcat读取gz压缩格式文件 • --readFilesIn:读入fastq文件的路径,首先是R2,然后是R1 • --soloOutFileNames:输出文件夹名 • --soloType CB_UMI_Simple:文库类型 • --soloCBwhitelist:所有使用的barcode序列 • --runThreadN:服务器使用的线程数
STARsolo的结果与cellranger count的结果相同, cellranger count也是调用STAR进行比对。STARsolo的运行速度比cellranger count快10倍。 • Barcodes.stats:对Barcode代码的统计。 • Features.stats:特征统计。 • Filtered:过滤后的表达矩阵,过滤方法与cellranger count 相同(用于下游分析,最主要)。 • Raw:未过滤的数据。 • Summary.csv:比对汇总。 • UMIperCellSorted.txt: UMI统计信息。
5.Salmon/Alevin:对FastQ文件进行基因定量¶
Salmon是一款不基于比对的定量软件,可用于转录组数据的快速定量,准确性高,其中Alevin功能 专门用于单细胞转录组定量【PMID: 28263959】。
安装
#要使用 conda 安装此包,请运行以下命令之一
conda install -c bioconda salmon
conda install -c bioconda/label/cf201901 salmon
time salmon index -t Homo_sapiens.GRCh38.cdna.all.fa.gz -i hg38_salmon