基于计算SNP熵值的混合样本鉴定方法与流程

基于计算snp熵值的混合样本鉴定方法
技术领域
1.本技术属于生信分析领域，具体涉及基于样本snp熵值的混合样本鉴定方法和系统。


背景技术：

2.单核苷酸多态性(singlenucleotidepolymorphism，snp)主要是指在基因组水平上由单个核苷酸的变异所引起的dna序列多态性。snp是最常见的遗传变异类型，一般表现为单个碱基的转换或颠换，也可能是碱基的插入或缺失，snp占人类基因组已知多态性的90％以上。snp多为双等位型标记,具有分布密度高,突变率低,位置不均匀等特点，具有较好的法医学和临床科学应用前景。
3.在法医和临床的dna检测过程中，经常会遇到多个体的混合样本，那么如何鉴定样本是否为混合样本？传统的鉴定混合样本的方法是通过某些位点包含3个等位基因来判断是混合样本，但该方法却至少存在以下缺点：
4.1)当样本的混合比率低于20％时，通常需要降低深度的判断标准来找到snp，判断标准不稳定，这样就要求检验人员有比较丰富经验；
5.2)为了最大化个体识别，通常选取的位点的次等位基因频率(maf)近似0.5，这样大多数个体都可能只有两个等位基因型。
6.有鉴于此，提出本技术，本技术提出一种直接从read出发，不需要判别snp分型即可判断混合样本的方法。


技术实现要素：

7.为解决上述技术问题，本技术提出如下具体技术方案：
8.本技术首先提供一种基于计算snp位点熵值的混合样本鉴定方法：包括如下步骤：
9.1)样本测序：测序样本文库构建，测序获得测序数据；
10.2)数据过滤：对测序数据进行序列过滤；
11.3)序列比对：对过滤后数据进行序列比对；
12.4)snp位点深度统计：对比对后序列进行snp位点碱基深度统计；
13.优选的，所述统计包括：snp染色体编号、在染色体上的位置、该位置总的覆盖度、在该位置测序到的a/c/g/t四种碱基的深度；
14.5)计算各snp位点熵值：利用snp位点碱基深度信息，计算snp位点各碱基的na(normalizedreadcount)、概率pa(probabilityofa)和熵值(entropyofonesnp)；
15.优选的，所述熵值越大，样本为混合样本可能性越高。
16.在一些方式中，所述步骤1)中，
17.所述文库构建过程中对每个样本加index；
18.所述测序数据根据index进行数据拆分。
19.在一些方式中，所述步骤2)中，
20.所述序列过滤为过滤低质量序列、短序列以及含n较多序列。
21.在一些方式中，所述步骤3)中，
22.所述序列比对采用bwa软件进行，得到比对的sam格式文件，具体的：第一步使用 bwa索引命令bwa index构建参考基因组索引；第二步，使用命令bwa-mem将序列比对到参考基因组，得到比对的sam格式文件；
23.在一些方式中，所述比对后进一步包括排序步骤：
24.将sam格式文件转换成bam格式文件，对bam文件进行排序，对排序后的bam文件建立索引。
25.在一些方式中，所述步骤4)中，
26.所述snp位点深度统计具体为：
27.准备snp位点的bed文件，所述bed文件内容为snp位点所在的染色体编号、在该染色体上的位置以及snp位点的rs编号，将排序后bam文件和snp位点bed文件作为输入，对各snp位点的碱基深度进行统计，统计文件内容包括snp染色体编号、在染色体上的位置、该位置总的覆盖度、在该位置测序到的a/c/g/t四种碱基的深度。
28.在一些方式中，所述步骤5)中，
29.所述na、pa和熵值的计算公式分别如下：
[0030][0031][0032]
entropy of one snp＝pa×
log(pa) pc×
log(pc) pg×
log(pc) p
t
×
log(p
t
)。
[0033]
在一些方式中，所述步骤5)中，
[0034]
所述计算还可包括：计算各样本snp位点的熵的平均值和熵的中位数，并标记样本的分组信息，得到关于样本、平均熵/中位数熵、样本分组的数据表格，对该数据进行可视化作图。
[0035]
另外，本技术还提供一种基于计算snp位点熵值的混合样本鉴定系统，所述模块用于执行上述任一所述方法的步骤。
[0036]
或者具体包括如下模块：
[0037]
1)样本测序模块：用于测序样本文库构建，测序获得测序数据；
[0038]
2)数据过滤模块：用于对测序数据进行序列过滤；
[0039]
3)序列比对模块：用于对过滤后数据进行序列比对；
[0040]
4)snp位点深度统计模块：用于对比对后序列进行snp位点碱基深度统计；
[0041]
优选的，所述统计包括：snp染色体编号、在染色体上的位置、该位置总的覆盖度、在该位置测序到的a/c/g/t四种碱基的深度；
[0042]
5)计算各snp位点熵值模块：用于利用snp位点碱基深度信息，计算snp位点各碱基的na(normalized read count)、概率pa(probability of a)和熵值(entropy ofone snp)；优选的，所述熵值越大，样本为混合样本可能性越高。
[0043]
在一些方式中，所述模块1)中，
[0044]
所述文库构建过程中对每个样本加index；
[0045]
所述测序数据根据index进行数据拆分。
[0046]
在一些方式中，所述模块2)中，
[0047]
所述序列过滤为过滤低质量序列、短序列以及含n较多序列。
[0048]
在一些方式中，所述模块3)中，
[0049]
所述序列比对采用bwa软件进行，得到比对的sam格式文件，具体的：第一步使用 bwa索引命令bwa index构建参考基因组索引；第二步，使用命令bwa-mem将序列比对到参考基因组，得到比对的sam格式文件；
[0050]
在一些方式中，所述比对后进一步包括排序步骤：
[0051]
将sam格式文件转换成bam格式文件，对bam文件进行排序，对排序后的bam文件建立索引。
[0052]
在一些方式中，所述模块4)中，
[0053]
所述snp位点深度统计具体为：
[0054]
准备snp位点的bed文件，所述bed文件内容为snp位点所在的染色体编号、在该染色体上的位置以及snp位点的rs编号，将排序后bam文件和snp位点bed文件作为输入，对各snp位点的碱基深度进行统计，统计文件内容包括snp染色体编号、在染色体上的位置、该位置总的覆盖度、在该位置测序到的a/c/g/t四种碱基的深度。
[0055]
在一些方式中，所述模块5)中，
[0056]
所述na、pa和熵值的计算公式分别如下：
[0057][0058][0059]
entropy of one snp＝pa×
log(pa) pc×
log(pc) pg×
log(pg) p
t
×
log(p
t
)。
[0060]
在一些方式中，所述步骤5)中，
[0061]
所述计算还可包括：计算各样本snp位点的熵的平均值和熵的中位数，并标记样本的分组信息，得到关于样本、平均熵/中位数熵、样本分组的数据表格，对该数据进行可视化作图。
[0062]
另外，本技术还提供一种计算机可读介质，其存储有计算机程序，所述计算机程序被处理器执行时，实现上述任一所述方法。
[0063]
另外，本技术还提供一种电子设备，包括处理器以及存储器，所述存储器上存储一条或多条可读指令，所述一条或多条可读指令被所述处理器执行时，实现上述任一所述方法。
[0064]
与现有技术相比，本技术至少具有如下优势：
[0065]
1)本技术直接根据snp位点的各碱基深度情况，即可计算其熵值，不需要分析 snp的分型。
[0066]
2)本技术不局限于样本的混合比例问题，当混合比例低于20％时，不需要用降低深度的判断标准来找snp。多次测试表明本技术的方法能够在1:19(5％)区分混合样本。
[0067]
3)本技术不局限于位点的maf值，maf接近于0.5也不受影响。本技术选用了 230个maf接近0.5的snp位点进行了多次测试，测试结果表明本技术的方法能够很好的区分混合
样本。
附图说明
[0068]
为了更清楚地说明本技术具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本技术的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0069]
图1、本技术流程图；
[0070]
图2、针对9947和9948不同比例混合，使用本技术的方法计算熵值，其中的横坐标代表不同比例混合的样本，前两个样本为纯合样本，后面是不同混合比例的样本，由于有生物学重复，每一种混合比例都有超过一个样本，纵坐标为计算的熵值的中位数。
[0071]
图3、针对样本a和样本a不同比例混合、样本b和样本b不同比例混合以及样本 a和样本b不同比例混合，使用本技术的方法计算熵值，其中的横坐标代表不同比例混合的样本，前两个样本为纯合样本，后面是不同混合比例的样本，由于有生物学重复，每一种混合比例都有超过一个样本，纵坐标为计算的熵值的中位数。
具体实施方式
[0072]
下面将结合附图对本技术的技术方案进行清楚、完整地描述，显然，所描述的实施例是本技术一部分实施例，而不是全部的实施例。基于本技术中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本技术保护的范围。
[0073]
以下术语或定义仅仅是为了帮助理解本技术而提供。这些定义不应被理解为具有小于本领域技术人员所理解的范围。
[0074]
除非在下文中另有定义，本技术具体实施方式中所用的所有技术术语和科学术语的含义意图与本领域技术人员通常所理解的相同。虽然相信以下术语对于本领域技术人员很好理解，但仍然阐述以下定义以更好地解释本技术。
[0075]
如本技术中所使用，术语“包括”、“包含”、“具有”、“含有”或“涉及”为包含性的(inclusive)或开放式的，且不排除其它未列举的元素或方法步骤。术语“由
…
组成”被认为是术语“包含”的优选实施方案。如果在下文中某一组被定义为包含至少一定数目的实施方案，这也应被理解为揭示了一个优选地仅由这些实施方案组成的组。
[0076]
在提及单数形式名词时使用的不定冠词或定冠词例如“一个”或“一种”，“所述”，包括该名词的复数形式。
[0077]
本技术中的术语“大约”、“大体”表示本领域技术人员能够理解的仍可保证论及特征的技术效果的准确度区间。该术语通常表示偏离指示数值的
±
10％，优选
±
5％。
[0078]
此外，说明书和权利要求书中的术语第一、第二、第三、(a)、(b)、(c)以及诸如此类，是用于区分相似的元素，不是描述顺序或时间次序必须的。应理解，如此应用的术语在适当的环境下可互换，并且本技术描述的实施方案能以不同于本技术描述或举例说明的其它顺序实施。
[0079]
本技术所述的基于计算snp位点熵值的区分或混合样本鉴定方法的流程如图1所示，大体包括了样本测序，数据过滤和比对，snp位点深度统计和计算各snp位点熵值等步
骤。
[0080]
示例性的，所述方法包括如下步骤：
[0081]
1)样本测序：测序样本文库构建，测序获得测序数据；2)数据过滤：对测序数据进行序列过滤；3)序列比对：对过滤后数据进行序列比对；4)snp位点深度统计：对比对后序列进行snp位点碱基深度统计；优选的，所述统计包括：snp染色体编号、在染色体上的位置、该位置总的覆盖度、在该位置测序到的a/c/g/t四种碱基的深度；5) 计算各snp位点熵值：利用snp位点碱基深度信息，计算snp位点各碱基的na(normali zed read count)、概率pa(probability of a)和熵值(entropy of one snp)；实践中，所述熵值越大，样本为混合样本可能性越高。
[0082]
在一些实施方式中，所述步骤1)中，所述文库构建过程中对每个样本加index；所述测序数据根据index进行数据拆分。
[0083]
在另一些实施方式中，所述步骤2)中，所述序列过滤为过滤低质量序列、短序列以及含n较多序列。
[0084]
在另一些实施方式中，所述步骤3)中，所述序列比对采用bwa软件进行，得到比对的sam格式文件，具体的：第一步使用bwa索引命令bwa index构建参考基因组索引；第二步，使用命令bwa-mem将序列比对到参考基因组，得到比对的sam格式文件；
[0085]
在另一些实施方式中，所述比对后进一步包括排序步骤：将sam格式文件转换成b am格式文件，对bam文件进行排序，对排序后的bam文件建立索引。
[0086]
在另一些实施方式中，所述步骤4)中，所述snp位点深度统计具体为：准备snp 位点的bed文件，所述bed文件内容为snp位点所在的染色体编号、在该染色体上的位置以及snp位点的rs编号，将排序后bam文件和snp位点bed文件作为输入，对各sn p位点的碱基深度进行统计，统计文件内容包括snp染色体编号、在染色体上的位置、该位置总的覆盖度、在该位置测序到的a/c/g/t四种碱基的深度。
[0087]
在一些另实施方式中，所述步骤5)中，所述na、pa和熵值的计算公式分别如下：
[0088][0089][0090]
entropy of one snp＝pa×
log(pa) pc×
log(pc) pg×
log(pg) p
t
×
log(p
t
)。
[0091]
在另一些实施方式中，所述步骤5)中，所述计算还可包括：计算各样本snp位点的熵的平均值和熵的中位数，并标记样本的分组信息，得到关于样本、平均熵/中位数熵、样本分组的数据表格，对该数据进行可视化作图。
[0092]
本技术方法对于测序数据不作过多要求，可适用于各种测序手段或测序仪器获得的测序数据，比如第二代测序或第三代测序，优选的为第二代测序。
[0093]
下面为具体的实施方法。
[0094]
实施例1本技术方法体系的构建和优化过程
[0095]
本技术整体分析的流程如图1所示。流程上，首先样本提取建库，测序获取测序数据，对测序数据进行拆分、过滤和比对，随后统计snp位点深度，最后计算各snp位点的熵值，具体如下：
[0096]
1、测序：样本提取，文库建库(文库构建过程中会给每个样本加上唯一的index，通过index来区分识别每一个样品)，上机进行高通量测序，测序数据下机最初为bcl 格式文件，准备samplesheet列表，列表中记录了样本与index间的对应关系，采用b cl2fastq软件根据样本的index进行数据拆分，得到各样本的fastq格式文件，即分析用原始数据。
[0097]
2、数据过滤：采用fastp软件对原始fastq数据做过滤，fastp软件能自动识别接头序列并进行裁剪，过滤低质量序列、太短的序列以及含n较多的序列，最终得到过滤后数据即clean data。
[0098]
3、序列比对及其优化：
[0099]
考虑到bwa软件是一种能够将差异度较小的序列比对到一个较大的参考基因组上的软件包，本实施例选用软件bwa做序列比对。bwa比对过程主要分为两步：第一步使用索引命令bwa index构建参考基因组的索引；第二步比对，但bwa有三种比对算法， bwa-backtrack是用来比对illumina的序列的，reads长度最长能到100bp；bwa-sw和 bwa-mem主要是用于比对长reads，支持的长度为70bp-1mbp，同时支持剪接性比对，但是bwa-mem比对运行更快，结果更加准确。因此本实施例使用bwa mem命令将序列比对到参考基因组上，得到比对的sam格式文件，其效果最优。
[0100]
而为了减少文件的存储，使用samtools view命令将sam格式转换成bam格式文件，bam文件是sam格式的二进制格式。接着用sambamba sort对bam文件进行排序。最后对排序后的bam文件用sambamba index命令建立索引，因为整个bam文件可能非常大，如果我们只关注很小的一段区域而将整个序列都读进内存是非常低效的，建立索引则方便针对性的提取特定区域。
[0101]
4、snp位点深度统计：首先准备关于snp位点的bed文件，文件内容为snp位点所在的染色体编号、在该染色体上的位置以及snp位点的rs编号，然后将排序后的bam 文件和snp位点的bed文件作为输入，用sambamba depth base命令对各snp位点的碱基深度进行统计，统计文件的主要内容为snp染色体编号、在染色体上的位置、该位置总的覆盖度、在该位置测序到的a/c/g/t四种碱基的深度。
[0102]
5、计算各snp位点的熵值：利用snp位点碱基深度信息，先计算该snp位点各碱基的normalized read count，再计算各碱基的概率，最后计算得到该snp位点的熵值，公式如下：
[0103][0104][0105]
entropy of one snp＝pa×
log(pa) pc×
log(pc) pg×
log(pg) p
t
×
log(p
t
)
[0106]
一个分组有多个生物学重复，计算各样本的snp位点的熵的平均值和熵的中位数，并标记样本的分组信息，得到关于样本、平均熵/中位数熵、样本分组的数据表格，使用 r语言的ggplot2包对该数据进行可视化作图，若熵值越大，样本为混合样本的可能性越高。
[0107]
实施例2本技术方法体系评估——使用标准品进行不同样本混合比例的测试
[0108]
一、测试样本准备：
[0109]
1、纯合样本：样本9947和样本9948，(样本9947和9948是法医学的标准品，其中样
本9947为女性样本，样本9948为男性样本)。
[0110]
2、混合样本：将两个纯合样本按不同比例混合得到不同比例的混合样本，分别为 1:19mix，1:14mix，1:9mix、1:5mix、1:2mix、1:1mix、2:1mix、5:1mix、 9:1mix、14:1mix，19:1mix，(例如1:9mix表示9947与9948样本按照1:9混合)。以上每个样品均需要至少做2-3个生物学重复。
[0111]
二、采用实施例1方法进行混合样本分析
[0112]
数据分析：样本经测序，得到原始下机bcl文件，用bcl2fastq软件做数据拆分，得到各样品的原始fastq数据文件；将原始数据用fastp过滤；过滤后数据与人类参考基因组进行比对，得到sam比对结果文件，将sam文件转换成二进制格式的bam文件，对bam文件进行排序并对排序后的bam文件建立索引；根据bam文件和对应的snp的b ed文件(共230个snp)，统计各snp位点的碱基深度，为了最大化个体识别的能力，选取的这230个snp的maf都接近0.5；根据碱基深度文件，先计算各样本的各snp位点的熵值，再计算每个样本的snp位点的平均熵值或者中位数熵值，纯样本和混合样本均有多个生物学重复，相同类型的样本标记为同一分组即index，例如：编号9-11样品均标记为1:9mix，得到一个关于熵值的矩阵，其内容为：样本编号、平均熵值/熵的中位数、分组index编号，分别作图(见图2)，如下所示：
[0113]
图2中横坐标为样本分组index，同一比例混合的样本在一个组中，图中的每一个点代表一个样本，纵坐标为样本的熵的中位数。熵值为0.07-0.25左右的时候可以都清晰的区分纯和和混合样本。
[0114]
由图2可知，纯合样本的熵值(中位数)明显低于不同比例混合样本的熵值，且样品混合比例越均衡，其熵值越高。熵值为0.47左右的时候仍可以区分分纯和和绝大部分混合样本。
[0115]
实施例3本技术方法体系评估——使用不同样本测试方法的可靠性
[0116]
一、测试样本准备：
[0117]
1、纯合样本：样本b1、样本b2；样本a1、样本a2。
[0118]
2、混合样本：将两个纯合样本按不同比例混合得到不同比例的混合样本，分别为1: 19mix，1:14mix，1:9mix、1:5mix、1:2mix、1:1mix、2:1mix、5:1mix、9:1 mix、14:1mix，19:1mix，(例如1:9mix表示样本1与样本2样本按照1:9混合)。
[0119]
二、一共进行了三组混合，包括样本a和样本a混合、样本b和样本b混合以及样本a和样本b混合。采用实施例1方法进行混合样本分析，结果如图3所示，纯合样本具有较小的熵值(小于0.1)，而混合样本从19:1(5％)到1:19(5％)都具有更大的熵值 (基本都大于0.3)，可见，混合比例远低于20％的时候本技术提出的方法也能够很好的区分混合样本。
[0120]
前述对本技术的具体示例性实施方案的描述是为了说明和例证的目的。这些描述并非想将本技术限定为所公开的精确形式，并且很显然，根据上述教导，可以进行很多改变和变化。对示例性实施例进行选择和描述的目的在于解释本技术的特定原理及其实际应用，从而使得本领域的技术人员能够实现并利用本技术的各种不同的示例性实施方案以及各种不同的选择和改变。本技术的范围意在由权利要求书及其等同形式所限定。

技术特征：
1.一种基于计算snp位点熵值的混合样本鉴定方法，其特征在于，包括如下步骤：1)样本测序：测序样本文库构建，测序获得测序数据；2)数据过滤：对测序数据进行序列过滤；3)序列比对：对过滤后数据进行序列比对；4)snp位点深度统计：对比对后序列进行snp位点碱基深度统计；优选的，所述统计包括：snp染色体编号、在染色体上的位置、该位置总的覆盖度、在该位置测序到的a/c/g/t四种碱基的深度；5)计算各snp位点熵值：利用snp位点碱基深度信息，计算snp位点各碱基的n
a
(normalized read count)、概率p
a
(probability of a)和熵值(entropy of one snp)；优选的，所述熵值越大，样本为混合样本可能性越高。2.根据权利要求1所述的鉴定方法，其特征在于，步骤1)中，所述文库构建过程中对每个样本加index；所述测序数据根据index进行数据拆分。3.根据权利要求1所述的鉴定方法，其特征在于，步骤2)中，所述序列过滤为过滤低质量序列、短序列以及含n较多序列。4.根据权利要求1所述的鉴定方法，其特征在于，步骤3)中，所述序列比对采用bwa软件进行，得到比对的sam格式文件，具体的：第一步使用bwa索引命令bwa index构建参考基因组索引；第二步，使用命令bwa-mem将序列比对到参考基因组，得到比对的sam格式文件；优选的，所述比对后进一步包括排序步骤：将sam格式文件转换成bam格式文件，对bam文件进行排序，对排序后的bam文件建立索引。5.根据权利要求1所述的鉴定方法，其特征在于，步骤4)中，所述snp位点深度统计具体为：准备snp位点的bed文件，所述bed文件内容为snp位点所在的染色体编号、在该染色体上的位置以及snp位点的rs编号，将排序后bam文件和snp位点bed文件作为输入，对各snp位点的碱基深度进行统计，统计文件内容包括snp染色体编号、在染色体上的位置、该位置总的覆盖度、在该位置测序到的a/c/g/t四种碱基的深度。6.根据权利要求1所述的鉴定方法，其特征在于，步骤5)中，所述n
a
、p
a
和熵值的计算公式分别如下：和熵值的计算公式分别如下：entropy of one snp＝p
a
×
log(p
a
) p
c
×
log(p
c
) p
g
×
log(p
g
) p
t
×
log(p
t
)。7.根据权利要求1所述的鉴定方法，其特征在于，步骤5)中，所述计算还可包括：计算各样本snp位点的熵的平均值和熵的中位数，并标记样本的分组信息，得到关于样本、平均熵/中位数熵、样本分组的数据表格，对该数据进行可视化作图。
8.一种基于计算snp位点熵值的混合样本鉴定系统，其特征在于，包括如下模块，所述模块用于执行权利要求1-7任一所述方法的步骤。9.一种计算机可读介质，其存储有计算机程序，所述计算机程序被处理器执行时，实现权利要求1-7任一所述方法。10.一种电子设备，包括处理器以及存储器，所述存储器上存储一条或多条可读指令，所述一条或多条可读指令被所述处理器执行时，实现权利要求1-7任一所述方法。

技术总结
本申请涉及生物信息学分析领域，具体提供一种基于样本SNP熵值的混合样本鉴定方法和系统，该方法和系统不需分析SNP分型，直接根据SNP位点各碱基深度情况计算熵值；本申请具有不局限于样本混合比例以及不局限于位点MAF值等优势。等优势。等优势。


技术研发人员：李梦 黄舒 郭茂平 申君毅 郭晋荣 胡欢 郑立 张奇 陈初光
受保护的技术使用者：北京阅微基因技术股份有限公司
技术研发日：2022.03.18
技术公布日：2022/5/25