DNA测序数据分析的分步指南

0

Javier Quilez Oliete博士,一个经验丰富的 生物信息学顾问 在Kolabtree上,提供了一份关于DNA测序数据分析的全面指南,包括用于读取数据的工具和软件。 

简介

脱氧核糖核酸(DNA)是携带大部分遗传信息的分子。 生物体的.在某些类型的病毒中,遗传信息是由核糖核酸(RNA)携带的)。 核苷酸(通常用字母A、C、G或T表示)是DNA分子的基本单位。从概念上讲。 DNA测序 是读取组成DNA分子的核苷酸的过程(例如,"GCAAACCAAT "是一个10个核苷酸的DNA字符串)。目前的测序技术产生了数百万个这样的DNA读数 在合理的时间内,以相对较低的成本进行测序。作为参考,人类基因组测序的成本--基因组是一个生物体内完整的DNA分子集--已经下降到 $100屏障 而且可以在几天内完成。这与第一项举措形成对比,即对 人类基因组该项目在十年内完成,耗资约1TP2-2.7亿美元。

这种高通量和低成本的DNA测序能力使越来越多的基于测序的方法和应用得到发展。例如,对疾病和健康个体的全基因组或其蛋白质编码区进行测序(两种方法分别称为全基因组和外显子组测序)可以提示致病的DNA改变。此外,对DNA转录的RNA进行测序--一种被称为RNA测序的技术--被用来量化基因的活性以及在不同条件下的变化(如未治疗与治疗)。另一方面,染色体构象捕获测序方法检测附近DNA分子之间的相互作用,从而帮助确定细胞内染色体的空间分布。

这些和其他DNA测序应用的共同点是产生了数千兆字节的数据集,包括数百万读数序列。因此,对高通量测序(HTS)实验的理解需要大量的数据分析能力。幸运的是,大多数HTS数据类型都有专门的计算和统计工具以及相对标准的分析工作流程。虽然一些(初始)分析步骤对大多数测序数据类型来说是共同的,但更多的下游分析将取决于数据的种类和/或分析的最终目标。下面我提供一个关于HTS数据分析的基本步骤的入门知识,并且我提到了流行的工具。 

下面的一些章节主要是对短读测序技术产生的数据的分析(主要是 伊利姆纳),因为这些技术在历史上一直主导着HTS市场。然而,产生较长读数的新技术(如 牛津纳米孔技术公司, 医学博士)正迅速得到发展。由于长读测序有一些特殊性(如较高的错误率),目前正在为分析这类数据开发专门的工具。 

原始读数的质量控制(QC

渴望的分析师将从FASTQ文件开始分析;而 FASTQ格式 长期以来一直是存储短读测序数据的标准。从本质上讲,FASTQ文件包含核苷酸序列和每个碱基的 对数百万条读数的调用质量。虽然文件的大小取决于实际的读数数量,但FASTQ文件通常是大的(在兆字节和千兆字节的数量级)和压缩的。值得注意的是,大多数使用FASTQ文件作为输入的工具可以处理压缩格式的文件,因此,为了节省磁盘空间,建议不要解压它们。作为惯例,我在这里将FASTQ文件等同于一个测序样本。

快速QC 可能是对原始读数进行质量控制的最流行的工具。它可以通过视觉界面或编程方式运行。对于不习惯使用命令行环境的用户来说,第一种选择可能更方便,而后者则提供了无可比拟的可扩展性和可重复性(想想看,为几十个文件手动运行该工具是多么的乏味和容易出错)。不管怎么说,FastQC的主要输出是一个 HTML文件 报告关于特定样本原始测序读数总体质量的关键汇总统计数据。逐一检查几十份FastQC报告是很乏味的,而且会使不同样品的比较变得复杂。因此,你可能想使用 多重控制它将来自FastQC的HTML报告(以及来自下游使用的其他工具,如适配体修剪、对齐)汇总到一个报告中。.

多重控制

质量控制信息的目的是让用户判断样品是否具有良好的质量,因此可以用于后续步骤,或者需要将其丢弃。遗憾的是,目前还没有一个基于FastQC指标的共识阈值来划分样品的质量好坏。我使用的方法如下。我希望所有经过相同程序(如DNA提取、文库制备)的样本都有类似的质量统计,并有大多数的 "通过 "标志。如果有些样品的质量低于平均水平,我仍然会在下游分析中使用它们,并记住这一点。另一方面,如果实验中的所有样品在多个指标中都系统地得到 "警告 "或 "失败 "标记(见 这个例子),我怀疑是实验中出了问题(如DNA质量不好,文库制备等),我建议重复实验。

读取修剪

原始读数的质量控制有助于识别有问题的样本,但它并不能提高读数的实际质量。为此,我们需要对读数进行修剪,以去除技术序列和低质量末端。

技术序列是实验过程中的遗留物(如测序适配器)。如果这样的序列与读数的真实序列相邻,比对(见下文)可能会将读数映射到基因组中的错误位置或降低给定比对的可信度。除了技术序列外,如果生物来源的序列在读数中高度存在,我们也可能要删除这些序列。例如,不理想的DNA制备程序可能会在样品中留下高比例的DNA转化的核糖体RNA(rRNA)。除非这种类型的核酸是测序实验的目标,否则保留来自rRNA的读数只会增加下游步骤的计算负担,并可能混淆结果。值得注意的是,如果技术序列、rRNA或其他污染物的水平非常高,这可能已经被质控部门强调过了,你可能要丢弃整个测序样品。

在短读测序中,DNA序列是一次确定一个核苷酸(技术上讲,每个测序周期确定一个核苷酸)。换句话说,测序周期的数量决定了读长度。HTS测序方法的一个已知问题是,随着测序周期的累积,确定核苷酸的准确性会下降。这反映在每碱基调用质量的整体下降,特别是在读数的末端。正如技术序列所发生的那样,试图对含有低质量末端的读数进行比对会导致错位或映射质量差。

为了去除技术/污染序列和低质量的末端,阅读修剪工具如 特里姆摩斯削减 存在并被广泛使用。从本质上讲,这些工具将去除技术序列(内部可用和/或由用户提供),并根据质量修剪读数,同时最大限度地增加读数长度。修剪后留下的太短的读数会被丢弃(过短的读数,如<36个核苷酸,会使比对步骤复杂化,因为这些读数可能会映射到基因组的多个位点)。你可能想看看修剪后幸存下来的读数的百分比,因为高比率的被丢弃的读数可能是数据质量不好的标志。 

最后,我通常在修剪过的读数上重新运行FastQC,以检查这个步骤是否有效,是否系统地改善了QC指标。

统一口径

除了例外情况(如 从头组装),对于大多数HTS数据类型和应用来说,比对(也被称为映射)通常是下一个步骤。读数比对包括确定读数序列在基因组中的位置(通常表示为染色体:起始端)。因此,在这一步,我们需要使用一个参考序列来对准/映射读数。

参考序列的选择将由多种因素决定。其一,被测序的DNA来自哪个物种。虽然具有高质量参考序列的物种数量正在增加,但对于一些研究较少的生物体来说,情况可能仍然不是这样的。在这些情况下,你可能想把读数与一个进化上接近的、有参考基因组的物种对齐。例如,由于没有土狼基因组的参考序列,我们可以使用与之密切相关的狗的基因组来进行读数比对。同样,我们可能仍然希望将我们的读数与存在更高质量参考序列的密切相关物种进行比对。例如,虽然长臂猿的基因组已经被 出版的在这种情况下,使用人类参考序列进行比对可能是有益的。

另一个需要考虑的因素是参考序列组合的版本,因为随着序列的更新和改进,新的版本也会发布。重要的是,一个给定的排列组合的坐标可以在不同的版本之间变化。例如,人类基因组的多个版本可以在 UCSC基因组浏览器.在任何物种中,我都强烈赞成在最新的汇编版本完全发布后迁移到该版本。这在过渡期间可能会造成一些麻烦,因为已经存在的结果将是相对于旧版本的,但从长远来看,这是有回报的。

此外,测序数据的类型也很重要。从DNA-seq、ChIP-seq或Hi-C协议产生的读数将与基因组参考序列对齐。另一方面,由于从DNA转录的RNA被进一步加工成mRNA(即去除内含子),许多RNA-seq读数将无法与基因组参考序列对齐。相反,当使用基因组序列作为参考时,我们需要将它们与转录组参考序列对齐,或者使用分裂意识对齐器(见下文)。与此相关的是对参考序列注释来源的选择,也就是具有基因、转录本、中心点等坐标的数据库。我通常使用 GENCODE注释 因为它结合了全面的基因注释和转录物序列。

已经开发了一长串的短读序列比对工具(见短读序列比对部分)。 这里).对它们的评论超出了本文的范围(关于这些工具背后的算法的细节可以找到 这里).根据我的经验,最受欢迎的是 Bowtie2, BWA, HISAT2, Minimap2, 星级帽子.我的建议是,你在选择对准器时要考虑一些关键因素,如HTS数据的类型。 我们需要的是一个能使我们的产品和应用与社区的接受程度、文件的质量和用户的数量相匹配的软件。例如,人们需要像STAR或Bowtie2这样的对准器,在将RNA-seq映射到基因组时能意识到外显子-外显子连接。 

大多数映射器的共同点是需要在实际比对之前为用作参考的序列编制索引。这个步骤可能很耗时,但对每个参考序列只需要做一次。大多数映射器将在SAM/BAM文件中存储比对,这些文件遵循的是 SAM/BAM格式 BAM文件是SAM文件的二进制版本)。在测序数据的分析中,比对是最需要计算和耗时的步骤之一,而且SAM/BAM文件很重(以千兆字节计)。因此,重要的是要确保你有必要的资源(见下面的最后一节),以便在合理的时间内运行比对并存储结果。同样,由于BAM文件的大小和二进制格式,避免用文本编辑器打开它们;而是使用Unix命令或专用工具,如 萨姆工具.

从排列组合来看

我想说的是,在对准之后没有一个明确的共同步骤,因为在这一点上,每个HTS数据类型和应用都可能有所不同。 

DNA-seq数据的一个常见的下游分析是变体调用,即识别基因组中相对于基因组参考和个体之间的不同位置。这种应用的一个流行的分析框架是 GATK 单核苷酸多态性(SNP)或小的插入/删除(indels)的检测。图2).由较大块的DNA组成的变体(也被称为结构变体)需要专门的调用方法(见 本条 以进行全面的比较)。)与对准器一样,我建议选择正确的工具,考虑关键因素,如变异的种类(SNP,indel或结构变异),社区的接受程度,文件的质量和用户的数量。

RNA-seq最常见的应用可能是基因表达定量。历史上,读数需要与参考序列对齐,然后将与特定基因或转录物对齐的读数作为代理来量化其表达水平。这种对齐+量化的方法是由以下工具完成的 袖扣, ǞǞǞ特征计数.然而,这种方法已被越来越多的新方法所超越,这些方法在软件中实现,如 卡利斯多三文鱼.从概念上讲,使用这种工具,读数的全部序列不需要与参考序列对齐。相反,我们只需要对准足够多的核苷酸,就可以确信一个读数来自于一个特定的转录本。简单地说,对齐+量化的方法被简化为一个步骤。这种方法被称为 "伪映射",大大增加了基因表达定量的速度。另一方面,请记住,伪映射不适合于需要完全对齐的应用(如从RNA-seq数据中调用变体)。

基于测序的应用中下游分析步骤和所需工具的差异的另一个例子是ChIP-seq。用这种技术产生的读数将被用于峰值调用,这包括检测基因组中读数明显过剩的区域,表明目标蛋白被结合的地方。有几个峰值调用器存在,并且 本出版物 调查它们。作为最后一个例子,我将提到Hi-C数据,其中的比对被用作确定相互作用矩阵的工具的输入,并从中确定基因组的三维特征。对本文范围之外的所有基于测序的检测方法进行评论(相对完整的清单见 本条).

在你开始之前...

本文的其余部分涉及到在分析HTS数据时可能没有严格考虑的步骤,而且在很大程度上被忽视的方面。与此相反,我认为,你有必要思考一下在 "HTS "中提出的问题。 表1 在你开始分析HTS数据(或任何种类的数据)之前,我曾就这些主题写过文章 这里这里.

表1

考虑一下吧 拟议的行动
你有分析所需的所有样品信息吗? 系统地收集实验的元数据
你能明确地识别你的样品吗? 建立一个系统,为每个样品分配一个独特的标识符。
数据和结果将在哪里? 数据的结构化和层次化组织
你是否能够无缝地处理多个样品? 代码的可扩展性、并行化、自动配置和模块化
你或其他人能够重现这一结果吗? 记录你的代码和程序!

 

如上所述,HTS的原始数据和分析过程中产生的一些文件是以千兆字节为单位的,因此,一个包括几十个样品的项目需要数千兆字节的存储空间,这并不罕见。此外,HTS数据分析中的一些步骤是计算密集型的(如排列)。然而,分析HTS数据所需的存储和计算基础设施是一个重要的考虑因素,它经常被忽视或没有被讨论。作为一个例子,作为最近分析的一部分,我们回顾了几十篇发表的论文,进行全表型关联分析(PheWAS)。现代PheWAS分析了100-1000个遗传变异和表型,这导致了重要的数据存储和计算能力。然而,在我们审查的论文中,几乎没有一篇对PheWAS分析所需的基础设施进行评论。毫不奇怪,我的建议是,你们要预先计划你们将面临的存储和计算要求,并与社区分享。

在分析DNA测序数据方面需要帮助?请联系 自由职业者生物信息学家基因组学专家 关于Kolabtree。 


Kolabtree帮助全球企业按需雇佣专家。我们的自由职业者已经帮助企业发表研究论文,开发产品,分析数据,以及更多。只需一分钟就可以告诉我们你需要做什么,并免费获得专家的报价。


分享。

关于作者

Javier Quilez Oliete博士是一名生物信息学家,在应用于基因组学的数据科学方面拥有10年以上的经验。他在管理和分析用SNP/甲基化微阵列和高通量测序(全基因组、靶向、RNA和ChIP测序、Hi-C)产生的大规模基因组数据集方面经验丰富。他还从事过全基因组关联研究(GWAS)、变异调用(SNPs、indels和较大的结构变异)、差异表达分析和ChIP测序峰值调用。他用Python、Bash、R、Perl编程,并使用IPython笔记本、GitHub(@jaquol)、Markdown和SQL数据库。

发表回复

值得信赖的自由职业者专家,随时为您的项目提供帮助


世界上最大的科学家自由职业平台  

不,谢谢,我现在不打算雇用。