CN114420212A - 一种大肠杆菌菌株鉴定方法和系统 - Google Patents
一种大肠杆菌菌株鉴定方法和系统 Download PDFInfo
- Publication number
- CN114420212A CN114420212A CN202210100336.XA CN202210100336A CN114420212A CN 114420212 A CN114420212 A CN 114420212A CN 202210100336 A CN202210100336 A CN 202210100336A CN 114420212 A CN114420212 A CN 114420212A
- Authority
- CN
- China
- Prior art keywords
- escherichia coli
- strain
- strains
- classification
- data
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 241000588724 Escherichia coli Species 0.000 title claims abstract description 239
- 238000000034 method Methods 0.000 title claims abstract description 58
- 108090000623 proteins and genes Proteins 0.000 claims abstract description 114
- 150000007523 nucleic acids Chemical group 0.000 claims abstract description 55
- 238000003908 quality control method Methods 0.000 claims abstract description 42
- 108091028043 Nucleic acid sequence Proteins 0.000 claims abstract description 37
- 238000012163 sequencing technique Methods 0.000 claims abstract description 34
- 241000894006 Bacteria Species 0.000 claims description 39
- 108091026890 Coding region Proteins 0.000 claims description 24
- 238000012545 processing Methods 0.000 claims description 10
- 238000009795 derivation Methods 0.000 claims description 6
- 238000003860 storage Methods 0.000 claims description 6
- 238000004422 calculation algorithm Methods 0.000 claims description 5
- 238000004891 communication Methods 0.000 claims description 3
- 238000003205 genotyping method Methods 0.000 abstract description 7
- 230000001580 bacterial effect Effects 0.000 abstract description 6
- 238000011160 research Methods 0.000 abstract description 6
- 238000004519 manufacturing process Methods 0.000 abstract description 5
- 239000003814 drug Substances 0.000 abstract description 2
- 238000005516 engineering process Methods 0.000 abstract 1
- 238000012360 testing method Methods 0.000 description 15
- 238000004458 analytical method Methods 0.000 description 11
- 239000012634 fragment Substances 0.000 description 8
- 238000010586 diagram Methods 0.000 description 6
- 238000012795 verification Methods 0.000 description 6
- 238000010276 construction Methods 0.000 description 5
- 238000012258 culturing Methods 0.000 description 5
- 238000001914 filtration Methods 0.000 description 5
- 238000007481 next generation sequencing Methods 0.000 description 4
- 102000004169 proteins and genes Human genes 0.000 description 4
- 238000004088 simulation Methods 0.000 description 4
- 230000009286 beneficial effect Effects 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 244000005700 microbiome Species 0.000 description 3
- 241000282414 Homo sapiens Species 0.000 description 2
- 238000012408 PCR amplification Methods 0.000 description 2
- OPTASPLRGRRNAP-UHFFFAOYSA-N cytosine Chemical compound NC=1C=CNC(=O)N=1 OPTASPLRGRRNAP-UHFFFAOYSA-N 0.000 description 2
- 238000007405 data analysis Methods 0.000 description 2
- 230000006870 function Effects 0.000 description 2
- UYTPUPDQBNUYGX-UHFFFAOYSA-N guanine Chemical compound O=C1NC(N)=NC2=C1N=CN2 UYTPUPDQBNUYGX-UHFFFAOYSA-N 0.000 description 2
- 238000013507 mapping Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000007781 pre-processing Methods 0.000 description 2
- 230000000717 retained effect Effects 0.000 description 2
- 241000894007 species Species 0.000 description 2
- 102000016911 Deoxyribonucleases Human genes 0.000 description 1
- 108010053770 Deoxyribonucleases Proteins 0.000 description 1
- 241001198387 Escherichia coli BL21(DE3) Species 0.000 description 1
- 241000901842 Escherichia coli W Species 0.000 description 1
- 241001465754 Metazoa Species 0.000 description 1
- 108700026244 Open Reading Frames Proteins 0.000 description 1
- 208000025174 PANDAS Diseases 0.000 description 1
- 208000021155 Paediatric autoimmune neuropsychiatric disorders associated with streptococcal infection Diseases 0.000 description 1
- 240000000220 Panda oleosa Species 0.000 description 1
- 235000016496 Panda oleosa Nutrition 0.000 description 1
- 108091005804 Peptidases Proteins 0.000 description 1
- 239000004365 Protease Substances 0.000 description 1
- 102100037486 Reverse transcriptase/ribonuclease H Human genes 0.000 description 1
- 241000258125 Strongylocentrotus Species 0.000 description 1
- 238000010367 cloning Methods 0.000 description 1
- 238000011109 contamination Methods 0.000 description 1
- 239000002537 cosmetic Substances 0.000 description 1
- 229940104302 cytosine Drugs 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 201000010099 disease Diseases 0.000 description 1
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 1
- 230000008030 elimination Effects 0.000 description 1
- 238000003379 elimination reaction Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 235000013305 food Nutrition 0.000 description 1
- 210000001035 gastrointestinal tract Anatomy 0.000 description 1
- 238000001502 gel electrophoresis Methods 0.000 description 1
- 238000012224 gene deletion Methods 0.000 description 1
- 230000002068 genetic effect Effects 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000000968 intestinal effect Effects 0.000 description 1
- 239000007788 liquid Substances 0.000 description 1
- 230000004807 localization Effects 0.000 description 1
- 210000004072 lung Anatomy 0.000 description 1
- 238000004949 mass spectrometry Methods 0.000 description 1
- 230000004060 metabolic process Effects 0.000 description 1
- 239000002207 metabolite Substances 0.000 description 1
- 210000000653 nervous system Anatomy 0.000 description 1
- 231100000252 nontoxic Toxicity 0.000 description 1
- 230000003000 nontoxic effect Effects 0.000 description 1
- 238000004806 packaging method and process Methods 0.000 description 1
- 238000012856 packing Methods 0.000 description 1
- 238000003909 pattern recognition Methods 0.000 description 1
- 239000013612 plasmid Substances 0.000 description 1
- 230000000241 respiratory effect Effects 0.000 description 1
- 238000012216 screening Methods 0.000 description 1
- 230000035945 sensitivity Effects 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- 239000003053 toxin Substances 0.000 description 1
- 231100000765 toxin Toxicity 0.000 description 1
- 108700012359 toxins Proteins 0.000 description 1
- 230000002485 urinary effect Effects 0.000 description 1
- 238000010200 validation analysis Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
- G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
- G16B40/00—ICT specially adapted for biostatistics; ICT specially adapted for bioinformatics-related machine learning or data mining, e.g. knowledge discovery or pattern finding
Landscapes
- Life Sciences & Earth Sciences (AREA)
- Physics & Mathematics (AREA)
- Health & Medical Sciences (AREA)
- Medical Informatics (AREA)
- Engineering & Computer Science (AREA)
- General Health & Medical Sciences (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Biophysics (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Biotechnology (AREA)
- Evolutionary Biology (AREA)
- Theoretical Computer Science (AREA)
- Artificial Intelligence (AREA)
- Chemical & Material Sciences (AREA)
- Analytical Chemistry (AREA)
- Proteomics, Peptides & Aminoacids (AREA)
- Bioethics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Data Mining & Analysis (AREA)
- Databases & Information Systems (AREA)
- Epidemiology (AREA)
- Evolutionary Computation (AREA)
- Public Health (AREA)
- Software Systems (AREA)
- Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
Abstract
本发明公开了一种大肠杆菌菌株鉴定方法和系统。该鉴定方法是一种用于细菌生物体(具体为大肠杆菌)的检测与识别的方法，基于待测大肠杆菌菌株二代测序获得的核酸序列数据，对核酸序列数据进行数据质控、数据比对和数据组装后，在构建的全基因组多位点序列分型数据库中比对查找到基因序列分型最接近的菌株，得到鉴定结果。本发明基于大肠杆菌全基因组多位点序列基因分型技术，提出了一种全新的大肠杆菌菌株鉴定方法和系统，能满足生物、医药、农业等多个领域生产和科研中更全面、复杂的大肠杆菌菌株鉴定需求。
Description
技术领域
本发明属于细菌生物体的检测与识别和测序领域，涉及一种大肠杆菌菌株鉴定方法和系统。
背景技术
大肠杆菌是一种革兰氏阴性的直杆菌。大肠杆菌是一种兼性厌氧微生物，能够进行呼吸代谢和发酵代谢。大肠杆菌分为多种菌株，一些菌株可以作为肠道微生物与人类形成有益的共生关系，而还有一些菌株进入人体后会产生毒素，引起肠道、泌尿系统、肺部和神经系统等部位的疾病。除此之外大肠杆菌也是重要的工程微生物，被广泛应用于化工、食品、生物医药、动物饲料和化妆品等多个领域。不同大肠杆菌工程菌株有截然不同的作用，例如DH5α菌株存在DNA酶缺陷，有利于保存质粒、克隆基因，但是该菌株容易降解蛋白质，不适合作为表达菌株。而BL21(DE3)菌株能够高效表达T7启动子驱动的外源基因，且存在蛋白酶基因缺失不容易降解蛋白质，适合用于非毒性蛋白质的表达。生产和科研中使用错误的大肠杆菌菌株会较大地影响效率，甚至会使结果偏离预期。对生产和科研中使用的大肠杆菌菌株进行菌株鉴定可以在质量控制环节中起到重要的作用。
目前大肠杆菌菌株鉴定的方法有：
培养分离法(见专利CN111235075A)，将样本在特定选择培养基上培养、挑选，然后对培养物进行特定基因的PCR扩增，根据凝胶电泳成像判断是否某种菌株阳性。
蛋白模式识别法(见专利CN109884160A)，对菌株进行培养，并对培养物进行质谱分析，根据结果中代谢物的模式识别特定类别的大肠杆菌菌株。
PCR法(见专利CN110982917A)，对菌株进行培养，获得可做PCR模板的菌液，使用特异性引物对特征序列片段进行PCR扩增，根据PCR产物条带情况鉴定是否某种菌株阳性。
这些大肠杆菌菌株鉴定方法的局限性在于：
1、这些鉴定方法均需要对大肠杆菌菌株进行培养，而菌株培养需要花费较多的时间，且容易在培养过程中产生污染。
2、这些鉴定方法都只能鉴定特定一类菌株，在生产和科研领域用到的菌株种类繁多。
3、只能鉴定菌株的一种或少数几种基因特征，而大肠杆菌菌株间的基因差异非常复杂，大肠杆菌核心基因组约为2000个基因，而泛基因组有18000多个基因(通常一种特定菌株包含4000多个基因)，这些方法难以覆盖复杂的鉴定需求。
4、这些鉴定方法都无法给出全面的基因序列分型信息。
发明内容
本发明基于大肠杆菌全基因组多位点序列基因序列分型，提出了一种全新的大肠杆菌菌株鉴定方法和系统，能满足生物、医药、农业等多个领域生产和科研中更全面、复杂的大肠杆菌菌株鉴定需求。
一方面，本发明公开了一种大肠杆菌菌株鉴定方法，根据待测大肠杆菌菌株二代测序的核酸序列数据，进行数据质控、数据比对和数据组装后，在构建的全基因组多位点序列分型数据库中比对查找到基因序列分型最接近的菌株，得到鉴定结果；
所述全基因组多位点序列分型数据库通过以下步骤构建获得：
S1、大肠杆菌菌株核酸序列获取：从NCBI获取大肠杆菌菌株核酸序列，得到fasta文件；
S2、大肠杆菌工程菌知识库建立：收集现有的大肠杆菌工程菌信息，建立所述大肠杆菌工程菌知识库；
S3、大肠杆菌菌株分类：使用mash程序分析所述步骤S1获得的fasta文件，计算所有目标菌株核酸序列两两之间的序列差异性；使用CL层次聚类算法将所有目标菌株根据核酸序列划分为N个大肠杆菌菌株分类；N取大于0的整数；
S4、参考基因组选取：对于每1个所述步骤S3获得的大肠杆菌菌株分类，计算分类中所有菌株与同类菌株的平均序列差异性；对分类中所有大肠杆菌菌株按与同类菌株的平均序列差异性从小到大排列，选取与分类内部所有菌株核酸序列平均序列差异性最小的菌株核酸序列作为该分类的参考基因组，从而得到大肠杆菌菌株参考基因组；所述大肠杆菌菌株参考基因组的数量为N个，对应N个大肠杆菌菌株分类；
S5、全基因组多位点序列分型靶基因选取：对于各个所述大肠杆菌菌株分类对应的所述大肠杆菌菌株参考基因组，从NCBI获取所有基因编码区序列；在同一个所述大肠杆菌菌株分类中，对于有相同序列的基因编码区，仅保留其中一个；将过滤后的基因作为该所述大肠杆菌菌株分类的全基因组多位点序列分型靶基因；
S6、菌株库去冗余：对于各个所述大肠杆菌菌株分类中的所有菌株，保留所述大肠杆菌工程菌知识库包含的菌株，对于所述大肠杆菌工程菌知识库之外的菌株进行去冗余处理：如果多个菌株之间序列差异性小于M，则仅保留与其他同类菌株平均序列差异性最小的菌株；M为0.00005-0.0005(即M具体值根据实践需求在0.00005-0.0005范围内调整。优选地，M＝0.0001)；
S7、构建得到所述全基因组多位点序列分型数据库：对于每1个所述大肠杆菌菌株分类，利用blat或blast将该分类中去冗余后的菌株一一比对到该分类的所述大肠杆菌菌株参考基因组的基因编码区上，得到各个菌株的全基因组多位点序列分型靶基因序列分型，构建完成所述全基因组多位点序列分型数据库。
在一些实施方案中，大肠杆菌工程菌知识库中含有128种大肠杆菌工程菌的信息。大肠杆菌工程菌信息包括NCBI编号、ATCC编号、菌株衍生关系等信息。
在一些实施方案中，菌株库去冗余后最终共保留1570个大肠杆菌菌株。
在一些实施方案中，包括以下步骤：
A1、数据质控；
A2、数据比对：通过数据比对得到样本在各个所述大肠杆菌菌株分类中的大肠杆菌菌株参考基因组的比对率、对大肠杆菌菌株参考基因组的基因组覆盖率和对大肠杆菌菌株参考基因组的基因组覆盖深度；
A3、数据组装；
A4、全基因组多位点序列分型靶基因检索：调用blat或blast程序将所述步骤A3拼接后的样本contigs比对到各个所述大肠杆菌菌株分类的大肠杆菌菌株参考基因组的基因编码区上，计算所述拼接后的样本contigs包含的各个分类的全基因组多位点序列分型靶基因数量、靶基因序列分型；
A5、样本大肠杆菌一级分类；在同一个样本中，依次按全基因组多位点序列分型靶基因数量、大肠杆菌菌株参考基因组的对比率和对大肠杆菌菌株参考基因组的基因组覆盖率对所述步骤A4比对的所述大肠杆菌菌株分类降序排列，取排名第一的大肠杆菌菌株分类为该样本所属的目标分类；
A6、样本大肠杆菌二级分类：在所述步骤A5找到的目标分类中，使用pyMLST程序的wgMLST流程线计算样本基因序列分型与所述目标分类的全基因组多位点序列分型数据库中各个菌株基因序列分型的差异情况，得到与所述样本基因序列分型相似度最高的菌株为该样本中大肠杆菌匹配的菌株种类。
在一些实施方案中，所述步骤A1具体为：数据质控：利用fastp软件对待测大肠杆菌菌株二代测序的核酸序列数据(比如待测大肠杆菌菌株二代测序下机数据，或者模拟二代测序软件ART生成的待测大肠杆菌菌株测序核酸序列数据)进行质控，去除低质量序列，得到质控后的二代测序数据。
在一些实施方案中，所述步骤A2具体为：数据比对：利用bowtie2软件将所述步骤A1得到的质控后的二代测序数据比对到各个所述大肠杆菌菌株分类的参考基因组上，计算样本数据在各个所述大肠杆菌菌株分类中的大肠杆菌菌株参考基因组的比对率、对大肠杆菌菌株参考基因组的基因组覆盖率和对大肠杆菌菌株参考基因组的基因组覆盖深度。
在一些实施方案中，所述步骤A3具体为：数据组装：基于Bruijn图原理的de novo数据组装，利用SPAdes软件将质控后的二代测序数据拼接为contigs(基因组长片段)。
在一些实施方案中，所述步骤A6还包括：再次用所述步骤A2数据比对将所述步骤A1数据质控后的样本数据比对到其匹配的菌株基因组序列上，从而获得精确的比对率、覆盖率、覆盖深度数据。
在一些实施方案中，还包括步骤A7、大肠杆菌工程菌信息注释：在所述大肠杆菌工程菌知识库中，检索样本中大肠杆菌匹配的菌株种类；对包含在所述大肠杆菌工程菌知识库中的菌株，对该菌株进行信息注释；所述信息注释包括NCBI编号、ATCC编号和/或菌株衍生关系。
在一些实施方案中，所述步骤S7具体为：构建得到所述全基因组多位点序列分型数据库：对于每1个所述大肠杆菌菌株分类，利用pyMLST程序的wgMLST流程线逐一构建每个分类的初始全基因组多位点序列分型数据库；利用blat或blast将该分类中去冗余后的菌株一一比对到该分类的所述大肠杆菌菌株参考基因组的基因编码区上，得到各个菌株的全基因组多位点序列分型靶基因序列分型，加入到对应分类的所述初始全基因组多位点序列分型数据库，构建完成所述全基因组多位点序列分型数据库。
在一些实施方案中，M＝0.0001，N＝10；所述大肠杆菌菌株参考基因组为GCF_014162235.1，GCF_001020945.2，GCF_000026265.1，GCF_002157245.1，GCF_004924275.1，GCF_001900735.1，GCF_008926085.1，GCF_001677475.2，GCF_009931435.1，GCF_013305705.1。
另一方面，本发明还公开了一种大肠杆菌菌株鉴定系统，包括如权利要求1所述的全基因组多位点序列分型数据库，以及处理模块、输入模块和显示模块；所述处理模块分别与所述全基因组多位点序列分型数据库、所述输入模块、所述显示模块通讯连接；
所述处理模块用于根据所述输入模块传入的待测大肠杆菌菌株二代测序的核酸序列数据，进行数据质控、数据比对和数据组装，并在所述全基因组多位点序列分型数据库中比对查找到序列最接近的菌株，得到鉴定结果，并传送给所述显示模块。
在一些实施方案中，所述处理模块中存储有可执行指令，所述可执行指令执行时实现如上所述的大肠杆菌菌株鉴定方法。
在一些实施方案中，所述大肠杆菌菌株参考基因组为GCF_014162235.1，GCF_001020945.2，GCF_000026265.1，GCF_002157245.1，GCF_004924275.1，GCF_001900735.1，GCF_008926085.1，GCF_001677475.2，GCF_009931435.1，GCF_013305705.1。
进一步地，M＝0.0001，N＝10。
第三方面，本发明公开了一种计算机可读存储介质，所述存储介质中存储有可执行指令，所述可执行指令执行时实现如上所述的大肠杆菌菌株鉴定方法。
第四方面，本发明公开了一种终端，所述终端包括：
存储器，用于存储可执行指令；
处理器，用于执行所述存储器中存储的可执行指令时，实现如上所述的大肠杆菌菌株鉴定方法。
本发明的有益效果有：
1、可鉴定多达128种工程大肠杆菌，及1570种NCBI大肠杆菌菌株。
2、基于多个大肠杆菌菌株参考基因组，分别鉴定每个大肠杆菌菌株参考基因组多达4000多个位点的基因序列分型。避免单一大肠杆菌菌株参考基因组造成的比对偏倚。
3、分析对象为通过二代测序或模拟二代测序软件ART获取的疑似大肠杆菌菌株样本的核酸序列数据，故无需对样本进行额外的纯化和培养。
4、结果包含全面的基因序列分型信息。
以下将结合附图对本发明的构思、具体结构及产生的技术效果作进一步说明，以充分地了解本发明的目的、特征和效果。
附图说明
图1是本发明所涉及的大肠杆菌菌株鉴定方法中的全基因组多位点序列分型数据库构建的流程框图。
图2是本发明所涉及的大肠杆菌菌株鉴定方法流程框图。
具体实施方式
为了使发明实现的技术手段、创造特征、达成目的和功效易于明白了解，下结合具体图示，进一步阐述本发明。但本发明不仅限于以下实施的案例。
须知，本说明书所附图式所绘示的结构、比例、大小等，均仅用以配合说明书所揭示的内容，以供熟悉此技术的人士了解与阅读，并非用以限定本发明可实施的限定条件，故不具技术上的实质意义，任何结构的修饰、比例关系的改变或大小的调整，在不影响本发明所能产生的功效及所能达成的目的下，均应仍落在本发明所揭示的技术内容得能涵盖的范围内。
实施例1
图1示出了本发明所涉及的大肠杆菌菌株鉴定方法中的全基因组多位点序列分型数据库构建的流程框架，具体阐述如下：
(1)大肠杆菌菌株核酸序列获取
通过NCBI获取已知的大肠杆菌菌株核酸序列，具体如下：
从NCBI官方网站下载refseq中的所有细菌核酸序列列表文件assembly_summary.txt(网址：https://ftp.ncbi.nlm.nih.gov/genomes/refseq/bacteria/assembly_summary.txt)，筛选版本最新(即version_status为latest)、拼接完整(即assembly_level为Complete Genome)、物种为大肠杆菌(即species_taxid为562)的所有核酸序列。
对于同一菌株的多个参考基因组(即infraspecific_name相同的多个核酸序列)，从列表文件assembly_summary.txt的ftp_path栏获取ftp数据路径，从ftp数据路径下载序列组装情况统计文件(即*_assembly_stats.txt)，获取各个核酸序列的组装质量标准N50(scaffold_N50)，将其按组装质量标准N50从大到小排序，只保留组装质量标准N50最大的一个核酸序列。
对于筛选后的所有核酸序列，从列表文件assembly_summary.txt的ftp_path栏获取ftp数据路径，从ftp数据路径下载序列对应的fasta文件(*_genomic.fna.gz)和文件指纹验证文件(md5checksums.txt)，使用文件指纹验证文件验证下载到fasta文件的完整性，重新下载验证失败的fasta文件，直到获取所有目标核酸序列的完整fasta文件。
(2)大肠杆菌工程菌知识库建立
收集已知的工程用大肠菌(共128种)，查阅文献资料收集NCBI编号、ATCC编号、菌株衍生关系等信息，建立工程用大肠菌知识库。
(3)菌株分类
使用mash程序(https://github.com/marbl/mash)分析上述步骤获得的fasta文件，计算所有目标核酸序列两两之间的序列差异性，具体使用方法举例如下：
mash sketch-p 24-l fa_list.txt-o build_prefix&&mash dist-p 24-tbuild_prefix.msh build_prefix.msh>build_prefix.dist
上述命令为使用mash软件计算多个fasta文件两两之间的遗传距离的具体方法，其中各个词段含义如下：
mash sketch：利用mash软件计算序列差异性的预处理指令；
-p 24：设置程序使用的CPU核数为24核；
-l fa_list.txt：将需要计算两两序列差异性的所有fasta文件本地路径逐行列在fa_list.txt文件中，输入mash程序进行分析；
-o build_prefix：输出文件前缀
mash dist：利用mash软件计算序列差异性的指令；
-t build_prefix.msh build_prefix.msh：将预处理步骤生成的文件输入mash程序，计算两两序列差异性；
>build_prefix.dist：将mash软件计算出的核酸序列两两之间的序列差异性数据输出到build_prefix.dist文件中。
基于上述步骤获得的所有目标核酸序列两两之间的序列差异性，使用CL层次聚类算法(complete linkage hierarchical clustering)将所有目标核酸序列划分为10类。
(4)参考基因组选取
对于每一个上述步骤获得的大肠杆菌菌株分类，计算分类中所有菌株与同类菌株的平均序列差异性，计算方法举例如下：
sum(inner_distance_list)/len(inner_distance_list)
该例通过python指令实现，具体代码解释如下：
inner_distance_list：分类中的某一菌株与同类菌株两两之间的序列差异性数值列表，从上一步mash程序结果中挑选出来，以python list变量形式存储；
sum：计算分类中的某一菌株与同类所有菌株两两之间的序列差异性数值总和；
len：计算分类中的某一菌株与同类菌株两两之间的序列差异性数值列表的长度，等于同类菌株个数。
下一步，对分类中所有大肠杆菌菌株按与同类菌株的平均序列差异性从小到大排列，选取与分类内部所有菌株核酸序列平均序列差异性最小的菌株核酸序列作为该分类的参考基因组。在本实施例中，得到的10个参考基因组为：GCF_014162235.1，GCF_001020945.2，GCF_000026265.1，GCF_002157245.1，GCF_004924275.1，GCF_001900735.1，GCF_008926085.1，GCF_001677475.2，GCF_009931435.1，GCF_013305705.1。
(5)全基因组多位点序列分型靶基因选取
对于各分类中选取的参考基因组序列，从列表文件assembly_summary.txt的ftp_path栏获取ftp数据路径，从ftp数据路径下载序列对应的编码区序列文件(*_cds_from_genomic.fna.gz)，使用文件指纹验证文件验证下载到编码区序列文件的完整性，重新下载验证失败的文件，直到获取所有参考基因组的编码区序列文件。
对于一个参考基因组的所有编码区序列，检查所有编码区序列重复情况，对于有相同序列的基因编码区，仅保留其中一个。将过滤后的所有编码区对应的基因作为该分类的全基因组多位点序列分型靶基因。
(6)菌株库去冗余
对于各个菌株类中的所有菌株，保留大肠杆菌工程菌知识库包含的所有菌株。对于工程菌知识库之外的菌株进行去冗余处理：基于之前步骤获得的所有目标核酸序列两两之间的序列差异性，使用CL层次聚类算法(complete linkage hierarchical clustering)将分类中菌株核酸序列以序列差异性0.0001为阈值进行再聚类。并按之前步骤提到的方法计算聚类中所有菌株与同聚类菌株的平均序列差异性，对于每个聚类，仅保留内部平均序列差异性最小的菌株作为代表菌株，从而过滤掉序列高度相似的冗余菌株。最终共保留1570个菌株。
(7)构建全基因组多位点序列分型数据库
对于(3)获得的各个菌株分类，利用pyMLST程序(https://github.com/bvalot/pyMLST)的wgMLST流程线逐一构建全基因组多位点序列分型数据库，并利用blat(类blast比对工具，与blast比有更快的比对速度，适合用于基因定位)将该分类中去冗余后的菌株一一比对到该分类的参考基因组基因编码区上，分析各个菌株的全基因组多位点序列分型靶基因序列分型，构建方法举例如下：
wgMLST create-f target.db reference_cds.fna&&wgMLST add target.dbstrain_genome.fna
其中各个词段含义如下：
wgMLST create：利用pyMLST程序基于参考基因组的编码区序列建立初始全基因组多位点序列分型数据库的指令；
-f：强制覆盖已有文件；
target.db：构建的全基因组多位点序列分型数据库存储路径；
reference_cds.fna：参考基因组编码区序列文件(解压后的*_cds_from_genomic.fna.gz)；
wgMLST add：调用blat将菌株比对到参考基因组基因编码区上，分析菌株靶基因序列分型的指令；
strain_genome.fna：特定菌株的核酸序列fasta文件。
实施例2
图2示出了本发明所涉及的大肠杆菌菌株鉴定方法的大体步骤，具体阐述如下：
流程(大体步骤)可以根据二代测序下机数据类型和需求不同进行灵活调整，本实施例以单端测序序列、分析使用24核CPU、纯菌株样本为例进行具体阐述。
(1)数据质控
利用fastp程序(https://github.com/OpenGene/fastp)对二代测序下机数据进行质控，去除低质量序列。具体分析方式举例如下：
fastp--thread 24--n_base_limit 3--in1 input.fastq--out1 trim.fastq
其中各个词段含义如下：
fastp：fastp质控程序；
--thread 24：设置程序使用的CPU核数为24核；
--n_base_limit 3：过滤掉N多于3个的reads；
--in1 input.fastq：输入下机核酸序列数据，路径为input.fastq；
--out1 trim.fastq：输出质控后的有效数据，路径为trim.fastq。
(2)数据比对
利用bowtie2程序(https://github.com/BenLangmead/bowtie2)将质控后的有效数据比对到10个大肠杆菌菌株分组的10个参考基因组上，计算各个参考基因组的比对率。并利用samtools软件计算各个参考基因组的覆盖率、覆盖深度。具体分析方式举例如下：
bowtie2-p 24--very-sensitive-x target.db.index-U trim.fastq-Soutput.sam&&\
samtools view-bS-@24output.sam|samtools sort-@24-o output.bam&&\
samtools coverage-o output.coverage_table.txt output.bam
其中各个词段含义如下：
bowtie2：比对程序；
-p 24：为比对程序配置使用的CPU核数为24核；
--very-sensitive：配置比对灵敏度为非常灵敏；
-x target.db.index：输入参考基因组(bowtie2 index形式)；
-U trim.fastq：将质控后的有效数据输入比对程序，与参考基因组进行比对；
-S output.sam：设置比对后的sam文件保存文件路径；
samtools view-bS：将sam文件转为bam文件的指令；
-@24：设置samtools程序使用的CPU核数为24核；
samtools sort：对bam文件进行排序的指令；
-o output.bam：设置转换并排序后的bam文件存储路径；
samtools coverage：计算参考基因组的覆盖率、覆盖深度的指令；
-o output.coverage_table.txt：将参考基因组的覆盖率、覆盖深度输出到output.coverage_table.txt文件。
(3)数据组装
利用spades程序(https://github.com/ablab/spades)基于Bruijn图原理将质控后的有效数据序列组装为contigs(基因组长片段)。具体方法和参数举例如下：
spades-s trim.fastq-o assembly_dir--isolate-t 24
其中各个词段含义如下：
spades：基因组组装程序
-s trim.fastq：将质控后的有效数据输入组装程序，基于Bruijn图原理进行denovo数据组装；
-o assembly_dir：输出组装中间文件及组装完成的contigs序列文件到assembly_dir文件夹；
--isolate：分析纯菌样本时配置spades以高纯度模式进行拼接，优化性能；
-t 24：为组装程序程序配置使用的CPU核数为24核。
(4)全基因组多位点序列分型靶基因检索
使用pyMLST程序(https://github.com/bvalot/pyMLST)的wgMLST流程线，调用blat(https://genome.ucsc.edu/cgi-bin/hgBlat？hgsid＝1223130093_U0R3ONzCRuXJCwwHRHJ3yIJhKwSP&command＝start)程序将拼接后的样本contigs比对到各个大肠杆菌菌株分类的参考基因组的基因编码区上，计算样本contigs包含的各个分类的全基因组多位点序列分型靶基因数量和靶基因序列分型。具体方法和参数举例如下：
wgmlst add target.db contigs.fa
其中各个词段含义如下：
wgMLST add：调用blat将样本contigs比对到参考基因组基因编码区上，分析样本contigs靶基因序列分型的指令；
target.db：已构建的全基因组多位点序列分型数据库；
contigs.fa：数据组装获得的样本contigs序列文件。
(5)样本大肠杆菌一级分类
综合比较样本数据在各个大肠杆菌菌株分类的参考基因组的比对率、对参考基因组的基因组覆盖率和覆盖的全基因组多位点序列分型靶基因数量。这三个数值越高，说明样本与对应分组菌株核酸序列相似程度越高。基于该原理找出样本大肠杆菌在10个分类中匹配的分类。具体方法举例如下：
该例通过python语言和pandas模块实现，具体代码解释如下：
identi_data：汇总后的各个参考基因组的比对率、对参考基因组的基因组覆盖率和覆盖的全基因组多位点序列分型靶基因数量等数据；
sort_values：对数据进行排序的指令；
by＝["sample_id","loci_number","mapping_rate","coverage"]：在同一个样本ID中依次按全基因组多位点序列分型靶基因数量、分类的参考基因组的比对率、对参考基因组的基因组覆盖率进行排序；
ascending＝[True,False,False,False]：不同样本按字符顺序，在同样本中，依次按全基因组多位点序列分型靶基因数量、分类的参考基因组的比对率、对参考基因组的基因组覆盖率降序排列；
inplace＝True：替换原表格
在同样本中，依次按全基因组多位点序列分型靶基因数量、分类的参考基因组的比对率、对参考基因组的基因组覆盖率降序排列指：将样本数据比对到各个大肠杆菌菌株分类的参考基因组上，比对分析之后，会产生比对率(mapping_rate)、覆盖率(coverage)和覆盖的靶基因数量(loci_number)三个数值，同一个样本在不同分类中这3个数据是不同的。所以可以排序。首先，按覆盖的靶基因数量由高到低排序(降序)，同一个样本在不同菌株分类中靶基因数量一般是不同的，如果都不同那排完这步算法就结束了，选出排序第一的为目标分类。但如果排在最前面的覆盖的靶基因数量这个值出现了相同，就对覆盖的靶基因数量相同的菌株分类进一步按比对率由高到低排序。同理，前两个值都相同的话再按覆盖率排序。排序后各个样本中排名第一的大肠杆菌菌株分类即为目标分类。
(6)样本大肠杆菌二级分类
在上步找到的样本大肠杆菌在10个分类中匹配的分类中，使用pyMLST程序(https://github.com/bvalot/pyMLST)的wgMLST流程线计算样本基因序列分型与分类数据库中各个菌株基因序列分型的差异情况。具体方法举例如下：
wgmlst distance-m 1-k-o distance.txt target.db
其中各个词段含义如下：
wgMLST distance：计算样本基因序列分型与分类数据库中各个菌株基因序列分型的差异情况的指令；
-m 1：过滤掉在分组所有菌株中均未覆盖的靶基因；
-k：过滤掉在所有菌株中基因序列分型一样的靶基因；
target.db：匹配到的分类的全基因组多位点序列分型数据库；
distance.txt：计算获得的差异结果输出文件。
最终与样本基因序列分型相似度最高的菌株即为样本中大肠杆菌匹配的菌株种类。最后由于前面进行步骤(2)时，是将质控后的样本数据比对到10个参考基因组中的，其得到的比对率、覆盖率、覆盖深度数据精度有限。故需要再次用步骤(2)同样的方法将质控后的样本数据比对到其匹配的菌株基因组序列上，从而获取精确的比对率、覆盖率、覆盖深度数据。
(7)大肠杆菌工程菌信息注释
在数据库构建步骤构建的大肠杆菌工程菌知识库中，检索样本中大肠杆菌匹配的菌株种类，如该菌株包含在工程菌知识库中，进一步对该菌株的NCBI编号、ATCC编号、菌株衍生关系等信息进行注释。
实施例3
利用模拟数据进行验证测试
一、数据生成
利用模拟二代测序软件ART(https://www.niehs.nih.gov/research/resources/software/biostatistics/art/index.cfm)生成菌株测序核酸序列数据，具体使用方法举例如下：
art_illumina-ss NS50-i GCF_016864475.1_ASM1686447v1_genomic.fna-l 75-f 100-o./GCF_016864475.1_ASM1686447v1
上述命令为使用软件“ART”生成模拟二代测序数据的具体方法，其中各个词段含义如下：
art_illumina：利用ART软件进行模拟illumina测序平台的二代测序数据生成
-ss NS50：模拟的设备型号为illumina NextSeq500 v2
-i GCF_016864475.1_ASM1686447v1_genomic.fna：对参考基因组GCF_016864475.1_ASM1686447v1_genomic.fna进行模拟测序
-l 75：生成序列读长为75bp
-f 100：生成平均覆盖深度为100x的测序数据
-o./GCF_016864475.1_ASM1686447v1：模拟生成的fastq文件输出路径前缀为./GCF_016864475.1_ASM1686447v1
以上是利用软件ART对NCBI上的大肠杆菌菌株参考基因组GCF_016864475.1进行模拟二代测序，生成测序核酸序列数据。除GCF_016864475.1外，本实施例中，还利用软件ART对NCBI上的大肠杆菌菌株参考基因组GCF_003367885.1、GCF_013167015.1、GCF_002899475.1、GCF_001276585.2进行模拟二代测序，生成测序核酸序列数据。
以上5个基因组对应的具体测试大肠杆菌菌株信息如表1所示。
表1、测试大肠杆菌菌株信息
二、分析测试
利用实施例1构建的全基因组多位点序列分型数据库，采用实施例2给出的数据分析方法(或称为大肠杆菌菌株鉴定方法)，对生成的模拟测序结果进行分析测试。具体数据分析步骤记录如下：
(1)数据质控
利用fastp程序(https://github.com/OpenGene/fastp)对生成的模拟测序结果(在其他实施例中也可以为真实样本的二代测序下机数据)进行质控，去除低质量序列(具体方法、参数同实施例2数据质控步骤)。数据质控结果见表2。
表2、数据质控结果表
| 样本编号 | NCBI基因组编号 | 总序列数 | 总数据量 | 数据Q30率 | GC比例 |
| 1 | GCF_013167015.1 | 6,080,300 | 456,022,500 | 93.44％ | 50.83％ |
| 2 | GCF_003367885.1 | 6,133,000 | 459,975,000 | 93.44％ | 50.78％ |
| 3 | GCF_016864475.1 | 6,329,600 | 474,720,000 | 93.45％ | 50.87％ |
| 4 | GCF_008868305.1 | 6,674,900 | 500,617,500 | 93.45％ | 50.82％ |
| 5 | GCF_008033295.1 | 6,631,500 | 497,362,500 | 93.45％ | 50.55％ |
表2为各示例样本(即测试大肠杆菌菌株)原始二代测序数据的质量结果数据，各列数据解释如下：
样本编号：与表1(测试大肠杆菌菌株信息表)对应的核酸序列数据样本编号；
NCBI基因组编号：示例样本对应的NCBI大肠杆菌基因组编号；
总序列数：各示例样本原始二代测序数据包含的序列总数；
总数据量：各示例样本原始二代测序数据包含的数据量总数(碱基总数)；
数据Q30率：各示例样本原始二代测序数据Phred数值大于30的碱基数量占总碱基数量的百分比；
GC比例：各示例样本原始二代测序数据中鸟嘌呤和胞嘧啶占总碱基数的比例。
(2)数据比对
利用bowtie2程序(https://github.com/BenLangmead/bowtie2)将质控后的有效数据比对到实施例1步骤(3)(菌株分类)和步骤(4)(参考基因组选取)获得的10个大肠杆菌菌株分组的10个参考基因组上，计算各个参考基因组的比对率。并利用samtools软件计算各个参考基因组的覆盖率、覆盖深度(具体方法、参数同实施例2数据比对步骤)。数据比对结果见表3。
表3、数据比对结果表
表3为各示例样本质控后的数据与各大肠杆菌菌株分类参考基因组比对结果数据，各列数据解释如下：
样本编号：与表1(测试大肠杆菌菌株信息表)对应的核酸序列数据样本编号；
大肠杆菌菌株分类：实施例1步骤(3)(菌株分类)获得的大肠杆菌菌株的10个大类；
菌株分类参考基因组：实施例1步骤(4)(参考基因组选取)获得的与大肠杆菌菌株大类对应的参考基因组，样本数据的比对目标；
比对率：示例样本质控后的数据与特定大类参考基因组比对获得的比对率，即能比对上参考基因组的序列数占总序列数的比例；
基因组覆盖率：示例样本质控后的数据与特定大类参考基因组比对获得的基因组覆盖率，即参考基因组上被样本序列覆盖的碱基长度占参考基因组总碱基长度的比例；
基因组覆盖深度：示例样本质控后的数据与特定大类参考基因组比对获得的基因组覆盖深度，即参考基因组上被样本序列覆盖的碱基平均被覆盖的序列数。
(3)数据组装
利用spades程序(https://github.com/ablab/spades)基于Bruijn图原理将质控后的有效数据序列组装为contigs(基因组长片段)(具体方法、参数同实施例2数据组装步骤)。数据组装结果见表4。
表4、数据组装结果表
| 样本编号 | 组装总长 | N50 | 样本数据组装率 |
| 1 | 4,478,932 | 54,733 | 0.9997 |
| 2 | 4,516,487 | 60,768 | 0.9998 |
| 3 | 4,657,799 | 75,344 | 0.9997 |
| 4 | 4,934,112 | 108,259 | 0.9996 |
| 5 | 4,859,161 | 57,767 | 0.9997 |
表4为各示例样本质控后的数据组装结果数据，各列数据解释如下：
样本编号：与表1(测试大肠杆菌菌株信息表)对应的核酸序列数据样本编号；
组装总长：各示例样本二代短片段组装为基因组长片段后的基因组长片段总长度(碱基数)；
N50：样本序列组装后会获得一系列不同长度的基因组长片段，将所有的基因组长片段按照从长到短进行排序，并逐一相加，当相加的长度达到组装总长的一半时，最后一个加上的基因组长片段的长度即为N50。为基因组组装的标准质控数据。
样本数据组装率：成功组装为基因组长片段的序列数占总序列数比例。
(4)全基因组多位点序列分型靶基因检索
使用pyMLST程序(https://github.com/bvalot/pyMLST)的wgMLST流程线，调用blat(https://genome.ucsc.edu/cgi-bin/hgBlat？hgsid＝1223130093_U0R3ONzCRuXJCwwHRHJ3yIJhKwSP&command＝start)程序将拼接后的样本contigs比对到各个大肠杆菌菌株分类的参考基因组的基因编码区上，计算样本contigs包含的各个分类的全基因组多位点序列分型靶基因数量和靶基因序列分型。(具体方法、参数同实施例2全基因组多位点序列分型靶基因检索步骤)。全基因组多位点序列分型靶基因检索结果见表5。
表5、全基因组多位点序列分型靶基因检索结果表
表5为各示例样本全基因组多位点序列分型靶基因检索结果数据，各列数据解释如下：
样本编号：与表1(测试大肠杆菌菌株信息表)对应的核酸序列数据样本编号；
大肠杆菌菌株分类：实施例1步骤(3)(菌株分类)获得的大肠杆菌菌株的10个大类；
菌株分类参考基因组：实施例1步骤(4)(参考基因组选取)获得的与大肠杆菌菌株大类对应的参考基因组，将样本contigs比对到其基因编码区上检索基因并比对基因序列分型；
覆盖靶基因数：样本contigs能比对到的参考基因组靶基因数。
(5)样本大肠杆菌一级分类
综合比较样本数据在各个大肠杆菌菌株分类的参考基因组的比对率、对参考基因组的基因组覆盖率和覆盖的全基因组多位点序列分型靶基因数量。这三个数值越高，说明样本与对应分组菌株核酸序列相似程度越高。基于该原理找出样本大肠杆菌在10个分类中匹配的分类(具体方法、参数同实施例2样本大肠杆菌一级分类步骤)。样本大肠杆菌一级分类结果见表6。
表6、样本大肠杆菌一级分类结果表
表6为各示例样本大肠杆菌一级分类结果数据，各列数据解释如下：
样本编号：与表1(测试大肠杆菌菌株信息表)对应的核酸序列数据样本编号；
最匹配的大肠杆菌菌株分类：大肠杆菌菌株的10个大类中与样本序列最相似的菌株分类；
比对率：示例样本质控后的数据与最匹配的大肠杆菌菌株分类参考基因组比对获得的比对率，即能比对上参考基因组的序列数占总序列数的比例；
基因组覆盖率：示例样本质控后的数据与最匹配的大肠杆菌菌株分类参考基因组比对获得的基因组覆盖率，即参考基因组上被样本序列覆盖的碱基长度占参考基因组总碱基长度的比例；
基因组覆盖深度：示例样本质控后的数据与最匹配的大肠杆菌菌株分类参考基因组比对获得的基因组覆盖深度，即参考基因组上被样本序列覆盖的碱基平均被覆盖的序列数；
覆盖靶基因数：样本contigs能比对到的最匹配的大肠杆菌菌株分类参考基因组基因数。
(6)样本大肠杆菌二级分类
在上步找到的样本大肠杆菌最匹配的大肠杆菌菌株分类中，使用pyMLST程序(https://github.com/bvalot/pyMLST)的wgMLST流程线计算样本基因序列分型与分类数据库中各个菌株基因序列分型的差异情况，最终与样本基因序列分型相似度最高的菌株即为样本中大肠杆菌匹配的菌株种类。(具体方法、参数同实施例2样本大肠杆菌二级分类步骤)。样本大肠杆菌二级分类结果见表7。
表7、样本大肠杆菌二级分类结果表
表7为各示例样本大肠杆菌二级分类结果数据，各列数据解释如下：
样本编号：与表1(测试大肠杆菌菌株信息表)对应的核酸序列数据样本编号；
最匹配的大肠杆菌菌株分类：大肠杆菌菌株的10个大类中与样本序列最相似的菌株分类；
最小靶基因序列分型差异数：样本contigs基因序列分型与一级分类得到的大肠杆菌菌株分类中基因序列分型最接近的菌株比较，基因序列分型不同的基因个数；
差异最小菌株：一级分类得到的大肠杆菌菌株分类中与样本contigs基因序列分型最接近的菌株名。
最后将质控后的样本数据比对到差异最小菌株基因组序列上，获取精确的比对率、覆盖率、覆盖深度数据，结果见表8。
表8、样本与差异最小菌株比对数据表
表8为样本与差异最小菌株比对数据，各列数据解释如下：
样本编号：与表1(测试大肠杆菌菌株信息表)对应的核酸序列数据样本编号；
差异最小菌株：一级分类得到的大肠杆菌菌株分类中与样本contigs基因序列分型最接近的菌株名；
比对率：示例样本质控后的数据与差异最小菌株参考基因组比对获得的比对率，即能比对上参考基因组的序列数占总序列数的比例；
基因组覆盖率：示例样本质控后的数据与差异最小菌株参考基因组比对获得的基因组覆盖率，即参考基因组上被样本序列覆盖的碱基长度占参考基因组总碱基长度的比例；
基因组覆盖深度：示例样本质控后的数据与差异最小菌株参考基因组比对获得的基因组覆盖深度，即参考基因组上被样本序列覆盖的碱基平均被覆盖的序列数。
(7)大肠杆菌工程菌信息注释
在实施例1步骤(2)(大肠杆菌工程菌知识库建立)构建的大肠杆菌工程菌知识库中，检索样本中大肠杆菌匹配的菌株种类，如该菌株包含在工程菌知识库中，进一步对该菌株的NCBI编号、ATCC编号、菌株衍生关系等信息进行注释。大肠杆菌工程菌信息注释结果见表9。
表9、大肠杆菌工程菌信息注释结果表
表9为大肠杆菌工程菌信息注释结果，各列数据解释如下：
样本编号：与表1(测试大肠杆菌菌株信息表)对应的核酸序列数据样本编号；
差异最小菌株：样本contigs基因序列分型与最匹配的大肠杆菌菌株分类中基因序列分型最接近的菌株名；
差异最小菌株ATCC编号：大肠杆菌工程菌知识库中差异最小菌株对应的ATCC编号；
差异最小菌株NCBI编号：大肠杆菌工程菌知识库中差异最小菌株对应的NCBI编号；
差异最小菌株衍生组：大肠杆菌工程菌知识库中差异最小菌株对应的菌株衍生组。
三、分析结果
分析结果表明本发明提供的大肠杆菌菌株鉴定方法正确地对大肠杆菌菌株进行了鉴定，获得了与实际菌株相符的鉴定结果，具体鉴定结果见表10和11。
表10、测试数据鉴定结果——菌株识别
| 样本编号 | 菌株实际名 | 本发明鉴定结果 | 鉴定是否与实际相符 |
| 1 | BL21(DE3) | Escherichia coli BL21(DE3) | 是 |
| 2 | C600 | Escherichia coli K-12C600 | 是 |
| 3 | Crooks | Escherichia coli Crooks | 是 |
| 4 | W | Escherichia coli W | 是 |
| 5 | J53 pMG223 | Escherichia coli K-12J53 | 是 |
表11、测试数据鉴定结果——工程菌信息注释
实施例4
被测大肠杆菌的测序核酸序列数据除了像实施例3那样来自模拟二代测序软件ART外，还可以来自真实的二代测序仪，即可以为二代测序下机数据。其他过程及分析同实施例3，效果也类似实施3，这里不再赘述。
实施例5
本实施例涉及大肠杆菌菌株鉴定系统，其包括实施例1的全基因组多位点序列分型数据库，以及处理模块、输入模块和显示模块。处理模块分别与全基因组多位点序列分型数据库、输入模块、显示模块通讯连接。处理模块用于根据输入模块传入的待测大肠杆菌菌株二代测序的原始数据，进行数据质控、数据比对和数据组装，并在全基因组多位点序列分型数据库中比对查找到序列最接近的菌株，得到鉴定结果，并传送给显示模块。
实施例6
本实施例涉及计算机可读存储介质。该存储介质中存储有可执行指令。该可执行指令执行时实现实施例2的大肠杆菌菌株鉴定方法。
实施例7
本实施例涉及终端。该终端包括：存储器，用于存储可执行指令；还包括处理器，用于执行该存储器中存储的可执行指令时，实现实施例2的大肠杆菌菌株鉴定方法。
以上详细描述了本发明的较佳具体实施例。应当理解，本领域的普通技术无需创造性劳动就可以根据本发明的构思作出诸多修改和变化。因此，凡本技术领域中技术人员依本发明的构思在现有技术的基础上通过逻辑分析、推理或者有限的实验可以得到的技术方案，皆应在由权利要求书所确定的保护范围内。
Claims (10)
1.一种大肠杆菌菌株鉴定方法，其特征在于，根据待测大肠杆菌菌株二代测序的核酸序列数据，进行数据质控、数据比对和数据组装后，在构建的全基因组多位点序列分型数据库中比对查找到基因序列分型最接近的菌株，得到鉴定结果；
所述全基因组多位点序列分型数据库通过以下步骤构建获得：
S1、大肠杆菌菌株核酸序列获取：从NCBI获取大肠杆菌菌株核酸序列，得到fasta文件；
S2、大肠杆菌工程菌知识库建立：收集现有的大肠杆菌工程菌信息，建立所述大肠杆菌工程菌知识库；
S3、大肠杆菌菌株分类：使用mash程序分析所述步骤S1获得的fasta文件，计算所有目标菌株核酸序列两两之间的序列差异性；使用CL层次聚类算法将所有目标菌株根据核酸序列差异性划分为N个大肠杆菌菌株分类；N取大于0的整数；
S4、参考基因组选取：对于每1个所述步骤S3获得的大肠杆菌菌株分类，计算分类中所有菌株与同类菌株的平均序列差异性；对分类中所有大肠杆菌菌株按与同类菌株的平均序列差异性从小到大排列，选取与分类内部所有菌株核酸序列平均序列差异性最小的菌株核酸序列作为该分类的参考基因组，从而得到大肠杆菌菌株参考基因组；所述大肠杆菌菌株参考基因组的数量为N个，对应N个大肠杆菌菌株分类；
S5、全基因组多位点序列分型靶基因选取：对于各个所述大肠杆菌菌株分类对应的所述大肠杆菌菌株参考基因组，从NCBI获取所有基因编码区序列；在同一个所述大肠杆菌菌株分类中，对于有相同序列的基因编码区，仅保留其中一个；将过滤后的基因作为该所述大肠杆菌菌株分类的全基因组多位点序列分型靶基因；
S6、菌株库去冗余：对于各个所述大肠杆菌菌株分类中的所有菌株，保留所述大肠杆菌工程菌知识库包含的菌株，对于所述大肠杆菌工程菌知识库之外的菌株进行去冗余处理：如果多个菌株之间序列差异性小于M，则仅保留与其他同类菌株平均序列差异性最小的菌株；M为0.00005-0.0005；
S7、构建得到所述全基因组多位点序列分型数据库：对于每1个所述大肠杆菌菌株分类，利用blat或blast将该分类中去冗余后的菌株一一比对到该分类的所述大肠杆菌菌株参考基因组的基因编码区上，得到各个菌株的全基因组多位点序列分型靶基因序列分型，构建完成所述全基因组多位点序列分型数据库。
2.如权利要求1所述的大肠杆菌菌株鉴定方法，其特征在于，包括以下步骤：
A1、数据质控；
A2、数据比对：通过数据比对得到样本在各个所述大肠杆菌菌株分类中的大肠杆菌菌株参考基因组的比对率、对大肠杆菌菌株参考基因组的基因组覆盖率和对大肠杆菌菌株参考基因组的基因组覆盖深度；
A3、数据组装；
A4、全基因组多位点序列分型靶基因检索：调用blat或blast程序将所述步骤A3拼接后的样本contigs比对到各个所述大肠杆菌菌株分类的大肠杆菌菌株参考基因组的基因编码区上，计算所述拼接后的样本contigs包含的各个分类的全基因组多位点序列分型靶基因数量、靶基因序列分型；
A5、样本大肠杆菌一级分类；在同一个样本中，依次按全基因组多位点序列分型靶基因数量、大肠杆菌菌株参考基因组的对比率和对大肠杆菌菌株参考基因组的基因组覆盖率对所述步骤A4比对的所述大肠杆菌菌株分类降序排列，取排名第一的大肠杆菌菌株分类为该样本所属的目标分类；
A6、样本大肠杆菌二级分类：在所述步骤A5找到的目标分类中，使用pyMLST程序的wgMLST流程线计算样本基因序列分型与所述目标分类的全基因组多位点序列分型数据库中各个菌株基因序列分型的差异情况，得到与所述样本基因序列分型相似度最高的菌株为该样本中大肠杆菌匹配的菌株种类。
3.如权利要求2所述的大肠杆菌菌株鉴定方法，其特征在于，所述步骤A6还包括：再次用所述步骤A2数据比对将所述步骤A1数据质控后的样本数据比对到其匹配的菌株基因组序列上，从而获得精确的比对率、覆盖率、覆盖深度数据。
4.如权利要求2所述的大肠杆菌菌株鉴定方法，其特征在于，还包括步骤A7、大肠杆菌工程菌信息注释：在所述大肠杆菌工程菌知识库中，检索样本中大肠杆菌匹配的菌株种类；对包含在所述大肠杆菌工程菌知识库中的菌株，对该菌株进行信息注释；所述信息注释包括NCBI编号、ATCC编号和/或菌株衍生关系。
5.如权利要求1所述的大肠杆菌菌株鉴定方法，其特征在于，所述步骤S7具体为：构建得到所述全基因组多位点序列分型数据库：对于每1个所述大肠杆菌菌株分类，利用pyMLST程序的wgMLST流程线逐一构建每个分类的初始全基因组多位点序列分型数据库；利用blat或blast将该分类中去冗余后的菌株一一比对到该分类的所述大肠杆菌菌株参考基因组的基因编码区上，得到各个菌株的全基因组多位点序列分型靶基因序列分型，加入到对应分类的所述初始全基因组多位点序列分型数据库，构建完成所述全基因组多位点序列分型数据库。
6.如权利要求1所述的大肠杆菌菌株鉴定方法，其特征在于，M＝0.0001，N＝10；所述大肠杆菌菌株参考基因组为GCF_014162235.1，GCF_001020945.2，GCF_000026265.1，GCF_002157245.1，GCF_004924275.1，GCF_001900735.1，GCF_008926085.1，GCF_001677475.2，GCF_009931435.1，GCF_013305705.1。
7.一种大肠杆菌菌株鉴定系统，其特征在于，包括如权利要求1所述的全基因组多位点序列分型数据库，以及处理模块、输入模块和显示模块；所述处理模块分别与所述全基因组多位点序列分型数据库、所述输入模块、所述显示模块通讯连接；
所述处理模块用于根据所述输入模块传入的待测大肠杆菌菌株二代测序获得的核酸序列数据，进行数据质控、数据比对和数据组装，并在所述全基因组多位点序列分型数据库中比对查找到基因序列分型最接近的菌株，得到鉴定结果，并传送给所述显示模块。
8.如权利要求7所述的大肠杆菌菌株鉴定系统，其特征在于，所述大肠杆菌菌株参考基因组为GCF_014162235.1，GCF_001020945.2，GCF_000026265.1，GCF_002157245.1，GCF_004924275.1，GCF_001900735.1，GCF_008926085.1，GCF_001677475.2，GCF_009931435.1，GCF_013305705.1。
9.一种计算机可读存储介质，其特征在于，所述存储介质中存储有可执行指令，所述可执行指令执行时实现如权利要求1-6任一项所述的大肠杆菌菌株鉴定方法。
10.一种终端，其特征在于，所述终端包括：
存储器，用于存储可执行指令；
处理器，用于执行所述存储器中存储的可执行指令时，实现如权利要求1-6任一项所述的大肠杆菌菌株鉴定方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210100336.XA CN114420212B (zh) | 2022-01-27 | 2022-01-27 | 一种大肠杆菌菌株鉴定方法和系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210100336.XA CN114420212B (zh) | 2022-01-27 | 2022-01-27 | 一种大肠杆菌菌株鉴定方法和系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114420212A true CN114420212A (zh) | 2022-04-29 |
| CN114420212B CN114420212B (zh) | 2022-10-21 |
Family
ID=81279797
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210100336.XA Active CN114420212B (zh) | 2022-01-27 | 2022-01-27 | 一种大肠杆菌菌株鉴定方法和系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114420212B (zh) |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115064215A (zh) * | 2022-08-18 | 2022-09-16 | 北京大学人民医院 | 一种通过相似度进行菌株溯源及属性鉴定的方法 |
| CN115083527A (zh) * | 2022-08-18 | 2022-09-20 | 北京大学人民医院 | 一种聚类泛基因组数据库构建方法 |
| CN117037912A (zh) * | 2023-09-13 | 2023-11-10 | 青岛极智医学检验实验室有限公司 | 一种泛基因组的构建方法、终端设备及存储介质 |
| CN117746980A (zh) * | 2023-12-18 | 2024-03-22 | 广州凯普医学检验所有限公司 | 一种流感病毒自动化快速分型方法、装置、设备及介质 |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2016124600A1 (en) * | 2015-02-02 | 2016-08-11 | Applied Maths | Method of typing nucleic acid or amino acid sequences based on sequence analysis |
| CN110714088A (zh) * | 2019-10-16 | 2020-01-21 | 北京出入境检验检疫局检验检疫技术中心 | 一种基于gMLST技术的沙门氏菌溯源分型的方法及应用 |
| WO2020055076A1 (ko) * | 2018-09-10 | 2020-03-19 | 주식회사 조앤김지노믹스 | 유산균 동정용 참조서열 제조방법 및 이를 이용한 유산균 동정방법 |
| CN111276185A (zh) * | 2020-02-18 | 2020-06-12 | 上海桑格信息技术有限公司 | 一种基于二代高通量测序的微生物鉴定分析系统及装置 |
| CN113373208A (zh) * | 2021-07-14 | 2021-09-10 | 上海序祯达生物科技有限公司 | 一种基于二代测序的人类白细胞抗原分型系统和方法 |
| CN114144843A (zh) * | 2019-07-12 | 2022-03-04 | 生物梅里埃公司 | 流行病学鉴定和监测细菌爆发的方法 |
-
2022
- 2022-01-27 CN CN202210100336.XA patent/CN114420212B/zh active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2016124600A1 (en) * | 2015-02-02 | 2016-08-11 | Applied Maths | Method of typing nucleic acid or amino acid sequences based on sequence analysis |
| WO2020055076A1 (ko) * | 2018-09-10 | 2020-03-19 | 주식회사 조앤김지노믹스 | 유산균 동정용 참조서열 제조방법 및 이를 이용한 유산균 동정방법 |
| CN114144843A (zh) * | 2019-07-12 | 2022-03-04 | 生物梅里埃公司 | 流行病学鉴定和监测细菌爆发的方法 |
| CN110714088A (zh) * | 2019-10-16 | 2020-01-21 | 北京出入境检验检疫局检验检疫技术中心 | 一种基于gMLST技术的沙门氏菌溯源分型的方法及应用 |
| CN111276185A (zh) * | 2020-02-18 | 2020-06-12 | 上海桑格信息技术有限公司 | 一种基于二代高通量测序的微生物鉴定分析系统及装置 |
| CN113373208A (zh) * | 2021-07-14 | 2021-09-10 | 上海序祯达生物科技有限公司 | 一种基于二代测序的人类白细胞抗原分型系统和方法 |
Non-Patent Citations (2)
| Title |
|---|
| D.BABENKO,ET AL.: "wgMLST as a standardized tool for assessing the quality of genome assembly data", 《INTERNATIONAL JOURNAL OF INFECTIOUS DISEASES》 * |
| 王玮玉: "致羔羊脑膜炎型大肠杆菌NMGCF-19菌株的分离鉴定及全基因组测序分析", 《CNKI硕士电子期刊》 * |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115064215A (zh) * | 2022-08-18 | 2022-09-16 | 北京大学人民医院 | 一种通过相似度进行菌株溯源及属性鉴定的方法 |
| CN115083527A (zh) * | 2022-08-18 | 2022-09-20 | 北京大学人民医院 | 一种聚类泛基因组数据库构建方法 |
| CN115064215B (zh) * | 2022-08-18 | 2023-10-24 | 北京大学人民医院 | 一种通过相似度进行菌株溯源及属性鉴定的方法 |
| CN117037912A (zh) * | 2023-09-13 | 2023-11-10 | 青岛极智医学检验实验室有限公司 | 一种泛基因组的构建方法、终端设备及存储介质 |
| CN117746980A (zh) * | 2023-12-18 | 2024-03-22 | 广州凯普医学检验所有限公司 | 一种流感病毒自动化快速分型方法、装置、设备及介质 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114420212B (zh) | 2022-10-21 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN114420212B (zh) | 一种大肠杆菌菌株鉴定方法和系统 | |
| Steinegger et al. | Protein-level assembly increases protein sequence recovery from metagenomic samples manyfold | |
| US10347365B2 (en) | Systems and methods for visualizing a pattern in a dataset | |
| Yan et al. | DeepTE: a computational method for de novo classification of transposons with convolutional neural network | |
| Meinicke | UProC: tools for ultra-fast protein domain classification | |
| Stranneheim et al. | Classification of DNA sequences using Bloom filters | |
| US20200294628A1 (en) | Creation or use of anchor-based data structures for sample-derived characteristic determination | |
| CN113744807B (zh) | 一种基于宏基因组学的病原微生物检测方法及装置 | |
| Dündar et al. | Introduction to differential gene expression analysis using RNA-seq | |
| CN115083521B (zh) | 一种单细胞转录组测序数据中肿瘤细胞类群的鉴定方法及系统 | |
| CN108710784A (zh) | 一种基因转录变异几率及变异方向的算法 | |
| CN114121160A (zh) | 一种检测样本中宏病毒组的方法和系统 | |
| CN115662516A (zh) | 一种基于二代测序技术的高通量预测噬菌体宿主的分析方法 | |
| EP2518656A1 (en) | Taxonomic classification system | |
| US7047137B1 (en) | Computer method and apparatus for uniform representation of genome sequences | |
| Utro et al. | Hierarchically labeled database indexing allows scalable characterization of microbiomes | |
| CN113380326B (zh) | 一种基于pam聚类算法的基因表达数据分析方法 | |
| CN114627964B (zh) | 一种基于多核学习预测增强子及其强度分类方法及分类设备 | |
| CN113744806B (zh) | 一种基于纳米孔测序仪的真菌测序数据鉴定方法 | |
| Hu et al. | A novel method for discovering local spatial clusters of genomic regions with functional relationships from DNA contact maps | |
| CN116153410B (zh) | 微生物基因组参考数据库及其构建方法和应用 | |
| Cai et al. | Application and research progress of machine learning in Bioinformatics | |
| CN114496089B (zh) | 一种病原微生物鉴定方法 | |
| CN113611355B (zh) | 基于氨基酸组成和蛋白质相互作用识别抗氧化蛋白方法 | |
| Neuhaus et al. | Transcription factor prediction using protein 3D structures |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| PE01 | Entry into force of the registration of the contract for pledge of patent right |
Denomination of invention: A method and system for identifying Escherichia coli strains Effective date of registration: 20231130 Granted publication date: 20221021 Pledgee: Industrial Bank Co.,Ltd. Shanghai Zhangyang Sub branch Pledgor: Shanghai xuzhenda Biotechnology Co.,Ltd. Registration number: Y2023310000791 |
|
| PE01 | Entry into force of the registration of the contract for pledge of patent right |