基于全基因组测序的法医单核苷酸多态性系谱推断研究*
2023-11-23赵雯婷方之晓徐景怡李彩霞
谢 全 桂 虚 赵雯婷 方之晓 徐景怡 刘 京 李彩霞**
（1）中国人民公安大学侦查学院，北京 100038；2）公安部物证鉴定中心，现场物证溯源技术国家工程实验室，北京市现场物证检验工程技术研究中心，北京 100038；3）四川省达州市公安局达川区分局，达州 635000）
单 核 苷 酸 多 态 性 （single nucleotide polymorphism，SNP）在人基因组中分布广泛，突变率低，遗传稳定性高，是法医遗传领域的第三代遗传标记，常用于族群地域推断、体貌表型预测、亲缘关系分析等［1］。其中，SNP在亲缘分析尤其是远距离亲缘分析方面具有独特的优势。第二代遗传标记短串联重复序列（short tandem repeat，STR）标记通常用于亲子鉴定、全同胞鉴定，但是对二级以上亲缘关系难以准确鉴定［2］。2018 年，美国“金州杀手案”的破获就是通过高密度SNP预测出嫌疑人的7级远亲关系，进而找到嫌疑人，该技术被称为法医SNP 系谱推断技术，并被《科学》（Science）杂志评为年度全球十大科学突破之一［3］。
法医SNP 系谱推断技术迅速成为近年来法医领域的研究热点，并为数百起冷案积案的侦破提供了关键线索［4］。本项目组也建立了相关技术并在命案积案侦破中发挥作用［5］。法医SNP 系谱推断技术是基于高密度SNP 数据，通过个体间的共祖片段（identity by descent，IBD）长度算法、状态一致性（identity by state，IBS）共享统计量等算法预测个体间的亲缘关系等级［6］。获得高密度SNP分型数据的方式主要为全基因组SNP 芯片技术和全基因组测序（whole genome sequencing，WGS）技术。其中，由于全基因组SNP芯片的价格优势，成为目前法医SNP 系谱推断最常用的数据获取方式［7］。SNP芯片是在一块基片表面固定序列已知的靶核苷酸的探针，用于特定SNP 位点组合的检测分型。国内外常用芯片包括GSA 芯片（65 万SNP）、CytoSNP 芯片（85 万SNP）、Affymetrix 定制芯片（20 万~260 万SNP）等，均 为Illumina、Affymetrix 等国外产品。而且，SNP 芯片技术对DNA 质量的要求较高，一般需200 ng 以上，250 pg时7级亲缘预测准确性由99.9%降至62.6%，对血样、精斑、唾液等样本检测效果好，但是微量降解等现场常见生物检材的检测能力差，模拟降解检材平均DNA 片段大小为150 bp 时，其准确率降至0［7］。WGS技术又称为大规模并行测序，区别于传统的Sanger（双脱氧法）测序，能够一次并行对大量核酸分子进行序列测定的技术。WGS 检测方案的一个优势在于，DNA 片段化处理是文库构建的第一步，增加了从降解的法医样本中获得SNP分型结果的可能性［8］。近期有研究发现，WGS 技术在微量降解等法医现场检材的检测能力优于全基因组SNP芯片技术［9］。
本项目组前期已构建起基于全基因组SNP 芯片技术的法医SNP 系谱推断技［5-6］，为建立国产平台检测能力，并进一步提升法医现场检材检测能力，本文使用华大MGISEQ-200RS 测序平台进行高深度的WGS，旨在探索从高深度WGS数据中获得准确的SNP 分型，并进行法医SNP 系谱推断的可行性和准确性研究。
1 材料与方法
1.1 样本来源及DNA提取
采集本实验室来自不同家系的静脉血样本A、B、C，标准品使用NA12878（具有北欧和西欧血统的犹他州女性）。所有样本在采集前均签署知情同意书，本研究通过了公安部物证鉴定中心伦理委员会审查（编号：2021-006）。静脉血均使用QIAamp Blood Midi Kit试剂盒（QIAGEN公司，德国）提取DNA，使用QubitTMdsDNA HS Assay Kit试剂盒在Invitrogen Qubit 4（Thermo Fisher 公司，美国）上进行DNA 定量，使用JY-SPDT 电泳仪（君意东方公司，中国）判断DNA完整度。
1.2 文库制备和测序
质检合格的DNA，使用MGIEasy酶切DNA文库制备试剂套装（华大智造公司，中国），对200 ng 起始样本进行酶切文库构建。利用MGISEQ-200RS高通量测序试剂盒将单链环状文库制备成DNA 纳米球（DNA nano ball，DNB），采用PE150 读长测序策略，进行平均深度30×WGS。所有样本均使用MGISEQ-200RS测序仪进行测序。
1.3 测序数据处理
利用FastQC 软件对测序仪下机数据的fastq 文件进行质检，使用Fastp 软件去除低质量片段和接头序列，生成Clean Fastq 文件。使用BWA 工具，将Clean Fastq 数据比对到人类参考基因组（hg19）上生成bam文件。用Sort工具，将bam文件按照染色体组型进行排序，生成sort.bam 文件。使用GATK软件Markduplicate工具标记重复的读段，选取唯一对比读段进行后续分析。通过GATK软件中BQSR工具用于矫正测序碱基的真实质量值，设置参考SNP 数据集作为真实变异集，这些数据集之外的位点，若测序数据与参考基因组存在差异，则认为是错误，以此来判定并矫正测序碱基的质量值。使用GATK的HaplotypeCaller工具，进行二倍体基因组的变异检测，识别出潜在变异区域。
使用Bcftools 软件从全基因组vcf 文件中提取感兴趣的位点，利用Python 3脚本过滤低质量位点进行质控，并删除Indel。使用Plink 软件将vcf 文件转换为二进制格式，与参考样本合并进行法医SNP系谱推断。
1.4 家系参考样本
家系参考样本来自项目组志愿者共计3 个家系，亲缘关系对应等级数量分布如表1，包含1~8级 亲 缘，234 份 样 本［6］，使 用Wegene GSA 芯 片（安澜智能公司，中国）获得645 199 个常染色体SNP 位点分型，参考刘京等［10］法医系谱亲缘关系等级进行划分。族群推断参考样本来自公开的千人基因组数据库。
Table 1 The quantity distribution of each kinship degree among samples
1.5 系谱推断
通过本项目组开发的亲缘关系预测系统，分别使用IBS和IBD算法对样本进行亲缘关系预测［5-6］。
IBS 算法通过计算亲缘关系系数∅和零IBD 共享统计量π0推断两个体间亲缘关系等级。零IBD共享统计量π0代表两个体在一个SNP 位点上共享同一祖先零个等位基因的概率，亲缘关系系数∅表示从两个体中随机抽取的两个等位基因来源于同一祖先的概率。根据系统预测的所有个体间∅和π0与推断等级标准范围对比，可进行个体间亲缘关系等级推断［6］。
IBD是指由于在遗传过程中，染色体会发生重组，两个有亲缘关系个体的基因组中有一段DNA来源于同一个祖先，通常用重组频率的测量单位厘摩（centimorgan，cM）衡量，IBD 片段越长，说明亲缘关系越近［11］。预测系统基于背景人群参考数据进行同源染色体分离，计算两两个体之间共祖片段IBD的cM总长度，进而推断个体间亲缘关系等级［12］。
法医系谱推断效能通过绝对准确率（accuracy，AC）、置信区间准确率（confidence interval accuracy， CIA）、 假 阴 性 率 （false negative，FN） 和假阳性率（false positive，FP）等指标评估。绝对准确率为预测亲缘等级与调查等级一致的关系对数占关系对数的比重，置信区间准确率为预测等级在调查等级的上下一级浮动的关系对数占等级数的比重，假阴性率为调查有亲缘关系被预测为无关的对数占亲缘对数的比重，假阳性率为调查无关被预测为有关的对数占亲缘对数的比重。绝对准确率和置信区间准确率越高，假阴性率和假阳性率越低，系谱推断体系越准确。
1.6 族群分析
以千人基因组数据库［13］样本为参考数据集，利用全基因组SNP 分型，使用Plink v1.9 软件进行PCA 分析，R 软件绘制PCA 图。提取实验室前期研究的27AISNP 体系位点，使用DAA 软件进行欧洲、非洲、东亚三大人群祖先来源推断［14］。
2 结 果
2.1 测序质量评估
从3 份中国北方人群的静脉血和标准品NA12878中获得的4份DNA，分为4次样本制备和测序，其中静脉血样本用于亲缘关系推断，标准品用于测序准确性评估。
测序芯片FCL PE150属于单线型泳道，满载输出约150 G fastq测序数据。表2中列出了芯片泳道的测序指标，其中每条泳道平均产生（642.65±25.63）百万次读数，平均芯片利用率为73.69%，平均Q30 比率为90.50%，过滤后有效点占比ESR的均值为73.69%。MGISEQ-200RS 平台与已发表的BGISEQ-500RS 平 台［15］和MGISEQ-2000RS 平台［16］测序质量对比，平台获得相似的Q30值。测序平均Q30 为90.50%，能够进行准确的基因组测序，是后续一致性和推断准确性分析的数据基础。
Table 2 Run metrics for 4 runs in MGISEQ-200RS sequencing platform
4 份测序样本总计获得2 612 M reads，比对到参考基因组hg19 的reads 为2 608 M，比对率为99.8%。平均比对质量为54.18，去除重复序列后样本深度在30×以上，基因组1×平均覆盖率为91.60%。GC含量均在40%左右，与人类的基因组GC含量一致，说明测序过程碱基的偏好性低，且读段覆盖更均匀［17］。所有样本的详细测序和SNP覆盖指标见表3。
Fig.1 Coverage rate of various samples in different depths
Table 3 Statistical table of bamQC parameters
全基因组覆盖率随深度阈值的变化如图1 所示，在1×至5×深度区间，随着深度阈值的提升覆盖率降低缓慢，但在大于5×区域下降快速。覆盖率随深度增加的变化趋势在平均深度处达到最大。以平均深度为阈值，各样本覆盖率接近50%，说明各位点的测序深度集中在平均测序深度附近，但不可避免存在一部分位点深度较低及未检测到的情况。不同深度与覆盖率的关系反应了建库片段回收的随机性与高通量测序的随机性。随着平均测序深度的增加，数据量越大，冗余越多，随机性会相应的降低。因此，平均深度30×的WGS，可获得大量的系谱推断SNP位点。
从4 份测序数据中提取Wegene GSA 常染色体645 199个SNP位点，获得1×及以上深度认为该位点被成功提取。NA12878、A、B、C 样本成功提取的SNP 位点数目分别为：644 897、645 029、644 882、645 014个。对SNP位点进行测序深度统计（图2），平均测序深度处SNP 位点数量达到峰值。
Fig.2 Quantity distribution of SNP at different depths
2.2 SNP位点分型一致性评估
从千人基因组数据库NA12878 的SNP 位点分型文件提取Wegene GSA常染色体645 199个位点，共获得621 128 个SNP 位点。从33.6×NA12878 标准品测序数据中提取该621 128个SNP位点，用于SNP 分型准确性评估，成功提取出621 127 个SNP位点。
33.6×未质控标准品SNP 位点与1000G 数据库［13］中NA12878的SNP位点进行一致性对比，一致率达99.62%，分型不一致位点2 365个。未质控标准品分型为杂合子而千人数据分型为纯合子的SNP 位点表现出位点分型的假阳性，假阳性位点589个。未质控标准品测序分型为纯合子而千人数据分型为杂合子的位点表现为杂合性丢失，为假阴性位点，共1 348 个。425 个SNP 位点表现为都为纯合位点且不一致，3个SNP位点表现为杂合不一致。分别统计假阳性、假阴性、纯合不一致和杂合不一致位点测序深度（depth，DP）及测序质量值（genotype quality，GQ）（图3）。
Fig.3 Na12878 and 1000G inconsistent SNP mass distribution
不一致位点中，假阴性与纯合不一致位点主要集中于测序质量值GQ较低的区域，为确保在较高检出率的前提下获得高一致率，以DP大于等于4，GQ大于等于10为阈值，对测序样本进行质控。质控样本的假阴性与纯合不一致位点数减少明显，假阳性位点数量却没有得到明显改善（图4）。假阳性位点存在的问题主要原因分为3点：a.文库构建过程中PCR 扩增错误；b.测序仪对位点荧光信号的错误判读；c.千人基因组数据该位点测序深度较低导致的SNP 位点杂合性丢失。因此，提高质控标准难以解决假阳性问题，并且可能造成大量分型正确位点的丢失。质控后的杂合不一致位点并没有被 过 滤 掉 ， 33.6×NA12878 的 rs6686181、rs72781591、rs9270230 分型为AT、GT、CT，位点测序深度分别为11×、26×和53×，1000G数据中rs6686181、rs72781591、rs9270230 分型AG、AG、CA。rs72781591 位点的26×测序中18 次检测到碱基G，2 次检测到碱基A，6 次检测到碱基T，rs9270230 位点出现了与rs72781591 一样的情况，23 次检测到碱基C，10 次检测到碱基A，20 次检测到碱基T。MGISEQ-200RS平台采用红绿双色荧光激发系统，碱基A能被红绿双色荧光激发，碱基T仅被绿光激发，碱基C仅被红光激发，碱基G不能被激发，可能为MGISEQ-200RS 的双色荧光激发系统在判读该位点时绿光信号偏强，从而导致碱基A、G被错误认为碱基T。
质控过滤后，纯合不一致和假阴性位点明显减少，不一致位点减少超过50%，一致率得到提升，但同时导致检出率的降低（表4）。提取率为深度大于等于1×的SNP位点个数在Wegene GSA芯片常染色SNP 位点集合的比重，检出率指通过深度大于等于4×、QG 值大于等于10 指控条件的SNP 位点个数在Wegene GSA芯片常染色SNP位点集合的比重。
Fig.4 Comparison of inconsistent loci in NA12878 after quality control
Table 4 Consistency comparison parameter statistics
从测序样本中提取Wegene GSA 常染色体位点，各样本进行GQ为10的质控过滤，质控后检出率和一致率如表5。由于A、B 样本仅有Wegene CGA 芯片分型数据做为对比。A、B 检出率为99.34% 和99.28%，提 取Wegene GSA 和Wegene CGA 芯片交集位点500 753 个，一致率分别为99.84% 和99.90%。C 检出率为99.36%，与其Wegene GSA 芯片分型进行对比，C 共有一致位点636 587个，一致率为99.62%。
Table 5 Samples call rate after quality control
2.3 SNP位点在染色体上分布情况
质控后，样本A、B、C 高深度测序数据均获得64 万以上的位点分型数据，位点在染色体上的分布如图5，结果显示各样本与GSA芯片位点分布完全一致，较为均匀的分布于22 对常染色体上，局部区域特征上端粒区域SNP 分布略高于染色体中部区域，其中所有样本6号染色体28Mb HLA区域SNP 分布密集。位点在染色体上分布情况，说明WGS不会导致染色体上某区域特异性检测问题，并且满足高密度SNP 系谱推断技术位点在人类基因组中均匀分布的要求。
2.4 亲缘关系预测
A、B、C均为平均深度30×测序数据，与家系参考样本Wegene GSA芯片数据合并，进行亲缘关系计算，将所有个体间预测的亲缘关系等级与实际调查的亲缘关系进行比较，评估亲缘关系预测准确性。
表6中展示了合并后数据进行IBS亲缘关系预测的准确性，本文把8 级及8 级以上亲缘定义为8th。从表中可以看出，1~4 级亲缘关系的CIA 为100%。由于IBS 算法的局限性，随着亲缘关系等级的增加，预测准确率也随之降低，5级置信区间准确率降为81.8%。5 级开始出现假阴性（FN），5级之后的亲缘关系绝对准确率（AC）明显下降。假阳性率（FP）为2.8%，将无关亲缘预测为8 级亲缘。30×测序数据使用IBS 算法能够预测四级及以内的亲缘关系。
Table 6 The evaluation of the accuracy of genetic relationship prediction by IBS
表7中展示了合并后的数据进行IBD亲缘关系预测的准确性。表中显示，IBD算法预测1~5级亲缘关系的置信区间准确率为100%，仅在预测8 级及以上关系对时出现假阴性，4级以上绝对准确率开始下降。假阳性率为2.8%，将无关亲缘预测为8级亲缘。因此，30×测序数据使用IBD 算法能够用于7 级及以内亲缘关系预测，将7 级及以下亲缘预测为有关系亲缘的准确率为100%。图6显示，30×测序数据使用IBD算法获得亲缘关系对的共享IBD长度与芯片数据基本一致。
Fig.6 Share IBD length distribution in kinship pairs
Table 7 The evaluation of the accuracy of genetic relationship prediction by IBD
经t检验，30×测序数据获得的IBD 片段长度与芯片数据对比，P值为0.93，远大于0.05，无显著差异。WGS数据可用于亲缘关系预测（图7）。
2.5 族群推断
图8显示千人数据库5大人群共2 504个个体与4个质控测序样本合并数据的主成分分析结果，左下为欧洲人群，左上为东亚人群，标准品NA12878 在参考数据中欧洲人群处聚集，趋于左下方。中国北方样本A、B、C 与东亚人群聚集，都趋向于左上方。主成分分析结果与调查结果一致，在全基因组层面能够区分出样本的族群来源。
Fig.7 t Test results
Fig.8 Results of whole genome principal component analysis
从WGS数据中提取4份样本的27AISNP分型，4 个样本的27SNP 位点与一代测序分型完全一致，准确率为100%。用该位点组合进行的族群推断结果如图9，主成分1 和主成分2 共解释了总方差的62.27%，与全基因组主成分分析结果一致。表明能够通过从WGS数据中提取特定位点集合的分型，准确预测出族群来源。
Fig.9 27AISNP inference ancestry results
3 讨 论
亲缘个体之间拥有相同的且来自共同祖先的基因片段称为IBD，随着减数分裂次数增加，非同源染色体重组次数增加，所以亲缘关系越远，IBD越短。法医SNP 系谱推断就是基于此原理进行亲缘关系的分析［9］。高密度SNP 数据是进行系谱推断的基础。在项目组前期研究中发现，针对微量降解检材，SNP芯片的检出率较低［10］。而WGS对法医中常见的微量降解检材具有优势，能获得较多的位点分型。国外有研究利用ThruPLEX®试剂盒针对仅50 pg 片段化DNA 构建文库，进而WGS 获得SNP进行系谱推断并破获案件［18］。针对高度降解的毛干等检材，WGS 也能检测出大量SNP 分型数据用于法医SNP系谱推断［19］。目前，国内尚无WGS技术用于法医系谱案件的研究报道。因此，本研究在国产MGISEQ-200RS 平台上构建WGS 法医SNP 系谱推断体系，并对其位点分型准确性和预测准确性进行探讨。
本研究通过标准品的测序数据与千人基因组数据对比，A、B、C 与芯片数据对比，探讨了高深度WGS 的一致性。以深度为4 质量值为10 质控过滤，在保证检出率的条件下明显降低假阴性和纯合不一致位点数量，标准品质控后检出率为99.21%，一致率达99.84%，不一致率低于Tillmar等［18］预测的30×深度不一致率0.8%。提取率高于Bode Technologies 研 讨 会 “Forensic Genealogy：Unlocking the Science of Genealogy”所研究的，血液WGS 起始量250~50 ng 提取率在97%~87%之间的范围，较GSA 芯片的检出率100%~95%低，但对极低起始量（1~0.25 ng）样本，WGS 提取率的稳定性远高于GSA 芯片［9］。对来自3 个不同家系的3份样本进行平均深度30×的WGS，与芯片分型数据对比，SNP 分型一致性达99.62%以上。有研究表明，高深度WGS 数据中假阴性位点集中在低深度位点部分，且低频突变（0.5%
目前，已有针对全基因组SNP 芯片数据预测亲缘关系准确性的系统研究。起始量200 ng 的Illumina GSA芯片分型数据IBD算法预测7级亲缘，准确性为74.6%［7］。Wgene GSA 芯片分型数据IBS算法预测4 级亲缘，绝对准确性为68.0%，置信区间准确率为98.1%［6］。在2020 年，国外开始出现WGS 获取积案中骨骼SNP 数据用于系谱推断的报道，获得平均深度32.3×共1 035 274个SNP，在第三方数据库中得到IBD 大于10 cM 的36 个潜在亲缘，但未能调查证实潜在亲缘的真实性［18］；WGS检测陈旧血迹获得60×数据共1 861 000 SNP 分型，在第三方数据库检索到共享IBD 为347 cM 的潜在亲缘，通过侦查锁定嫌疑人，利用STR 图谱认定为作案人［21］。但除了上述WGS系谱推断的个案报道，国内外尚无在真实家系中研究不同等级亲缘推断准确性的系统研究。由于系谱推断的结果将成为指导侦查方向的重要依据，为避免造成方向错误和浪费警力，非常有必要针对该技术进行系统研究，评估预测结果准确性。通过本研究发现，使用IBS算法，可用于1~4 级亲缘关系的预测，4 级置信区间准确率达100%。4 级亲缘推断准确性较已有研究［6］稍高且无假阴性，可能是由于4 级关系对数量较少，预测不准确的小概率事件并没有发生。使用IBD 算法，能够对1~7 级亲缘进行预测，7 级及以内亲缘预测为有亲缘的准确率达100%，假阴性为0%。IBD 算法对预测8 级亲缘也具有一定的价值，即使在IBD片段长度不能达到预测分级的标准时，也可作为案件侦查的一个排查方向。经t检验，30×测序结果与芯片预测结果无显著差异，可用于法医系谱推断。
本研究还探究了使用WGS 用于族群地域推断的可行性，全基因组SNP 位点和提取AISNP 体系位点推断结果均与调查结果一致，即WGS 法医SNP 系谱推断体系可用于DNA 供者的族群来源推断。
4 结 论
基于国产测序平台的WGS 获取的SNP 位点组合与全基因组SNP 芯片的分型一致率高，可与全基因组SNP 芯片数据合并分析，进行准确的系谱推断。此外，本研究使用的测序仪器、耗材和试剂均为国产，数据安全性高，易于在法医实验室普及。下一步将继续优化体系性能，探索并建立不同测序深度的测序技术对微量降解检材的检测能力，以及不同等级亲缘关系的预测性能。
