一种低深度测序群体基因型填充计算内存优化方法技术

本发明专利技术公开了一种低深度测序群体基因型填充计算内存优化方法，其包括：对单核苷酸多态性位点SNP进行分块，依据分块来设定检查点值；依据设定的检查点值按分块进行计算前向辅助变量和后向辅助变量，并计算单倍型观测条件概率。本发明专利技术具有原理简单、可有效降低计算内存消耗、提升基因型填充效率等优点。提升基因型填充效率等优点。提升基因型填充效率等优点。

【技术实现步骤摘要】
一种低深度测序群体基因型填充计算内存优化方法


[0001]本专利技术主要涉及到基因测序及生物信息分析
，特指一种低深度测序群体基因型填充计算内存优化方法。

技术介绍

[0002]随着测序成本的降低，大规模群体低深度全基因组测序逐渐成为群体遗传研究中不可或缺的技术手段。通过对大量样本进行低深度测序分型，可以以相同甚至更低成本获得比传统芯片分型和简化基因组测序更密集的分子标记，从而提高GWAS定位的功效。然而，由于测序深度非常低，该方法在技术上不可避免地造成了大量单核苷酸多态性位点(SNP)的基因型缺失，进而影响了下游遗传图谱构建、QTL定位和群体GWAS研究的准确性，因此可见从群体层面提升基因型的完整度具有十分重要的意义。
[0003]群体基因型填充的基本原理是：通过已知的单倍型参考数据集或群体大量样本之间的共享单倍型信息，采用统计学方法对样本中的缺失基因型进行推断。
[0004]现有通过已知的单倍型参考数据进行基因型填充的步骤为：
[0005]1)对目标样本进行单倍型定相(phasing)，从而得到染色体区域内的单倍型信息；
[0006]2)通过参考单倍型集合对缺失基因型进行填充。
[0007]上述方法需要预建参考单倍型数据集，这对于一些单倍型图谱研究不是十分充分的物种来说，该数据集通常无法获取，进而无法使用该方法。另外，当目标群体与构建参考单倍型群体存在差异时，基因型填充的准确性也会受到影响。2016年，有从业者Davies等研究出不依赖参考单倍型而借助大量个体的测序数据中共享读段信息对祖先单倍型进行建模的STITCH方法，进而实现了对每个个体缺失基因型的推断，解决了无参考单倍型数据的物种基因型填充问题。
[0008]另有STITCH方法通过一阶隐马尔可夫模型(HMM)对群体中的祖源单倍型进行建模，将染色体视作K个未知祖源单倍型的嵌合体(对于远交群体而言，可以认为这些单倍型集合捕获了一个区域内的所有不同单倍型的完整信息，因此K可能很大)。该方法通过EM算法迭代更新HMM模型的参数，E步对每个样本计算其祖源单倍型观测概率，M步通过多样本单倍型的观测概率对祖源单倍型频率和其他HMM参数进行更新。从计算上来讲，该算法单次迭代的时间复杂度与样本数N、SNP数目T和祖源单倍型数目K2呈线性关系，其计算瓶颈在于：E步对每个样本计算其祖源单倍型观测概率。
[0009]E步对每个样本计算祖源单倍型观测概率可以归类为HMM模型三个基本问题中的评估问题，即给定模型λ＝(Λ，B，π)和观测序列O＝o1，o2，...，o
T
，计算条件概率p{O|λ}。如果对该问题直接进行概率运算，其时间复杂度为{K2}
T
，显然是不可行的，解决该问题常见的方法是采用前向
‑
后向算法，该前向
‑
后向算法通过引入辅助变量，有效降低了时间复杂度，其计算步骤如下：
[0010]1、计算前向辅助变量α
t
(i)：
[0011]前向辅助变量α
t
(i)定义为部分观测序列o1，o2，...o
t
在状态i中止的概率，数学上
表示为α
t
(i)＝p{o1，o2，...o
t
，q
t
＝i|λ}，则有：
[0012][0013]其中α1(i)＝π
i
b
i
(o1)，1≤i≤K2，b
i
(o
t+1
)为发射概率，α
ij
为隐状态转移概率，π
i
为状态i的初始概率，j表示第j种祖源单倍型组合，t表示第t个SNP位点，应用此递推公式，可以计算出α
T
(i)，1≤i≤K2；
[0014]2、计算后向辅助变量β
t
(i)：
[0015]后向辅助变量β
t
(i)定义为部分观测序列o
t+1
，o
t+2
，...o
T
以状态i为起始的概率，数学上表述为p{o
t+1
，o
t+2
，...o
T
|q
t
＝i，λ}，同理计算β
t
(i)的递推公式为：
[0016][0017]其中β
T
(i)＝α
T
(i)，1≤i≤K2，b
j
(o
t+1
)为发射概率，α
ij
为隐状态转移概率，j表示第j种祖源单倍型组合，t表示第t个SNP位点。
[0018]则有α
t
(i)β
t
(i)＝p{O，q
t
＝i|λ}，1≤i≤K2，1≤t≤T，于是有：
[0019][0020]3、计算上述概率过程中引入的辅助变量α
t
(i)和β
t
(i)虽然可以有效降低计算的时间复杂度，但是对计算机内存的要求非常高。通过简单计算可知，对于一个包含3百万SNP位点，祖源单倍型数目为30的染色体，单个α
t
(i)矩阵的存储将耗费20G内存，单个样本在E步消耗的内存将达到60GB以上，这使得该方法在实际应用时由于内存受限无法并行处理多个样本，进而无法充分利用众核服务器的CPU计算资源。
[0021]由以上背景可知，单个样品辅助变量矩阵的内存消耗与SNP数目T和祖源单倍型数目K2的乘积成正比，考虑到单倍型区域的局部性，有从业者提出一种降低内存的方法，其是将SNP集合划分为互不包含的连续的N个子集，每个子集含有相同的SNP数目；然后分别对每个SNP子集应用上述算法，该方法没有考虑到计算辅助变量矩阵时前后列之间的数据依赖关系，进而影响了基因型填充的准确性。
[0022]另有从业者提出的一种改进的方式是在划分SNP集合时，允许相邻两个集合重叠固定长度的染色体区间，这虽然在一定程度上改善了上一方法的缺陷，但是由于染色体不同区域，不同群体的连锁不平衡程度不同，人为选取的固定长度区间仍然不可避免会带来填充精度的损失。

技术实现思路

[0023]本专利技术要解决的技术问题就在于：针对现有技术存在的技术问题，本专利技术提供一种原理简单、可有效降低计算内存消耗、提升基因型填充效率的低深度测序群体基因型填充计算内存优化方法。
[0024]为解决上述技术问题，本专利技术采用以下技术方案：
[0025]一种低深度测序群体基因型填充计算内存优化方法，其包括：
[0026]对单核苷酸多态性位点SNP进行分块，依据分块来设定检查点值；
[0027]依据设定的检查点值按分块进行计算前向辅助变量和后向辅助变量，并计算单倍型观测条件概率。
[0028]作为本专利技术方法的进一步改进：利用第一轮前向辅助变量和后向辅助变量算法来进行检查点的确定；在检查点值的基本文档来自技高网...

【技术保护点】

【技术特征摘要】
1.一种低深度测序群体基因型填充计算内存优化方法，其特征在于，包括：对单核苷酸多态性位点SNP进行分块，依据分块来设定检查点值；依据设定的检查点值按分块进行计算前向辅助变量和后向辅助变量，并计算单倍型观测条件概率。2.根据权利要求1所述的低深度测序群体基因型填充计算内存优化方法，其特征在于，利用第一轮前向辅助变量和后向辅助变量算法来进行检查点的确定；在检查点值的基础上利用第二轮前向辅助变量和后向辅助变量算法计算块内的前向辅助变量和后向辅助变量。3.根据权利要求2所述的低深度测序群体基因型填充计算内存优化方法，其特征在于，在依据设定的检查点值按分块进行计算前向辅助变量和后向辅助变量过程中采用分块并行化的处理方式。4.根据权利要求2
‑
3中任意一项所述的低深度测序群体基因型填充计算内存优化方法，其特征在于，对单核苷酸多态性位点SNP进行分块的流程包括：记单核苷酸多态性位点SNP个数为T，祖源单倍型个数为K，则单个前向辅助变量矩阵的尺寸为T*K2；记分块大小为M，则分块个数为存储检查点所需向量尺寸为单个分块的前向辅助变量的尺寸为M*K2。5.根据权利要求4所述的低深度测序群体基因型填充计算内存优化方法，其特征在于，当将检查点辅助变量值全部放入内存，单个分块的前向辅助变量消耗内存为6.根据权利要求4所述的低深度测序群体基因型填充计算内存优化方法，其特征在于，当时，峰值内存为最小，其中M
*
为最优分块大小；当将辅助变量检查点向量存储至磁盘，根据主机内存调整SNP分块大小以适配主机内存容量。7.根据权利要求4所述的低深度测序群体基因型填充计算内存优化方法，其特征在于，利用第一轮算法来进行前向辅助变量的检查点值的确定流程包括：递推计算前向辅助变量α
t
(i)；前向辅助变量α
t
(i)定义为部分观测序列o1，o2，...o
t
在状态i中止的概率，数学上表示为α
t
(i)＝p{o1，o2，...o
t
，q
t
＝i|λ}，则有：其中，α1(i)＝π
i
b
i
(o1)，1≤i≤K2，b
i
(o
t+1
)为发射概率，α
ij
为隐状态转移概率，π
i
为状态i的初始概率，j表示第j种祖源单倍型组合，t表示第t个SNP位点，应用此递推公式，计算出α
T
(i)，1≤i≤K2；获取检查...

【专利技术属性】
技术研发人员：蒋艳凰，马丑贤，王振国，杨仁武，毛海波，黄立磊，冯博伦，
申请(专利权)人：人和未来生物科技长沙有限公司，
类型：发明
国别省市：