CN116597893B - 预测耐药基因-病原微生物归属的方法 - Google Patents
预测耐药基因-病原微生物归属的方法 Download PDFInfo
- Publication number
- CN116597893B CN116597893B CN202310702713.1A CN202310702713A CN116597893B CN 116597893 B CN116597893 B CN 116597893B CN 202310702713 A CN202310702713 A CN 202310702713A CN 116597893 B CN116597893 B CN 116597893B
- Authority
- CN
- China
- Prior art keywords
- genome
- drug
- kmer
- species
- drug resistance
- 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.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 39
- 244000005700 microbiome Species 0.000 title claims description 20
- 238000004902 predicting drug resistance Methods 0.000 title claims description 8
- 108090000623 proteins and genes Proteins 0.000 claims abstract description 207
- 206010059866 Drug resistance Diseases 0.000 claims abstract description 118
- 229940079593 drug Drugs 0.000 claims description 105
- 239000003814 drug Substances 0.000 claims description 105
- 241000894007 species Species 0.000 claims description 89
- 244000052769 pathogen Species 0.000 claims description 73
- 230000001717 pathogenic effect Effects 0.000 claims description 60
- 238000001514 detection method Methods 0.000 claims description 36
- 244000052616 bacterial pathogen Species 0.000 claims description 33
- 238000010276 construction Methods 0.000 claims description 21
- 238000012216 screening Methods 0.000 claims description 19
- 230000001580 bacterial effect Effects 0.000 claims description 18
- 238000001914 filtration Methods 0.000 claims description 18
- 238000004458 analytical method Methods 0.000 claims description 12
- 241000894006 Bacteria Species 0.000 claims description 9
- 238000004590 computer program Methods 0.000 claims description 6
- 244000000010 microbial pathogen Species 0.000 abstract description 5
- 238000012163 sequencing technique Methods 0.000 description 24
- 241000588747 Klebsiella pneumoniae Species 0.000 description 20
- 238000004364 calculation method Methods 0.000 description 10
- 238000012360 testing method Methods 0.000 description 10
- 208000015181 infectious disease Diseases 0.000 description 6
- 230000035945 sensitivity Effects 0.000 description 6
- 241000588724 Escherichia coli Species 0.000 description 5
- 230000002458 infectious effect Effects 0.000 description 5
- 238000005516 engineering process Methods 0.000 description 4
- 238000011160 research Methods 0.000 description 4
- 238000004088 simulation Methods 0.000 description 4
- 108700039887 Essential Genes Proteins 0.000 description 3
- 206010034133 Pathogen resistance Diseases 0.000 description 3
- 241000588626 Acinetobacter baumannii Species 0.000 description 2
- 101100026178 Caenorhabditis elegans egl-3 gene Proteins 0.000 description 2
- 230000003115 biocidal effect Effects 0.000 description 2
- 210000004369 blood Anatomy 0.000 description 2
- 239000008280 blood Substances 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000011156 evaluation Methods 0.000 description 2
- 238000012268 genome sequencing Methods 0.000 description 2
- 238000007481 next generation sequencing Methods 0.000 description 2
- 108020004707 nucleic acids Proteins 0.000 description 2
- 150000007523 nucleic acids Chemical class 0.000 description 2
- 102000039446 nucleic acids Human genes 0.000 description 2
- 101100456369 Aquifex aeolicus (strain VF5) mdh1 gene Proteins 0.000 description 1
- 208000003322 Coinfection Diseases 0.000 description 1
- 101100447530 Dictyostelium discoideum gpi gene Proteins 0.000 description 1
- 241000588697 Enterobacter cloacae Species 0.000 description 1
- 101150039774 GAPA1 gene Proteins 0.000 description 1
- 101000710013 Homo sapiens Reversion-inducing cysteine-rich protein with Kazal motifs Proteins 0.000 description 1
- 101000661807 Homo sapiens Suppressor of tumorigenicity 14 protein Proteins 0.000 description 1
- 101150058595 MDH gene Proteins 0.000 description 1
- 101100384788 Methanocaldococcus jannaschii (strain ATCC 43067 / DSM 2661 / JAL-1 / JCM 10045 / NBRC 100440) comC gene Proteins 0.000 description 1
- 101100038261 Methanococcus vannielii (strain ATCC 35089 / DSM 1224 / JCM 13029 / OCM 148 / SB) rpo2C gene Proteins 0.000 description 1
- 101100020705 Mycoplasma gallisepticum (strain R(low / passage 15 / clone 2)) ldh gene Proteins 0.000 description 1
- 206010036790 Productive cough Diseases 0.000 description 1
- 108010049395 Prokaryotic Initiation Factor-2 Proteins 0.000 description 1
- 241000589517 Pseudomonas aeruginosa Species 0.000 description 1
- 101100282114 Pseudomonas aeruginosa (strain UCBPP-PA14) gap2 gene Proteins 0.000 description 1
- 101100290490 Rattus norvegicus Mdh1 gene Proteins 0.000 description 1
- -1 ST147 Proteins 0.000 description 1
- 241000191967 Staphylococcus aureus Species 0.000 description 1
- 101100120969 Streptomyces coelicolor (strain ATCC BAA-471 / A3(2) / M145) pgi1 gene Proteins 0.000 description 1
- WKDDRNSBRWANNC-UHFFFAOYSA-N Thienamycin Natural products C1C(SCCN)=C(C(O)=O)N2C(=O)C(C(O)C)C21 WKDDRNSBRWANNC-UHFFFAOYSA-N 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000003766 bioinformatics method Methods 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 210000001175 cerebrospinal fluid Anatomy 0.000 description 1
- 239000003153 chemical reaction reagent Substances 0.000 description 1
- 229940121657 clinical drug Drugs 0.000 description 1
- 238000007418 data mining Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000003745 diagnosis Methods 0.000 description 1
- 201000010099 disease Diseases 0.000 description 1
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 1
- 238000001647 drug administration Methods 0.000 description 1
- 210000003608 fece Anatomy 0.000 description 1
- 238000005111 flow chemistry technique Methods 0.000 description 1
- 101150073818 gap gene Proteins 0.000 description 1
- 101150091570 gapA gene Proteins 0.000 description 1
- ZSKVGTPCRGIANV-ZXFLCMHBSA-N imipenem Chemical compound C1C(SCC\N=C\N)=C(C(O)=O)N2C(=O)[C@H]([C@H](O)C)[C@H]21 ZSKVGTPCRGIANV-ZXFLCMHBSA-N 0.000 description 1
- 229960002182 imipenem Drugs 0.000 description 1
- 108090000237 interleukin-24 Proteins 0.000 description 1
- 239000007788 liquid Substances 0.000 description 1
- 210000004072 lung Anatomy 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000005065 mining Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 101150053253 pgi gene Proteins 0.000 description 1
- 101150058164 phoE gene Proteins 0.000 description 1
- 238000013468 resource allocation Methods 0.000 description 1
- 230000000717 retained effect Effects 0.000 description 1
- 101150085857 rpo2 gene Proteins 0.000 description 1
- 101150090202 rpoB gene Proteins 0.000 description 1
- 238000002864 sequence alignment Methods 0.000 description 1
- 239000002689 soil Substances 0.000 description 1
- 210000003802 sputum Anatomy 0.000 description 1
- 208000024794 sputum Diseases 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000001502 supplementing effect Effects 0.000 description 1
- 101150080758 tonB gene Proteins 0.000 description 1
- 238000012549 training Methods 0.000 description 1
- 238000010200 validation analysis Methods 0.000 description 1
- PJVWKTKQMONHTI-UHFFFAOYSA-N warfarin Chemical compound OC=1C2=CC=CC=C2OC(=O)C=1C(CC(=O)C)C1=CC=CC=C1 PJVWKTKQMONHTI-UHFFFAOYSA-N 0.000 description 1
- 229960005080 warfarin Drugs 0.000 description 1
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
- G16B15/00—ICT specially adapted for analysing two-dimensional or three-dimensional molecular structures, e.g. structural or functional relations or structure alignment
- G16B15/30—Drug targeting using structural data; Docking or binding prediction
-
- 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
- G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
- G16B20/20—Allele or variant detection, e.g. single nucleotide polymorphism [SNP] detection
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
Landscapes
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Health & Medical Sciences (AREA)
- Physics & Mathematics (AREA)
- Chemical & Material Sciences (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Medical Informatics (AREA)
- General Health & Medical Sciences (AREA)
- Theoretical Computer Science (AREA)
- Biophysics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Biotechnology (AREA)
- Evolutionary Biology (AREA)
- Pharmacology & Pharmacy (AREA)
- Medicinal Chemistry (AREA)
- Crystallography & Structural Chemistry (AREA)
- Analytical Chemistry (AREA)
- Genetics & Genomics (AREA)
- Molecular Biology (AREA)
- Proteomics, Peptides & Aminoacids (AREA)
- Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
Abstract
本申请属于生物信息学领域，具体涉及一种预测耐药基因的病原微生物归属的方法，相比于传统方法，本申请准确度更高、可靠性更好。
Description
技术领域
本申请属于生物信息学技术领域，具体涉及一种预测耐药基因-病原微生物归属的方法。
技术背景
宏基因组测序(mNGS)是近十年发展的新兴的病原检测手段，相较于传统临床微生物的检测方法，mNGS不需要筛选得到环境中所有的病原体的纯培养物，就可以从少量的样本中一次性读取所有的病原体核酸，能够覆盖更广范围的病原微生物。随着测序技术发展以及测序成本降低，将mNGS测序技术用于临床标本的检测以预测临床病人感染病原的细菌耐药性相对于传统的细菌培养，可以快速指导疾病的针对性用药。目前，基于临床标本宏基因组测序(尤其是针对BALF、血液等高宿主核酸含量样本)用于感染病原病原检测的研究，很多是还关注在对病原菌种的鉴定性能上，至于进一步预测病原菌的耐药表型性能如何，相关研究较少。近期我们在journal of clinical microbiology上发表了一项基于mNGS预测鲍曼不动杆菌耐药表型的研究，结果证实了mNGS用于预测病原菌耐药表型的可靠性，但模型的性能也会受到耐药基因-物种归属准确性的较大影响。关于耐药基因-物种归属，针对粪便、土壤等微生物含量高的标本，有研究采用Contig/Reads bining等方法可以得到单菌基因组草图，以此同时可解决基因-物种来源归属问题，但是因为临床样本(如肺泡灌洗液、痰液、血液、脑脊液等)受宿主基因组污染较严重，往往微生物含量较低(往往低于5％)，进而常规测序数据量下(20M reads或5G测序数据量)，得到的目标病原微生物reads数量相对较少(细菌基因组覆盖很多不超过1X)，采用Contig/Reads bining等方式实现目标基因的物种来源归属不现实，需要开发一种方法或技术，专门针对临床标本mNGS可实现目标基因-病原微生物物种的准确归属。
不同于细菌纯培养分离株，临床标本中往往包含各种不同的细菌，且丰度也各不一样。在现有常规的20M左右reads数的测序量下，临床标本测得的细菌基因组覆盖很多不超过1X。在此测序数据量情况下，临床标本mNGS如何能够实现耐药病原菌的准确鉴定，除了明确耐药基因型与表型的一致性以及实现耐药基因的准确鉴定外，准确预测检出耐药基因的病原菌物种来源也尤为重要。针对宏基因组测序应用于临床样本进行耐药菌检测过程中因耐药基因-病原物种归属错误导致药敏结果预测不准的问题，本申请尝试提出一种新的策略准确鉴定耐药基因-病原微生物物种归属的方法。有鉴于此，提出本申请。
发明内容
为解决上述技术问题，本申请创新性的提出了一种基于与目标耐药基因相关联的物种特异性kmer序列检出的方法来鉴定耐药基因-病原微生物物种归属，并进一步提出了联合多策略来准确鉴定耐药基因-病原微生物物种归属的方法，包括：策略i)根据耐药基因的物种注释推断可能归属的物种；策略ii)根据耐药基因和物种检出的序列情况推算出来的耐药基因拷贝数是否落在正常范围内；策略iii)与目标耐药基因相关联的物种特异性kmer序列是否被同时检出。
具体的，本申请提出如下技术方案：
本申请首先提供一种用于预测耐药基因-病原微生物归属模型的构建方法，以及相应的预测耐药基因-病原微生物归属的方法，所述方法包括如下步骤：
构建与目标耐药基因相关联的物种特异性kmer序列库，将待预测的病原耐药基因序列与物种特异性kmer序列库进行比对，判断该耐药基因的物种来源；如果在序列库中检测到病原耐药基因相关联的物种特异性kmer，那么确认双方物种归属关系，否则拒绝双方物种归属关系。
进一步的，所述物种特异性kmer序列库的构建方法如下：
基于公共数据库获得目标病原菌株基因组kmer，筛选得到与各目标耐药基因阳性一致性高的序列；通过数据库比对注释，挑选出特异性比对到目标病原物种的kmer序列，以此作为物种特异性kmer序列库。
进一步的，所述公共数据库包括NCBI数据库；
进一步的，所述与各目标耐药基因阳性一致性高的序列为与各目标耐药基因阳性一致性PPV>0.90的序列，所述与各目标耐药基因阳性一致性PPV计算公式为：同时携带目标耐药基因和kmer的菌株数/携带目标耐药基因的所有菌株数。
更进一步的，所述物种特异性kmer序列库的构建方法具体如下：
a)针对特定目标病原菌，从NCBI genome数据库下载所有或足够量的病原菌株基因组；将所有基因组打碎为kmer序列，通过筛选得到与目标耐药基因的阳性检出一致性PPV>0.90的kmer序列；
b)将筛选出的kmer序列在NCBI数据库进行比对注释，保留仅特异性比对到该病原的kmer序列，以此作为物种特异性kmer序列库。
在一些优选的方式中，所述步骤a)具体包括如下步骤：
i.针对特定目标病原菌，从NCBI genome数据库下载所有或足够量的病原菌株基因组，对下载到的目标病原菌所有菌株基因组进行耐药基因检测分析，得到各菌携带的耐药基因信息，统计各重要耐药基因在所有菌株基因组上的发生率；
ii.计算目标病原基因组上的各区域发生率，筛选出较高发生率区域：根据细菌多位点序列分型MLST分型，在发生频次较高的多个MLST分型中各选择1个代表株，以多个代表株做为参考株，将剩余菌株基因组与参考株基因组比对，统计各个参考株的基因组上各区域发生率；以步骤i中统计得到的重要耐药基因发生率为阈值，筛选出高于阈值的各参考株基因组上的所有区域，作为较高发生率的区域；
iii.目标病原基因组序列kmer化，过滤掉位于基因组较高发生率区域上的kmer：针对每个菌株，将其基因组打碎为kmer序列，与步骤ii得到的基因组较高发生率区域进行比对和过滤，过滤掉位于基因组较高发生率区域上的kmer；
iv.针对过滤后剩余的各kmer序列，计算其在所有菌株中的发生率以及相对于目标耐药基因的菌株发生一致率PPV，挑选出PPV在0.90以上的kmer序列。
本申请还提供一种用于综合预测耐药基因-病原微生物归属模型的构建方法，包括如下步骤：
1)基于耐药基因的物种注释，推断可能的病原物种来源；
2)依据耐药基因和病原检出序列情况，评价耐药基因的拷贝数是否正常，保留拷贝数正常的基因-物种归属关系；
3)上述任一所述的构建方法。
进一步的，所述步骤1)具体为：将耐药基因与NCBI nt库进行比对，获得各耐药基因的所有物种注释信息，并与CARD耐药数据库信息整合，构建耐药基因的潜在来源物种信息库，
进一步的，所述步骤2)中，所述耐药基因的拷贝数是否正常是通过来自于NCBIgenome数据库的相应病原所有基因组进行耐药基因检测分析统计得到；
所述耐药基因的拷贝数的计算公式如下：
其中，CoverageARG为耐药基因的覆盖度，Cover_depthARG为耐药基因的覆盖深度，
Coveragegenpme为基因组的覆盖度，Cover_depEhgenome为基因组的覆盖深度。
更进一步的，所述步骤2)具体为：
基于NCBI genome数据库中的目标病原菌所有菌株基因组，将各菌株基因组contig序列直接与耐药基因数据库进行比对，保留identity>90％、coverage>60％、与NCBINDARO数据库注释记录比较测试性能大于0.90的结果；对于每个contig对齐区域，选择hit值最高的区域，同时计算检测到的耐药基因拷贝数(目标耐药基因的覆盖范围*覆盖范围的深度)/(目标基因组的覆盖范围/基因组覆盖区域的深度)，并统计目标病原菌各耐药基因的携带发生频次情况；然后，根据临床标本实际检出的耐药基因和病原菌检出reads情况，计算耐药基因的拷贝数，通过与以上统计得到的目标病原菌耐药基因发生频数范围进行比较，查看是否正常来判断物种耐药归属关系。
本申请还提供一种预测耐药基因-病原微生物归属的模型，其特征在于，包括如下模块：
1)kmer序列库构建模块：用于构建与目标耐药基因相关联的物种特异性kmer序列库；
2)序列比对模块：用于将病原耐药基因序列与物种特异性kmer序列库进行比对，判断耐药基因的物种来源；如果在序列库中检测到病原耐药基因相关联的物种特异性kmer，那么确认双方的物种归属关系，否则拒绝物种归属关系。
进一步的，所述kmer序列库序列库模块的构建方法如下：
基于公共数据库获得目标病原菌株基因组kmer，筛选得到与各目标耐药基因阳性一致性高的序列；通过数据库比对注释，挑选出特异性比对到目标病原物种的kmer序列，以此作为物种特异性kmer序列库。
进一步的所述公共数据库包括NCBI数据库；所述阳性一致性高的序列为基因阳性一致性PPV>0.90的序列，所述PPV计算公式为：同时携带目标耐药基因和kmer的菌株数/携带标耐药基因的所有菌株数。
更进一步的，所述物种特异性kmer序列库的构建方法具体如下：
a)针对特定目标病原菌，从NCBI genome数据库下载所有或足够量的病原菌株基因组；将所有基因组打碎为kmer序列，通过筛选得到与目标耐药基因的阳性检出一致性PPV>0.90的kmer序列；
b)将筛选出的kmer序列在NCBI数据库进行比对注释，保留仅特异性比对到该病原的kmer序列，以此作为物种特异性kmer序列库。
本申请还提供一种综合预测耐药基因-病原微生物归属的模型，包括如下模块：
1)物种注释推断模块：基于耐药基因的物种注释，推断可能的病原物种来源；
2)拷贝数评价模块：依据耐药基因和病原检出序列情况，评价耐药基因的拷贝数是否正常，保留拷贝数正常的基因-物种归属关系；
3)kmer库评价模块：内容同上述预测耐药基因-病原微生物归属模型。
进一步的，所述模块1)具体为：将耐药基因与NCBI nt库进行比对，获得各耐药基因的所有物种注释信息，并与CARD耐药数据库信息整合，构建耐药基因的潜在来源物种信息库，据此推断可能的病原物种来源。
进一步的，所述模块2)中，所述耐药基因的拷贝数是否正常是通过来自于NCBIgenome数据库的相应病原所有基因组进行耐药基因检测分析统计得到；所述耐药基因拷贝数的计算公式如下：
其中，CoverageARG为耐药基因的覆盖度，Cover_depthARG为耐药基因的覆盖深度，Coveragegenpme为基因组的覆盖度，Cover_depthgenpme为基因组的覆盖深度。
本申请还提供一种电子设备，包括：处理器和存储器；所述处理器和存储器相连，其中，所述存储器用于存储计算机程序，所述处理器用于调用所述计算机程序，以执行上述任一项所述的方法。
本申请还提供一种计算机存储介质，所述计算机存储介质存储有计算机程序，所述计算机程序包括程序指令，所述程序指令当被处理器执行时，执行上述任一项所述的方法。
本申请有益技术效果：
1)本申请针对mNGS耐药预测过程中耐药基因-病原物种归属不准问题，提供一套有效的解决方案，确保了检出耐药基因可以准确地匹配到感染病原中，从而可以进一步地准确预测病原菌的耐药表型结果，为临床标本mNGS应用于临床耐药菌的感染诊断提供了有力支撑。
2)本申请创新性地开发了一种基于基因组kmer序列库，筛选与目标耐药基因相关联的物种特异性序列方法，该方法克服了需求极高服务器计算资源配置的要求。一方面，构建筛选出来的kmer序列可以用于提升基于mNGS准确鉴定耐药基因的物种来源问题；另一方面，也为其他基因组数据挖掘研究(如耐药特征挖掘)提供了一种方法流程。
3)为解决计算量较大导致需要较高要求的服务器资源等问题，本申请还开发了一种更加节约计算资源的kmer序列库构建方法和相应的预测物种归属的方法。
4)本申请还定义一种基于mNGS测序计算基因拷贝数的方法，该方法可以用于有效解决临床标本存在单一病原感染菌时(尤其目标感染病原具有明显的丰度优势时)的耐药基因-物种归属问题。
附图说明
图1、本申请实施例1全流程技术路线图；
图2、本申请优选节约计算资源的物种归属流程图；
图3、肺炎克雷伯菌参考基因组kmer覆盖频率图；
图4、基于测序reads模拟测试验证基因-物种归属的结果比较；
图5、基于肛试子标本验证基因-物种归属结果图。
具体实施方式
下面将结合实施例对本申请的实施方案进行详细描述，但是本领域技术人员将会理解，下列实施例仅用于说明本申请，而不应视为限制本申请的范围。实施例中未注明具体条件者，按照常规条件或制造商建议的条件进行。所用试剂或仪器未注明生产厂商者，均为可以通过市场购买获得的常规产品。
部分术语定义
除非在下文中另有定义，本申请具体实施方式中所用的所有技术术语和科学术语的含义意图与本领域技术人员通常所理解的相同。虽然相信以下术语对于本领域技术人员很好理解，但仍然阐述以下定义以更好地解释本申请。
如本申请中所使用，术语“包括”、“包含”、“具有”、“含有”或“涉及”为包含性的(inclusive)或开放式的，且不排除其它未列举的元素或方法步骤。术语“由…组成”被认为是术语“包含”的优选实施方案。如果在下文中某一组被定义为包含至少一定数目的实施方案，这也应被理解为揭示了一个优选地仅由这些实施方案组成的组。
在提及单数形式名词时使用的不定冠词或定冠词例如“一个”或“一种”，“所述”，包括该名词的复数形式。
本申请中的术语“大约”表示本领域技术人员能够理解的仍可保证论及特征的技术效果的准确度区间。该术语通常表示偏离指示数值的±10％，优选±5％。
此外，说明书和权利要求书中的术语第一、第二、第三、(a)、(b)、(c)以及诸如此类，是用于区分相似的元素，不是描述顺序或时间次序必须的。应理解，如此应用的术语在适当的环境下可互换，并且本申请描述的实施方案能以不同于本申请描述或举例说明的其它顺序实施。
本申请的用于预测耐药基因-病原微生物归属模型的构建方法，以及相应的预测耐药基因-病原微生物归属的方法，包括如下步骤：
通过构建与目标耐药基因相关联的物种特异性kmer序列库，将病原耐药基因序列与物种特异性kmer序列库进行比对，判断耐药基因的物种来源；在一些实施方式中，如果在序列库中检测到病原耐药基因相关联的物种特异性kmer，那么确认双方的物种归属关系，否则拒绝物种归属关系。
所述物种特异性kmer序列库的构建方法可以如下：
a)基于公共数据库获得目标病原菌株基因组kmer，筛选得到与各目标耐药基因阳性一致性高的序列；b)通过数据库比对注释，挑选出特异性比对到目标病原物种的kmer序列，以此作为物种特异性kmer序列库。
在一些实施方式中，所述公共数据库优选为NCBI数据库；所述阳性一致性高的序列为基因阳性一致性PPV>0.90的序列，PPV计算公式为：同时携带目标耐药基因和kmer的菌株数/携带目标耐药基因的所有菌株数。
在一些实施方式中，所述物种特异性kmer序列库的构建方法具体如下：
a)针对特定目标病原菌，从NCBI genome数据库下载所有或足够量的病原菌株基因组；将所有基因组打碎为kmer序列，通过筛选得到与目标耐药基因的阳性检出一致性PPV>0.90的kmer序列；
b)将筛选出的kmer序列在NCBI数据库进行比对注释，保留仅特异性比对到该病原的kmer序列，以此作为物种特异性kmer序列库。
考虑到实践中筛选与目标耐药基因相关联的kemr序列时，往往计算量较大，导致需要较高要求的服务器资源。因此在一些特定的实施例中，通过如下方法构建物种特异性kmer序列库会更加节约计算资源，具体的，在步骤a)方面，具体步骤如下：
i.针对特定目标病原菌，从NCBI genome数据库下载所有或足够量的病原菌株基因组，对下载到的目标病原菌所有菌株基因组进行耐药基因检测分析，得到各菌携带的耐药基因信息，统计各重要耐药基因在所有菌株基因组上的发生率；
ii.计算目标病原基因组上的各区域发生率，筛选出较高发生率区域：根据细菌多位点序列分型MLST分型，在发生频次较高的多个MLST分型中各选择1个代表株，以多个代表株做为参考株，将剩余菌株基因组与参考株基因组比对，统计各个参考株的基因组上各区域发生率；以步骤i中统计得到的重要耐药基因发生率为阈值，筛选出高于阈值的各参考株基因组上的所有区域，作为较高发生率的区域；
iii.目标病原基因组序列kmer化，过滤掉位于基因组较高发生率区域上的kmer：针对每个菌株，将其基因组打碎为kmer序列，与步骤ii得到的基因组较高发生率区域进行比对和过滤，过滤掉位于基因组较高发生率区域上的kmer；
iv.针对过滤后剩余的各kmer序列，计算其在所有菌株中的发生率以及相对于目标耐药基因的菌株发生一致率PPV，挑选出PPV在0.90以上的kmer序列。
本申请的进一步涉及的综合预测耐药基因-病原微生物归属模型的构建方法，包括如下步骤或三种策略：
1)基于耐药基因的物种注释，推断可能的病原物种来源；
2)依据耐药基因和病原检出序列情况，评价耐药基因的拷贝数是否正常，保留拷贝数正常的基因-物种归属关系；
3)上述kmer序列库比对的方法。
在一些实施方式中，所述步骤1)具体为：将耐药基因与NCBI nt库进行比对，获得各耐药基因的所有物种注释信息，并与CARD耐药数据库信息整合，构建耐药基因的潜在来源物种信息库，
在一些实施方式中，所述步骤2)中，所述耐药基因的拷贝数是否正常是通过来自于NCBI genome数据库的相应病原所有基因组进行耐药基因检测分析统计得到；
所述耐药基因的拷贝数的计算公式如下：
其中，CoverageARG为耐药基因的覆盖度，Cover_depthARG为耐药基因的覆盖深度，
Coveragegenpme为基因组的覆盖度，Cover_depthgenpme为基因组的覆盖深度。
在一些更为具体的实施例中，所述步骤2)具体为：
基于NCBI genome数据库中的目标病原菌所有菌株基因组，将各菌株基因组contig序列直接与耐药基因数据库进行比对，保留identity>90％、coverage>60％、与NCBINDARO数据库注释记录比较测试性能大于0.90的结果；对于每个contig对齐区域，选择hit值最高的区域，同时计算检测到的耐药基因拷贝数(目标耐药基因的覆盖范围*覆盖范围的深度)/(目标基因组的覆盖范围/基因组覆盖区域的深度)，并统计目标病原菌各耐药基因的携带发生频次情况；然后，根据临床标本实际检出的耐药基因和病原菌检出reads情况，计算耐药基因的拷贝数，通过与以上统计得到的目标病原菌耐药基因发生频数范围进行比较，查看是否正常来判断物种耐药归属关系。
下面结合具体实施例来阐述本申请。
实施例1本申请方法建立
图1是本申请全流程方法学技术路线图，详细各个步骤的描述如下：
1.基于临床标本mNGS检测鉴定病原和耐药基因
1.1病原菌鉴定
针对临床常见的病原细菌(如肺炎克雷伯菌、大肠埃希菌、鲍曼不动杆菌、铜绿假单胞菌、阴沟肠杆菌、金黄色葡萄球菌等)，从NCBI genome数据库上搜索下载目标病原菌参考基因组，以此作为目标病原菌种鉴定的参比序列库。使用minimap2软件(v2.17，参数：-xsr-a--secondary＝no-L)，将NGS测序reads序列(如Illumina SE75测序reads或华大MGISE50测序reads)与以上设定的目标病原参考基因组序列库进行比对，然后计算比对上目标病原菌种的检出序列数、基因组覆盖度和覆盖深度。
1.2耐药基因鉴定
使用blastn(版本2.9.0+)软件将NGS测序reads(如Illumina测序或华大MGI测序reads)与耐药基因数据库(如CARD库)进行比对(参数设置：-evalue 1e-5-outfmt 6)，先进行过滤仅保留identity高于90％的hit，然后挑选每条reads序列的最高score的hit即besthit为最终hit，如果最高socre的hit存在多个相同值的(即多重比对)，则对这多个hit采用LCA算法进行该read序列的最终注释(即针对单条read序列，由于多重比对导致注释不到基因型，进而向更高一层级注释为基因家族水平)，然后统计得到样本中耐药基因检出的特异性reads数和多重比对reads数以及耐药基因所属家族的特异性reads数，并计算各检出基因的覆盖度和覆盖区域深度。
2.基于临床标本测序，实现耐药基因-病原物种的正确归属
联合使用3种策略，实现检出基因正确地归属给目标病原物种,具体如下：
策略I：基于耐药基因的物种注释，推断可能的病原物种来源。
将耐药基因与NCBI nt库进行比对和统计注释，获得各基因序列的物种注释结果后，并将其与耐药基因数据库CARD本身对各基因的物种来源信息进行合并，以此作为各耐药基因所有可能的物种来源信息。具体地，首先下载NCBI nt参考数据库到本地服务器(ascp-i～/.aspera/connect/etc/asperaweb_id_dsa.openssh-l 100M-k 1-T anonftp@ftp.ncbi.nlm.nih.gov:/blast/db/FASTA/nt.gz./)，然后对下载的NCBI NT数据库进行格式化建库(makeblastdb-dbtype nucl-in nt.fa-out nt.fa)，最后使用blastn将耐药基因与上述建库成功的NT库进行比对(blastn-query test.fa-db nt.fa-num_threads 48-max_target_seqs 5-outfmt 6)，保留Identity>90％和subject_coverage>90％的结果,统计获得各个耐药基因的所有物种注释信息，并且与CARD耐药数据库信息整合，得到耐药基因潜在的来源物种信息库。
策略II：依据耐药基因和病原检出序列情况，计算并查看耐药基因的拷贝数是否正常，保留拷贝数正常的基因-物种归属关系。
首先，耐药基因的拷贝数的正常范围是通过来自于NCBI genome数据库的相应病原所有基因组进行耐药基因检测分析统计得到的。具体地，基于NCBI genome数据库中的目标病原菌所有菌株基因组，将各菌株基因组contig序列直接与耐药基因数据库进行比对，使用BLASTN软件将组装好的基因组contig序列与ARG参考数据库进行比对(Version:ncbi-blast-2.9.0+,参数:-evalue 1e-5-outfmt 0-num_alignments 10000)，检测抗生素耐药基因是否存在SNP/Indels，仅保留identity>90％、coverage>60％、与NCBI NDARO数据库注释记录比较，测试性能大于0.90的结果。对于每个contig对齐区域，选择hit值最高的区域(hit_score/subject_coverage)，同时，使用公式计算检测到的耐药基因拷贝数：(目标耐药基因的覆盖范围*覆盖范围的深度)/(目标基因组的覆盖范围/基因组覆盖区域的深度)，以检测各菌株的耐药基因携带情况，并统计目标病原菌各耐药基因的携带发生频次情况。然后，按照如下定义的公式，根据临床标本实际检出的耐药基因和病原菌检出reads情况，计算耐药基因的拷贝数，通过与以上统计得到的目标病原菌耐药基因发生频数范围进行比较，查看是否正常来判断物种耐药归属关系。
公式中，耐药基因拷贝数(copy nARG)＝【耐药基因的覆盖度(CoverageARG)*耐药基因的覆盖深度(Cover depthARG)】/【基因组的覆盖度(Coverage genome)*基因组的覆盖深度(Cover depth genome)】。
策略III：构建与目标耐药基因相关联的物种特异性kmer序列库，通过检出有无来判断基因-物种归属关系。
针对某个病原，从NCBI genome数据库下载得到所有菌株基因组，然后将各基因组打碎为一定长度的kmer序列(长度为30bp)，筛选得到与目标耐药基因的阳性检出一致性较高(PPV>0.90)的kmer序列。最后将筛选出的kmer序列进行NCBI nt和NCBI Refseq genome数据库进行比对注释，并过滤保留仅特异性比对上该病原的kmer序列，以此作为物种特异性kmer序列库。通过检测与目标耐药基因相关联的物种特异性kmer的有无，来推断相应耐药基因的物种来源归属关系。同样地，可依此获得其他各病原菌的物种特异性kmer序列库。
具体地，在针对某个病原菌，筛选与目标耐药基因相关联的kemr序列时，往往计算量较大导致需要较高要求的服务器资源。这里，我们设计了开发了一种更加节约计算资源的可行性技术方案，技术路线如图2，详细如下如下步骤1-6。
步骤1病原微生物基因组下载：
针对某个病原菌(如肺炎克雷伯菌)，根据NCBI genome数据库的所有基因组组装信息表(wget https://ftp.ncbi.nlm.nih.gov/genomes/ASSEMBLY_REPORTS/assembly_summary_refseq.txt或wget https://ftp.ncbi.nlm.nih.gov/genomes/ASSEMBLY_REPORTS/assembly_summary_genbank.txt)，搜索该病原菌条目信息，找到对应的Assembly号和下载地址，然后下载所有菌株基因组数据。
步骤2耐药基因检测
对下载到的目标病原菌所有菌株基因组进行耐药基因检测分析，得到各菌携带的耐药基因信息，并统计得到各重要耐药基因在所有菌株基因组上的发生率。基于基因组contig比对检测耐药基因分析过程描述如下：使用BLASTN软件将组装好的基因组contig序列与ARG参考数据库进行比对(Version:ncbi-blast-2.9.0+,参数:-evalue 1e-5-outfmt0-num_alignments 10000)，检测抗生素耐药基因是否存在SNP/Indels，仅保留identity>90％、coverage>60％、与NCBI NDARO数据库注释记录比较，测试性能大于0.90的结果。对于每个contig对齐区域，选择hit值最高的区域(hit_score/subject_coverage)，同时，使用公式计算检测到的耐药基因拷贝数：(目标耐药基因的覆盖范围*覆盖范围的深度)/(目标基因组的覆盖范围/基因组覆盖区域的深度)
步骤3计算目标病原基因组上的各区域发生率，并筛选出较高发生率的区域
考虑到目标病原种群多态性，其整个泛基因组相对于单个菌株基因组可能会大很多，那么仅选用一个菌株作为参考基因组，并依此计算得到的目标病原基因组上各区域发生率可能无法包含其他附属区域的信息。因此我们根据MLST分型，在发生频次较高的5-10MLST分型里各选择出1个代表株，这样最终以这多个代表株做为参考株，然后统计计算出各个参考株的基因组上各区域发生率。然后，我们以步骤2中统计得到的各重要耐药基因发生率最高的值为阈值，来筛选出高于这一阈值的各参考株基因组上发生频率的所有区域，作为较高发生率的区域，以此用作后续的各菌株基因组kmer过滤使用。具体过程如下：首先对收集到的病原微生物基因组计算细菌多位点序列分型(Multilocussequencetyping，MLST)，从pubMLST数据库下载目标病原的管家基因信息，使用blast软件(参数：blastn-evalue 1e-5-num_threads 6-outfmt 0-num_alignments 10000)将菌株基因组contig与目标病原管家基因(如肺克的gapA、infB、mdh、pgi、phoE、rpoB以及tonB)进行比对和注释，获得各菌株ST分型，按照ST分型统计每个ST分型的菌株数目，并进行从大到小排序，从每个ST型菌株中挑选一株组装质量较高的基因组作为参考基因组。接着，将剩余的所有其他菌株基因组使用blastn比对(补充过滤阈值)到各个参考基因组，分别统计覆盖深度后，根据步骤2耐药基因检测得到的重要耐药基因的发生率最高的作为阈值，筛选出目标病原各参考基因组上较高发生率的区域。
步骤4目标病原基因组序列kmer化，并过滤掉位于基因组较高发生率区域上的kmer
针对每个菌株，使用jellyfish软件将其基因组按照1碱基为步长滑窗，打碎为kmer序列(长度为30bp)，然后与步骤3得到的基因组较高发生率区域使用blastn进行比对和过滤。
步骤5相对于目标耐药基因，计算各kmer的一致发生率PPV，并挑选较高PPV的kmer作为与目标基因相关联的kmer序列
针对过滤后剩余的各kmer序列，计算其在所有菌株中的发生率以及相对于目标耐药基因的菌株发生一致率(即PPV，公式为：同时携带目标耐药基因和kmer的菌株数/携带目标耐药基因的所有菌株数)，并挑选出PPV在0.90以上的kmer序列。
步骤6挑选物种特异性kmer序列
将步骤5获得PPV在0.90以上的kmer进行NCBI NT和NCBI RefSeq genome数据库比对，过滤掉多重比对上非本病原菌种的kmer，保留特异性比对上目标病原的kmer，即为与目标基因相关联的物种特异性kmer序列。
综上，联合使用3种策略进行耐药基因的物种归属过程为：例如，针对检出的某一个耐药基因，如果当前检出的病原属于该耐药基因的可能来源病原，则进一步计算耐药基因相对于此病原的基因拷贝数，并查看计算得到的拷贝数据是否位于正常范围内(通过训练集所有菌株的耐药基因检测结果统计得到)，如果满足上述条件，则将该基因归属给当前病原，否则拒绝当前的基因-物种归属关系。如果满足上述条件的病原存在多个，则根据是否检出与耐药基因相关联的物种特异性kmer序列，进行最终来源物种的判定。
3.预测目标病原-抗生素药敏结果
定义阴阳性判断指标Score(公式如下)，结合通过模拟测试确定的报告规则和阈值，进行检出目标病细菌的药敏结果预测。具体地，根据以下公式计算Score指标值，然后与阈值进行比较。当Score值大于或等于阈值，预测为耐药，当Score值小于阈值且目标病原测得的基因组覆盖较高(如大于50％)时预测为敏感，否则预测为”/”，即未知。
式中，arg_Wi表示相应基因型的权重系数，genefamily_Wi表示相应基因家族的权重系数。当检出基因型(genetype，即gt)且基因型权重系数>0，以基因型权重系数计算；当检出基因型，但基因型权重系数为0或无权重系数，则以基因家族权重系数计算(基因家族无权重系数时记为0)。
实施例2、与耐药基因相关联的目标病原菌种特异性kmer序列库构建
本实施例以针对肺炎克雷伯菌为例，对本申请方法进行阐述。
步骤1：针对肺炎克雷伯菌，从NCBI genome数据下载菌株基因组。基于NCBI数据库下载46354株肺炎克雷伯属的病原基因组数据(wget https://ftp.ncbi.nlm.nih.gov/genomes/ASSEMBLY_REPORTS/assembly_summary_genbank.t xt)。
步骤2：查找出肺炎克雷伯菌主要亚群的基因组较高发生率区域。将下载的肺炎克雷伯基因组数据进行MLST分型(使用blast比对管家基因，统计管家基因的数量，参数：blastn-evalue 1e-5-num_threads 6-outfmt 0-num_alignments 10000)，最后选出ST群体最大的几个ST型(ST11、ST14、ST147、ST15、ST16、ST258)，每个ST型群体挑选一株组装质量较高的基因组作为参考基因组，之后将剩余的菌株基因组使用blastn分别比对到这6个参考基因组(Blastn-evalue 1e-5-num_threads 6-outfmt 0-num_alignments 10000)。
步骤3：统计每个碱基位置被其余非参考菌株基因组覆盖的深度，最后根据此深度确定肺炎克雷伯菌属的基因组较高发生率区域(为每个碱基的位置进行排序编号，菌株基因组每覆盖该位置一次，该位置计数加一，依次循环，直到所有的kmer计算完毕，选择比对频率>23960作为基因组较高发生率区域；比对频率阈值23960的确定条件：对下载的肺炎克雷伯菌株基因组进行耐药基因检测，对所有的key耐药基因的物种发生率进行统计排序，选择排名靠前的sul1-like基因的物种发生频率55％作为阈值，通过公式确定：肺炎克雷伯菌株数*55％)，如下图3所示。
步骤4：将各菌株基因组打碎为kmer序列(长度为30bp)，然后进行与以上基因组较高发生率区域比对和过滤，仅保留位于非基因组较高发生率区域的kmer。将所有的病原物种基因组使用jellyfish软件(V2.3.0)按照1碱基为步长滑窗打碎为kmer(参数：dump-c-t-U 1000)，通过kmer与基因组较高发生率区域比对，计算。使用blast软件将所有的kmer(570w个)比对到基因组较高发生率区域(参数blastn-evalue 1e-5-num_threads 6-outfmt 0-num_alignments 10000)，过滤完全比对的kmer，将剩余的kmer作为后续使用的kmer(150w个)，对上述kmer计算在46354株肺炎克雷伯属菌株中的发生率(公式：每个kmer在所有肺炎克雷伯基因组组出现的次数/46354*100％)。同时对上述的肺炎克雷伯属菌株进行耐药检测，统计key耐药基因在46354株肺炎克雷伯属菌株中的发生情况(公式：每个key耐药基因在所有肺炎克雷伯基因组组出现的次数/46354*100％)，计算得到每个kmer的PPV，最后选取PPV>0.90的kmer，使用blast比对到NT库后(参数：blastn-task megablast-evalue 1e-05-perc_identity 70-qcov_hsp_perc 70-outfmt 6-out AAC3_IIe.0.m8-num_threads 8)过滤掉其他物种的kmer，完成建库。
实施例3、混合感染测序reads模拟测试验证基因-物种归属的性能
步骤1：基于从NCBI genome数据库下载的肺炎克雷伯菌和大肠埃希菌菌株基因组耐药基因检测结果，挑选出携带KPC-2肺炎克雷伯菌株基因组(不含NDM-5)和携带NDM-5的大肠埃希菌株基因组(不含KPC-2)各30例。
步骤2：基于步骤1下载的基因组，使用ART软件(Version2.5.8，参数-ss NS50-l75-f 5-nf 0-rs 1)分别模拟0.01x、0.02x、0.03x、0.05x、0.1x、0.2x、0.3x、0.4x、0.5x、0.6x、0.7x、0.8x、0.9x、1x、2x、3x、4x、5x、10x等不同梯度数据量的短测序reads(IlluminaSE50)，然后将模拟的的肺炎克雷伯菌基因组测序reads和大肠埃希菌基因组测序reads，按照不同数据量任意两两混合形成模拟的混合感染的样本。
步骤3：将步骤2得到的混合模拟样本，进行本申请的耐药基因检测和基因-物种归属分析，评估耐药基因物种归属的准确性，并对比仅策略I/II联合与策略I/II/III联合的性能差异。
结果显示(如图4)，仅联合使用策略I与策略II，大多数模拟的混合样本的基因-物种归属准确率都比较高，准确率大多分布在70％-100％，但是存在一些准确率较低的组合，尤其是混合的两个病原菌的含量相差较小时(如在8倍差异内)；而将策略III也联合使用时，基因-物种归属的准确率明显得到提高，绝大多数样本的准确率达到分布在95％-100％的区间范围。
实施例3、临床14例肛试子标本mNGS测序验证基因-物种归属的性能
步骤1：收集到临床14例大肠埃希菌和肺炎克雷伯菌培养阳性的肛试子样本，然后进行mNGS测序，进行病原鉴定和药敏预测。
步骤2：基于步骤1的mNGS测序数据，进行生物信息学分析，完成病原鉴定，耐药基因检测以及药敏结果预测，然后与培养药敏结果对比。
结果显示(见图5)，针对亚胺培南，相对于基于培养的AST结果，mNGS-based药敏预测(基因-物种模块联合使用了本申请3种策略)，阳性预测一致性和阴性预测一致性分别为100％(7/7)、100％(4/4)。
上述结果表明，在mNGS-based药敏预测过程中，本申请的方法可以非常准确地实现耐药基因的检测及其来源物种的归属。
最后应说明的是：以上各实施例仅用以说明本申请的技术方案，而非对其限制；尽管参照前述各实施例对本申请进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本申请各实施例技术方案的范围。
Claims (8)
1.一种用于预测耐药基因-病原微生物归属模型的构建方法，其特征在于，所述方法包括如下步骤：
构建与目标耐药基因相关联的物种特异性kmer序列库，将待预测的病原耐药基因序列与物种特异性kmer序列库进行比对，判断该耐药基因的物种来源；
所述物种特异性kmer序列库的构建方法如下：
基于公共数据库获得目标病原菌株基因组kmer，筛选得到与各目标耐药基因阳性一致性高的序列；通过数据库比对注释，挑选出特异性比对到目标病原物种的kmer序列，以此作为物种特异性kmer序列库；
所述公共数据库包括NCBI数据库；所述与各目标耐药基因阳性一致性高的序列为与各目标耐药基因阳性一致性PPV>0.90的序列；所述与各目标耐药基因阳性一致性PPV＝同时携带目标耐药基因和kmer的菌株数/携带目标耐药基因的所有菌株数。
2.根据权利要求1所述的构建方法，其特征在于，所述物种特异性kmer序列库的构建方法具体如下：
a)针对特定目标病原菌，从NCBI genome数据库下载所有或足够量的病原菌株基因组；将所有基因组打碎为kmer序列，通过筛选得到与目标耐药基因的阳性检出一致性PPV>0.90的kmer序列；
b)将筛选出的kmer序列在NCBI数据库进行比对注释，保留仅特异性比对到该病原的kmer序列，以此作为物种特异性kmer序列库。
3.根据权利要求2所述的构建方法，其特征在于，所述步骤a)具体包括如下步骤：
i.针对特定目标病原菌，从NCBI genome数据库下载所有或足够量的病原菌株基因组，对下载到的目标病原菌所有菌株基因组进行耐药基因检测分析，得到各菌携带的耐药基因信息，统计各重要耐药基因在所有菌株基因组上的发生率；
ii.计算目标病原基因组上的各区域发生率，筛选出较高发生率区域：根据细菌多位点序列分型MLST分型，在发生频次较高的多个MLST分型中各选择1个代表株，以多个代表株做为参考株，将剩余菌株基因组与参考株基因组比对，统计各个参考株的基因组上各区域发生率；以步骤i中统计得到的重要耐药基因发生率为阈值，筛选出高于阈值的各参考株基因组上的所有区域，作为较高发生率的区域；
iii.目标病原基因组序列kmer化，过滤掉位于基因组较高发生率区域上的kmer：针对每个菌株，将其基因组打碎为kmer序列，与步骤ii得到的基因组较高发生率区域进行比对和过滤，过滤掉位于基因组较高发生率区域上的kmer；
iv.针对过滤后剩余的各kmer序列，计算其在所有菌株中的发生率以及相对于目标耐药基因的菌株发生一致率PPV，挑选出PPV在0.90以上的kmer序列。
4.一种用于综合预测耐药基因-病原微生物归属模型的构建方法，其特征在于，包括如下步骤：
1)基于耐药基因的物种注释，推断可能的病原物种来源；
2)依据耐药基因和病原检出序列情况，评价耐药基因的拷贝数是否正常，保留拷贝数正常的基因-物种归属关系；
3)权利要求1-3任一所述的构建方法。
5.根据权利要求4所述的构建方法，其特征在于，所述步骤1)具体为：将耐药基因与NCBI nt库进行比对，获得各耐药基因的所有物种注释信息，并与CARD耐药数据库信息整合，构建耐药基因的潜在来源物种信息库，据此推断可能的病原物种来源。
6.根据权利要求5所述的构建方法，其特征在于，所述步骤2)中，所述耐药基因的拷贝数是否正常是通过来自于NCBI genome数据库的相应病原所有基因组进行耐药基因检测分析统计得到；所述耐药基因的拷贝数的计算公式如下：
其中，ConerageARG为耐药基因的覆盖度，Cover_depthARG为耐药基因的覆盖深度，Coverage1enome为基因组的覆盖度，Cover_depth1enome为基因组的覆盖深度。
7.一种电子设备，其特征在于，包括：处理器和存储器；所述处理器和存储器相连，其中，所述存储器用于存储计算机程序，所述处理器用于调用所述计算机程序，以执行如权利要求1-6任一项所述的方法。
8.一种计算机存储介质，其特征在于，所述计算机存储介质存储有计算机程序，所述计算机程序包括程序指令，所述程序指令当被处理器执行时，执行如权利要求1-6任一项所述的方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310702713.1A CN116597893B (zh) | 2023-06-14 | 2023-06-14 | 预测耐药基因-病原微生物归属的方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310702713.1A CN116597893B (zh) | 2023-06-14 | 2023-06-14 | 预测耐药基因-病原微生物归属的方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116597893A CN116597893A (zh) | 2023-08-15 |
| CN116597893B true CN116597893B (zh) | 2023-12-15 |
Family
ID=87611699
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310702713.1A Active CN116597893B (zh) | 2023-06-14 | 2023-06-14 | 预测耐药基因-病原微生物归属的方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116597893B (zh) |
Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103384828A (zh) * | 2010-12-23 | 2013-11-06 | 雀巢产品技术援助有限公司 | 使用基于抗体的阵列进行恶性癌症疗法的药物选择 |
| CN106566874A (zh) * | 2016-08-24 | 2017-04-19 | 首都医科大学附属北京友谊医院 | 检测肺炎支原体耐药突变基因的特异性引物组及检测试剂盒 |
| CN110021365A (zh) * | 2018-06-22 | 2019-07-16 | 深圳市达仁基因科技有限公司 | 确定检测靶点的方法、装置、计算机设备和存储介质 |
| CN112530519A (zh) * | 2020-12-14 | 2021-03-19 | 广东美格基因科技有限公司 | 一种检测样本中微生物和耐药基因的方法和系统 |
| CN112863601A (zh) * | 2021-01-15 | 2021-05-28 | 广州微远基因科技有限公司 | 病原微生物耐药基因归属模型及其建立方法和应用 |
| CN114067912A (zh) * | 2021-11-23 | 2022-02-18 | 天津金匙医学科技有限公司 | 基于机器学习筛选细菌耐药表型相关重要特征基因的方法 |
| CN114333987A (zh) * | 2021-12-30 | 2022-04-12 | 天津金匙医学科技有限公司 | 一种基于宏基因组测序的预测耐药表型的数据分析方法 |
| CN115188411A (zh) * | 2022-06-29 | 2022-10-14 | 南方科技大学 | 一种基于单分子自适应采样测序的病原微生物和耐药基因的检测方法及系统 |
| CN116189919A (zh) * | 2023-04-20 | 2023-05-30 | 杭州迪安医学检验中心有限公司 | 一种微生物药敏的计算机分析方法、分析系统及其应用 |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6924094B1 (en) * | 1996-02-08 | 2005-08-02 | Affymetrix, Inc. | Chip-based species identification and phenotypic characterization of microorganisms |
| US10095831B2 (en) * | 2016-02-03 | 2018-10-09 | Verinata Health, Inc. | Using cell-free DNA fragment size to determine copy number variations |
| WO2019113577A1 (en) * | 2017-12-10 | 2019-06-13 | Yan Wang | A Multiplexed Method for Detecting DNA Mutations and Copy Number Variations |
| KR102349921B1 (ko) * | 2018-09-05 | 2022-01-12 | 주식회사 천랩 | 시료 미생물의 동정 및 분류 방법 |
-
2023
- 2023-06-14 CN CN202310702713.1A patent/CN116597893B/zh active Active
Patent Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103384828A (zh) * | 2010-12-23 | 2013-11-06 | 雀巢产品技术援助有限公司 | 使用基于抗体的阵列进行恶性癌症疗法的药物选择 |
| CN106566874A (zh) * | 2016-08-24 | 2017-04-19 | 首都医科大学附属北京友谊医院 | 检测肺炎支原体耐药突变基因的特异性引物组及检测试剂盒 |
| CN110021365A (zh) * | 2018-06-22 | 2019-07-16 | 深圳市达仁基因科技有限公司 | 确定检测靶点的方法、装置、计算机设备和存储介质 |
| CN112530519A (zh) * | 2020-12-14 | 2021-03-19 | 广东美格基因科技有限公司 | 一种检测样本中微生物和耐药基因的方法和系统 |
| CN112863601A (zh) * | 2021-01-15 | 2021-05-28 | 广州微远基因科技有限公司 | 病原微生物耐药基因归属模型及其建立方法和应用 |
| CN114067912A (zh) * | 2021-11-23 | 2022-02-18 | 天津金匙医学科技有限公司 | 基于机器学习筛选细菌耐药表型相关重要特征基因的方法 |
| CN114333987A (zh) * | 2021-12-30 | 2022-04-12 | 天津金匙医学科技有限公司 | 一种基于宏基因组测序的预测耐药表型的数据分析方法 |
| CN115188411A (zh) * | 2022-06-29 | 2022-10-14 | 南方科技大学 | 一种基于单分子自适应采样测序的病原微生物和耐药基因的检测方法及系统 |
| CN116189919A (zh) * | 2023-04-20 | 2023-05-30 | 杭州迪安医学检验中心有限公司 | 一种微生物药敏的计算机分析方法、分析系统及其应用 |
Non-Patent Citations (4)
| Title |
|---|
| Genomic Features Associated with the Degree of Phenotypic Resistance to Carbapenems in Carbapenem-Resistant Klebsiella pneumoniae;Zackery P. Bulman 等;《mSystems》;第06卷(第05期);1-17 * |
| Novel Clinical mNGS-Based Machine Learning Model for Rapid Antimicrobial Susceptibility Testing of Acinetobacter baumannii;Xuejiao Hu 等;《Journal of Clinical Microbiology》;第61卷(第05期);1-15 * |
| 沙门氏菌泛耐药基因组分析及耐药表型预测;何逸尘;《中国优秀硕士学位论文全文数据库 基础科学辑》(第(2022)01期);A006-975 * |
| 表面增强拉曼光谱用于病原菌种属鉴定及耐药性分析研究;陈雪萍;《中国博士学位论文全文数据库 医药卫生科技辑》(第(2020)03期);E060-66 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116597893A (zh) | 2023-08-15 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Anand et al. | Next generation sequencing of pooled samples: guideline for variants’ filtering | |
| CN113689912B (zh) | 基于宏基因组测序的微生物对比结果校正的方法和系统 | |
| CN114333987B (zh) | 一种基于宏基因组测序的预测耐药表型的数据分析方法 | |
| Sekizuka et al. | TGS-TB: total genotyping solution for Mycobacterium tuberculosis using short-read whole-genome sequencing | |
| CN113160882B (zh) | 一种基于三代测序的病原微生物宏基因组检测方法 | |
| JP6240210B2 (ja) | 標的シーケンシングリードの正確かつ迅速なマッピング | |
| Hall et al. | Evaluation of Nanopore sequencing for Mycobacterium tuberculosis drug susceptibility testing and outbreak investigation: a genomic analysis | |
| Mallawaarachchi et al. | Genomic diagnostics in polycystic kidney disease: an assessment of real-world use of whole-genome sequencing | |
| Shi et al. | Fast and accurate metagenotyping of the human gut microbiome with GT-Pro | |
| Shuch et al. | Estimation of the carrier frequency of fumarate hydratase alterations and implications for kidney cancer risk in hereditary leiomyomatosis and renal cancer | |
| US20230135480A1 (en) | Molecular technology for detecting a genome sequence in a bacterial genome | |
| CN115083521A (zh) | 一种单细胞转录组测序数据中肿瘤细胞类群的鉴定方法及系统 | |
| Sobih et al. | Metaflow: Metagenomic profiling based on whole-genome coverage analysis with min-cost flows | |
| Liu et al. | De novo identification of microbial contaminants in low microbial biomass microbiomes with Squeegee | |
| Panoutsopoulou et al. | Quality control of common and rare variants | |
| WO2019046804A1 (en) | IDENTIFICATION OF FALSE POSITIVE VARIANTS USING A MODEL OF IMPORTANCE | |
| CN116597893B (zh) | 预测耐药基因-病原微生物归属的方法 | |
| Oh et al. | Identification of signature gene set as highly accurate determination of metabolic dysfunction-associated steatotic liver disease progression | |
| US20150032385A1 (en) | Methods of Analyzing Massively Parallel Sequencing Data | |
| Mauguen et al. | Estimating the probability of clonal relatedness of pairs of tumors in cancer patients | |
| Abdelrazik et al. | Benchmarking of Antimicrobial Resistance Gene Detection Tools in Assembled Bacterial Whole Genomes | |
| Hall et al. | Nanopore sequencing for Mycobacterium tuberculosis drug susceptibility testing and outbreak investigation | |
| Aldawiri et al. | A Novel Approach for Mapping Ambiguous Sequences of Transcriptomes | |
| CN116825182B (zh) | 一种基于基因组ORFs筛选细菌耐药特征的方法及应用 | |
| Kelmemi et al. | Determining the genome-wide kinship coefficient seems unhelpful in distinguishing consanguineous couples with a high versus low risk for adverse reproductive outcome |
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 | ||
| CP03 | Change of name, title or address | ||
| CP03 | Change of name, title or address |
Address after: Room 110, 1st Floor, Building 3, No. 2 East Binhe Road, Youanmenwai, Fengtai District, Beijing, 100050 Patentee after: Beijing Jinshao Medical Laboratory Co.,Ltd. Country or region after: China Patentee after: Tianjin JinKe Medical Technology Co.,Ltd. Patentee after: Tianjin Huazhinuo Technology Co.,Ltd. Address before: Room 110, 1st Floor, Building 3, No. 2 East Binhe Road, Youanmenwai, Fengtai District, Beijing, 100050 Patentee before: Beijing Jinshao Medical Laboratory Co.,Ltd. Country or region before: China Patentee before: Tianjin JinKe Medical Technology Co.,Ltd. Patentee before: Jinshi Zhizao (Tianjin) Medical Technology Co.,Ltd. |