CN116564405B - 一种基于平均无序度的基因组测序突变位点过滤方法 - Google Patents

Abstract

本申请属于生信分析技术领域，具体涉及一种基于平均无序度的基因组测序突变位点过滤算法，可用于突变位点假阳过滤。

Description

一种基于平均无序度的基因组测序突变位点过滤方法

技术领域

本申请属于生物信息学技术领域，具体涉及一种基于平均无序度的基因组测序突变位点过滤方法及应用。

技术背景

基因组测序技术目前已经在肿瘤诊断项目中广泛开展。通过对癌组织进行基因组测序，诊断者可以明确癌细胞的突变类型，进行针对性用药。但测序中发现的突变存在非生物来源，如仪器误差、文库制备污染、扩增阶段错配等，属于假性突变，需要加以过滤。

假性突变通常具有一定的特征，即支持突变的序列上存在大量无序突变。目前业界广泛采用支持突变序列上的错配数进行过滤，如Vardict软件的NM指标通过计算支持突变的序列上的平均错配数(含插入/缺失)，用于突变位点的过滤。但是此指标只能说明支持突变序列上错配多少，并不能代表是否混乱；另外对于多核苷酸多态性位点(MNP),由于该类突变是由相邻的一组单核苷算多态位点(SNP)组成，错配数会增加导致NM统计量升高，会将真实的突变误判成假阳。

鉴于此，提出本申请。

发明内容

为解决上述技术问题，本申请通过生物信息学分析研究，建立一套基于平均无序度的基因组测序突变位点过滤方法，该方法可用于突变位点假阳过滤。

具体的，本申请提出如下技术方案：

本申请首先提供一种基于平均无序度的基因组测序突变位点过滤方法，所述方法是通过平均无序度指标对突变位点进行判定。

进一步的，所述平均无序度指标是基于信息熵构建。

进一步的，所述方法包括如下步骤，

1)基于测序数据获得突变位点相关信息；

2)建立突变位点附近的碱基数量矩阵；

3)计算碱基数量矩阵中碱基出现频率，利用信息熵计算位点无序度：

4)统计碱基数量矩阵中突变位点覆盖深度，计算位点的平均无序度；

5)基于无序度阈值进行突变过滤。

进一步的，所述步骤1)中，所述突变位点相关信息包括突变位点所在染色体、位置坐标、突变后的碱基与支持突变的序列数量。

进一步的，所述步骤2)中，所述碱基数量矩阵具有4行和2n+2列；每列分别表示指定观测范围内的基因组坐标位点，其中第n列表示突变位点所在位置，n的取值是使得2n+2接近序列比对文件中的序列平均长度；每行分别表示每个基因组位点上4种碱基被覆盖的次数。

进一步的，所述步骤3)中，所述计算为分别计算碱基数量矩阵的每一列4种碱基出现的频率P_A、P_C、P_G、P_T，利用信息熵计算位点无序度DSOi：

DSOi＝-P_Alog₂(P_A)-P_Tlog₂(P_T)-_Clog₂(P_C)-P_Glog₂(P_G)。

进一步的，所述步骤4)中，所述统计为对碱基数量矩阵的每列分别统计4种碱基出现次数的总和，得到该位点覆盖深度；

进一步的，优选的，对覆盖深度大于等于矩阵第n列突变位点覆盖深度50％的位点进行计数；将覆盖深度不足突变位点50％的位点的无序度DSOpos置为0。

进一步的，所述步骤4)中，所述根据覆盖深度计算位点的平均无序度具体为：对突变位点附近覆盖深度大于等于突变位点本身覆盖深度一半的位点，计算其平均无序度DSO，所述DSO计算如下：

进一步的，所述步骤5)具体为：对实际样本进行审核，确定无序度阈值(阈值优选为1.0)，大于阈值的位点被认为支持突变的序列混乱，用于后续突变过滤。

本申请还提供一种平均无序度指标在测序突变位点假阳性过滤中的应用，

优选的，所述平均无序度指标是基于信息熵构建；

更优选的，所述应用具体包括如下步骤：

1)基于测序数据获得突变位点相关信息；

2)建立突变位点附近的碱基数量矩阵；

3)计算碱基数量矩阵的碱基出现频率，利用信息熵计算位点无序度：

4)统计碱基数量矩阵中突变位点的覆盖深度，计算位点的平均无序度；

5)基于无序度阈值进行突变过滤。

进一步的，所述步骤1)中，所述突变位点相关信息包括突变位点所在染色体、位置坐标、突变后的碱基与支持突变的序列数量。

进一步的，所述步骤2)中，所述碱基数量矩阵具有4行和2n+2列；每列分别表示指定观测范围内的基因组坐标位点，其中第n列表示突变位点所在位置，n的取值是使得2n+2接近序列比对文件中的序列平均长度；每行分别表示每个基因组位点上4种碱基被覆盖的次数。

进一步的，所述步骤3)中，所述计算为分别计算碱基数量矩阵的每一列4种碱基出现的频率P_A、P_C、P_G、P_T，利用信息熵计算位点无序度DSOi：

DSOi＝-P_Alog₂(P_A)-P_Tlog₂(P_T)-_Clog₂(P_C)-P_Glog₂(P_G)。

进一步的，所述步骤4)中，所述统计为对碱基数量矩阵的每列分别统计4种碱基出现次数的总和，得到该位点覆盖深度；

进一步的，优选的，对覆盖深度大于等于矩阵第n列突变位点覆盖深度50％的位点进行计数，得到计数值；将覆盖深度不足突变位点50％的位点的无序度DSOpos置为0。

进一步的，所述步骤4)中，所述根据覆盖深度计算位点的平均无序度具体为：对突变位点附近覆盖深度大于等于突变位点本身覆盖深度一半的位点，计算其平均无序度DSO，所述DSO计算如下：

进一步的，所述步骤5)具体为：对实际样本进行审核，确定无序度阈值(阈值优选为1.0)，大于阈值的位点被认为支持突变的序列混乱，用于后续突变过滤。

本申请还提供一种计算机存储介质，所述计算机存储介质存储有计算机程序，所述计算机程序包括程序指令，所述程序指令当被处理器执行时，实现上述任一项所述方法。

本申请还提供一种电子设备，所述电子设备包括存储器和处理器，所述存储器和处理器相连，所述存储器用于存储计算机程序，所述处理器用于调用所述计算机程序，实现上述任一项所述方法。

本申请有益技术效果：

本申请通过探索发现，基于信息熵构建的平均无序度指标，可用于描述位点附近的序列混乱程度，可对突变位点的真实性进行准确判定，可对假阳性位点进行有效过滤。

附图说明

图1.本申请方法的流程图；

图2.不同VD阈值下DSO在样本中的密度分布曲线；

图3.平均错配数较高的有序突变形态示例，基因组数据可视化软件IGV(Integrative Genomics Viewer)界面截图。图中每条横向的条块表示一条测序序列，条块的一端为箭头状，表示序列测序的方向(与箭头指向的方向一致)；灰色条块是正常序列，红色条块表示序列的插入片段长度过大，其他颜色条块表示双端测序中与本序列成对的另一条序列被比对到其他染色体上。条块中出现的字母表示在本序列的相应位置相对基因组出现突变，无字母显示的区域表示序列与基因组碱基一致。字母颜色的深浅代表了碱基的质量分数，越深则该碱基质量越高。部分附图由于放缩比例过小的缘故，字母无法看清，IGV会自动将碱基字母调整为相应的色块(比如A为绿色、C为蓝色、G为黄色、T为红色、N为灰色)以便查看。

图4.平均错配数较低的无序突变形态示例，基因组数据可视化软件IGV(Integrative Genomics Viewer)界面截图。图中每条横向的条块表示一条测序序列，条块的一端为箭头状，表示序列测序的方向(与箭头指向的方向一致)；灰色条块是正常序列，红色条块表示序列的插入片段长度过大，其他颜色条块表示双端测序中与本序列成对的另一条序列被比对到其他染色体上。条块中出现的字母表示在本序列的相应位置相对基因组出现突变，无字母显示的区域表示序列与基因组碱基一致。字母颜色的深浅代表了碱基的质量分数，越深则该碱基质量越高。部分附图由于放缩比例过小的缘故，字母无法看清，IGV会自动将碱基字母调整为相应的色块(比如A为绿色、C为蓝色、G为黄色、T为红色、N为灰色)以便查看。

图5.实际样本中的无序突变，基因组数据可视化软件IGV(Integrative GenomicsViewer)界面截图。图中每条横向的条块表示一条测序序列，条块的一端为箭头状，表示序列测序的方向(与箭头指向的方向一致)；灰色条块是正常序列，红色条块表示序列的插入片段长度过大，其他颜色条块表示双端测序中与本序列成对的另一条序列被比对到其他染色体上。条块中出现的字母表示在本序列的相应位置相对基因组出现突变，无字母显示的区域表示序列与基因组碱基一致。字母颜色的深浅代表了碱基的质量分数，越深则该碱基质量越高。部分附图由于放缩比例过小的缘故，字母无法看清，IGV会自动将碱基字母调整为相应的色块(比如A为绿色、C为蓝色、G为黄色、T为红色、N为灰色)以便查看。

图6.实际样本中的有序突变，基因组数据可视化软件IGV(Integrative GenomicsViewer)界面截图，图中每条横向的条块表示一条测序序列，条块的一端为箭头状，表示序列测序的方向(与箭头指向的方向一致)；灰色条块是正常序列，红色条块表示序列的插入片段长度过大，其他颜色条块表示双端测序中与本序列成对的另一条序列被比对到其他染色体上。条块中出现的字母表示在本序列的相应位置相对基因组出现突变，无字母显示的区域表示序列与基因组碱基一致。字母颜色的深浅代表了碱基的质量分数，越深则该碱基质量越高。部分附图由于放缩比例过小的缘故，字母无法看清，IGV会自动将碱基字母调整为相应的色块(比如A为绿色、C为蓝色、G为黄色、T为红色、N为灰色)以便查看。

具体实施方式

下面将结合实施例对本申请的实施方案进行详细描述，但是本领域技术人员将会理解，下列实施例仅用于说明本申请，而不应视为限制本申请的范围。实施例中未注明具体条件者，按照常规条件或制造商建议的条件进行。所用试剂或仪器未注明生产厂商者，均为可以通过市场购买获得的常规产品。

部分术语定义，除非在下文中另有定义，本申请具体实施方式中所用的所有技术术语和科学术语的含义意图与本领域技术人员通常所理解的相同。虽然相信以下术语对于本领域技术人员很好理解，但仍然阐述以下定义以更好地解释本申请。

本申请中的术语“大约”表示本领域技术人员能够理解的仍可保证论及特征的技术效果的准确度区间。该术语通常表示偏离指示数值的±10％，优选±5％。

如本申请中所使用，术语“包括”、“包含”、“具有”、“含有”或“涉及”为包含性的(inclusive)或开放式的，且不排除其它未列举的元素或方法步骤。术语“由…组成”被认为是术语“包含”的优选实施方案。如果在下文中某一组被定义为包含至少一定数目的实施方案，这也应被理解为揭示了一个优选地仅由这些实施方案组成的组。

此外，说明书和权利要求书中的术语第一、第二、第三、(a)、(b)、(c)以及诸如此类，是用于区分相似的元素，不是描述顺序或时间次序必须的。应理解，如此应用的术语在适当的环境下可互换，并且本申请描述的实施方案能以不同于本申请描述或举例说明的其它顺序实施。

本申请的基于平均无序度的基因组测序突变位点过滤方法是通过平均无序度指标对突变位点进行判定。

在一些实施方式中，所述平均无序度指标是基于信息熵构建。

在一些实施方式中，所述方法包括如下步骤，

1)基于测序数据获得突变位点相关信息；

2)建立突变位点附近的碱基数量矩阵；

3)计算碱基数量矩阵的碱基出现频率，利用信息熵计算位点无序度：

4)统计碱基数量矩阵中突变位点的覆盖深度；

5)根据覆盖深度计算位点的平均无序度；

6)基于无序度阈值进行突变过滤。

在一些实施方式中，所述步骤1)中，所述突变位点相关信息包括突变位点所在染色体、位置坐标、突变后的碱基与支持突变的序列数量。

在一些实施方式中，所述步骤2)中，所述碱基数量矩阵具有4行和2n+2列；每列分别表示指定观测范围内的基因组坐标位点，其中第n列表示突变位点所在位置，n的取值是使得2n+2接近序列比对文件中的序列平均长度；每行分别表示每个基因组位点上4种碱基被覆盖的次数。

在一些实施方式中，所述步骤3)中，所述计算为分别计算碱基数量矩阵的每一列4种碱基出现的频率P_A、P_C、P_G、P_T，利用信息熵计算位点无序度DSOi：

DSOi＝-P_Alog₂(P_A)-P_Tlog₂(P_T)-P_Clog₂(P_C)-P_Glog₂(P_G)。

在一些实施方式中，所述步骤4)中，所述统计为对碱基数量矩阵的每列分别统计4种碱基出现次数的总和，得到该位点覆盖深度；

在一些优选的实施方式中，对覆盖深度大于等于矩阵第n列突变位点覆盖深度50％的位点进行计数；将覆盖深度不足突变位点50％的位点的无序度DSOpos置为0。

在一些实施方式中，所述步骤5)中，所述根据覆盖深度计算位点的平均无序度具体为：对突变位点附近覆盖深度大于等于突变位点本身覆盖深度一半的位点，计算其平均无序度DSO，所述DSO计算如下：

在一些实施方式中，所述步骤6)具体为：对实际样本进行审核，确定无序度阈值(阈值优选为1.0)，大于阈值的位点被认为支持突变的序列混乱，用于后续突变过滤。

下面结合具体实施例来阐述本申请。

实施例1本申请方法体系建立

本申请通过生信分析探索建立一套基于平均无序度的基因组测序突变位点过滤算法，具体包括了如下步骤：

算法主要逻辑描述如下：

1.阅读来自上游突变位点分析程序输出的变异分析结果文件中记录的突变位点，确认突变位点所在的染色体名chr、位置坐标pos、突变后的碱基alt与支持突变的序列数量vd。当多线程运行时，每个线程自队列中取出变异分析结果分别处理。

2.建立突变位点附近的碱基数量矩阵，矩阵具有4个行和2n+2个列，每列分别表示指定观测范围内的基因组坐标位点，其中第n列表示突变位点所在的位置，n的值应当使得2n+2接近序列比对文件中的序列平均长度；每行分别表示每个基因组位点上4种碱基被覆盖到的次数。

3.对步骤2中所述碱基数量矩阵的每一列分别计算4种碱基出现的频率P_A、P_C、P_G、P_T，然后利用信息熵的定义计算位点无序度：

DSOi＝-P_Alog₂(P_A)-P_Tlog₂(P_T)-P_Clog₂(P_C)-P_Glog₂(P_G)

4.为了消除覆盖深度不足对统计结果造成的影响，对步骤2中所述碱基数量矩阵的每列分别计算4种碱基出现次数的总和，即该位点的覆盖深度，对覆盖深度大于等于矩阵第n列突变位点覆盖深度50％的位点进行计数，得到计数值COUNT；同时将覆盖深度不足突变位点50％的位点的无序度DSOpos置为0。

5.计算突变位点附近覆盖深度大于等于突变位点本身覆盖深度一半的位点的平均无序度：

6.存储时将该值放大50倍以避免过多的小数位数，默认值为0。

7.对实际样本进行审核，确定无序度的阈值(在本发明中设定为1.0)，大于阈值的位点被认为支持突变的序列混乱，用于后续的突变过滤。

实施例2DSO阈值划定

由于DSO是一个统计学量，支持其计算的原始数据(即支持突变的read条数VD)过少会导致统计受到异常值的强烈影响而出现不可靠的问题，故而需要在实际样本中寻找使得DSO计算相对稳定的VD阈值。

本申请通过设定不同的VD阈值，对某样本中的全部突变计算平均无序度DSO，并调查DSO在整个样本中的分布情况(只调查非零的DSO值)。

随着VD阈值的提升，DSO的分布应当逐渐趋于稳定。为了确保统计量的稳定性，本申请以2为间隔，等间隔加大VD阈值，选取DSO分布密度峰值与前值小于0.01后的第一个VD阈值作为实际应用的VD阈值。

由图2，VD阈值被定为8。

根据VD＝8时DSO的分布情况，本申请选取密度分布曲线最高点之后的第一个整数作为DSO的阈值划定标准。由此DSO的阈值被设定为1，DSO高于1的位点被认定为序列过乱。

实施例3本申请DSO与常规指标NM的效果比较

本实施例基于实际病例样本，测试无序度指标DSO与常规错配数指标NM对突变过滤的效果。突变的真假一般基于IGV截图确认，其中突变位点周围的混乱程度是判定突变真假其中一个指标，常用的混乱程度描述指标NM划定阈值为3，NM>＝3认为是混乱假位点，NM<3为有序真实突变，测试表明本申请的DSO方法在多种突变类型中表现优于NM。。

1)如图3所示，IGV复核为一个真实突变，而NM指标计算出为3(NM＝3)，被判定成为一个混乱假阳位点；而基于本申请的DSO指标计算值为0(DSO＝0)，小于域值1，表明是有序的真阳位点，这与IGV复核结果一致；

2)如图4所示，IGV复核为一个假阳突变，而NM指标计算值为1.78(NM＝1.78)，被判定为有序真突变；而基于DSO指标计算值为1.662(DSO＝1.662)，大于域值1，因此是混乱假阳突变，这与IGV复核结果一致。

实施例4本申请突变位点无序度分辨能力评估

本实施例在实际病例样本中测试了共计34个报出突变。其中，图5为测试中3个突变位点DSO指标判定为混乱假阳位点举例，图6为测试中3个DSO指标判定为有序高质量真阳位点举例。可见，本申请DSO指标对于突变位点无序度分辨结果与IGV截图复核结果一致。

最后应说明的是：以上各实施例仅用以说明本申请的技术方案，而非对其限制；尽管参照前述各实施例对本申请进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本申请各实施例技术方案的范围。

Claims (8)

1.一种基因组测序突变位点的过滤方法，其特征在于，所述方法通过平均无序度指标对突变位点进行判定；所述平均无序度指标是基于信息熵构建；

所述方法具体包括如下步骤：

1)基于测序数据获得突变位点相关信息；

2)建立突变位点附近的碱基数量矩阵；

3)计算碱基数量矩阵中碱基出现频率，基于信息熵计算位点无序度：

4)统计碱基数量矩阵中突变位点的覆盖深度，计算位点的平均无序度；

5)基于无序度阈值进行突变过滤。

2.根据权利要求1所述的过滤方法，其特征在于，所述步骤1)中，所述突变位点相关信息包括突变位点所在染色体、位置坐标、突变后的碱基与支持突变的序列数量。

3.根据权利要求1所述的过滤方法，其特征在于，所述步骤2)中，所述碱基数量矩阵具有4行和2n+2列；每列分别表示指定观测范围内的基因组坐标位点，其中第n列表示突变位点所在位置，n的取值是使得2n+2接近序列比对文件中的序列平均长度；每行分别表示每个基因组位点上4种碱基被覆盖的次数。

4.根据权利要求1所述的过滤方法，其特征在于，所述步骤3)中，所述计算为分别计算碱基数量矩阵的每一列4种碱基出现的频率P_A、P_C、P_G、P_T，利用信息熵计算位点无序度DSOi：

DSOi＝-P_Alog₂(P_A)-P_Tlog₂(P_T)-P_Clog₂(P_C)-P_Glog₂(P_G)。

5.根据权利要求1所述的过滤方法，其特征在于，所述步骤4)中，所述统计为对碱基数量矩阵的每列分别统计4种碱基出现次数的总和，得到该位点覆盖深度。

6.根据权利要求1所述的过滤方法，其特征在于，所述步骤4)中，根据覆盖深度计算位点的平均无序度具体为：对突变位点附近覆盖深度大于等于突变位点本身覆盖深度一半的位点，计算其平均无序度DSO，所述DSO计算如下：

7.根据权利要求1所述的过滤方法，其特征在于，所述步骤5)具体为：对实际样本进行审核，确定无序度阈值，大于阈值的位点被认为支持突变的序列混乱，用于后续突变过滤。

8.基于平均无序度指标在测序突变位点假阳性过滤中的应用；所述平均无序度指标是基于信息熵构建，所述应用具体包括如下步骤：

1)基于测序数据获得突变位点相关信息；

2)建立突变位点附近的碱基数量矩阵；

3)计算碱基数量矩阵中碱基出现频率，基于信息熵计算位点无序度：

4)统计碱基数量矩阵中突变位点的覆盖深度，计算位点的平均无序度；

5)基于无序度阈值进行突变过滤。