全基因组测序临床生物信息学流程

Partners HealthCare 个性化医疗全基因组测序的生物信息学工作流程

摘要

精准医疗的有效实施需要充分了解每位患者的遗传构成，以更好地治疗其表现出的症状或减轻疾病的发作。理想情况下，这包括每个个体完整基因组的序列信息。在Partners HealthCare 个性化医疗，我们已开发出一种用于全基因组测序（WGS）的临床流程，适用于健康个体以及患病个体。在本文中，我们将介绍我们的生物信息学策略，以高效地处理和交付基因组数据给遗传学家进行临床解读。我们描述了从FASTQ到最终用于临床审查的变异列表的数据处理过程。我们还将讨论验证该工作流程的方法以及运行全基因组测序的成本影响。

关键词 ：临床测序；全基因组测序；下一代测序；生物信息学；验证；精准医疗

1. 引言

精准医疗正日益成为医学研究的重点 [1]。为了实现个性化临床诊疗所需的分辨率，人们越来越关注患者基因组的更高分辨率。下一代测序（NGS）提供了一种经济高效的方法，可在碱基对分辨率下对已知的疾病基因进行靶向测序[2]。此外，外显子组测序的出现使得导致孟德尔疾病的新基因得以快速发现。尽管基因面板和外显子组测序在将基因组结果快速、低成本地反馈给患者方面已证明有效，但这些技术受限于我们当前对外显子组的认识，而这种认识会随时间变化。此外，靶向捕获的使用可能引入数据偏差，包括PCR重复、覆盖深度差异以及难以扩增的目标区域的检测失败 [3]。

测序成本、数据处理与维护以及数据分析的复杂性等实际问题是实验室在考虑开展新的下一代测序（NGS）项目时的重要考量因素。由于数据量庞大，这些问题在全基因组测序（WGS）中更为突出，长期以来一直是临床实验室采用全基因组测序（WGS）的进入壁垒。尽管全基因组测序（WGS）能够检测整个基因组，但临床解读通常仍集中于仅覆盖基因组的3%（即，外显子组数据、药物基因组学风险变异以及与复杂疾病风险相关的单核苷酸变异）[4–7]。因此，由于计算负荷加重、需分析的变异数量增加以及数据存档规模更大，全基因组测序服务在临床环境中的应用可能因成本上升和周转时间延长而被忽视。然而，测序和存储成本的持续下降现在使实验室能够考虑开展基因组测序。全基因组测序，特别是无PCR全基因组测序，还可减少因编码序列或靶向区域变化、新基因发现或新参考基因组发布而需要重复测序的频率。在启动全基因组测序项目时，必须权衡成本、周转时间、准确性与完整性。本文中，我们描述了在临床环境中支持全基因组测序生物信息学工作流程所采用的方法及面临的挑战。

2. 结果

2.1. 生物信息学验证

我们的流程表现稳健，不同的输入点以及同一样本的不同运行均返回完全相同的变异。利用 HapMap样本NA12878和先前的Sanger确认区域，我们确定了深度质量分数(QD) ≥ 4和费舍尔链偏差(FS) ≤ 30的变异阈值，以达到敏感性和特异性的最佳平衡。该样本已被我们实验室以及其他研究团队充分表征[8,9]。在总共425个已确认的变异中，基因组测序检测到了全部425个变异，敏感性为100%（95% CI：SNVs为99.1%–100%，插入缺失突变为79.6%–100%；表1a）。此外，有四个可能的参考序列错误位点（这些位置仅发现过替代等位基因）被基因组测序正确地鉴定为替代等位基因的纯合子，其中三个在此前的正交检测中存在错误的基因型判定。除这些真阳性变异外，还检出了21个假阳性(FP)变异，包括20个替换和一个插入缺失突变（表1b）。基于前述最优QD和FS阈值进行变异过滤后，仅剩一个FP变异。

表1. 全基因组测序的敏感性和特异性。利用HapMap样本NA12878，我们先前通过桑格测序在~700 kb序列范围内的195个基因中识别并确认了425个变异。

(a) 特异性
| 变异类型 | 假阳性（阈值前） | 假阳性（阈值后） |
| :— | :— | :— |
| SNVs | 20 | 1 |
| 插入缺失突变 | 1 | 0 |

(b) 敏感性
| 变异类型 | # FN | 敏感性 | 95% 置信区间 |
| :— | :— | :— | :— |
| SNVs | 410 | 0 | 100% 99.1%–100% |
| 插入缺失突变 | 15 | 0 | 100% 79.6%–100% |

(c) 与1000基因组数据的一致性
| 变异类型 | 1K基因组变异 | 存在于 NGS检测 | %存在于 NGS检测 | 存在于NGS检测中与匹配的基因型 | %存在于NGS中检测与匹配的基因型 |
| :— | :— | :— | :— | :— | :— |
| SNVs | 2,762,933 | 2,735,592 | 99.01% | 2,730,826 | 98.84% |
| 插入缺失突变 | 3,27,474 | 299,300 | 91.39% | 285,401 | 87.15% |
| 总计 | 3,090,407 | 3,034,892 | 98.20% | 3,016,227 | 97.60% |

我们将通过本流程检测到的NA12878变异与在相似覆盖度下通过1000基因组计划检测到的变异进行了全基因组比较，结果一致（表1c）。在我们流程检测出的270万个单核苷酸变异和28.5万个插入缺失突变中，发现其中98.8%的单核苷酸变异也在1000基因组数据集中被检出。此外，总共检出的变异中有97.6%的基因型一致。对NA12878也进行了注释和变异过滤。随机选取了50多个变异进行人工检查，以确保其得到了正确的注释和过滤。

2.2. 已知的低覆盖度区域

我们从Illumina的CLIA实验室接收到的所有基因组均具有至少30倍覆盖度，平均覆盖度为43倍。然而，这种覆盖度在基因组的编码序列中存在差异，并影响基因组的临床相关区域和临床未知区域。基因水平上可检测碱基百分比的列表见补充表S1。总计，在ClinVar中被临床实验官认定为具有至少五个致病性或可能致病性变异的1381个基因中，有94个基因的编码区覆盖度低于<90%。覆盖度指标最差的20个基因中包含许多具有高患病率和临床相关性的基因（表2）。在我们的临床流程中，针对基于指征基因列表的驱动分析，我们会报告所分析基因的覆盖度，并突出显示存在覆盖度问题的基因。这为判断潜在的假阴性结果提供了有用的指导，包括那些由于难以用下一代测序技术进行测序的区域所导致的结果。

表2. 临床相关性差的前20个覆盖不佳基因。临床相关性定义为在提交实验室或工作组报告的基因中，ClinVar数据库中至少有五个致病性或可能致病性变异。

Gene	#临床显著变异	%可检测	疾病	疾病患病率
STRC	8	20	感音神经性耳聋	常见
ADAMTSL2	5	32	胶样物理发育不良	Rare
CYP21A2	13	44	先天性肾上腺皮质增生	常见
ARX	19	45	X连锁婴儿痉挛综合征	Rare
MECP2	250	53	雷特综合征	常见
GJB1	16	53	夏科‐马里‐图斯病	常见
ABCD1	33	57	X连锁肾上腺脑白质营养不良	中等
EMD	11	57	埃默里‐德赖富斯肌营养不良	中等
G6PD	16	58	葡萄糖‐6‐磷酸脱氢酶缺乏症	常见
GATA1	12	60	红细胞生成异常性贫血和血小板减少症	Rare
AVPR2	15	62	肾性尿崩症	Rare
EDA	37	63	少汗性外胚层发育不良	中等
SLC16A2	11	63	艾伦‐赫恩敦‐达德利综合征	Rare
FLNA	42	64	耳腭指综合征	Rare
EBP	24	64	X连锁点状软骨发育不良	Rare
RPGR	17	64	视网膜色素变性	常见
TAZ	17	64	巴思综合征	Rare
IDS	16	64	亨特综合征	中等
FGD1	8	64	阿尔斯科格‐斯科特综合征	Rare
GPR143	6	65	眼白化病	中等

2.3. 成本分析与可扩展性

我们全基因组测序的生物信息学成本主要在于数据存储，而非计算处理。作为临床工作流程的一部分，FASTQ、BAM和VCF文件会定期归档到一个复制的辅助存储站点。我们还存储MD5校验和，以确保在传输和存储过程中的数据一致性。这些磁盘主要用于存储数据，对这些文件的任何访问都必须先将其复制回更快的存储位置（例如，主存储）。根据所用CPU时间、集群寿命以及生物信息学分析师的每小时成本进行计算，我们估计处理一个基因组的计算时间为77美元，所需人工操作时间为128美元。然而，每个基因组大约需要300 GB的存储来保存包含未比对和已比对读段以及最终变异检测结果的文件（相应的FASTQ、BAM和VCF文件）。在表3中，我们估计每个基因组每年的存储成本约为40至55美元，具体取决于该基因组是存储在主存储设备上还是辅助（或归档）磁盘上；我们实验室处理并存储一个基因组一年的生物信息学流程总成本约为245美元。测序数据的磁盘存储除测序本身和解释性组件外，还会产生显著成本，并且根据实验室需求和政策，数据在活动存储或深度存储上的可用性可能有所不同。该流程基于高性能计算（HPC）集群构建，因此在特定时间内处理多个基因组的能力在很大程度上取决于集群规模。在当前基础设施下，我们不会同时运行超过10个基因组。该流程在集群规模方面具有可扩展性，但我们尚未探索将此流程扩展到同时处理数百或数千个基因组。

表3. 全基因组测序数据存储的成本分析。主存储假设为无复制的活动存储，具有高输入/输出（I/O）能力。辅助存储假设为复制的深度存储，具有低输入/输出（I/O）能力。处理一个基因组以及在主存储和辅助存储上保留数据一年的成本为~$245。

存储类型	基因组/月（美元）	基因组/年（美元）	基因组/5年（$）
主要	4.42	53.04	265.20
次要	3.48	41.76	208.80
总计	7.90	94.80	474.00
## 3. 讨论
我们上文所述的经过验证的生物信息学工作流程已被用于处理数百个基因组。自最初实施该流程以来，我们已多次更新其功能并重新验证该流程的各个组件。我们流程的模块化特性使得在有新的或更新的注释资源可用时，只需重复进行注释、上传和过滤流程。类似地，也可以从Oracle数据库中的持久性变异数据中轻松地重新查询过滤结果。

我们还将我们的NA12878结果与美国国家标准与技术研究院（NIST）主导的基因组瓶中联盟（GIAB）数据集进行了比较[9]。另一个类似的工具称为遗传检测参考材料协调计划（GeT-RM），它允许用户查询同一数据集的不同比对结果，以确认您的流程发现的目标变异是否也被其他流程检测到。该数据集包含了由我们实验室生成的变异数据。这两个资源在验证临床基因组测序方面具有极大的价值。

基因组测序的局限性是一个重要的考虑因素。尽管通常被称为“全基因组”，但当前技术仍无法准确测定其中许多区域。这些区域对当前标准的临床全基因组测序分析不可及，包括具有高同源性的区域（读段无法唯一比对）、参考基因组存在错误的区域、具有多个参考单倍型的区域，以及超出测序读段长度的串联重复区域。这些区域可能与临床相关区域重叠，若不谨慎处理，可能导致假阴性结果。在表2中，我们标示了在我们的WGS检测中始终覆盖度最低的一部分区域。目前管理这些区域的技术包括：针对重复扩增疾病使用重复引物PCR，针对存在同源性问题的基因使用长片段PCR，以及利用基于等位基因的条形码技术将读段锚定到唯一的邻近区域。比对至GRCh38参考基因组可能有助于处理具有多个参考等位基因的区域，但这些区域如何整合到临床工作流程中仍有待明确。最终，长读长测序技术以及可能的从头组装方法的应用，将能够实现对这些困难区域更准确的比对和变异检测。

随着基因组测序变得越来越常规，该流程将不断进行进一步的改进。在我们的实验室中，这些改进主要集中在为变异过滤提供一个用户界面，以便我们的临床工作人员进行实时过滤，从全基因组测序数据中识别拷贝数和结构变异，以及升级流程以与GRCh38参考基因组对齐。这项工作需要一支具备广泛技能的生物信息学家团队，能够直接与临床工作人员和信息技术团队协作。找到这种协作模式可能是创建临床基因组测序项目最具挑战性的方面之一。

4. 方法与材料

4.1. 生物信息学分析流程

生物信息学分析流程以相当自动化的方式处理由Illumina临床服务实验室（美国加利福尼亚州圣地亚哥）交付的硬盘中的数据，直至生成用于临床解读的变异报告文件。当收到这些硬盘后，会在实验室中进行编号登记，并上传至我们的高性能计算（HPC）环境。硬盘中提供的BAM文件包含由因美纳开发的流程使用CASAVA（或近期数据交付中使用的ISAAC）比对后的读段。这些文件同时包含已比对和未比对的读段，因此将BAM文件转换回FASTQ格式时，可将比对信息与读段序列分离，而不会造成测序信息的丢失。

从FASTQ文件到VCF文件的数据处理主要采用各软件包推荐的参数 [10–13]。使用bwa 0.6.1-r104对测序读段进行与hg19参考基因组的双端比对。通过samtools 0.1.18对已比对的读段进行排序并去除PCR重复。使用Genome Analysis Toolkit (GATK) 2.2.5以及布罗德研究所（美国马萨诸塞州剑桥市）GATK开发团队提出的最佳实践流程中的建议，进行局部插入缺失重比对、碱基质量重校准、使用UnifiedGenotyper进行变异检测以及变异位点重校准。整个流程如图1所示。为了更高的计算效率，我们经过验证的流程包括通过将文件分割为较小部分并在集群中均匀分布来实现并行化。

该方法已在我们将此流程应用于MedSeq项目（一项评估基因组测序在临床实践中影响的随机临床试验）时进行了描述 [14]以及应用于发现一个导致新型隐性综合征的推定致病位点，该综合征表现为骨骼畸形和恶性淋巴增殖性疾病 [15]。此前，其他临床实验室在使用下一代测序技术时也成功应用了类似的流程[4,6,16–18]。

连续的变异注释通过一组脚本完成，这些脚本独立地对数据集进行广泛的信息注释，包括：（1）来自Alamut（Interactive Biosoftware，法国鲁昂）的转录本和基因注释；（2）来自变异效应预测器（VEP）的基因注释；（3）来自1000基因组计划、ClinVar和外显子组测序项目（ESP）的变异注释；以及（4）由我们实验室在GeneInsight（美国马萨诸塞州剑桥市）中维护的临床解读[19–22]。

作为一家临床实验室，我们的重点是那些已知或可能与疾病相关的基因组区域。因此，我们将基因组上检测到的变异注释限制在我们的目标区域内。目标区域包含所有映射到hg19的基因的RefSeq编码序列、药物基因组学变异，以及来自NHGRI GWAS目录的关联或风险变异[23]。这些区域两侧各扩展50个碱基对（bp）。在进行变异注释之前，VCF文件已被限制为此目标区域。此外，在过滤和覆盖度分析过程中，我们进一步聚焦于最可能包含可报告的致病性变异的区域。因此，我们创建了临床关注区域（ROI）文件，其中包含每个基因所有外显子的编码区，分别用于过滤和覆盖度分析的±15 bp或±2 bp。

样本数据和变异注释被上传到Oracle SQL数据库中。变异过滤流程旨在查询该数据库，并将感兴趣的变异以Excel电子表格形式返回，以便进行仔细的临床审查。变异过滤是一个不断演进的过程，过滤规格根据遗传病的常见程度、外显率和预测的遗传模式而有所不同。表4列出了用于生成最终报告的一些常见过滤方法。

通常，最终的变异报告是通过布尔逻辑组合上述多个过滤器生成的。整个生物信息学流程的概述，从比对和变异检测到最终报告生成，如图2所示。

表4。示例过滤方法。这些过滤器使用布尔逻辑应用于每个个体，以生成最终的过滤后变异列表。

筛选器名称	参数	描述
频率	X（例如，0.01或0.05）	保留ESP中具有频率的变异 1000基因组≤ X
功能丧失		保留可能暗示基因功能丧失的变异，包括以下序列本体注释的变异：序列本体关键词：移码_变异, 终止_获得, 终止_丢失, 剪接_受体_变异, 起始_密码子_变异, 剪接_供体_变异。
基因列表	基因列表（在 HGNC 中）命名法	基因过滤基于选择位于其内的变异特定基因。我们检查一个变异是否被注释有兴趣的临床区域内的基因符号
已报告致病性		选择被分类为致病的或可能致病的变异在变异数据库（包括ClinVar）中为致病的变异
GeneInsight		选择被分类为致病的或可能致病的变异在我们的内部GeneInsight数据库中为致病的
复合型杂合子		如果至少有两个功能丧失和错义变异，则选择它们基因中可能影响两个等位基因功能的变异

4.2. 生物信息学验证

从科里尔细胞库（美国新泽西州卡姆登）获取并在因美纳平台上测序的NA12878全基因组测序数据集被用于验证该流程的有效性。该流程在同一数据集上运行了两次，以测试其稳健性，证明其能够对给定数据集持续提供一致的结果。在本次测试中，流程从两个不同的上游输入点启动，分别是因美纳比对后的BAM文件和FASTQ文件，并对测试生成的最终变异输出文件进行了比较。另一部分验证来源于对变异检测流程敏感性和特异性的评估。这些指标是通过使用已评估NA12878基因组水平变异的不同标准数据集进行估算的。首先将变异检测结果与实验室内部生成的测序数据进行比较。

NA12878已被许多实验室作为验证标准，我们已对来自心肌病、听力相关疾病的195个基因中的约700 kb序列进行了测序，利用桑格测序、基于阵列的测序和靶向下一代测序等多种正交技术，进行与呼吸系统、努南综合征和马凡综合征相关的检测[2,24,25]。因此，我们仅针对这一特定目标区域评估了全基因组测序流程的敏感性和特异性。根据GATK早期版本的建议，我们调整了所调用变异的QD和FS过滤参数，以寻找最优阈值。此外，我们还检查了我们的全基因组数据与1000基因组计划高覆盖度全基因组测序数据集（>30×平均覆盖度）中NA12878报告的变异的一致性。[19]通过每个VCF文件中是否存在相应变异（变异一致率）以及每个变异所调用基因型的一致性（基因型一致率）来评估变异。

4.3. 低覆盖度区域的特征分析

我们选取了在Illumina临床服务实验室进行测序，并于2014财年通过我们的生物信息学分析流程处理的15个基因组，以生成RefSeq基因编码序列（外加2个碱基对的缓冲区域以覆盖剪接位点）的覆盖度列表。我们对GATK 2.2.5进行了定制，运行CallableLoci以确定每个编码碱基的可检测性，其中“可检测”定义为：一个碱基至少具有8X的高质量读段覆盖度，且该区域总读段中低比对质量读段的比例低于10%。该覆盖度列表中每个基因的可检测碱基对数目是基于15个基因组计算的平均值。我们的流程未对线粒体染色体进行变异检测，因此线粒体基因未包含在该列表中。

接下来，我们希望评估在覆盖度列表中发现的低覆盖区域的临床相关性。为此，我们选择了ClinVar数据库中的已报道变异集，这是一个由医学遗传学界提供的关于变异及其临床意义的公共档案[26]。我们筛选了ClinVar数据集XML版本（2015年8月发布）中“临床意义”字段标注为“致病的”和“可能致病的”变异。ClinVar对变异有多种分类方式，我们将具有唯一MeasureSet ID的变异定义为一个独特变异。为了识别当前在临床测序实验室中检测的临床相关基因，我们重点关注由实验室和工作组提交的变异。为此，我们排除了来自大型研究级数据库（包括OMIM）的变异分类。对于具有冲突临床意义分类的MeasureSet ID（例如，被一家实验室报告为致病，而被另一家实验室报告为良性），我们也予以排除。然而，若某变异被一家实验室判定为致病，而在另一家实验室中被判定为意义未明，则视为一致，不作为冲突变异排除。

4.4. 成本分析与可扩展性

全基因组测序（WGS）比当前其他测序方案需要更多的计算资源和存储空间。由于临床测序要求数据必须保存一定时间，因此会产生额外的成本。计算成本是根据购买计算节点的价格，并结合完成一次基因组运行所需的CPU小时数以及硬件使用寿命（Partners HealthCare保证五年）进行估算的。人工操作时间成本是根据我们每小时收取的生物信息学分析费用，乘以成功启动流程、将变异结果交付给遗传学家以及归档数据所需的时间来估算的。我们还基于当前内部提供的价格计算了存储成本。用于运行基因组测序流程的集群资源由多个生物信息学流程共享。目前，我们拥有10个计算节点，每个节点配备128 GB内存、16核英特尔至强2.60 GHz处理器，运行CentOS 6.5和IBM LSF作业调度器。

5. 结论

我们已在Partners HealthCare 个性化医疗的实验室中描述了临床全基因组测序的实施。我们的实施方案包括专门的脚本，用于执行以下任务：（1）比对与变异检测；（2）变异注释；（3）将数据导入SQL数据库；以及（4）变异过滤。这些组件均已通过临床验证，并可独立更新。全基因组测序已被证明能够可靠地识别遗传变异，但主要挑战仍在于对这些变异的临床解读。这对于绝大多数非蛋白质编码变异而言尤为如此，而这些变异在外显子组测序中大多被忽略。随着非编码区公共知识库（如ENCODE数据和GWAS信号）的不断扩大，我们可以开始建立包含更大比例非编码区的过滤方法，使其在临床环境中逐渐具备应用价值，最终实现完整基因组测序的全部效用。