CN115019886B - 一种宏基因组绝对定量实验全流程的数字孪生方法 - Google Patents

Abstract

本发明公开了一种宏基因组绝对定量实验全流程的数字孪生方法，包括三个模块：构建样本库、构建内标库、宏基因组分析流程。本发明通过计算机模拟的方式展示了宏基因组绝对定量的全流程，使新手能够更快更方便地理解宏基因组绝对定量的流程，并可以对比内标法宏基因组绝对定量各因素之间的影响。

Description

一种宏基因组绝对定量实验全流程的数字孪生方法

技术领域

本发明属于计算机模拟技术领域，具体涉及一种宏基因组绝对定量实验全流程的数字孪生方法。

背景技术

宏基因组学是为了解决传统方法分离纯化复杂微生物群落难的痛点，从基因组水平来研究微生物群落。将一个微生物群落的DNA全部提取出来，随机打断、测序，返回大量的短DNA片段。通过宏基因组学技术不仅可以得到微生物群落中的物种结构，还可以进行功能分析，重构一些难以培养的微生物的基因组。

目前宏基因组研究中得到微生物含量的方式主要有两种，分别是相对定量和绝对定量，但前者只能看到微生物的相对变化情况，并不能提供物种丰度变化的程度和方向的信息，忽略了微生物总体变化情况。绝对定量分析则是计量样本每种微生物的拷贝数或者数目，从而实现绝对定量，相对于常规相对定量扩增子测序，绝对定量分析能反映样本每种微生物的真实数量和组间样本的真实差异，因此绝对定量分析相对于相对定量分析更能反映样本细菌群落的真实变化，是进行微生态研究的首选。所以，微生物的绝对定量是非常必要的。

而且在内标法测宏基因组绝对定量中，内标加入的时间、多少都会影响到绝对定量的结果，本发明通过提供在提取前加入内标和在提取后加入内标两种加入内标的方式以及三种不同梯度浓度的内标来对比以上因素对内标法绝对定量的影响。

此外，现实中的宏基因组测序需要购买菌株、测序、计算分析，每一步的费用都很高，但对于计算机模拟来说，菌株DNA可以通过已有库模拟生成，测序也可以通过加入误差因子的方式来模拟现实结果，最后通过结果分析出微生物的绝对定量。计算机模拟的方式成本几乎为0，而且每一步的结果都是可观测的，便于新人学习和理解宏基因组绝对定量流程，对实际的宏基因组实验也有一定的实际意义。

发明内容

为了克服现有技术的不足，本发明提供了一种宏基因组绝对定量实验全流程的数字孪生方法，包括三个模块：构建样本库、构建内标库、宏基因组分析流程。本发明通过计算机模拟的方式展示了宏基因组绝对定量的全流程，使新手能够更快更方便地理解宏基因组绝对定量的流程，并可以对比内标法宏基因组绝对定量各因素之间的影响。

本发明解决其技术问题所采用的技术方案包括如下模块：

一种宏基因组绝对定量实验全流程的数字孪生方法，包括如下模块：

模块1：构建样本库；

设定菌群的相对定量信息，包括每种菌的名称代号、各种菌之间的比例关系和DNA序列的总数量；各种菌的比例和为1；

通过现有的微生物信息DNA库，根据各种菌之间的比例关系和DNA序列的总数量，确定整个菌群的绝对定量信息，模拟用户需要的菌群环境；

DNA序列的总数量乘以菌株的相对定量为该菌株的绝对定量，即需要从现有微生物信息DNA库中提取的菌株DNA的数量；最终将所有菌株提取完毕，生成用户需要的菌群环境，即样本库；

模块2：构建内标库；

内标库的生成有两种方式：一种是用户自定义内标，需要用户自定义输入菌群的相对定量信息和DNA序列的总数量，并且用户设定的菌株需要满足与样本库中相似度低于80％的条件；另一种是系统自动生成内标，不需要用户额外输入，系统自动根据样本库挑选相似度低于80％的菌群，生成内标库；内标库生成的总量为样本库的 0.1％、1％、10％三种内标浓度供用户选择；

采用两种方式插入内标：一种是在步骤1中提取菌株DNA前插入，加入内标后一起提取，内标会和DNA一起提取，这种插入方式只在后续测序过程产生误差；第二种方式是在步骤1中提取菌株DNA后插入内标DNA，由于DNA会在提取过程中产生误差，而内标DNA没有这部分误差，此时会在提取和测序中产生两段误差；

模块3：宏基因组分析流程；

本模块分为四个步骤；

步骤1：提取样本；

在提取样本时，本模块会随机生成95％-100％之间保留两位小数的百分数，作为本次提取的回收率；样本提取后将DNA序列随机切断成长度为300bp-500bp之间的片段，随机切断后的DNA库作为待测序文库以作测序使用；

步骤2：测序；

对待测序文库中的DNA序列进行模拟测序，测序是将生物信息中的“ACGT”四种碱基转换为电子信息“ACGT”四种字母，在本步骤中，直接将步骤1中的DNA数据作为测序的结果；

采用梯度方式模拟测序误差，具体实现如下：测序长度为150bp，前100bp的误差设置为一个固定的值0.01％，随着测序的后移，误差呈梯度上升的，在最后达到峰值；梯度上升过程是非线性上升的，波动范围为0.5％-2％；

步骤3：计算；

通过下式计算微生物群落的绝对定量的归一化因子：

其中，n为插入的内标DNA总数，此项为用户输入；C_s,i为第i个内标DNA的浓度，此项在生成内标DNA库时得到；Z_s,i为第i个内标DNA的reads数，此项在生成内标 DNA库时得到；L_s,i为第i个内标DNA的的碱基长度，此项为固定长度150bp，如果在提取前加入内标则是打断后的长度；

步骤4：分析对比；

归一化因子乘以未知样本DNA的Z_t/L_t，即得到未知样本DNA的浓度；Z_t、L_t分别表示样本DNA t的reads数和碱基长度，都是在测序后得到；然后用未知样本DNA 的相对定量乘样本总量得到未知样本DNA的绝对定量，将结果与模块1中生成DNA 库时的绝对定量进行对比，观察内标法绝对定量的精确度；

通过对比两种内标插入方式和三种内标浓度来比较内标法的效果。

本发明的有益效果如下：

本发明提供了一种对宏基因组绝对定量实验全流程的数字孪生方法，它能够通过现有数据库中的DNA信息，按照用户输入的相对定量信息来生成用户需要的DNA库，作为样本DNA库或者内标DNA库；并完成与现实实验相对应的流程：提取、插入、测序、计算和分析。通过计算机模拟的方式展示了宏基因组绝对定量的全流程，使新手能够更快更方便地理解宏基因组绝对定量的流程。并对内标法绝对定量实验进行了模拟，如：在提取DNA前和提取DNA后分别插入，并测试两种插入方式对内标法准确度的影响；以一种梯度浓度的方式来插入内标，观察何种浓度的内标法效果更佳，具有一定的现实意义。

附图说明

图1为本发明方法总体流程图。

图2为本发明实施例的软件结构图。

具体实施方式

下面结合附图和实施例对本发明进一步说明。

本发明的目的在于个性化生成用户需要的菌群环境、提供较低成本的模拟测序流程以及对比内标法宏基因组绝对定量方案优劣，通过计算机程序模拟宏基因组中对微生物群落进行绝对定量实验的流程，总体由三个主要模块来实现本发明的功能。

一种宏基因组绝对定量实验全流程的数字孪生方法，包括如下模块：

模块1：构建样本库；

设定菌群的相对定量信息(包括的每种菌的名称代号和他们之间的比例，要求比例和为1)和总体DNA序列的数量；

通过现有的微生物信息DNA库，模拟出用户需要的菌群环境，各种菌之间的比例确定，总的DNA序列数也确定，这样就确定了整个菌群的绝对定量信息，这些信息可以作为模块3中的分析绝对定量的标准答案，以此作为评分标准来判定何种内标法的准确度更优。

总DNA序列数量乘以某菌株的相对定量，即为该菌株的绝对定量，也就是需要从现有DNA库中提取的菌株DNA的数量，最终将所有菌株提取完毕，生成用户需要的菌群环境。

本发明此模块的优势在于，用户可以根据自己的需求来自定义样本库，本发明使用现有DNA库来模拟生成满足用户需求的菌群环境，由于采用的都是真实的DNA数据，所以更加真实可靠；且在当下菌株价格居高不下的环境中，节省了大量的成本。

模块2：构建内标库；

内标库的生成有两种方式，一种是用户自定义内标，需要用户自定义输入菌群的相对定量信息和DNA序列的总数量，并且用户设定的菌株需要满足与样本库中相似度低于80％的条件。另一种是系统自动生成内标，这种方法不需要用户额外输入，系统会自动根据样本库挑选相似度低于80％的菌群，生成内标库。内标库生成的的总量为样本DNA库的0.1％、1％、10％，三种不同梯度供用户选择。

本发明采用两种方式插入内标，一种是在DNA提取前插入，加入内标后一起提取，内标会和样本DNA一起提取，该情况中只在后续测序过程产生误差；第二种方式是在提取DNA后插入内标DNA，由于样本DNA会在提取过程中有一部分误差，而内标DNA没有这部分误差，此时会在提取和测序中产生两段误差。

模块3：宏基因组分析流程；

本模块分为四个步骤；

步骤1：提取样本；

在实际实验的提取过程中，尤其是在试剂盒裂解、纯化和洗脱的过程中，会导致DNA的消耗，业界会对提取的DNA进行质控，以回收率来判定DNA提取的程度；一般业界中较好的回收率都在95％-100％这个区间内，为了与现实更接近，本发明也采用回收率95％-100％这个范围，在提取时，本模块会随机生成95％-100％之间保留两位小数的百分数，作为本次提取的回收率。提取后需要将DNA序列随机切断成长度为300bp-500bp之间的片段，随机打断后的DNA库以作测序使用；

步骤2：测序；

对待测序文库中的DNA序列进行模拟测序，测序的本质就是将生物信息中的“ACGT”四种碱基转换为电子信息“ACGT”四种字母，由于本发明中从始至终都是电子模拟的，所以不需要实际实验中的荧光读取等。直接步骤1中的数据作为测序的结果，但实际实验中的测序是有少量误差的。

在测序过程中，测序仪通过检测释放的荧光信号强度达到实时测定DNA序列的目的。然而由于化学反应的特殊性，以及一些技术原因，通常检测的光信号存在抖动起伏，进而导致了测序误差。而这种误差情况多出现在测序的后半部分，且越靠后，误差的概率越大。基于此，本发明需要采用梯度的方式模拟这种误差，具体实现如下：测序长度为150bp，前100bp的误差很低，设置为一个固定的值(0.01％)，然后随着测序的后移，越后面出现误差的概率越大，这个误差是呈梯度上升的，在最后达到峰值，但又不是线性上升的，会有一个参数波动范围(0.5％-2％)，但宏观上还是梯度上升，这样模拟的更加真实。

步骤3：计算；

通过下式计算微生物群落的绝对定量的归一化因子：

其中，n为插入的内标DNA总数，此项为用户输入；C_s,i为第i个内标DNA的浓度，此项在生成内标DNA库时得到；Z_s,i为第i个内标DNA的reads数，此项在生成内标DNA库时得到；L_s,i为第i个内标DNA的的碱基长度，此项为固定长度150bp，如果在提取前加入内标则是打断后的长度；

步骤4：分析对比；

归一化因子乘于未知样本DNA的Z_t/L_t，即可得到未知样本DNA的浓度；Z_t、L_t分别表示样本DNA t的reads数和碱基长度，都是在测序后得到；然后用未知样本DNA 的相对定量乘样本总量得到未知样本DNA的绝对定量，将结果与模块1中生成DNA 库时的绝对定量进行对比，观察内标法绝对定量的精确度。通过对比两种内标插入方式和三种内标浓度来比较内标法的何种方式效果更佳。

具体实施例：

本实施例使用Python程序语言作为编程语言，图1中的构建样本库模块中，需要用户输入各种菌之间的相对定量和DNA序列的总数目。现有的基因基准库中有大量的物种DNA数据，每种都有自己的编号，用户需要输入编号，以及编号所代表的物种的相对丰度。如：Al001 21.32％,Ah003 12.25％。物种编号后空格隔开，后面的百分数是该物种所占的比例，即相对定量。每个物种之间的相对定量信息之间用逗号隔开。然后本发明将从基因基准库中随机选取满足用户要求的DNA序列，且不会随机到重复序列。提取出的DNA序列，会以字符串的形式存储在一个txt文件中。程序会一直循环生成这些DNA序列，直到数量达到了用户输入的数量。且由于用户输入了相对定量信息和总的DNA数量，所以可以得到构建出的DNA库的绝对定量信息，该信息可作为模块3的判定标准。本发明模块1默认构建的是样本DNA库，如果用户需要构建内标库，则需要通过模块3调用模块1来构建。所有生成的样本DNA序列都存储在SampleDNA文件夹中。从已有数据库中提取DNA来建库会使数据更接近现实，有较高的现实意义。

图1中的构建内标库模块中，内标库的生成有两种方式，一种是用户自定义内标，需要用户自定义输入菌群的相对定量信息和DNA序列的总数量，并且用户设定的菌株需要满足与样本库中相似度低于80％的条件，此处的相似度对比是通过BLAST对比，如果用户输入的信息不符合该条件，则会提示用户其输入的哪种菌不符合。另一种是系统自动生成内标，这种方法不需要用户额外输入，系统会自动根据样本库挑选相似度低于80％的菌群，生成内标库。内标库有三种梯度浓度供用户选择，分别是0.1％、 1％、10％相对样本浓度。内标库建立完成后，用户可选择两种插入模式，一种是在样本提取前插入，一种是在样本提取后测序前插入；

图1中的模拟宏基因组分析流程模块中，我们从模块1中已经得到了原始的DNA 信息，现在需要先提取出来并制备文库，然后才能测序。首先是提取阶段，由于本发明提取DNA只需要直接读取即可，不存在损耗，而实际实验中提取过程中的试剂盒裂解、DNA纯化和洗脱过程中会使DNA有不同程度的损耗。所以本发明为了弥补与真实条件之间的差异，会采用回收率为95％-100％这个闭区间之间的随机保留两位小数的百分数。具体实施则是使用python中的random库进行随机化操作，且本发明中所有的随机化操作均是采用该库。随后则是对已经制备好的文库进行测序，本发明会将序列长度较大的序列随机打断，最终使得95％的DNA序列都已经在500bp以下即可。测序的目的是将现实中DNA的核苷酸转化为可视的数据，但本发明属于数字孪生，本来就是数据化的模拟，所以不需要进行测序，直接就可以得到数据，这样与真实的测序操作相比缺少了测序产生的误差，进而在后续的计算中会与真实实验中的结果有一些差距。为了减少与真实实验的差距，本发明引入了一种梯度的随机误差参数，通过实际测序过程中测序精确度递减这一现象，来模拟测序中的误差，使结果更加真实。具体则是通过函数sequencingDeviation()来完成，该函数会读取已经打断后的文库，并根据DNA序列的长度来决定如何模拟误差，测序的读长是固定的150bp，前100bp 的误差很低，设置为一个固定的值(0.01％)，主要误差将集中在后50bp上。整体误差趋势呈一种波动上升的趋势，设定总体趋势为方程y＝kx+c，y就是某一位碱基产生突变、丢失的概率(突变与丢失的概率比为3：7)，k为0.8，c为0.01；且波动参数会在0.5％-2％闭区间内随机生成。在计算中，通过前面步骤插入的内标DNA数据可以知道内标基因i的浓度为C_s,i，内标基因i的reads数为Z_s,i，内标基因i的碱基长度： L_s,i。然后计算出绝对定量中的归一化因子：归一化因子与其他未知基因的Z_t/L_t相乘即可得到基因t的相对浓度，乘于样本总量即可得到基因t的绝对丰度。最后与模块1中的标准答案进行分析对比，观察内标法两种插入方式、三种梯度浓度中，哪种效果最佳。

Claims (1)

1.一种宏基因组绝对定量实验全流程的数字孪生方法，其特征在于，包括如下模块：

模块1：构建样本库；

设定菌群的相对定量信息，包括每种菌的名称代号、各种菌之间的比例关系和DNA序列的总数量；各种菌的比例和为1；

通过现有的微生物信息DNA库，根据各种菌之间的比例关系和DNA序列的总数量，确定整个菌群的绝对定量信息，模拟用户需要的菌群环境；

DNA序列的总数量乘以菌株的相对定量为该菌株的绝对定量，即需要从现有微生物信息DNA库中提取的菌株DNA的数量；最终将所有菌株提取完毕，生成用户需要的菌群环境，即样本库；

模块2：构建内标库；

内标库的生成有两种方式：一种是用户自定义内标，需要用户自定义输入菌群的相对定量信息和DNA序列的总数量，并且用户设定的菌株需要满足与样本库中相似度低于80％的条件；另一种是系统自动生成内标，不需要用户额外输入，系统自动根据样本库挑选相似度低于80％的菌群，生成内标库；内标库生成的总量为样本库的0.1％、1％、10％三种内标浓度供用户选择；

采用两种方式插入内标：一种是在步骤1中提取菌株DNA前插入，加入内标后一起提取，内标会和DNA一起提取，这种插入方式只在后续测序过程产生误差；第二种方式是在步骤1中提取菌株DNA后插入内标DNA，由于DNA会在提取过程中产生误差，而内标DNA没有这部分误差，此时会在提取和测序中产生两段误差；

模块3：宏基因组分析流程；

本模块分为四个步骤；

步骤1：提取样本；

在提取样本时，本模块会随机生成95％-100％之间保留两位小数的百分数，作为本次提取的回收率；样本提取后将DNA序列随机切断成长度为300bp-500bp之间的片段，随机切断后的DNA库作为待测序文库以作测序使用；

步骤2：测序；

对待测序文库中的DNA序列进行模拟测序，测序是将生物信息中的“ACGT”四种碱基转换为电子信息“ACGT”四种字母，在本步骤中，直接将步骤1中的DNA数据作为测序的结果；

采用梯度方式模拟测序误差，具体实现如下：测序长度为150bp，前100bp的误差设置为一个固定的值0.01％，随着测序的后移，误差呈梯度上升的，在最后达到峰值；梯度上升过程是非线性上升的，波动范围为0.5％-2％；

步骤3：计算；

通过下式计算微生物群落的绝对定量的归一化因子：

其中，n为插入的内标DNA总数，此项为用户输入；C_s,i为第i个内标DNA的浓度，此项在生成内标DNA库时得到；Z_s,i为第i个内标DNA的reads数，此项在生成内标DNA库时得到；L_s,i为第i个内标DNA的的碱基长度，此项为固定长度150bp，如果在提取前加入内标则是打断后的长度；

步骤4：分析对比；

归一化因子乘以未知样本DNA的Z_t/L_t，即得到未知样本DNA的浓度；Z_t、L_t分别表示样本DNAt的reads数和碱基长度，都是在测序后得到；然后用未知样本DNA的相对定量乘样本总量得到未知样本DNA的绝对定量，将结果与模块1中生成DNA库时的绝对定量进行对比，观察内标法绝对定量的精确度；

通过对比两种内标插入方式和三种内标浓度来比较内标法的效果。