论文部分内容阅读
DNA序列的插入缺失标记（insertion and deletion），又被统称为indel，是指基因组上发生一定长度的核苷酸序列的增加或删减，是一种重要的遗传变异类型。由于indel在染色体上的分布广泛且密集、遗传模式稳定、多态性强，在法医学鉴定，遗传病诊断，肿瘤检测等领域有着广阔的应用前景。同时，下一代测序技术（Next Generation Sequencing，NGS）的迅猛发展也推动着indel检测率的提高。然而，全基因组或外显子组测序，对于大样本、小靶点的检测序是成本高昂的测序策略。基于多重PCR的靶向DNA富集是一种经济、快速、准确的测序建库策略，在大样本群体的indel检测中有着极大的潜力。本实验室在2017年开发了使用钝发夹结构优化的PCR引物来提升多重PCR靶向测序建库效率，该方法能使靶向测序深度均一化同时降低引物二聚体的发生。在数据分析工具上，目前仍未出现针对多重PCR靶向测序数据开发的indel识别软件或分析流程。在NGS数据中识别遗传变异位点的过程被称为calling。Indel calling相比于SNP calling的难度要大，因为插入/缺失的存在，本身就干扰比对排序过程，这种干扰会影响Indel calling本身的准确性。因此，当前主流indel识别软件在靶向测序数据集中识别indel性能的评估，以及自动化indel识别流程的建立，亟待研究。
目的：本文使用两组人类样本的多重PCR靶向重测序数据集（每组数据集的样本数>1000），旨在评估目前主流的indel识别软件在测试数据集中识别indel变异的性能，并构建且优化indel自动化分析流程。
方法：第一，本实验室搭建了高性能计算机（High Performance Computer，HPC）用于分析和存储海量的高通量测序数据。用于搭建HPC的硬件主要包括Intel公司的Xeon E5-2620处理器，和Supermicro公司的X10DRH主板；测试的软件分为比对软件和indel识别软件，比对软件包括BWA_MEM和Minimap2，indel识别软件包括Pindel、VarScan2、SAMtools、Sentieon Haplotyper（SH）、Sentieon Genotyper（SG）、以及Genome Analysis ToolKit（GATK）工具包中的UnifiedGenotyper（GATK_UG）和HaplotypeCaller（GATK_HC）。第二，本研究使用了两组人类样本的多重PCR靶向测序数据集，总数据量约151Gb，样本数分别为1248和2496个、样本平均reads数分别为301818和288332、重测序位点平均深度达2800以上，原始数据集经过切除接头和质量控制预处理数据量。第三，软件评估标准的纳入。主要包括运行时间、准确度、精准度、召回率与F1-measure值等指标，同时技术性重复样本被用于评估indel识别软件的一致性。第四，使用整合基因组浏览器（Integrative Genomics Viewer，IGV）手动审查样本基因型。为了客观地评价变异识别软件的分类性能，通过利用IGV软件对测序比对文件的可视化，来手动审查每个样本的indel基因型。第五，模拟变异数据集的构造与使用。使用自定义的Python程序来生成模拟indel数据集，该模拟数据集被视为相对真集，用来统计各indel识别软件对真实数据的真、假阳性率。第六，GATK关键参数调试。通过打开/关闭GATK3.X的downsampling参数，探究了该参数对识别indel基因型的影响。第七，分析流程的搭建与优化。使用Python语言将用于分析测序数据的各个组件进行封装，构建完整的分析流程。分析流程主要分为测序数据分析和变异数据报告系统两大部分。测序数据分析部分包括了测序数据预处理（如切接头、质控）、reads比对、过滤和indel识别；变异数据报告系统包括解析记录变异信息的Variant Calling Format(VCF)文件、统计测序深度、基因型数据标准化、以及计算最小等位基因频率（Minimum Allele Frequency，MAF）等。
结果：（1）运行时间评估发现，2款比对软件和9款indel识别软件组成的18个不同的pipeline的运行时间差异较大。就比对软件而言，Minimap2比BWA_MEM在reads比对过程中消耗的时间更少，使用Minimap2能为整个NGS数据分析在序列比对步骤中节省最多约70.85%的时间；就indel识别软件而言，运行时间最长是Pindel，最短的是GATK_UG3.7。（2）利用技术性重复样本评估各indel识别软件的一致性，结果发现不同的indel识别软件存在较大差异。其中，一致性表现最好的依次是GATK_HC3.3、GATK_HC3.7、Sentieon-Haplotyper和GATK_HC4.0。同时，还发现GATK_HC3.X的默认变异等位基因频率（Variant Allele Frequency，VAF）阈值小于0.2，会造成判定杂合子基因型的准确度降低。（3）不同软件在准确度、精准度与召回率评估上，存在较大差异。准确度最高的前三个pipeline分别是Minimap2-SH、Minimap2-HC3.7和BWA_MEM-SH，精准度最高的前三个pipeline分别是Minimap2-UG3.3、BWA_MEM-UG3.3和Minimap2-UG3.7，召回率表现最好的前三个pipeline分别是Minimap2-HC3.7、Minimap2-SH和Minimap2-HC3.3。综合性能较好的indel识别软件是HC3.3、HC3.7和Sentieon-Haplotyper（4）对比研究发现，GATK关键参数downsampling严重影响运行时间与基因型结果。关闭downsampling参数后，使用GATK_HC的pipeline运行时间显著提升，BWA_MEM-HC3.3、Minimap2-HC3.3、BWA_MEM-HC3.7和Minimap2-HC3.7的运行时间分别提升了96.0%、119.4%、236.2%和277.7%。同时还发现，关闭downsampling参数可更精准的提供VCF中记录的深度值，亦有部分样本的基因型随之改善。（5）优化并建立了自动化indel检测分析流程。使用Minimap2作为比对软件，GATK_HC3.3和GATK_HC3.7作为变异识别软件建立了自动化indel检测分析流程，并对downsampling参数及GATK_HC3.X的默认VAF阈值进行了优化，使用Python语言对各软件的命令行进行封装，同时开发了变异信息报告系统用于解析pipeline生成的VCF，从中提取等位基因、深度值、计算MAF以供信息解读和存储，最后将该分析流程和报告系统上传至Github。
结论：本研究首次系统性评估了多重PCR靶向测序数据中indel识别软件的性能，包括运行时间、准确度、精准度和召回率等指标。最终发现Minimap2可以在reads比对过程中比BWA_MEM节省最多70.85%的运行时间，同时GATK_HC3.3、GATK_HC3.7和Sentieon-Haplotyper在基于PCR建库的靶向测序数据中识别indel的表现优于其他变异识别软件。此外，还对GATK_HC3.X中downsampling参数和默认的VAF阈值进行了优化，使GATK的应用更适合处理靶向测序数据。基于这些结果，本研究首次开发了经过优化的基于多重PCR靶向测序数据中indel的分析流程。该分析流程主要优点包括分析速度快，中间操作简单，错误回溯方便，数据报告可读性强等优点，适用于多种不同物种的大规模群体并行靶向测序研究。该流程的开发和优化能够极大推动使用靶向测序数据检测indel这一新型遗传变异标记在分子遗传育种、孟德尔遗传疾病诊断、基因组学功能研究等领域的更深入应用。
目的：本文使用两组人类样本的多重PCR靶向重测序数据集（每组数据集的样本数>1000），旨在评估目前主流的indel识别软件在测试数据集中识别indel变异的性能，并构建且优化indel自动化分析流程。
方法：第一，本实验室搭建了高性能计算机（High Performance Computer，HPC）用于分析和存储海量的高通量测序数据。用于搭建HPC的硬件主要包括Intel公司的Xeon E5-2620处理器，和Supermicro公司的X10DRH主板；测试的软件分为比对软件和indel识别软件，比对软件包括BWA_MEM和Minimap2，indel识别软件包括Pindel、VarScan2、SAMtools、Sentieon Haplotyper（SH）、Sentieon Genotyper（SG）、以及Genome Analysis ToolKit（GATK）工具包中的UnifiedGenotyper（GATK_UG）和HaplotypeCaller（GATK_HC）。第二，本研究使用了两组人类样本的多重PCR靶向测序数据集，总数据量约151Gb，样本数分别为1248和2496个、样本平均reads数分别为301818和288332、重测序位点平均深度达2800以上，原始数据集经过切除接头和质量控制预处理数据量。第三，软件评估标准的纳入。主要包括运行时间、准确度、精准度、召回率与F1-measure值等指标，同时技术性重复样本被用于评估indel识别软件的一致性。第四，使用整合基因组浏览器（Integrative Genomics Viewer，IGV）手动审查样本基因型。为了客观地评价变异识别软件的分类性能，通过利用IGV软件对测序比对文件的可视化，来手动审查每个样本的indel基因型。第五，模拟变异数据集的构造与使用。使用自定义的Python程序来生成模拟indel数据集，该模拟数据集被视为相对真集，用来统计各indel识别软件对真实数据的真、假阳性率。第六，GATK关键参数调试。通过打开/关闭GATK3.X的downsampling参数，探究了该参数对识别indel基因型的影响。第七，分析流程的搭建与优化。使用Python语言将用于分析测序数据的各个组件进行封装，构建完整的分析流程。分析流程主要分为测序数据分析和变异数据报告系统两大部分。测序数据分析部分包括了测序数据预处理（如切接头、质控）、reads比对、过滤和indel识别；变异数据报告系统包括解析记录变异信息的Variant Calling Format(VCF)文件、统计测序深度、基因型数据标准化、以及计算最小等位基因频率（Minimum Allele Frequency，MAF）等。
结果：（1）运行时间评估发现，2款比对软件和9款indel识别软件组成的18个不同的pipeline的运行时间差异较大。就比对软件而言，Minimap2比BWA_MEM在reads比对过程中消耗的时间更少，使用Minimap2能为整个NGS数据分析在序列比对步骤中节省最多约70.85%的时间；就indel识别软件而言，运行时间最长是Pindel，最短的是GATK_UG3.7。（2）利用技术性重复样本评估各indel识别软件的一致性，结果发现不同的indel识别软件存在较大差异。其中，一致性表现最好的依次是GATK_HC3.3、GATK_HC3.7、Sentieon-Haplotyper和GATK_HC4.0。同时，还发现GATK_HC3.X的默认变异等位基因频率（Variant Allele Frequency，VAF）阈值小于0.2，会造成判定杂合子基因型的准确度降低。（3）不同软件在准确度、精准度与召回率评估上，存在较大差异。准确度最高的前三个pipeline分别是Minimap2-SH、Minimap2-HC3.7和BWA_MEM-SH，精准度最高的前三个pipeline分别是Minimap2-UG3.3、BWA_MEM-UG3.3和Minimap2-UG3.7，召回率表现最好的前三个pipeline分别是Minimap2-HC3.7、Minimap2-SH和Minimap2-HC3.3。综合性能较好的indel识别软件是HC3.3、HC3.7和Sentieon-Haplotyper（4）对比研究发现，GATK关键参数downsampling严重影响运行时间与基因型结果。关闭downsampling参数后，使用GATK_HC的pipeline运行时间显著提升，BWA_MEM-HC3.3、Minimap2-HC3.3、BWA_MEM-HC3.7和Minimap2-HC3.7的运行时间分别提升了96.0%、119.4%、236.2%和277.7%。同时还发现，关闭downsampling参数可更精准的提供VCF中记录的深度值，亦有部分样本的基因型随之改善。（5）优化并建立了自动化indel检测分析流程。使用Minimap2作为比对软件，GATK_HC3.3和GATK_HC3.7作为变异识别软件建立了自动化indel检测分析流程，并对downsampling参数及GATK_HC3.X的默认VAF阈值进行了优化，使用Python语言对各软件的命令行进行封装，同时开发了变异信息报告系统用于解析pipeline生成的VCF，从中提取等位基因、深度值、计算MAF以供信息解读和存储，最后将该分析流程和报告系统上传至Github。
结论：本研究首次系统性评估了多重PCR靶向测序数据中indel识别软件的性能，包括运行时间、准确度、精准度和召回率等指标。最终发现Minimap2可以在reads比对过程中比BWA_MEM节省最多70.85%的运行时间，同时GATK_HC3.3、GATK_HC3.7和Sentieon-Haplotyper在基于PCR建库的靶向测序数据中识别indel的表现优于其他变异识别软件。此外，还对GATK_HC3.X中downsampling参数和默认的VAF阈值进行了优化，使GATK的应用更适合处理靶向测序数据。基于这些结果，本研究首次开发了经过优化的基于多重PCR靶向测序数据中indel的分析流程。该分析流程主要优点包括分析速度快，中间操作简单，错误回溯方便，数据报告可读性强等优点，适用于多种不同物种的大规模群体并行靶向测序研究。该流程的开发和优化能够极大推动使用靶向测序数据检测indel这一新型遗传变异标记在分子遗传育种、孟德尔遗传疾病诊断、基因组学功能研究等领域的更深入应用。