DUGMO：利用纯细菌样本的高通量测序数据检测未知转基因生物的工具

摘要

背景

欧洲共同体对转基因生物的传播和使用采取了非常严格的政策。事实上，“无转基因”标签允许的最大污染阈值为0.9%。近年来，检测到进口了未描述的转基因生物。它们的序列没有描述，因此无法通过常规方法（如PCR）检测到。

结果

我们基于Illumina配对测序数据开发了DUGMO，这是一种生物信息学管道，用于检测转基因细菌，包括未知的转基因细菌。该方法目前主要用于检测纯细菌样本中可能含有部分转基因的转基因细菌。在初步步骤中，编码序列（CDS）通过两个连续的BLASTN与宿主全基因组进行比对，并使用相关的调谐参数来区分属于野生型基因组（wgCDS）的CDS和潜在的GM编码序列（pgmCDS）。然后，根据基因组词汇的差异，计算wgCDS和每个pgmCDS之间的Bray-Curtis距离。最后，基于Bray-Curtis距离和GC内容等六个变量，采用两种机器学习方法，即随机森林和广义线性模型，以目标为真正的GM CDS。对GM进行的测试枯草芽孢杆菌显示25个阳性的CDSs对应于氯霉素抗性基因和插入质粒的CDSs。在一个野生类型上枯草芽孢杆菌，未检测到假阳性序列。

结论

DUGMO在高通量测序数据中检测外源CDS、截断、融合或高度突变的野生CDS，并被证明在检测GM序列方面是有效的，但它也可以用于鉴定最近的水平基因转移。

背景

转基因生物（GMO）是一种“通过交配和/或自然重组以非自然方式改变遗传物质的生物”（欧洲委员会，2001年）[1]. 我们称宿主基因组为野生的原始基因组，对其进行基因改造以构建GMO。在这篇文章中，我们重点关注外源基因整合到宿主基因组中产生的转基因生物。欧洲立法要求，当转基因食品中的成分含量超过0.9%时，应报告其在食品中的存在[2]. 该法规要求所有转基因生物或来自转基因生物的产品必须向欧洲当局申报，其核酸序列的修改必须记录在适当的数据库中，并且必须提供基于标准分子技术（例如PCR）的检测方法。GMOMETHODS公司[三]，尤金[4]和GMO-Amplicon[5]数据库列出了欧洲市场上授权的转基因生物及其相关检测方法。此外，Morisset等人，2014年[6]开发了一个名为GMOseek的软件程序，以优化实验室中的常规转基因检测，从而帮助在分析的所有阶段做出决策。最近，随着高通量全基因组测序（WGS）的发展[1,7,8,9,10]基于比对相似性的测序数据集生物信息学分析[1,7,8]已开发用于描述引用的GMO元素。

Willems等人，2016年[11]提出了一个基于下一代测序（NGS）的统计框架，以预测检测插入序列所需的读取次数，从而确认其整合到宿主基因组中。只要所需的插入顺序已知，所有这些工具都有助于常规实验室测试中的检测工作。然而，也会产生未描述的转基因生物。对于这些基因，插入宿主基因组的序列是完全未知的（GMO设计者没有提供）。因此，用上述方法很难检测到它们。我们必须区分3种未知的GM类型，即仅含有抑制基因的GM、利用CRISPR/TALEN技术进行精确突变的GM，以及最终由外源CDS插入、CDS截断或与外源CDS融合组成的GM。为了检测这第三类，也是最常见的未知GM，我们提出了一种新的生物信息学工具，称为DUGMO，目前专门用于分析细菌基因组。

每个物种都有自己的词汇表，可以支持统计分析。研究最多的词汇属性是密码子用法[12]，但在统计方法中使用较长的单词来解读具有生物功能的基序[13]，检测水平基因转移[14]或鉴定物种特异性病毒序列[15]例如。迄今为止，没有使用这些类型属性的应用程序用于GMO检测。DUGMO填补了这一空白。

DUGMO机器学习步骤使用两组CDS。第一组对应于宿主基因组的CDS，具有极其精确的词汇定义。第二组对应于具有非常广泛词汇的已知GMO CDS，仅略述了可能的多样性。为了检测GM插入，DUGMO将把预测集中词汇与宿主基因组词汇不同的所有CDS报告为GM。

DUGMO检测基因改造（GM）细菌未知外源插入物的序列。它使用高通量测序数据，这些数据之前通过专用管道和宿主全基因组上的两个BLASTN比对进行了清理和分类（第2.2节）。在清理和分类步骤之后，清理已知GMO插入物的数据库，以删除该物种的野生CDS（第2.2节）。所提出的工具基于与自动学习步骤（第2.4节）相关的距离的统计计算（第2.3节）。由于删除或替换而产生的微小修改而产生的转基因生物将不予考虑。CRISPR-CAS9等方法产生的转基因生物[16]只有当它们结合了至少一种外源性CDS或涉及基因截短或融合时才能被检测到。检测到染色体整合和质粒插入。

结果

数据集集合

转基因菌株枯草芽孢杆菌已在德国巴伐利亚健康和食品安全局（LGL）的Hi-seq 1500系统中进行测序[17,18]. 此数据集在下面称为“Data1”。野生型配对Hi-seq 2500测序数据枯草芽孢杆菌所使用的9407菌株可在NCBI序列读取档案（SRA）上获得，登录号为SRR8935610。此数据集在下文中称为“Data2”。三株转基因菌株的配对Mi-seq测序数据大肠杆菌[19]已在NCBI SRA上使用，注册号为SRR9304542、SRR930.4539和SRR9304 540。这三个基因组来自同一菌株，但每个基因组都有一个不同的质粒，包含第E部分来自三种不同物种之一的基因：根癌农杆菌,结核分枝杆菌和化脓性链球菌在撰写本文时，使用的NGS数据集是唯一公开可用的明确描述为转基因细菌的数据集。两个参考基因组，枯草芽孢杆菌街道168（NC_000964.3）和大肠杆菌街道K-12子街道。还使用了MG1655（U00096.3）。构建了两个全基因组，一个用于枯草芽孢杆菌由37株菌组成，其中一株用于大肠杆菌由45株菌组成。泛基因组是人工构建的，涵盖了所有分析物种，可在DUGMO git库的“testdata”目录中找到。对于枯草芽孢杆菌，我们提取了一个已知的全基因组[20]并用NCBI上提供的各种环境中的野生菌株对其进行富集（附加文件中给出的登录号列表1，第2节）。对于大肠杆菌，5个主要物种亚群[21]考虑并丰富了NCBI上各种环境中更新的完整PacBio基因组和质粒（附加文件中给出的登录号列表1，第3节）。创建了一个由JRC提供的转基因编码序列（CDS）数据库，并以文献中的细菌转基因生物插入CDS作为补充（附加文件中规定的列表2，附加文件中给出了相关文献参考三). 该数据库目前不包括本文分析的细菌的转基因CDS。DUGMO需要单个细菌物种的成对Illumina测序数据；它不处理宏基因组或多基因组数据。本文提供的转基因数据库中包含的JRC序列是通过对包括专利和可用的全植物基因组在内的公共核苷酸序列数据库的PCR模拟筛选获得的[22].

清理原始数据、宿主泛基因组和已知GMO数据库的管道

本节的目的是描述三组数据的产生：宿主基因组CDS、潜在的GM CDS和已知的GMO CDS。为了创建这些数据集，实施了一条管道，通过删除属于参考基因组的所有序列来清除可疑转基因细菌的原始数据。此步骤经常用于处理高通量测序数据（例如，参考Baron等人2018[23]去除染色体读数进行组装大肠杆菌质粒）。该管道将CDS分为两类：潜在转基因插入物的CDS和属于野生型基因组的CDS。需要几个步骤（图1). 首先，使用Shovill 1.0.4对潜在GM细菌的原始测序数据进行两次组合[24]使用修剪选项，一次在整个基因组上，一次集中在质粒上（SPAdes 3.13.1的质粒选项[25]). 拆除覆盖深度小于3的脚手架；该阈值旨在消除细菌基因组的潜在污染，例如由另一个测序文库或细菌所在基质造成的污染。这两个装配结果用Prokka 1.12注释[26]使用use-genus选项预测CDS。合并这些CDS并过滤出重复序列。同时，使用Trimmomatic软件0.39对潜在转基因细菌的原始测序数据进行“修剪”[27]，并按照Bowtie2 2.3.4.3校准[28]参考基因组。然后，使用BWA MEM 0.7.17将非对齐读取（包括任何GMO插入）与之前Prokka CDS中预测的CDS对齐[29]. 使用Samtools 1.9从与BWA对齐的读取中推断出涵盖CDS的共识[30]. 在这些步骤的最后，本BWA比对共识中存在的CDS以及包含潜在GM插入序列的CDS被保留为潜在GM CDS。其余CDS作为与宿主基因组相关的CDS保存。

图1

原始数据清洗管道图，允许对潜在GM基因组的CDS进行初步分类

全基因组将同一物种中存在的所有基因组合在一起。两个BLASTN 2.6.0[31]清洁管道产生的潜在GM CDS在潜在GM细菌的全基因组上进行比对，一个在CDS上，另一个在整个基因组上（附加文件1，图S1）。这些BLASTN比对的目的是进行二次排序（第一次是通过清洁管道进行的），以便从潜在的GM CDS中分离与宿主基因组相关的CDS。全基因组CDS上的第一个BLASTN删除了潜在GM插入物中存在的属于宿主基因组的CDS（有关参数，请参阅附加文件1，第4.1节）。全基因组上的第二个BLASTN旨在删除Prokka预测的CDS，但在我们的全基因组中对应于非编码序列（有关参数，请参阅附加文件1，第4.2节）。将至少在两个BLASTN步骤中的一个步骤中对齐的CDS添加到与清洁管道末端定义的宿主基因组相关的CDS中。不结盟CDS作为潜在的GM CDS保存。

在这两个排序步骤之后，对已知GMO插入的数据库进行过滤，以定义用于距离计算的一组已知GM CDS（第2.3节）。只有当CDS在不同于其自然存在的基因组的环境中考虑时，CDS才被认定为GM。当目的是确定细菌或生物体是否为GMO时，使用BLASTN比对从我们的GMO CDS数据库中消除疑似菌株物种基因组中自然存在的所有CDS（有关参数，请参阅附加文件1，第4.3节）。此外，在此BLASTN期间对齐的CDS部分也会从数据库中的CDS中删除，从而生成一个GMO CDS的过滤库。我们计算了每个宿主基因组CDS和所有宿主基因组CDSs集合的三个距离（第2.3节中的P L3M1、P L4M2和F L9M7）的中位数，并将它们分别称为medL3M1、medL4M2和medL9M7。最后，丢弃已知GMO数据库中的CDS，该数据库与所有宿主基因组CDS集的距离已验证（P L3M1<medL3M 1和P L4M2<med L4M 2和F L9M7<medL 9M7）。

经过两个排序步骤和数据库过滤，我们生成了三组不相交的CDS，其中每个CDS分别与所有宿主基因组CDS进行比较：（a）GMO CDS的过滤库，（b）宿主基因组CDSs，和（c）潜在的GM CDS。

宿主基因组词汇和bray-Curtis（BC）距离

一般信息

为了评估每个潜在的GM CDS是否距离宿主基因组的CDS较远，提出了一种距离计算方法，即距离越大，CDS成为GM插入物的概率越高。我们的假设是，引入未知GMO的插入CDS应该与GMO来源的基因组具有不同的“词汇”。事实上，基因组有自己的词汇表，由单词组成。单词是具有预定义长度的短核苷酸序列，例如“ATGCCT”。词汇表，表示为M（M），在我们的方法中使用的可以是给定长度的所有单词的集合我或这些单词的子集，特定于宿主基因组。在第二种情况下，最初使用R'MES软件进行单词选择[32]搜索给定序列中的异常单词，即相对于给定的马尔可夫模型，明显过度或不足的单词。对于给定长度的每个单词，R'MES计算例外性分数，该分数衡量其观察到的和预期出现次数之间差异的显著性。为此，出现次数的分布近似于高斯分布，很适合频繁出现的单词。

布雷-柯蒂斯（BC）距离

比较多个备选距离后（附加文件1，第5节），选择了最相关的Bray-Curtis差异。这个相似度指数并不能验证三角形不等式，因此不应将其称为距离。然而，为了便于理解，本文将使用术语“Bray-Curtis距离”。布雷-柯蒂斯距离最初是一种生态距离[33]评估两个给定样本在可用物种丰度方面的差异。当两个比较数据集没有公共点时，该距离等于1；当它们的组成相同时，该间距等于0。我们的目的是比较一个序列S公司和一组序列H（H）根据他们的构词法，通过以下方程式

$$BC\left（S，H\right）=1-\frac{2{\sum}_{m\ in m}\left[f\left$$

(1)

哪里M（M）是所选单词的长度集我（如果没有选择，M（M）是所有4^我可能的长度词我字母{A，C，G，T}），m是一个单词M（M）,（f）是设置为的函数F类（对于频率；在以下等式中定义(2))或P（P）（用于比例；定义如下公式(三))，与图案的数量有关米在特别的序列中。实际上，S公司通常是CDS（是否为潜在GM），并且H（H）是推导出的主机CDS集，导致这两个序列集的累积长度存在显著差异，然后在S公司和H（H）为了避免这个问题，Bray-Curtis距离是根据归一化计数计算的，并考虑了两种归一化，称为频率与比例，如下所述。

BC频率距离

这个（f）等式中的函数(1)设置为

$$f\左（m，S\右）=f\左（m，S\左）：=\frac{C\左（米，S\右侧）}{sum_{w\在m}C\左$$

(2)

哪里C类(米,S公司)表示单词的出现次数或计数米按顺序S公司，以及w个每个单词都在里面吗M（M）换句话说，f（m，S）是m在S中的归一化频率。对于Bray-Curtis距离的计算，只考虑CDS密码子第三位置的核苷酸，用CDS表示_三附加文件中给出了密码子第三位置串联的理由1，第6节。CDS_三set允许我们考虑单词的大小n个包含大小为3的单词的有用信息n个在整个CDS中。

BC距离（按比例）

这个（f）等式中的函数(1)设置为

$$f\左（m，S\右）={P} _小时\left（m，S\ right）：=C\ left（m，S\ right）\times\ frac{\sum_{w\ in m}C\ lert（w，H\right）}{\sum_{w\ inM}C\left（w，S\ reght）}$$

(3)

换句话说，f（米，秒）是计数的标准化米在里面S公司相当于考虑米在一个与H（H）注意，根据此定义，当S=高度，我们得到f（m，H）=C（m，H）.

初步结果

在初始步骤中，以频率和比例计算BC距离，单词长度为我从3到9，最大阶马尔可夫模型k=l-2以及所选单词的不同比例。这些参数组合在我们的数据集Data1和Data2上进行了测试，其中已知GM细菌状态（第2.1节），以便从距离方面更好地将宿主CDS与GM CDS分开。在下面描述的机器学习步骤中，选择了三种不同的组合来共同使用：

• （P L3M1）：1阶马尔可夫模型下所有3个字母单词的BC比例。

• （P L4M2）：在2阶马尔可夫模型下，所有4个字母单词的BC比例。

• （F L9M7）：在7阶马尔可夫模型下，10%最具代表性的9个字母单词的BC频率。

之所以选择P L3M1，是因为此计算表征了CDS的三个字母单词，这是特定于宿主基因组词汇的，包括密码子用法。F L9M7使用长的过度表示单词（由于密码子第三个位置的串联），使在比较CDS中发现的单词非常特定（因为长度），并适合宿主基因组密码子的用法（只考虑CDS 27个字母中的第三个字母）。P L4M2允许我们通过使用不同于密码子用法的较小且因此更频繁的单词来表征CDS，以便在CDS上获得更准确的单词分布。

这些Bray-Curtis距离计算（P L3M1、P L4M2和F L9M7）是在清洁管道末端（第2.2节）获得的三组CDS中的每一组CDS上进行的，如图所示2，以便准备机器学习步骤（下一段）。注意，只使用长度大于或等于27个核苷酸的CDS，因为在计算Bray-Curtis距离的频率时考虑了单词大小的最小长度。

图2

用于获取机器学习步骤中使用的学习和预测数据的比例（L3M1、L4M2）和频率（L9M7）的Bray-Curtis距离

用机器学习方法预测GM镶块及其选择

方法

我们的目标是瞄准真正的通用插入式CDS。为此，在机器学习过程中应用了大量分类方法（R 3.5.1中的插入符号6.0_81包[34])，并选择了在预测GM CDS方面最有效的方法。使用的12种分类方法是：4种参数化方法：广义线性模型（logit）、逐步线性判别分析（stepLDA）、逐步二次判别分析（stepQDA）和偏最小二乘法（plsda）；8种非参数方法：神经网络（nnet）、径向基函数核支持向量机（svmRadial）、K-最近邻（knn），以及5种基于决策树的方法：分类树（Rpart）、随机森林（RF）、袋装分类（Treebag）、极端梯度提升树（xgboost）、，监督分类算法C5.0。在此过程中，使用了K折交叉验证程序。它包括将整个数据集拆分为K个子集，即校准集（由K-1个子集组成）和验证集（剩余的K个子集）。第一步，使用校准集选择模型参数及其拟合能力。第二步，在预测集上验证该模型，以获得其预测能力。该校准和预测过程重复了K次。每种机器学习方法的参数均按照10倍的交叉验证程序进行优化。在两倍交叉验证的基础上，对这些方法的性能进行了比较，其中数据之前是分层和集中的。值得注意的是，分层抽样在抽样数据中提供了GM CDS的均匀分布。

比较方法的标准

枯草芽孢杆菌用GM（Data1）选择最佳预测方法，提出四个性能标准：假阴性率、特异性、敏感性[35]，和假阳性率。在这四个标准中，最关键的一个是，该方法不应将GM插入CDS预测为宿主基因组CDS（即无假阴性）。还计算了每种方法变量的重要性。对于每种方法，计算ROC曲线和六个可用变量中每个变量使用的重要性（附加文件1，第7节）。

数据

学习数据是宿主基因组的CDS和筛选出的已知GMO CDS的数据库（图。2a和b）。从这些学习数据中删除编码RNA的CDS。预测数据为候选GM CDS（图。2c） 与泛基因组不一致（第2.2节）。相关变量为L4M2、L3M1比例的Bray-Curtis距离和L9M7频率、CDS长度、平均值A类_H（H）（S）R'MES在L4M2和L9M7宿主基因组中提供的异常分数（参见公式。4)，CDS中4个和9个字母单词的每个核苷酸的计数密度（计数总和C（米，秒）所有选定单词的米除以CDS的长度S公司)CDS GC含量。每个核苷酸的计数密度是一种测量CDS中宿主异常词比例的方法。例外性得分的平均值是衡量CDS中单词的例外性的一个指标，考虑到宿主CDS的组成

$$ {A} _小时（S） =\ frac{\sum_{m\ in m}K\ left（m，H\ right）\ times C\ left$$

(4)

哪里K（米，高）是一个单词的R'MES异常性得分米在主机CDS中H。

结果

结果如图所示三，从运行50次的12种分类方法中获得，并应用于学习数据（图。2a和b）枯草芽孢杆菌转基因细菌（数据1）。这些过程冗余确保了结果的稳健性，有些方法是非确定性的。

图3

使用的学习数据对50个仿真进行了12种分类方法的结果枯草芽孢杆菌转基因细菌（数据1）。校准和预测结果来自双重交叉验证程序

RF、C5.0、Treebag和xgboost方法在校准和预测方面给出了最佳结果。为了进行校准，RF、Treebag和C5.0方法的灵敏度和特异性为1，假阳性率和假阴性率为0。xgboost法的敏感性和特异性为0.99，假阳性率和假阴性率为0.01。为了进行预测，RF和C5.0的敏感性和特异性为0.98，阳性率和假阴性率为0.01。Treebag和xgboost方法的敏感性和特异性为0.97，假阳性率和假阴性率为0.015。对于logit方法，预测数据和校准数据具有相同的结果。logit方法的灵敏度为0.95，特异性为0.94，假阳性率为0.05，假阴性率为0.04。

根据我们的性能标准（即假阴性率、特异性、敏感性和假阳性率），在自动学习步骤中保留了RF和logit方法阳性结果的非冗余结合。然后，推导出程序的全局敏感性和特异性（附加文件1第7.1节）。使用Random Forest方法预测CDS的最重要变量是GC含量和Bray-Curtis距离（比例L4M2），使用logit方法预测频率L9M7和比例L3M1的Bray-Curitis距离（附加文件1第7.2节）。根据其变量重要性，它们具有互补的选择标准，并且不属于同一方法家族。

DUGMO在实际数据中的应用

DUGMO是使用来自枯草芽孢杆菌（即数据1、数据2）。验证于大肠杆菌基因组。表1总结了DUGMO在这两个方面的结果枯草芽孢杆菌基因组（即Data1、Data2）大肠杆菌基因组。GM上的GM CDS枯草芽孢杆菌DUGMO正确预测的（Data1）是位于质粒（核黄素）上的主要基因和产生氯霉素抗性的基因。将RF和logit方法的结果结合起来，为准确识别潜在的GM插入物提供了附加值。野生型自动学习步骤结束时枯草芽孢杆菌（数据2），没有CDS被预测为GM插入（参见GM行枯草芽孢杆菌和表中的机器学习列1). 与总经理大肠杆菌携带的基因A.肿瘤,结核分枝杆菌或化脓性链球菌，结果分别包含5个、4个和5个预测为GM插入的CDS。然而，在以下结果中发现了一个假阳性大肠杆菌携带的基因化脓链球菌。找到的CDS对应于带注释的aslA公司基因存在于大肠杆菌在我们的菌株中被强烈截断（984nt而不是泛基因组中的1650nt或更多）。如果将来自宿主基因组中代表性较差的蛋白质家族的高度截断蛋白质并入，则该序列可被视为GM插入。如果提交的未知基因组实际上是转基因细菌，则DUGMO可以靶向插入的CDSs。

表1两种药物的DUGMO结果枯草芽孢杆菌基因组（即Data1、Data2）和三个大肠杆菌基因组

为了完成我们工具的评估，DUGMO在另外6种野生型细菌上进行了测试：空肠弯曲菌,乳酸乳球菌,单核细胞增生李斯特菌,结核分枝杆菌,鼠伤寒沙门氏菌和金黄色葡萄球菌检查结果中几乎没有假阳性。然后，基于这六种带有七个外源基因的细菌（包括来自这些外源细菌的五个基因、一个人类基因和一个水稻基因）的合成转基因样品，对DUGMO进行了45次额外测试。每个基因插入参考基因组中的两个CDS之间。然后，ART 2.5.8[36]使用art_ilumina progam对每个修改后的基因组使用，参数为：-ss HS25，−p，−l 150，−f 11，−m 200，−s 10，以创建合成illumina测序数据。用户可以使用这些测试中使用的所有全基因组进行转基因分析。提供了全基因组，结果详见附加文件1，第8节。

在对野生型细菌的6项测试中，我们只获得了一项假阳性，即结核分枝杆菌我们强烈怀疑来自最近水平基因转移的基因。事实上，它没有巨细胞[37]与NCBI wgs数据库中发现的所有放线菌基因组的94.7%相比，一致性更好（2020年5月1日）。在45次合成GM测试中，我们只得到了一次假阳性，这可能是由于金黄色葡萄球菌pangenome和以环境和相似基因组特征插入的C.空肠基因（详见附加文件1第8.6节）。

然后，我们用具体的合成数据评估了该方法的局限性。我们评估了该方法对双齿类优化的鲁棒性。为此枯草芽孢杆菌核黄素合成酶基因编码，优化用于大肠杆菌使用双齿类动物，被随机插入大肠杆菌然后利用ART和前面描述的参数从该基因组生成Illumina测序数据。该程序运行了10次。在所有情况下，DUGMO检测到枯草芽孢杆菌基因作为GM插入物，从而证明该方法对二齿类优化不敏感，因此密码子优化。

为了评估使DUGMO检测到野生型基因所需的替代比例，我们在野生型基因中人为引入替代枯草芽孢杆菌两个野生基因独立测试，一个短（417核苷酸序列）和一个长序列（1317核苷酸）。野生型Surfactin长基因和CadI短基因枯草芽孢杆菌基因组经历了n个随机替换，避免DUGMO很容易检测到的突变事件，例如起始密码子、终止密码子中的替换，以及序列中引入早期终止密码子的替换。每个修改过的基因都被替换成野生型枯草芽孢杆菌组件。最后，使用ART生成Illumina测序数据[36]如上所述的软件。这个过程重复了10次。DUGMO结果表明，超过9%的突变，突变CDS被检测为GM。

讨论

我们提出了一种检测未描述的转基因细菌的方法。DUGMO旨在评估纯化细菌后是否存在细菌转基因，如Paracchini等人[18]. 一旦DUGMO确认可疑的转基因，就可以使用针对DUGMO识别序列的PCR对食品或环境样品中的转基因进行常规检测。DUGMO利用pangenome和基因组物种特性来提供机器学习，并将转基因CDSs与宿主CDSs区分开来。泛基因组的稳健性在很大程度上依赖于其完整性的实现，也依赖于将其提交给公共数据库时对各自物种的准确基因组归属。物种归属通常是从16S RNA推导出来的，导致归属中经常出现错误，进而可能导致混合物种全基因组[38]. 由于宿主CDS排序步骤中BLAST匹配不一致，在全基因组中引入合成菌株也可能导致统计特性有偏差（第2.3节）。这突出表明迫切需要专门用于评估全基因组完整性的工具[39].

最近，Berbers等人[40]描述了未经授权的GM的修订枯草芽孢杆菌菌株2014–3557，具有一个pGMrib质粒和53kb染色体插入，Paracchini等人首次描述（数据1）[18]带有四个质粒pGMBsub01–04。为了清楚起见，在以下结果中，我们提供了Paracchini等人使用的括号内的基因名称枯草芽孢杆菌包含一个整合质粒pGMrib，设计用于生产维生素B2（核黄素）和几个染色体修饰。首先猫将与氯霉素抗性相关的基因插入染色体记录基因.其次，在短链蛋白A基因。在这个集成中里德阿特和部分里贝达操纵子与肋骨的操纵子解淀粉芽孢杆菌此外，在53 kb的插入中发现了多个β-内酰胺酶和卡那霉素耐药基因拷贝。最后，存在博来霉素抗性基因。

两者都是猫和记录DUGMO检测到染色体部分为GM CDS。一个枯草芽孢杆菌基因被机器学习步骤丢弃，即染色体RNA聚合酶δ亚基。DUGMO区分GM解淀粉芽孢杆菌scpB野生染色体插入的CDS枯草杆菌scpBpGMrib质粒。pGMrib质粒（附加文件中的质粒可视化1，部分9）包括两个部分。第一个有一个纯里氏枯草杆菌DUGMO识别并丢弃的操纵子，包括sipS公司(麻风杆菌)和GAY71_RS22375(gerPA公司). 第二种携带对红霉素产生耐药性的杆菌基因。在第二部分中，DUGMO将九个基因识别为GM(GAY71_RS22270(代表S),仓3/GAY71_RS22275(贝塔)，soj公司(三角洲)，泽塔，GAY71_RS22265(警察)，GAY71_RS22365(ORF磅/平方英寸)，GAY71_RS22350(ORFφ)，GAY71_RS22355部分(ORF卡)和GAY71_RS22370(hyprORF-t型))，而其他十个则被丢弃(GAY71_RS22250(埃塔),GAY71_RS22255(θ）,GAY71_RS22260(奥塔),顶部B/GAY71_RS22280(伽马射线),GAY71_RS22290/GAY71_RS22330（Ω）,GAY71_RS22295/GAY71_RS22325（ε）,erm（B）,GAY71_RS22315号具有从核苷酸14107至14138的额外部分(电流互感器)，GAY71_RS22340(陶加玛)和GAY71_RS22345(小西隆)). 其中电流互感器和陶加玛基因在全基因组的BLASTN步骤中被丢弃；有趣的是，这些基因的核苷酸序列与枯草芽孢杆菌并且被包括在主机CDS学习数据中。大多数其他被丢弃的CDS与肠球菌基因最匹配，这些基因在系统发育上与革兰氏+杆菌基因接近枯草芽孢杆菌部分pGMrib（Paracchini等人的pGMBsub04）被认为来源于专利的pMX45质粒枯草芽孢杆菌提到了起源和基因工程方法[18]. 然而，一些质粒具有非常广泛的宿主范围，并不能总是断言其来源。在这种情况下，这十个基因的遗传距离，如果它们不是枯草芽孢杆菌出身，不允许他们歧视。有趣的是，DUGMO发现了两次未在pGMrib（pGMBsub04部分）中描述的额外CDS，一次在质粒上，一次是在染色体上。该CDS位于pGMrib中的2652至2822位置，包括GAY71_RS22250（pGMBsub04质粒中为1至171），对应于粪肠球菌质粒（登录号AP018546.1）。此外，pGMrib（pGMBsub03部分）有两个额外的非枯草芽孢杆菌基因，部分GAY71_RS22445(代表B)和的缩写版本泰特（L）(四氯化碳)在DUGMO结果中发现GM。53kb的染色体插入（pGMBsub01部分）携带不同来源的基因(解淀粉芽孢杆菌和金黄色葡萄球菌). 这个金黄色葡萄球菌基因GAY71_RS12530(混音器（bleR）)和aadD1(堪萨斯州)和肋骨D,里巴,肋骨T和scpB型的CDSB.淀粉液化菌（ribDEATH）操纵子被DUGMO正确检测为GM CDS。β-内酰胺酶TEM-116（a）基因mpR/放大器，染色体上有6个位点,pGMrib中的一个），肋骨和短链蛋白ACDS在BLASTN参数方面与全基因组CDS一致，并包含在宿主CDS学习数据中，显示了CDS分类对密切相关物种的限制。这可能也解释了未检测到代表和里布（ribE）.Prokka注释提供了两个额外的假设CDS，匹配的位置为4759–5103和10112–10201B.淀粉液化菌（ribDEATH）操纵子（pGMBsub01部分）（位点GAY71_RS12515,GAY71_RS12560,GAY71_RS12605,GAY71_RS12650,GAY71_RS12725,GAY71_RS12770和GAY71_RS12480,GAY71_RS12690分别位于染色体上）。第二种是通用汽车，第一种不是。

DUGMO，关于大肠杆菌菌株携带的基因化脓性链球菌，发现为假阳性CDS，对应于强截断阿尔卑斯CDS（芳基硫酸酯酶蛋白编码），其全长版本自然存在于大肠杆菌基因组。在分析期间，全基因组CDS步骤上的BLASTN验证了匹配CDS的长度变化不超过15%。对于其他两个大肠杆菌基因组，携带A.肿瘤或结核分枝杆菌，的阿尔卑斯CDS是组件中的全长，DUGMO未检测到。被截断的阿尔卑斯CDS输入大肠杆菌携带的基因化脓链球菌在一个新的contig的末尾被定位，15%的缺失部分在另一个contig中被发现。因此，DUGMO能够检测到可能包含在野生CDS中长时间删除或插入的CDS。该样本数据的平均覆盖深度为41，而其他两个样本的覆盖深度平均值大肠杆菌菌株分别为63株和68株。此示例表明，由于CDS截断，平均重叠深度为41的程序集是部分的，并导致误报。另一方面，覆盖深度平均值大于60的数据集合在我们的数据集中不存在基因截断问题。大多数细菌汇编程序建议覆盖深度为80[41]和100（N50显示随着覆盖深度的增加而增加，直到细菌基因组的值为100[42]). 覆盖深度为60可能是在不截断CDS的情况下获得完整组件的可接受下限。这些观察结果强调了需要获得良好的装配，从而获得足够的覆盖深度。在未来的版本中，DUGMO将能够接受程序集，包括长读程序集，从而消除由于程序集截断而导致的误报情况。

Paracchini等人描述的转基因细菌分析过程[18]需要一个多月的工作（个人通信），而DUGMO需要三个小时，使用64个内核i7 CPU的10个线程GB RAM。数据处理的速度是该工具的一个主要优势，此外，除了验证结果中检测到的有限可能的转基因基因集外，终端用户不需要事先具备转基因专家知识。用户必须具备Linux命令行技能和足够的生物学知识才能构建泛基因组。此外，DUGMO具有持续改进的潜力，最终用户可以在学习数据中添加新的已确认GMO CDS，以及在全基因组中添加新野生型基因组和CDS。

结论

所提出的DUGMO工具将高通量测序数据清理管道与pangenomes上的BLASTN比对以及与机器学习方法、随机森林和广义线性模型的组合相关的不同Bray-Curtis距离计算相结合，被选为对转基因细菌数据最具预测性的工具。该工具需要可靠的物种特异性全基因组和参考基因组，即疑似转基因物种的成员。在GM上测试枯草芽孢杆菌和三位总经理大肠杆菌，DUGMO能够检测外源性或截断或融合的GM CDS，并生成很少或没有假阳性和假阴性。该工具是基于野生型基因组的选定特性，对未知单一转基因细菌可检测性概念的证明。由于DUGMO使用与CDS第三位相关的统计特性（特定于密码子的使用），因此它不会只检测到tRNA或rRNA插入物的转基因细菌。从概念上讲，DUGMO是为了找到不使用宿主基因组词汇的CDS。DUGMO无法区分水平转移基因和GM基因。DUGMO的使用可能会转向在野生型细菌中发现最近水平基因转移获得的基因。这种DUGMO的使用没有经过测试。我们计划通过调整距离计算参数，将DUGMO扩展到植物和动物转基因检测。最后，DUGMO首次能够快速、系统地检测可疑样本中未知的转基因细菌。