通过靶向纳米孔测序对恶性疟原虫疟疾进行灵活且经济有效的基因组监测

多重 PCR panel 的开发

NOMADS8 和 NOMADS16 面板是使用 beta 版本的 multiply（称为 pf-multiply）生成的，可从以下位置获取： https://github.com/JasonAHendry/pf-multiply（设计文件为 游牧者8 和 游牧者16）。 两者都使用 BED (*.bed) 文件来描述 mdr1 扩增子。 首先使用以下命令生成 NOMADS8：

蟒蛇multiply.py -d设计/pf-nomads8-mdr1part.ini

NOMADS16是使用pf-multiply的augment命令创建的。 NOMADS8 多重 PCR 引物以等摩尔量组合，总引物浓度为 0.6 μM； 在 R9.4.1 Flongle Flow Cell 上使用模拟样品进行初始测序后，根据观察到的扩增子丰度粗略调整引物浓度（加倍或减半），保持总引物浓度固定。 对 NOMADS16 重复相同的过程； 即进行一轮引物浓度调整。

创建恶性疟原虫和人类 DNA 的模拟样本

我们订购了实验室菌株 3D7、Dd2、GB4 和 HB3 以及柬埔寨田间衍生菌株 IPC 5202 (kelch13 R539T) 的恶性疟原虫基因组 DNA； IPC 4912（kelch13 I543T）、IPC 3445（kelch13 C580Y）； 和 IPC 3663 (kelch13 WT)63 来自 BEI 资源（www.beiresources.org）。 为了创建 10,000 p/μl 体外 DNA 混合物，我们将来自 36 个 HapMap 细胞系的池中的 25ng/μL 人类基因组 DNA 稀释至 0.25ng/μl69。 然后将 DNA 混合物以不同的数量和比例组合，以复制不同比例或 COI 的混合感染，和/或在额外的人类基因组 DNA 中连续稀释，以复制较低的寄生虫血症感染（补充表 2）。 为了验证 hrp2/3 缺失，如前所述，将寄生虫系 3D7 (NF54)、Dd2 和 HB3 在从 DRK Blutspendienst Nord-Ost gemeinnützige Gmb 获得的商业红细胞中以 5% 血细胞比容进行培养70。 使用 Qiagen 血液和组织试剂盒在用 0.15% 皂苷裂解的寄生虫沉淀上提取所有品系的基因组 DNA。 如上所述，将提取的 DNA 与人类基因组 DNA（Roche，11691112001）结合以产生 10,000 p/μl 库存。 较低的寄生虫血症菌株是通过将 10,000 p/μl 库存稀释到人类基因组 DNA 中 2 倍连续稀释而产生的。

现场样本采集

赞比亚的样本来自两项研究。 第一个样本是根据赞比亚国家卫生研究局根据卫生部实验室实验室质量改进研究 (NHRA000004/16/11/2021) 授予的道德豁免收集的。 在赞比亚西部省卡奥马一家诊所就诊的有症状患者被诊断为 RDT，同时还收集了显微镜载玻片和 DBS（在 Whatmann No3 滤纸上）。 所有样本均已去识别化，并且没有记录人口统计或临床数据。 按照制造商的说明，使用 Qiagen QIAamp 试剂盒从 DBS 中提取大量 DNA。 通过光学显微镜对薄膜血玻片定量寄生虫血症。 第二组是在赞比亚卫生部通过 ERES Converge IRB 下的国家疟疾消除中心进行的治疗功效研究 (TES) 期间收集的，该研究对赞比亚选定地点的蒿甲醚-本芴醇、青蒿琥酯-阿莫地喹和双氢青蒿素-哌喹进行了治疗功效测试。 TES 研究是常规进行的，旨在评估三种用于治疗单纯性疟疾的 ACT 抗疟药物的疗效。 到赞比亚西北部索卢韦齐区一家诊所就诊的有症状患者通过 RDT 被诊断为疟疾，同时收集了显微镜载玻片和 DBS。 对每个阳性 RDT 的寄生虫血症进行量化，任何寄生虫血症达到 1000 个寄生虫/μl 或更高的患者都可以选择通过书面同意程序参加该研究。 收集其他临床信息，例如发烧状态和人口统计数据（即年龄、身高、体重）。 患者的家庭住址被记录下来用于研究随访目的，但在进行任何分析之前 DBS 的身份信息已被取消。

实验室方案和测序

完整的实验室方案，包括材料和引物序列，可在protocols.io 上在线获取（https://www.protocols.io/） 作为 “恶性疟原虫疟疾的具有成本效益的靶向纳米孔测序”。 简而言之，对于每个样品，使用 10–40ng 提取的基因组 DNA 作为 50μl sWGA 反应中的模板51，但减少了 phi29 DNA 聚合酶 (NEB #M0269S) 以最大限度地降低成本。 之后，2μl sWGA 产品用作 25μl 多重 PCR 的模板，使用 KAPA HiFi 聚合酶（罗氏 #KK2101）和 NOMADS8 或 NOMADS16 引物池。 使用 0.5X 比例的 NEBNext 样品纯化珠 (NEB #E7103) 清洁多重 PCR 产物，并在 15 µl 无核酸酶水中洗脱。 使用 Qubit (ThermoFisher #Q32854) 对 DNA 洗脱物进行定量，并将 100 至 600 ng 的 DNA 进行条形码和测序。 我们使用 Oxford Nanopore Technologies (ONT) 天然条形码连接测序试剂盒（带有 EXP-NBD104 的 SQK-LSK109、用于 R9.4.1 流动池的 EXPNBD114；用于 R10.4.1 流动池的 SQK-NBD114.96）的独特条形码与每个样品连接一种改进的一锅式条形码协议46。 然后在接头连接和测序之前汇集样本，我们遵循 ONT 协议建议。

生物信息学管道

对于使用 R9.4.1 Flow Cells 的实验，MinKNOW 软件生成的 FAST5 文件使用 Guppy (v5.0.11) 进行碱基调用，最小质量分数阈值为 8。对于 Flongle 实验，我们使用超准确 (SUP) 碱基调用模型对于所有其他实验，我们使用高精度 (HAC) 碱基识别模型。 对于使用 R10.4.1 Flow Cells 的实验，POD5 文件使用 dorado（v0.34；v0.34； https://github.com/nanoporetech/dorado）使用超精确（SUP）模型。 然后使用 Guppy (v5.0.11) 对 FASTQ 文件进行多路分解，无需设置 ––require_both_ends 标志，即使用单端多路分解。 解复用的 FASTQ 文件被映射到从 PlasmoDB 下载的恶性疟原虫 3D7 参考基因组的第 52 版71 （https://plasmodb.org）使用小地图272 (v2.24-r1122) 和 -ax-ont 参数设置。 在生成的 BAM (*.bam) 文件中，使用 samtools 识别无法映射到 3D7 参考的读取73 (v1.16)，使用命令 samtools view -f 0x904。 这些未映射的读数使用 samtools fastq 转换回 FASTQ 文件，然后重新映射到从 NCBI 基因组数据库下载的 GRCh38 人类参考基因组 (https://www.ncbi.nlm.nih.gov/genome）并随后从下游分析中排除。 来自目标的读数被定义为与从 PlasmoDB 下载的恶性疟原虫 3D7 参考基因组第 52 版的基因特征格式 (GFF) (*.gff) 中定义的编码序列重叠的读数71。 使用 Clair3 的奇点图像对映射到 3D7 参考基因组的读数进行变体调用47 （v1.0.4； https://github.com/HKU-BAL/Clair3）在二倍体模式下，设置了标志 ––platform=’ont’ ––include_all_ctgs ––enable_phasing。 对于图 1 中的 SNP 调用分析 4，我们通过设置 ––var_pct_full=1.0 和 –ref_pct_full=1.0 将所有变体发送到对齐模型。

SNP 调用准确性分析

下采样读取

我们使用 samtools 将映射到 3D7 参考的读数划分为与每个目标基因重叠的读数73 (v1.16)，从而为我们的每个目标生成 BAM 文件。 对于每个目标 BAM，我们使用 samtools view 命令和 -s/––subsample 标志对读取进行下采样，以实现所需的读取数量。 对 Dd2 和 HB3 样品重复此过程； 对于每个目标，下采样至 10、20、30、40、50、60、70、80、90 和 100 个读数，每个读数下采样 10 次。 因此，对于每个给定的深度和目标，我们为 Dd2 和 Hb3 生成了 10 个随机下采样的 BAM 文件； 总共 20 个重复。 在图中。 4b “全部”类别是通过连接以这种方式为 NOMADS8 面板中的所有目标（不包括 msp2）生成的 BAM 文件而创建的。

创建一组真实的变体

Dd2 和 HB3 已使用 Pacific Bioscience Sequencing SMRT 技术进行深度测序并组装，所得 FASTA (*.fasta) 序列可在 PlasmoDB 上使用71。 为了识别这些组件中相对于 3D7 参考基因组的变异，我们从 FASTA 文件中模拟了高质量（Phred 60）无错的计算机读取，并使用 minimap2 (v2.24-r1122) 将它们映射到 3D7 参考基因组，然后使用 bcftools 识别变体74 (v1.16) 编译和调用命令。 特别是，我们基于 Dd2 和 HB3 的 GFF 文件，通过提取跨越目标 +/-4kbp 的 FASTA 序列，为 NOMADS8 面板中的每个目标模拟了 60 个无错误读取（一半正向链和一半反向链）。

分层变异调用比较。 我们使用了hap.py这个工具49（https://github.com/Illumina/hap.py）从 Docker 镜像 (jmcdani20/hap.py:v0.3.12) 中计算不同目标区域的变体调用准确性，并与上述真实变体集进行比较。 为了将准确度测量分析限制在编码序列上，我们将 3D7 GFF 子集化为仅 CDS 特征，识别与目标相关的行，并将染色体、起始位置和终止位置输出为 BED (*.bed) 文件。 然后我们使用 hap.py 的 ––stratifications 标志来计算这些间隔上的度量。 我们使用 hap.py 生成的带注释的 VCF (*.vcf) 文件在 Python 中生成跨目标的假阳性率和真阳性率的位置图。

msp2 读取分析

计算编码序列长度

使用 minimap2 将读数映射到 3D7 参考基因组后72 (v2.24-r1122)，我们使用bedtools提取了与msp2 (PF3D7_0206800)编码序列完全重叠的reads74 (v2.31.0) 使用 intersect -F 1.0 命令。 从生成的 BAM 文件中，我们通过仅保留每个读取在间隔内对齐的部分，将这些读取修剪到 msp2 编码测序的范围 [273689, 274507] 2 号染色体 (Pf3D7_02_v3)； 如果插入缺失两侧的碱基在间隔内对齐，则插入缺失将被保留。 异常短的修剪读数（<–400bp）被作为可能的假象而被删除。 修剪的读数用于创建长度分布图。 它们是使用 minimap2 独立绘制的72 (v2.24-r1122)，发布从 PlasmoDB 下载的 3D7、Dd2、GB4 和 HB3 的 52 个参考基因组71。 我们让 minimap2 输出一个 PAF (*.paf) 文件，并通过将第 10 列（比对匹配数）除以第 11 列（总比对长度）来计算映射比对的身份。

全局成对比对

我们实现了 Needleman-Wunsch 算法的带状版本来计算修剪的 msp2 读段对之间的全局比对分数。 我们对评分模型进行参数化，以便分数反映两个观察到的读数源自相同的基础单倍型序列的对数概率； 即，所有比对差异都是由测序错误引起的。 假设插入缺失率为 5%，这与观察到的错误率大致一致，我们使用线性间隙分数 ({log }_{10}(0.05))。 对于替代分数，我们考虑了 Guppy 生成的碱基质量分数，如下所示。 将 x 和 y 定义为匹配中观察到的两个碱基，它们由相同单倍型碱基 h 生成的可能性为

$$P(x,y| h,{p}_{x},{p}_{y})=left{begin{array}{ll}(1-{p}_{x}) (1-{p}_{y})+frac{{p}_{x}{p}_{y}}{3}hfillquad &{{{{{{{rm{if} }}}}}}},x=y\ {p}_{x}(1-{p}_{y})+{p}_{y}(1-{p}_{x} )+frac{2{p}_{x}{p}_{y}}{3}quad &{{{{{{{rm{if}}}}}}}},x ne y\ quad end{array}right.$$

(1)

其中 ({p}_{x}=1{0}^{frac{{Q}_{x}}{-10}}) 和 ({p}_{y}=1{0} ^{frac{{Q}_{y}}{-10}})，其中 Qx 和 Qy 是 x 和 y 的 Phred 缩放碱基质量得分。 然后，替换分数计算为 log10(P(x, y–h,–px,–py))。 对于本研究中的所有比对，使用以全局比对矩阵对角线为中心的 80bp 带宽。 使用 SciPy 中的 scipy.cluster.hierachy.linkage 函数对结果分数进行层次聚类75 (v1.4.1) Python 包。

hrp2/3缺失检测的统计模型

数据和符号

我们首先描述我们的模型使用的数据和符号（表 2）。 扩增子测序数据表示为正整数 X 的二维矩阵，它保存质量过滤、解复用和映射后的读数计数。 矩阵 X 的行代表样本，其索引为 i ‾ {1,≫2,≫.≫.≫.≫,≫n}，列代表目标基因，其索引为 j – {1, – 2, –. –. –., – m}; 每个元素 xij 代表样本 i 映射到目标基因 j 的读数数量。 我们定义一个相应的二元矩阵 C，其中每个元素 cij 表示给定的目标基因 j 在样本 i 中存在（cij = 1）或删除（cij = 0）。 我们定义一个大小为 n 的向量 a，使得 ai ˆˆ [0, 1]，代表测序文库中每个样本的相对丰度。 最后，我们定义两个标量参数：读取错误分类率，ϵ ∈ [0, 1]，它表示从样本 k 导出的读数对另一个样本的读数计数（即 xi = k,j）的贡献率，无论是由于污染还是在解复用过程中错误的样本分配； 以及读取计数分散项 (nu in {mathbb{R}}+)，下面给出了精确的数学定义。

表2 hrp2/3缺失检测模型的符号模型

读取计数矩阵的每一列 xj=–{x1j,–x2j,–.–.–xnj} 包含给定目标基因的读取计数j 所有 n 个样本。 我们使用狄利克雷多项分布对 xj 进行建模：

$$P({{{{{{{{bf{x}}}}}}}}_{j};{N}_{j},{{{{{alpha}}}}} _ {j})=frac{{{Gamma }}(mathop{sum }nolimits_{i=1}^{n}{alpha }_{ij}){{Gamma }}({ N }_{j}+1)}{{{Gamma }}({N}_{j}+mathop{sum }limits_{i=1}^{n}{alpha }_{ij } )}mathop{prod }limits_{i=1}^{n}frac{{{Gamma }}({x}_{ij}+{alpha }_{ij})}{{ { 伽玛 }}({alpha }_{ij}){{伽玛 }}({x}_{ij}+1)},$$

(2)

其中 ({N}_{j}=mathop{sum }nolimits_{i=1}^{n}{x}_{ij}) 是所有样本中目标 j 的总读取计数； 且 αj=α{α1j,α2j,α.α.α.,αnj} 是复合向量目标基因 j 的参数 αij。 这些 αij 确定给定目标基因的每个样本的预期读取数量，并分三个步骤计算。 首先，我们使用样本 ai 的相对丰度及其目标基因 j 的删除状态 cij 来计算应源自样本 i 的目标基因 j 生成的读数的预期比例：

$${p}_{ij}=frac{{c}_{ij}{a}_{i}}{mathop{sum }limits_{k=1}^{n}{c}_ {kj}{a}_{k}}.$$

(3)

请注意，如果样本 i 中删除了目标基因 j，则 pij 要么等于 0，要么等于 0。 或测序文库中样本 i 的比例，但仅相对于存在该基因的样本。 在生成这些读数的过程中，这些预期比例会因读数错误分类和样本污染而改变，从而在最终数据中反映一组不同的预期比例 Ï€ij。 在这里，我们假设读取错误分类以固定速率 ϵ 发生，并且在样本中均匀发生，从而得出表达式：

$${pi}_{ij}={p}_{ij}(1-epsilon )+(1-{p}_{ij})epsilon 。$$

(4)

对于给定的 j，pij 和 Ï€ij 之和均为 1。 最后，我们将狄利克雷多项式参数化为 α ij = α β ij 。 其效果是，预期读取计数 xij 等于目标基因 j 的读取总数 Nj 乘以删除状态和误差调整样本比例 Ï€ij 的乘积：

$$E[{x}_{ij}]={N}_{j}frac{{alpha }_{ij}}{mathop{sum }nolimits_{k=1}^{n}{alpha }_{kj}}={N }_{j}frac{nu {pi }_{ij}}{mathop{sum }nolimits_{k=1}^{n}nu {pi }_{kj}}={ N}_{j}{pi }_{ij}.$$

(5)

xij 的方差等于：

$$V[{x}_{ij}]={N}_{j}{pi }_{ij}(1-{pi }_{ij})left(frac{{N}_{j}+nu }{1+nu }右），$$

(6)

随着 δ 向零收缩而增加，或随着 δ 向无穷大增长而渐近地接近二项式分布的方差； β 控制相对于二项式分布的读数计数过度分散。 总之，αij 包含有关目标基因的删除状态、文库中每个样本的相对丰度、测序运行中的错误分类率以及样本之间的读数计数过度分散量的信息。

推理

我们的目标是使用读数计数中的所有显着信息来推断给定样本中是否存在（cij = 1）或不存在（cij = 0）感兴趣的目标基因矩阵 X。我们用贝叶斯公式来处理这个问题：将如上所述的狄利克雷多项分布视为似然并计算 cj=”(c1j,”c2j,” 上的后验概率。 �.�.�,�cnj) 为：

$$P({{{{{{{{bf{c}}}}}}}}}_{j}| }_{j};epsilon,v,{{{{{{{bf {a}}}}}}}})propto P({{{{{{{{bf{x}}} }}}}}}_{j}|{{{{{{{bf {c}}}}}}}}_{j};epsilon,v,{{{{{{{ bf{a}}}}}}})P({{{{{{{{ bf{c}}}}}}}}}_{j}).$$

(7)

每个 cij 的先验的自然选择将是伯努利分布，cij ~ Bern(θ)，其中 1 — — θ 给出预期的删除概率。 在这里，为了简单起见，我们选择了一个相当于 α=0.5 的统一先验，尽管原则上这可以根据样本来源区域中目标基因 j 的缺失流行率的先前知识进行调整。集。

同样为了简单起见，我们选择将 ϵ、v 和 a 视为固定参数，并使用点估计来拟合它们。 令 δ 为包含测序运行中包含的所有阴性对照样本的索引的集合，使得 δδ 代表阴性对照的数量。 我们首先通过采用所有这些负控制的 xij 的经验平均值来对错误分类率进行简单的点估计：

$$epsilon=frac{1}{| delta | m}mathop{sum}limits_{iin delta }mathop{sum }limits_{j=1}^{m}{x}_{ij}/{N}_{j}。 $$

(8)

这利用了以下事实：E[xij]/Nj = ϵ 对于 ai = 0 的情况，根据阴性对照的定义，这是正确的。 接下来我们计算 a 和 δ 参数的点估计。 为了进行这些估计，我们定义了一个子集 Ï• � {1, 2, . . . , � m} 表示目标基因的索引没有已知的删除。 在 NOMADS16 panel 的背景下，这包括十个靶基因，不包括 hrp2、hrp2, up.、hrp2, down.、hrp3、hrp3, up. 和 hrp3, down。 狄利克雷或狄利克雷多项式参数的闭合形式最大似然估计不存在76，因此我们使用集合 Ï• 中目标基因的 xij/Nj 经验平均值来估计 ai：

$${a}_{i}=frac{1}{| phi | }mathop{sum}limits_{jin phi }{x}_{ij}/{N}_{j}.$$

(9)

然后，跟随 Minka (2012)76，我们使用以下方法估计 δ：

$$log (nu )=frac{1}{| phi | -1}mathop{sum}limits_{jin (| phi | -1)}log left(frac{{a}_{i}(1-{a}_{i}) {var({x}_{ij}/{N}_{j})}-1right).$$

(10)

通过 δ、δ 和 a 的点估计，我们使用马尔可夫链蒙特卡罗 (MCMC) 计算 cj 的后验分布。 对于每个目标基因 j，我们运行独立的 Metropolis-Hastings 算法来计算后验 cj。 我们用 cij=1 来初始化 MCMC，对于所有 i––{1,–2,–.–.–.,–n} 。 在每次迭代中，我们通过从 i 中均匀选择来提出一个新的 ({{{{{{{{bf{c}}}}}}}}}_{j}^{{prime} })，然后通过计算({c}_{ij}^{{prime} }=1-{c}_{ij})来切换对应cij的删除状态。 由于这个提案是对称的，黑斯廷斯比率为1，我们接受更新的概率为：

$$分钟左[1,frac{P({{{{{{{{bf{x}}}}}}}}}_{j};{N}_{j},{{{{{{{{bf{c}}}}}}}}}_{j}^{{prime} },epsilon,v,{{{{{{{bf{a}}}}}}}})P({{{{{{{{bf{c}}}}}}}}}_{j}^{{prime} })}{P({{{{{{{{bf{x}}}}}}}}}_{j};{N}_{j},{{{{{{{{bf{c}}}}}}}}}_{j},epsilon,v,{{{{{{{bf{a}}}}}}}})P({{{{{{{{bf{c}}}}}}}}}_{j})}right]$$

(11)

我们总共对每个目标基因进行了 10,000 次 MCMC 迭代，并丢弃前 500 次作为老化。 最后，给定样本 i 携带目标基因 j 缺失的后验概率等于 cij = 0 的迭代分数。

报告摘要

有关研究设计的更多信息，请参阅 自然投资组合报告摘要 链接到这篇文章。

1708026211
2024-02-15 19:30:48
#通过靶向纳米孔测序对恶性疟原虫疟疾进行灵活且经济有效的基因组监测