CN110600083A - 基于无拼接组装wgs数据的醋酸钙—鲍曼不动杆菌复合群鉴定方法 - Google Patents

Abstract

本发明公开了基于无拼接组装WGS数据的醋酸钙—鲍曼不动杆菌复合群鉴定方法。本发明提供了一种基于无拼接组装WGS数据对待鉴定菌进行种属鉴定的方法。本发明的基本原理是通过建立一个完整的菌种基因组指纹特征数据库，然后将待鉴定菌株的WGS测序reads直接打散为片段序列，通过和特征指纹数据库的各菌种进行比较打分，从而实现复合群菌种的鉴定。本发明方法不需要对测序reads进行组装，因而非常简单快捷，利用了全基因组的信息，另一方面由于本发明构建了包含2279种的细菌指纹特征数据库，不仅可以用于鉴定醋酸钙—鲍曼不动杆菌复合群中菌种，同样可适用于其它复合群菌种或其它菌种的鉴定。

Description

基于无拼接组装WGS数据的醋酸钙—鲍曼不动杆菌复合群鉴 定方法

技术领域

本发明涉及生物技术领域，具体涉及基于无拼接组装WGS数据的醋酸钙—鲍曼不动杆菌复合群鉴定方法。

背景技术

不动杆菌属(Acinetobacter genus)是一种革兰氏阴性菌，目前不动杆菌属包含55个菌种。该属病原菌在医院环境中分布很广且可以长期存活，在ICU病房中尤其流行，极易造成危重患者的感染，常导致菌血症，肺炎，脑膜炎，尿路感染，手术部位感染等。其中的鲍曼不动杆菌(Acinetobacter baumannii)是院内感染最主要，最流行的致病菌。不动杆菌属的菌种遗传上高度相似难以鉴定，将其鉴定到种的水平仍然难以把握，其中最难区分的是醋酸钙-鲍曼不动杆菌复合群(Acinetobacter calcoaceticus-Acinetobacterbaumannii(ACB)complex)，该复合群主要包括，鲍曼不动杆菌(Acinetobacterbaumannii)，醋酸钙不动杆菌(Acinetobacter calcoaceticus)，皮特不动杆菌(Acinetobacter pittii)，医院不动杆菌(Acinetobacter nosocomialis)4个菌种，这几个菌种表型和遗传上都非常相似，难以鉴定区分，临床上经常只将结果报为醋酸钙-鲍曼不动杆菌复合群。其中除醋酸钙不动杆菌分布在环境中外，其它几个菌种都会导致人感染，尽管鲍曼不动杆菌被认为是不动杆菌属中最流行和最致命的成员，但皮特不动杆菌和医院不动杆菌也会引起严重的侵入性疾病。

鲍曼不动杆菌、醋酸钙不动杆菌、皮特不动杆菌和医院不动杆菌虽然遗传上和表型相似，但流行病学特征却大相径庭，醋酸钙不动杆菌主要存在于环境标本，皮特不动杆菌主要存在于皮肤表面和环境标本，医院不动杆菌主要存在于临床标本，而鲍曼不动杆菌则为院内感染的主要病原体之一。随着近年来不动杆菌对常用抗生素的耐药率呈增加趋势，甚至出现了多重耐药和泛耐药菌株，更加引起临床医生和微生物研究人员的高度关注。很多研究发现，醋酸钙-鲍曼不动杆菌复合群中各菌种的耐药特性存在较大差异，并且在临床表现和治疗上都有差别。菌种鉴定不精确的缺陷,导致目前很多有关鲍曼不动杆菌的耐药及流行病学资料存在片面性，严重影响了对鲍曼不动杆菌耐药现状及临床分布状况的认识。

临床上菌种鉴定主要依靠表型差异，由于不动杆菌属中菌种相似度极高，传统的生化鉴定方法在不动杆菌种属鉴定中存在局限性，而醋酸钙-鲍曼不动杆菌复合群更无法用基于表型的检测方法区分。目前VITEK-2等微生物自动分析仪凭借高效、便捷的优点成为医院临床进行菌种分析最常用的技术方法，但对于不动杆菌属的精确鉴定存在一定的缺陷，VITEK 2只能鉴定不动杆菌属中有限的几种，对生化特性相近的鲍曼不动杆菌、醋酸钙不动杆菌、皮特不动杆菌和医院不动杆菌则难以区分，一般只能鉴定为醋酸钙-鲍曼不动杆菌复合菌。

不动杆菌属鉴定除了生化方法外，也发展了一些分子生物学的方法，研究者发展出利用PCR方法扩增gyrB基因，或者多重PCR扩增16S-23S rRNA基因的基因间区以及看家基金gyrB和recA基因来鉴定区分不动杆菌属菌种，近些年来随测序的兴起，又出现了基于测序的鉴定方法，比如基于16S rRNA基因、rpoB基因、gyrB基因和recA基因测序鉴定法。这些分子方法都是只是利用了菌株的部分基因序列，由于使用的信息单一，此类方法对遗传上相似的菌种鉴定存在缺陷，虽然可以区分出不动杆菌属的某些菌种，但是仍无法精确有效区分醋酸钙-鲍曼不动杆菌复合群内的菌种。

另外，基质辅助激光解吸电离飞行时间质谱(Matrix-Assisted LaserDesorption Ionization-Time of Flight Mass Spectrometry，MALDI-TOF MS)是目前临床实验室中应用最为广泛的质谱技术之一，在微生物检测与鉴定中有着广泛和大量的应用。MALDI-TOF MS与生化表型分析和分子实验方法相比具有快速可靠、简单经济的优点，仍存在一些不足，方法的准确性依赖具备的数据处理分析软件和质谱数据库中菌种的信息量，多项研究指出使用MALDI-TOF MS的标准分析流程鉴定醋酸钙-鲍曼不动杆菌复合菌群，会出现错误的鉴定结果，如将医院不动杆菌确定为鲍曼不动杆菌，这主要由于目前国内菌种数据库中质谱分析图谱量较少，从而影响了其准确性。

目前，全基因组测序(Whole Genome Sequencing,WGS)技术的快速发展，使得我们可以利用全基因组测序数据对微生物进行鉴定。开发使用菌株的包括非编码区在内的基因组全部信息进行微生物鉴定和分析的方法会具有极高的分辨率。基于WGS全基因组测序数据的鉴定方法将能够精确区分醋酸钙-鲍曼不动杆菌复合群的菌种以及不动杆菌属的其它菌种。最新研究已表明醋酸钙-鲍曼不动杆菌复合群几个菌种在对抗生素的敏感性，致病性，及临床表现会显示出截然不同的特征，能够准确鉴定出鲍曼不动杆菌，并区分出醋酸钙-鲍曼不动杆菌复合群其它菌种对于院内感染的治疗、预后及监测有重要意义。

发明内容

本发明的目的是提供一种基于无拼接组装WGS数据的醋酸钙—鲍曼不动杆菌复合群鉴定方法。

第一方面，本发明要求保护一种基于无拼接组装WGS数据对待鉴定菌进行种属鉴定的方法。

本发明所要求保护的基于无拼接组装WGS数据对待鉴定菌进行种属鉴定的方法，可包括如下步骤：

(A)从菌基因组数据库中获取全部已测序的菌全基因组数据，建立菌株基因组信息和菌株分类学信息之间的关系；

(B)按照如下步骤构建菌种指纹特征数据库：

(b1)针对步骤(A)获得的每株菌的全基因组核酸序列，均以核酸序列长度为k(即k个bp长度)、步长为1(即1个bp长度)对全基因组核酸序列进行切分(长度为L的基因组可划分出L-k+1个核酸片段)，切分后序列重复的片段只保留一个，将所得片段集合称为集合A；每个菌株得到一个所述集合A；组成所述集合A的片段称为片段A；其中，k和L均为大于1的正整数。

(b2)将四种碱基的字母表示转换为00、01、10、11这种数字存储，从而使得每个所述片段A都能转换为一个2k位的数字，即0至2^2k-1之间的一个数；

这样存储只占用k/4个bytes(2k个bit)，而使用原始的字符(用A、T、C、G表示4种碱基)存储，一个所述片段A要占用k个byte，压缩比为4倍，且通过这样，在后续的计算上可以将字符串的比较转换为数字的查询，速度上则显著提升。

(b3)遍历所有0至2^2k-1之间的数字，对所有菌株的所述集合A中的所有所述片段A均一一做记录，记录各所述片段A所对应的2k位数字以及对应的种属信息，筛选出仅记录1个菌种的所述片段A对应的2k位数字，按照菌种分类存储，即得所述菌种指纹特征数据库。

进一步地，步骤(b3)可通过用python或其它语言编写计算机程序来实现。遍历所有0至2^2k-1之间的数字(即切分产生所有可能的所述片段A)，并用数字作为key建立字典；然后对所有菌株的所述集合A中的所有所述片段A(数字存储)在字典里做记录，发现一个数字(即所述片段A)后，则在相应所述片段A的字典里记录代表相应菌种的TaxID，在遍历过程中在如果在其它某菌种发现该数字(即该片段A)，则在字典里直接删除key(也就是去除该片段A)。直到所有菌种的全部菌株处理完成。最终筛选出仅记录1个菌种的所述片段A对应的2k位数字，按照菌种分类存储，即得所述菌种指纹特征数据库。

(C)将待鉴定菌株的全基因组测序所得每一条测序读长均按照步骤(b1)的方式进行切分，得到所述待鉴定菌株的所述集合A；按照步骤(b2)的方式将所述待鉴定测菌株的所述集合A中的所有所述片段A均转换为2k位数字；然后与步骤(b3)所得所述菌种指纹特征数据库进行比对，根据比对结果计算Score分值，得分最高的菌种视为所述待鉴定菌株所属菌种；

所述待鉴定菌株用于表示所述片段A的2k位数字与所述菌种指纹特征数据库中某一菌种A用于表示所述片段A的2k位数字，两者交集的绝对数量越多，则对于所述待鉴定菌株而言，所述菌种A的得分越高；

所述待鉴定菌株用于表示所述片段A的2k位数字与所述菌种指纹特征数据库中某一菌种A用于表示所述片段A的2k位数字，两者交集占所述菌种指纹特征数据库中所述菌种A用于表示所述片段A的2k位数字的数量比例越高，则对于所述待鉴定菌株而言，所述菌种A的得分越高。

实际应用中，综合“所述待鉴定菌株用于表示所述片段A的2k位数字与所述菌种指纹特征数据库中某一菌种A用于表示所述片段A的2k位数字，两者交集的绝对数量”以及“所述待鉴定菌株用于表示所述片段A的2k位数字与所述菌种指纹特征数据库中某一菌种A用于表示所述片段A的2k位数字，两者交集占所述菌种指纹特征数据库中所述菌种A用于表示所述片段A的2k位数字的数量比例”两方面因素通过加权求和确定Score分值的得分。

具体的，步骤(C)中，所述Score分值可按照如下公式计算获得：

Score＝α*normalization(N)+(1-α)*normalization(P)

N＝card(S_x∩S_i)

M＝card(S_i)

其中，S_x代表所述待鉴定测菌株用于表示所述片段A的2k位数字的集合；S_i代表所述菌种指纹特征数据库中某一菌种用于表示所述片段A的2k位数字的集合；α是加权系数；card表示求集合里的元素个数；N代表所述待鉴定测菌株用于表示所述片段A的2k位数字的集合与所述菌种指纹特征数据库中某一菌种用于表示所述片段A的2k位数字的集合的交集中的元素个数；M代表所述菌种指纹特征数据库中某一菌种用于表示所述片段A的2k位数字的集合中的元素个数。

根据我们的测试，参数α＝0.48，鉴定结果最好。N，P需要归一化处理为：

N_min，N_max，P_min，P_max是与所述菌种指纹特征数据库中全部菌种计算后得到的N，P的最小值和最大值。

在所述方法中，步骤(A)可按照如下进行：从NCBI的Genome数据库获取全部基因组状态Status为Complete Genome的菌全基因组序列数据；从NCBI的Taxonomy数据库获取生物分类(taxonomy)元数据信息，根据TaxID建立菌株基因组信息和菌株分类学信息之间的关系。

在步骤(b1)中，所述k应具有足够高的特异度，又综合考虑占用存储空间，选择长度16在程序设计时方便存储为一个32bit整数。所以，在本发明的具体实施方式中，所述k具体为16。

在步骤(b1)中，对全基因组核酸序列进行切分具体可使用python或其它语言编写脚本程序完成。

在所述步骤(A)中，所述菌可为细菌。

在所述方法中，所述待鉴定菌可为细菌。进一步地，所述细菌可为鲍曼不动杆菌。更进一步地，所述鲍曼不动杆菌可为隶属于醋酸钙-鲍曼不动杆菌复合群中的鲍曼不动杆菌。

第二方面，本发明要求保护一种构建菌种指纹特征数据库的方法。

本发明所要求保护的构建菌种指纹特征数据库的方法，可包括前文第一方面所述方法的步骤(A)和步骤(B)。

第三方面，本发明要求保护利用第二方面所述方法构建得到的菌种指纹特征数据库。

第四方面，本发明要求保护第三方面所述菌种指纹特征数据库在基于无拼接组装WGS数据对待鉴定菌进行种属鉴定中的应用。

其中，所述待鉴定菌可为细菌。进一步地，所述细菌可为鲍曼不动杆菌。更进一步地，所述鲍曼不动杆菌可为隶属于醋酸钙-鲍曼不动杆菌复合群中的鲍曼不动杆菌，具体可为鲍曼不动杆菌(Acinetobacter baumannii)、醋酸钙不动杆菌(Acinetobactercalcoaceticus)、皮特不动杆菌(Acinetobacter pittii)或医院不动杆菌(Acinetobacternosocomialis)。

在本发明的具体实施方式中，所述鲍曼不动杆菌除了上述4种隶属于醋酸钙-鲍曼不动杆菌复合群中的鲍曼不动杆菌外，还可为溶血不动杆菌(Acinetobacter_haemolyticus)、约氏不动杆菌(Acinetobacter_johnsonii)、琼氏不动杆菌(Acinetobacter_junii)、嗜油不动杆菌(Acinetobacter_oleivorans)、申氏不动杆菌(Acinetobacter_schindleri)或土壤不动杆菌(Acinetobacter_soli)。

在本发明的具体实施方式中，所述细菌除了上述鲍曼不动杆菌外，还可为其他种属的细菌，如图3中所示的至少74种细菌中任一种。

如已确定所述待鉴定菌为醋酸钙-鲍曼不动杆菌复合群菌种，则步骤(C)只需要计算数据库中醋酸钙-鲍曼不动杆菌复合群内各菌种的评分Score值进行排序即可。

本发明所提供的方法不需要对测序reads进行组装，因而非常简单快捷，利用了全基因组的信息，另一方面由于本发明构建了包含2279种的细菌指纹特征数据库，不仅可以用于鉴定醋酸钙—鲍曼不动杆菌复合群中菌种，同样可适用于其它复合群菌种或其它菌种的鉴定。

附图说明

图1为本发明鉴定方法原理及流程图。

图2为醋酸钙-鲍曼不动杆菌复合群菌种鉴定准确率与测序数据mapping比例变化关系。

图3为其它74种常见菌种的鉴定准确率。

具体实施方式

下述实施例中所用的数据、工具等，如无特殊说明，均可从商业途径得到。

本发明设计了一种直接使用全基因组WGS测序reads而不经过组装拼接对醋酸钙—鲍曼不动杆菌复合群菌种进行鉴定的方法。

本发明的基本原理是通过建立一个完整的菌种基因组指纹特征数据库，然后将待鉴定菌株的WGS测序reads直接打散为片段序列，通过和特征指纹数据库的各菌种进行比较打分，从而实现复合群菌种的鉴定。

本发明鉴定方法原理及流程图如图1所示。

本发明设计的方法具体包括下列步骤：

1、获取全部已测序的细菌全基因组数据

从NCBI(National Center for Biotechnology Information)Genome数据库获取全部可用的细菌全基因组序列数据，首先根据ftp.ncbi.nlm.nih.gov/genomes/GENOME_REPORTS/prokaryotes.txt文件得到测序数据的元信息，根据信息筛选出基因组状态Status为Complete Genome的细菌菌株，包含7571菌株，包含2279个菌种。并获取菌株全基因组核酸序列数据。

从NCBI的Taxonomy数据库获取生物分类(taxonomy)元数据信息，链接：ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdump.tar.gz，该数据提供了物种的分类学信息，可以获得菌株的种属及整个分类学信息。根据TaxID建立菌株基因组信息和菌株分类学信息之间的关系。

2、构建各菌种指纹特征片段数据库

对获取的菌株全基因组核酸序列，采取下述步骤得到每个菌种的指纹特征片段，从而构建涵盖能获取到全基因组数据的全部细菌菌种的指纹片段数据库：

(1)对获取到的每株细菌的全基因组核酸序列进行片段化，假定某株菌基因组核酸序列长度为L，选取碱基长度为kbp的核酸片段作为特征片段，以步长为1bp对全基因组核酸序列进行切分，长度为L的基因组可划分出L-k+1个核酸片段，将这样的片段集合称为集合A，即所有长度为k的子串，每个菌株得到一个所述集合A，组成所述集合A的片段称为片段A，切分后对于每株菌重复的所述片段A只保留1个。其中，k和L均为大于1的正整数。

(2)为了实现指纹片段和数据库的快速比对，本发明将碱基的字母表示转换为数字存储，4种碱基分别表示如下：A：00，C：01，G：10，T：11，这样每个kbp长度碱基片段都转换为一个2k位的数字，即0至2^2k-1之间的一个数，这样存储只占用k/4个bytes(2k个bit)，而使用原始的字符存储(用A、T、C、G表示4种碱基)，一个k-mer片段要占用k个byte，压缩比为4倍，且通过这样，在后续的计算上可以将字符串的比较转换为数字的查询，速度上则显著提升。

(3)求得每个菌种具有代表特征的指纹片段，按照如下步骤：

①根据(2)处理后所有所述片段A均为0至2^2k-1之间的一个数，以这个数作为所述片段A的key建立字典，遍历每个菌种里每一个菌株的切分得到的所有所述片段A(数字存储)，发现一个数字(即所述片段A)后，则在相应所述片段A的字典里记录代表该菌种的TaxID，在遍历过程中在如果在其它某菌种发现该所述片段A，则在字典里直接删除整个key，也就是去除这个所述片段A。直到所有菌种的全部菌株处理完成。

②遍历处理字典里剩余的所有数字/所述片段A(0至2^2k-1之间)，筛选出那些仅记录1个菌种的所述片段A(数字存储)，按照每个菌种分类存放，即得到每个菌种的指纹特征片段集合(数字存储)，形成涵盖超过2000多个菌种的指纹特征片段数据库。

3、对待鉴定菌株全基因组测序数据进行片段化

对待鉴定菌株WGS测序数据fastq文件中的每一条测序read，按照step＝1的滑动窗口将其切分成长度为k的碱基片段，这样每条测序read被处理成的碱基片段个数为：read长度-k+1，将所有读长处理完毕后，去除掉重复的碱基片段得到无冗余的kbp长的所述集合A，根据构建指纹数据库步骤中(2)所描述方法将每一个长度为kbp的所述片段A转换为数字存储，得到待鉴定菌株基因组的所有所述集合A。

4、与菌种指纹数据库对比

将待鉴定菌株的全基因组片段集合与构建的菌种指纹数据库进行对比，根据对比结果判定待鉴定菌所属菌种，具体方法如下：根据待鉴定菌株的全部所述集合A和指纹数据库中的每一个菌种的指纹片段集合计算Score分值，然后对结果进行排名，得分最高的菌种认为是待测菌株所属菌种，即鉴定结果。

打分函数的原理是考虑这个两个因素：1、待测菌株与数据库中菌种指纹片段交集的绝对数量，2、交集占某菌种指纹数据库中片段数的比例。我们综合这两方面因素通过加权求和确定得分，打分函数Score设计为：

Score＝α*normalization(N)+(1-α)*normalization(P)

N＝card(S_x∩S_i)

M＝card(S_i)

其中，S_x代表所述待鉴定测菌株用于表示所述片段A的2k位数字的集合；S_i代表所述菌种指纹特征数据库中某一菌种用于表示所述片段A的2k位数字的集合；α是加权系数；card表示求集合里的元素个数；N代表所述待鉴定测菌株用于表示所述片段A的2k位数字的集合与所述菌种指纹特征数据库中某一菌种用于表示所述片段A的2k位数字的集合的交集中的元素个数；M代表所述菌种指纹特征数据库中某一菌种用于表示所述片段A的2k位数字的集合中的元素个数。

根据我们的测试，参数α＝0.48，鉴定结果最好。N，P需要归一化处理为：

N_min，N_max，P_min，P_max是与所述菌种指纹特征数据库中全部菌种计算后得到的N，P的最小值和最大值。

下述实施例中所使用的实验方法如无特殊说明，均为常规方法。

下述实施例中所用的材料、试剂等，如无特殊说明，均可从商业途径得到。

下面的实施例对本发明中的步骤进行详细描述，实施例针对的是可公开下载的测序数据进行鉴定，显然方法适用于以任何方式获取的菌株WGS测序数据。本发明实施通过编写python计算机程序实现。

实施例1、基于WGS的的醋酸钙—鲍曼不动杆菌复合群鉴定方法

1、获取全部已测序菌株的基因组核酸序列数据

从NCBI(National Center for Biotechnology Information)Genome数据库获取全部可用的细菌全基因组序列数据，首先根据ftp.ncbi.nlm.nih.gov/genomes/GENOME_REPORTS/prokaryotes.txt文件得到测序数据的元信息，其中包含物种名称Organism/Name，分类学编号TaxID，基因组组装状态Status，基因组数据的唯一标识AssemblyAccession，及数据路径FTP Path等。TaxID可以追溯菌株的种属和分类学信息，Status表示菌株的基因组状态，根据基因组序列拼接完成度，共分为四种：Complete Genome，Chromosome，Scaffold，Contig。Assembly Accession，FTP Path里面有个由9个数字组成的ID，作为菌株基因组数据的唯一的标识。

本发明选择Status为Complete Genome的菌株基因组核酸序列数据构建数据库，包含7571菌株，包含2279个菌种。FTP Patth提供了基因组数据的下载链接。*_genomic.fna.gz文件是菌株的基因组核酸序列数据，数据示例链接：ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/009/085/GCA_000009085.1_ASM908v1/GCA_000009085.1_ASM908v1_genomic.fna.gz，其中GCA_000009085.1中的9位数字000009085作为基因组核酸序列数据的ID。我们将作为研究对象的7571个菌株的基因组核酸序列数据全部下载到本地以备构建数据库使用。

从NCBI的Taxonomy数据库获取生物分类(taxonomy)元数据信息，链接ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdump.tar.gz，根据该文件和TaxID我们建立菌株基因组核酸序列数据和菌株分类学信息之间的关系。

2、构建菌种指纹特征数据库

我们对每一菌株的基因组核酸序列进行处理，按照前述方法使用python语言编程将全基因组核酸序列切分为长度为16bp的片段(考虑到片段越长，会有足够的特异度，但太长会占用太多存储空间，16bp长度的片段转换后程序存储的时候正好是1个32bit整数)，得到菌株无冗余的16bp片段集合后，按照二进制编码：A：00，C：01，G：10，T：11，将每一个长度为16bp的片段转换为一个32位二进制数，即0至2³²-1之间的一个数。

按照前述方法用python语言编写计算机程序，遍历所有0至2³²-1的数字，即我们切分产生所有可能的16bp片段，并用该数字作为key建立字典。然后对每个菌种里每一个菌株的切分得到的所有16bp片段(数字存储)，在字典里做记录，处理完所有菌种的菌株后，筛选出那些仅记录1个菌种的16bp片段(数字存储)，按照每个菌种分类存储，即得到包含2279个菌种的指纹特征数据库。

3、获得醋酸钙-鲍曼不动杆菌复合群二代测序原始数据

从NCBI的SRA数据库中搜索所有不动杆菌属二代测序数据，筛选出测序数据量在100MB到3000MB之间的全部原始WGS测序reads，对于超过100个记录的种，我们随机抽样100株菌数据，而对少于100个记录的菌种，我们使用全部菌株。

WGS测序数据量大，推荐使用ascp工具(https://download.asperasoft.com/download/sw/connect/3.8.0/ibm-aspera-connect-3.8.0.158555-linux-g2.12-64.tar.gz)进行下载，其稳定性和速度均远远快于普通http及ftp下载方式，参考命令ascp-i～/.aspera/connect/etc/asperaweb_id_dsa.openssh-k 1-Tv-l 200m download_url target_dir。

SRA数据库以sra格式存储二代测序数据，需要将其转化为常用的fastq格式，以备下一步使用，利用sra toolkit(https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/current/)工具中的fastq-dump工具，命令：fastq-dump--split-3$sra_file--gzip-Ooutput_dir，生成fastq文件后根据输出信息清除那些异常的reads文件，最终得到315株复合群菌种测序reads数据。

4、利用WGS测序reads生成待鉴定菌株基因组核酸片段

对fastq文件中的每一条测序read，按照step＝1的滑动窗口将其切分成长度为16的碱基片段，这样每条测序read被处理成的碱基片段个数为：read长度-16+1＝read长度-15，将所有read处理完毕后，去除掉重复的碱基片段得到无冗余的16bp长片段集合。这一步同样利用python编写程序把reads序列切割成16bp片段，并按照前文所述编码转换为数字储存。

5、与指纹数据库进行对比得到菌种鉴定结果

按照前文方法中的公式计算待鉴定菌株与数据库中全部菌种的评分Score值，我们需要计算：(1)待测菌株切分得到16bp片段(数字存储)与数据库中每个菌种指纹片段(数字存储)交集的绝对数量，(2)交集占该菌种指纹数据库中指纹片段数的比例。根据(1)(2)的值结合发明前文步骤4中公式计算所有菌种的Score值进行排序，排名第一的菌种就是待鉴定菌株所属菌种。这步中，如已确定待鉴定菌为醋酸钙-鲍曼不动杆菌复合群菌种，则只需要计算数据库中复合群内菌种的评分Score值进行排序即可。

6、鉴定结果分析

(1)醋酸钙-鲍曼不动杆菌复合群中菌种鉴定

本发明对鲍曼不动杆菌(Acinetobacter baumannii)、醋酸钙不动杆菌(Acinetobacter calcoaceticus)、皮特不动杆菌(Acinetobacter pittii)和医院不动杆菌(Acinetobacter nosocomialis)4个菌种进行鉴定。对某一个待鉴定菌株而言，根据本方法鉴定得到的菌种名称与NCBI Taxonomy数据库得到的菌种名称一致的，视为鉴定正确，否则视为鉴定错误。

图2显示了菌种鉴定准确率与原始测序reads能映射到菌种参考基因组上比例mapping ratio的关系，利用bwa(http://bio-bwa.sourceforge.net/)工具将测序reads比对(mapping)到菌种的参考基因组上，根据工具计算得到mapping到参考基因组上reads占全部reads的比例，即mapping ratio。可以看出整个复合群的鉴定准确率很高，即使mapping ratio很低时，也超过90％，而当mapping ratio超过80％时，在测试数据上准确率达到100％。分菌种来看，当mapping ratio小于50％时，鲍曼不动杆菌(Acinetobacterbaumannii)、医院不动杆菌(Acinetobacter nosocomialis)和皮特不动杆菌(Acinetobacter pittii)鉴定准确率较高，而醋酸钙不动杆菌(Acinetobactercalcoaceticus)准确率较低；但当mapping ratio超过80％后，4个菌种的鉴定准确率在测试数据集上都达到100％。

由上述结果可以看出，当测序reads质量较高时，本方法对于醋酸钙-鲍曼不动杆菌复合群中菌种能够以很高的准确率有效区分鉴定。

(2)鲍曼不动杆菌其它菌种鉴定结果

为了测试本方法对于不动杆菌属其它菌种的鉴定能力，获取了其它几个菌种的测序reads进行鉴定，包括还可为溶血不动杆菌(Acinetobacter_haemolyticus)、约氏不动杆菌(Acinetobacter_johnsonii)、琼氏不动杆菌(Acinetobacter_junii)、嗜油不动杆菌(Acinetobacter_oleivorans)、申氏不动杆菌(Acinetobacter_schindleri)或土壤不动杆菌(Acinetobacter_soli)共6个菌种,鉴定结果如表1所示，显示本方法对于不动杆菌属的其它菌种鉴定同样适用。

表1不动杆菌属其它6个菌种鉴定情况

3、对于其它常见菌种的鉴定结果

对于其它菌种，获取测序数据后测试结果如图3，结果表明对其它74种菌种鉴定准确率同样非常高，一些菌种甚至可以达到100％。有些菌种由于数据质量问题导致鉴定准确率有所下降，但仍超过95％。根据醋酸钙-鲍曼不动杆菌复合群的例子可以推断当数据质量提升时，这些菌种的鉴定准确率也会提高。

从实施例可以看出，用本发明的方法可以准确鉴定醋酸钙-鲍曼不动杆菌复合群中菌种，当数据质量好时，结果非常可靠准确；此外可以鉴别不动杆菌属的其它菌种以及用于其它菌种的鉴定。

Claims (10)

1.一种基于无拼接组装WGS数据对待鉴定菌进行种属鉴定的方法，包括如下步骤：

(A)从菌基因组数据库中获取全部已测序的菌全基因组数据，建立菌株基因组信息和菌株分类学信息之间的关系；

(B)按照如下步骤构建菌种指纹特征数据库：

(b1)针对步骤(A)获得的每株菌的全基因组核酸序列，均以核酸序列长度为k、步长为1对全基因组核酸序列进行切分，切分后序列重复的片段只保留一个，将所得片段集合称为集合A；每个菌株得到一个所述集合A；组成所述集合A的片段称为片段A；

(b2)将四种碱基的字母表示分别转换为00、01、10、11这种数字存储，从而使得每个所述片段A都能转换为一个2k位的数字，即0至2^2k-1之间的一个数；

(b3)遍历所有0至2^2k-1之间的数字，对所有菌株的所述集合A中的所有所述片段A均一一做记录，记录各所述片段A所对应的2k位数字以及对应的种属信息，筛选出仅记录1个菌种的所述片段A对应的2k位数字，按照菌种分类存储，即得所述菌种指纹特征数据库。

(C)将待鉴定菌株的全基因组测序所得每一条测序读长均按照步骤(b1)的方式进行切分，得到所述待鉴定菌株的所述集合A；按照步骤(b2)的方式将所述待鉴定测菌株的所述集合A中的所有所述片段A均转换为2k位数字；然后与步骤(b3)所得所述菌种指纹特征数据库进行比对，根据比对结果计算S_core分值，得分最高的菌种视为所述待鉴定菌株所属菌种；

所述待鉴定菌株用于表示所述片段A的2k位数字与所述菌种指纹特征数据库中某一菌种A用于表示所述片段A的2k位数字，两者交集的绝对数量越多，则对于所述待鉴定菌株而言，所述菌种A的得分越高；

所述待鉴定菌株用于表示所述片段A的2k位数字与所述菌种指纹特征数据库中某一菌种A用于表示所述片段A的2k位数字，两者交集占所述菌种指纹特征数据库中所述菌种A用于表示所述片段A的2k位数字的数量比例越高，则对于所述待鉴定菌株而言，所述菌种A的得分越高。

2.根据权利要求1所述的方法，其特征在于：步骤(C)中，所述S_core分值按照如下公式计算获得：

Score＝α*normalization(N)+(1-α)*normalization(P)

N＝card(S_x∩S_i)

M＝card(S_i)

其中，S_x代表所述待鉴定测菌株用于表示所述片段A的2k位数字的集合；S_i代表所述菌种指纹特征数据库中某一菌种用于表示所述片段A的2k位数字的集合；α是加权系数；card表示求集合里的元素个数；N代表所述待鉴定测菌株用于表示所述片段A的2k位数字的集合与所述菌种指纹特征数据库中某一菌种用于表示所述片段A的2k位数字的集合的交集中的元素个数；M代表所述菌种指纹特征数据库中某一菌种用于表示所述片段A的2k位数字的集合中的元素个数；

N，P需要归一化处理为：

N_min，N_max，P_min，P_max是与所述菌种指纹特征数据库中全部菌种计算后得到的N，P的最小值和最大值；

进一步地，所述α为0.48。

3.根据权利要求1或2所述的方法，其特征在于：在所述方法中，步骤(A)按照如下进行：从NCBI的Genome数据库获取全部基因组状态Status为Complete Genome的菌全基因组序列数据；从NCBI的Taxonomy数据库获取生物分类元数据信息，根据TaxID建立菌株基因组信息和菌株分类学信息之间的关系。

4.根据权利要求1-3中任一所述的方法，其特征在于：步骤(A)中，所述菌为细菌。

5.根据权利要求1-3中任一所述的方法，其特征在于：所述待鉴定菌为细菌。

6.根据权利要求5所述的方法，其特征在于：所述细菌为鲍曼不动杆菌。

7.根据权利要求6所述的方法，其特征在于：所述鲍曼不动杆菌为隶属于醋酸钙-鲍曼不动杆菌复合群中的鲍曼不动杆菌。

8.一种构建菌种指纹特征数据库的方法，包括权利要求1-7中所述方法的步骤(A)和步骤(B)。

9.利用权利要求8所述方法构建得到的菌种指纹特征数据库。

10.权利要求9所述菌种指纹特征数据库在基于无拼接组装WGS数据对待鉴定菌进行种属鉴定中的应用。