非核心型耐药基因在肺炎克雷伯菌药敏预测中的应用的制作方法

非核心型耐药基因在肺炎克雷伯菌药敏预测中的应用
1.本技术为中国专利申请cn202111400540.5的分案申请。
技术领域
2.本技术涉及基因测序技术领域，具体涉及一种基于机器学习技术筛选细菌耐药表型相关重要特征基因的方法。
技术背景
3.全基因组关联研究(genome-wide association study,gwas)是一种从基因组水平筛选与某表型(phenotype)显著相关的遗传变异，进而阐明表型遗传机制的方法。相较于传统的分子遗传学方法，gwas并不对表型产生的遗传机制做任何假设，而是直接从表型出发，设置合理的对照组，通过大样本的数据统计分析找到与表型关联的遗传变异。在人类复杂疾病研究中，gwas已取得丰硕成果，极大增进了人们对复杂表型的认识。同样，gwas也可以用于细菌研究，可为宿主适应性、耐药、毒力等复杂表型的遗传机制探索提供新思路。
4.针对细菌耐药性研究，鉴于耐药产生的复杂性和目前已了解的多样化的耐药性产生机制，与耐药表型相关的遗传数据可概括为snp、插入缺失(indel)、非核心基因获得缺失及k-mer等特征。目前已发表的细菌gwas研究大都集中于对核心基因组层面的单核苷酸多态性(snp)、插入/缺失(indel)或k-mer特征的筛选，如结核杆菌、金黄色葡萄球菌、肺炎链球菌、鲍曼不动杆菌等，且这些不同菌的研究所用到的软件算法各异，导致所得出的结果存在不同程度的偏移(假阳性)。实际上，由于长期的抗生素广泛使用和滥用，临床现今面临更多的是细菌获得性耐药导致的药物耐药率逐年升高加剧问题。细菌获得性耐药，一方面是由于种群本身存在某亚群携带耐药基因导致耐药外，另外较常见的是由耐药基因的水平转移介导产生的耐药(如质粒的水平转移)，因此除了与耐药产生相关的snp等变异特征外，筛选找到与耐药表型相关的重要非核心型耐药基因也是十分重要的工作。再就是，不同药物其耐药性的产生往往也具有多机制途径的特点。尽管，携带单个耐药基因并表达会导致细菌产生耐药表现，如kpc-2基因会导致碳青霉烯耐药，但也存在单一机制或单个耐药基因表达却不足以直接造成耐药表现的情况，所以有必要在筛选耐药相关特征的同时，能够对每个特征量化其对耐药产生的贡献，即对这些特征基因进行重要性排序。
5.有鉴于此，提出本技术。


技术实现要素：

6.本技术要解决的核心问题是寻求一种筛选耐药表型相关重要特征基因(非核心型耐药基因)的方法。为解决上述问题，本技术提出如下技术方案：
7.本技术首先提供一种基于机器学习技术筛选细菌耐药表型相关重要特征基因的方法，所述方法包括如下步骤：
8.步骤1)、获取目标细菌菌株基因组数据，同时收集对应药敏测试结果数据；
9.步骤2)、基于细菌基因组的contig序列进行耐药数据库的比对注释；
10.步骤3)、针对目标药物，进行基因型与耐药表型数据关联分析，筛选与耐药产生相关的重要特征基因，并计算重要特征基因权重系数；
11.优选的，所述重要特征基因为非核心耐药基因。
12.进一步的，所述方法还包括：
13.步骤4)、roc分析评估基于筛选出的重要基因预测药敏结果的模型性能。
14.进一步的，所述roc分析如下：基于步骤3)得到的重要基因权重系数的矩阵，定义并计算score值，并以此作为阴阳性判读指标，绘制roc曲线以及确定cut off值，用验证集样本进行验证和评估模型性能；所述其中arg_w

表示检出相应基因的权重系数值。
15.进一步的，所述步骤1)中所述菌株基因组个数》＝100，菌株来源涵盖各种亚型，耐药株与敏感株的菌株数目比例均衡；
16.在一些优选方式中，所述获取从公共数据库搜索下载已发表的目标基因组序列，或者通过现行收集临床培养鉴定出的细菌菌株进行测序组装获得；
17.在一些更优选的方式中，所述从公共数据库搜索下载为：从ncbi ndaro数据库和patric数据库平台搜集记录有药敏测试结果的细菌菌株信息，整理表型数据，并根据基因组组装id号从ncbi genome数据库或根据patric id从patric数据库批量下载基因组数据。
18.进一步的，所述步骤2)中所述比对注释为：将contig序列与card耐药基因参考序列库进行比对，过滤identity和覆盖度低的hit(优选的，先过滤identity小于90％或参比基因覆盖度subject coverage小于90％的hit)再将每条contig比对上的区域选择best hit作为该contig区域的最终比对结果，并添加耐药基因的注释信息。
19.进一步的，所述步骤3)中所述关联分析采用拉索回归模型进行关联分析。
20.更进一步的，所述步骤3)中拉索回归模型关联分析方法具体为：以基因检出分布矩阵和抗生素药敏测试结果数据矩阵作为输入，使用glmnet程序包进行基因型和耐药表型数据的关联分析，并做k(优选k＝5～15)重交叉验证，筛选得到与耐药表型相关的重要特征基因，计算重要特征基因的权重系数；
21.进一步的，所述重要特征基因具体为：根据不同个数特征基因下的模型cv错误率、auc变化曲线，选择在cv错误率最低且此时模型auc值相对稳定处所对应的数目基因为重要特征基因。
22.进一步的，步骤3)中进一步可包括人工召回，所述人工召回为：对与耐药表型ppv较高(优选的，ppv》＝0.8)的基因进行人工召回，并基于以上得到的重要基因其权重系数值计算得到召回基因的权重系数。
23.进一步的，本技术所述细菌包括但不限于大肠埃希式菌、肺炎克雷伯式菌、鲍曼不动杆菌、铜绿假单胞菌、阴沟肠杆菌复合群、金黄色葡萄球菌、屎肠球菌、粪肠球菌、肺炎链球菌、化脓性链球菌、流感嗜血杆菌、表皮葡萄球菌；优选的为肺炎克雷伯菌。
24.进一步的，本技术所述耐药表型包括但不限于耐受碳青霉烯类、头孢菌素类、青霉素类、β内酰胺类抗生素抑制剂、氨基糖苷类、磺胺类、四环素类、喹诺酮类、糖肽类、恶唑烷酮类、多粘菌素药物的表型；优选的，所述耐药表型为耐受碳青霉烯药物的表型。
25.本技术还提供一种基于机器学习技术筛选细菌耐药表型相关重要特征基因的模型，所述模型包括如下模块：
26.模块1)、用于获取目标细菌菌株基因组数据，同时收集对应药敏测试结果数据；
27.模块2)、用于基于细菌基因组的contig序列进行耐药数据库的比对注释；
28.模块3)、用于针对目标药物，进行基因型与耐药表型数据关联分析，筛选与耐药产生相关的重要特征基因，并计算重要特征基因权重系数；
29.优选的，所述重要特征基因为非核心型耐药基因。
30.进一步的，所述模块还包括：
31.模块4)、roc分析评估基于筛选出的重要基因预测药敏结果的模型性能。
32.进一步的，上述各模块分别执行上述方法步骤，因而模块的进一步限定同上述方法中各步骤的限定。
33.本技术还提供一种基因kpc-3、kpc-2、oxa-181、oxa-232和ndm-1作为非核心型耐药基因在肺炎克雷伯菌辅助药敏预测中的应用；
34.进一步的，所述药敏预测包括耐药预测和敏感预测，优选为敏感预测；
35.进一步的，所述药敏针对碳青霉烯类药物，例如亚胺培南或美罗培南。
36.当为亚胺培南时，所述基因还包括基因kpc-14、ndm-5、ndm-7、oxa-65、oxa-83、oxa-24、adc-214、adc-56、vim-27、vim-1、vim-19、veb-1中的一个或多个或全部；
37.当为美罗培南时，所述基因还包括kpc-4、oxa-163、ndm-7、oxa-65、imp-4、adc-214、ges-19、ges-26、ges-20、vim-27、vim-1、vim-19、vim-4、oxa-24中的一个或多个或全部。
38.本发明还提供一种针对非核心型耐药基因kpc-3、kpc-2、oxa-181、oxa-232和ndm-1检测试剂在制备肺炎克雷伯菌辅助药敏预测试剂盒中的应用；
39.进一步的，所述药敏预测包括耐药预测和敏感预测，优选为敏感预测；
40.进一步的，所述药敏针对碳青霉烯类药物，例如亚胺培南或美罗培南。
41.当为亚胺培南时，所述基因还包括基因kpc-14、ndm-5、ndm-7、oxa-65、oxa-83、oxa-24、adc-214、adc-56、vim-27、vim-1、vim-19、veb-1中的一个或多个或全部；
42.当为美罗培南时，所述基因还包括kpc-4、oxa-163、ndm-7、oxa-65、imp-4、adc-214、ges-19、ges-26、ges-20、vim-27、vim-1、vim-19、vim-4、oxa-24中的一个或多个或全部。
43.本发明还提供一种肺炎克雷伯菌辅助药敏检测试剂盒，包含对针对作为非核心型耐药基因kpc-3、kpc-2、oxa-181、oxa-232和ndm-1的检测试剂；
44.优选的：
45.当针对亚胺培南药物时，所述基因还包括基因kpc-14、ndm-5、ndm-7、oxa-65、oxa-83、oxa-24、adc-214、adc-56、vim-27、vim-1、vim-19、veb-1中的一个或多个或全部；
46.当针对美罗培南药物时，所述基因还包括kpc-4、oxa-163、ndm-7、oxa-65、imp-4、adc-214、ges-19、ges-26、ges-20、vim-27、vim-1、vim-19、vim-4、oxa-24中的一个或多个或全部。
47.本技术有益技术效果
48.1、本技术拓展了机器学习技术在细菌耐药性研究方向的应用，主要针对与耐药表型产生相关的非核心型基因获得缺失特征进行关联分析，找出造成耐药表型贡献度高的重要耐药基因，同时计算得到这些耐药基因相应的权重系数，便于后续进一步向临床耐药检
测的转化应用研究。
49.2、本技术涉及到的基因型数据主要为非核心型基因获得缺失矩阵(0-1矩阵)，同时可以拓展用于snp、插入缺失等特征的筛选挖掘。表型数据为药敏测试结果，可涵盖基本所有常见的药物类型，如碳青霉烯类、头孢菌素类、氨基糖苷类、磺胺类、四环素类、喹诺酮类等。
50.3、本技术在检测耐药基因时，直接基于基因组contig序列与公共的card耐药数据库比对方式，进行耐药基因检测和注释，绕过基因预测以及基于预测得到的cds序列比对检测耐药基因环节，这样避开了基因预测过程可能引入的偏差。
51.4、本技术在进行基因型和表型数据关联分析时，通过对比三种机器学习算法(glm广义线性模型、拉索回归模型、随机森林)的差异，最终选择效果最优的拉索回归模型算法，其在筛选得到与耐药表型相关的重要基因的同时，可以对这些重要基因进行排序，计算得到各耐药基因的权重系数。
附图说明
52.图1、本技术技术总体技术流程；
53.图2、针对肺炎克雷伯菌，各抗生素药物模型性能测试，关联分析模型比较结果；
54.图3、不同个数特征基因下的模型cv错误率、auc变化曲线；
55.图4、roc曲线验证筛选出的耐药亚胺培南的特征基因的可靠性；
56.图5、roc曲线验证筛选出的针对美罗培南特征基因的可靠性。
具体实施方式
57.下面将结合实施例对本技术的实施方案进行详细描述，但是本领域技术人员将会理解，下列实施例仅用于说明本技术，而不应视为限制本技术的范围。实施例中未注明具体条件者，按照常规条件或制造商建议的条件进行。所用试剂或仪器未注明生产厂商者，均为可以通过市场购买获得的常规产品。
58.部分术语定义
59.除非在下文中另有定义，本技术具体实施方式中所用的所有技术术语和科学术语的含义意图与本领域技术人员通常所理解的相同。虽然相信以下术语对于本领域技术人员很好理解，但仍然阐述以下定义以更好地解释本技术。
60.如本技术中所使用，术语“包括”、“包含”、“具有”、“含有”或“涉及”为包含性的(inclusive)或开放式的，且不排除其它未列举的元素或方法步骤。术语“由
…
组成”被认为是术语“包含”的优选实施方案。如果在下文中某一组被定义为包含至少一定数目的实施方案，这也应被理解为揭示了一个优选地仅由这些实施方案组成的组。
61.在提及单数形式名词时使用的不定冠词或定冠词例如“一个”或“一种”，“所述”，包括该名词的复数形式。
62.本技术中的术语“大约”表示本领域技术人员能够理解的仍可保证论及特征的技术效果的准确度区间。该术语通常表示偏离指示数值的
±
10％，优选
±
5％。
63.此外，说明书和权利要求书中的术语第一、第二、第三、(a)、(b)、(c)以及诸如此类，是用于区分相似的元素，不是描述顺序或时间次序必须的。应理解，如此应用的术语在
适当的环境下可互换，并且本技术描述的实施方案能以不同于本技术描述或举例说明的其它顺序实施。
64.本技术所述的“耐药”又称抗药性，系指微生物、寄生虫以及肿瘤细胞对于药物作用的耐受性，耐药性一旦产生，药物作用就明显下降。本技术优选的指体内感染为细菌对于抗生素类药物的耐药性。
65.本技术所述的“耐药表型”，通常是指为生物呈现的耐药特征称为耐药表型(resistant phenotype)，而拥有的耐药基因称为耐药基因型(resistant genotype)。
66.本技术所述的“非核心基因”是这样一类基因，针对某细菌种群，仅在部分菌株中存在的基因，与之相对应的是核心基因，即所有菌株中都存在的基因。本技术方法检测的耐药基因主要是针对这种非核心基因。
67.以下术语或定义仅仅是为了帮助理解本技术而提供。这些定义不应被理解为具有小于本领域技术人员所理解的范围。
68.本技术所述的基于机器学习技术筛选细菌耐药表型相关重要特征基因(非核心型耐药基因。)的方法，大体包括如下思路：基于bgwas思想搜集公共平台上目标细菌基因组或者现行收集测序组装后得到的大样本量菌株基因组数据及其对应的抗生素药物药敏测试结果，使用机器学习方法进行基因型与耐药表型两者之间的关联分析，以筛选出与耐药表型相关的重要特征基因，同时计算得到这些重要特征基因的权重系数，最后使用roc分析确定找出来的各药物相关耐药基因的可靠性。具体可参见附图1。
69.可以理解，本技术是拓展了机器学习技术在细菌耐药性研究方向的应用，主要针对与耐药表型产生相关的非核心基因获得缺失特征进行关联分析，找出造成耐药表型贡献度高的重要耐药基因，同时计算得到这些耐药基因相应的权重系数，便于后续进一步向临床耐药检测的转化应用研究。
70.本技术方法对于耐药细菌的种类不做限制，其具有普适性，例如大肠埃希式菌、肺炎克雷伯式菌、鲍曼不动杆菌、铜绿假单胞菌、阴沟肠杆菌复合群、金黄色葡萄球菌、屎肠球菌、粪肠球菌、肺炎链球菌、化脓性链球菌、流感嗜血杆菌、表皮葡萄球菌；优选的为肺炎克雷伯菌等；而且本技术方法还可以拓展用于snp、插入缺失等特征的筛选挖掘等；另外，本技术的表型数据为药敏测试结果，可涵盖基本所有常见的药物类型，如但不限于耐受碳青霉烯类、头孢菌素类、青霉素类、β内酰胺类抗生素抑制剂、氨基糖苷类、磺胺类、四环素类、喹诺酮类、糖肽类、恶唑烷酮类、多粘菌素药物的表型。
71.本技术方法中在菌株基因组收集的要求大体如下：菌株基因组个数至少100，且菌株来源广泛多样涵盖各种亚型，耐药株与敏感株的菌株数目比例尽量均衡。
72.在一些实施方式中，细菌菌株基因组的获取可以从公共数据库搜索下载得到已发表的目标基因组序列，或者通过现行收集临床培养鉴定出的细菌菌株进行测序组装获得，同时收集对应菌株的药敏测试结果信息。
73.在一些具体的实施方式中，公共基因组数据的搜索和下载：从ncbi ndaro数据库(https://www.ncbi.nlm.nih.gov/pathogens/isolates)和patric数据库平台(https://patricbrc.org)搜集记录有药敏测试结果的细菌菌株信息，整理表型数据，并根据基因组组装id号从ncbi genome数据库或根据patric id从patric数据库批量下载基因组数据。
74.在一些实施方式中，本技术方法基于基因组contig序列进行card耐药库的比对和
耐药基因(arg)注释统计。具体的，采用blastn(版本2.9.0 )将contig序列与card耐药基因参考序列库(v3.1.0)进行比对，先过滤掉identity小于90％或参比基因覆盖度小于90％的hit，然后将每条contig比对上的区域选择best hit(first hit)作为该contig区域的最终比对结果，并添加耐药基因的注释信息。
75.在一些实施方式中，针对各抗生素药物，进行基因型与耐药表型数据的关联分析，以筛选找出与耐药产生相关的重要特征基因。
76.在一些实施方式中，使用拉索回归模型进行关联分析，筛选得出与耐药表型相关的关键耐药基因，并计算得到这些基因的权重系数。
77.在一些具体的实施方式中，拉索回归分析方法如下：以基因检出分布矩阵(0-1型矩阵，0表示未检出耐药基因，1表示有检出该耐药基因)和抗生素药敏测试结果数据矩阵作为输入，使用glmnet程序包进行基因型和耐药表型数据的关联分析，并做诸如10重交叉验证，筛选得到与耐药表型相关的重要特征基因。优选的，所述重要特征基因具体为：根据不同个数特征基因下的模型cv错误率、auc变化曲线，选择在cv错误率最低且此时模型auc值相对稳定处所对应的数目基因为重要特征基因。
78.在一些实施方式中，为避免由于模型筛选可能漏掉的低频次发生的重要耐药基因，需要对与耐药表型ppv较高(》＝0.8)的基因进行人工召回，并基于以上得到的重要基因其权重系数值，同时计算得到召回基因的权重系数等。
79.在一些实施方式中，最后通过绘制roc曲线，确定基于当前筛选出来的重要耐药基因进行药敏预测的模型性能。
80.在一些优选的实施方式中，所述roc分析如下：基于步骤3)得到的重要基因权重系数的矩阵，定义并计算score值，并以此作为阴阳性判读指标，绘制roc曲线以及确定cut off值，用验证集样本进行验证和评估模型性能；所述其中arg_w

表示检出相应基因的权重系数值。
81.基于上述方法的描述，本技术还可以包含执行上述方法的模型、装置或存储介质等，这些都是本技术的保护范围。
82.另外，通过本技术方法，鉴定出一系列基因在不同药物耐药中的应用。
83.例如，在一些实施例中，发现基因kpc-3、kpc-2、oxa-181、oxa-232和ndm-1可作为非核心型耐药基因在肺炎克雷伯菌辅助药敏预测中的应用，比如耐药预测和敏感预测；当为亚胺培南时，所述基因还包括基因kpc-14、ndm-5、ndm-7、oxa-65、oxa-83、oxa-24、adc-214、adc-56、vim-27、vim-1、vim-19、veb-1中的一个或多个或全部；当为美罗培南时，所述基因还包括kpc-4、oxa-163、ndm-7、oxa-65、imp-4、adc-214、ges-19、ges-26、ges-20、vim-27、vim-1、vim-19、vim-4、oxa-24中的一个或多个或全部。
84.下面结合具体实施例来阐述本技术。
85.实施例1本技术设计优化
86.如本技术

背景技术：
部分所提及，现有细菌耐药筛选研究大都集中于对核心基因组层面的单核苷酸多态性(snp)、插入/缺失(indel)或k-mer特征的筛选，但除了这些特征外，筛选找到与耐药表型相关的重要非核心耐药基因也是十分重要的工作。
87.图1为本技术的设计思路，基于该思路，本实施例设计优化耐药筛选的具体方法。以模型选择为例展示本技术的确立过程。
88.在进行基因型与药敏结果数据的关联分析时，本技术对比使用了glm广义线性模型(r语言glm()、stepaic()命令)、拉索回归模型(r语言glmnet程序包)、随机森林(r语言randomforest程序包)三种分类模型方法的效果差异，具体参见图2，由图2可知，针对肺炎克雷伯菌，各抗生素药物模型性能。针对每个抗生素，有两个值，tr为训练集auc值，va为验证集auc值。auc值高低在单元格中通过不同色泽深度显示。
89.除了头孢吡肟(cpm)模型性能稍低些外，相对于逻辑回归和随机森林，其他抗生素的拉索模型都具有相类似或更高的性能表型(auc值)。并且，整体上，氨基糖苷类、β内酰胺类抗生素模型auc值较高，可以达到0.9以上。
90.因此，选择使用拉索回归模型方法，进行筛选得出与耐药表型相关的关键耐药基因，并计算得到这些基因的权重系数。
91.通过优化设计，最终确立本技术的方法步骤如下：
92.1、细菌菌株基因组的获取可以从公共数据库搜索下载得到已发表的目标基因组序列，或者通过现行收集临床培养鉴定出的细菌菌株进行测序组装获得，同时收集对应菌株的药敏测试结果信息。公共基因组数据的搜索和下载：从ncbi ndaro数据库(https://www.ncbi.nlm.nih.gov/pathogens/isolates)和patric数据库平台(https://patricbrc.org)搜集记录有药敏测试结果的细菌菌株信息，整理表型数据，并根据基因组组装id号从ncbi genome数据库或根据patric id从patric数据库批量下载基因组数据。本步骤要求：菌株基因组个数至少100，且菌株来源广泛多样涵盖各种亚型，耐药株与敏感株的菌株数目比例尽量均衡。
93.2、基于基因组contig序列进行card耐药库的比对和耐药基因(arg)注释统计。采用blastn(版本2.9.0 )将contig序列与card耐药基因参考序列库(v3.1.0)进行比对，先过滤掉identity小于90％或参比基因覆盖度小于90％的hit，然后将每条contig比对上的区域选择best hit(first hit)作为该contig区域的最终比对结果，并添加耐药基因的注释信息。
94.3、针对各抗生素药物，进行基因型与药敏结果数据的关联分析，以筛选找出与耐药产生相关的重要特征基因。使用拉索回归模型(r语言glmnet程序包)进行筛选得出与耐药表型相关的关键耐药基因，并计算得到这些基因的权重系数。
95.拉索回归分析方法为：以基因检出分布矩阵(0-1型矩阵，0表示未检出耐药基因，1表示有检出该耐药基因)和抗生素药敏测试结果数据矩阵作为输入，使用glmnet程序包进行基因型和耐药表型数据的关联分析，并做k(比如10重)重交叉验证，筛选得到与耐药表型相关的重要特征基因。
96.4、绘制roc曲线，确定基于当前筛选出来的重要耐药基因进行药敏预测的模型性能。
97.实施例2肺炎克雷伯菌耐碳青霉烯、头孢菌素相关中药特征基因的筛选和验证
98.步骤1、从公共数据库库，搜索下载肺炎克雷伯菌株基因组及其对应的抗生素药敏测试结果数据。
99.从ncbi ndaro数据库下载：打开网址https://www.ncbi.nlm.nih.gov/pathogens/isolates，在搜索栏输入”klebsiella pneumoniae”检索肺炎克雷伯菌信息，然后在matched isolates子窗口,点击”choose columns”选择”ast pheotypes”将此列信息
展示出来，接着下载整个窗口的表格数据，整理具有药敏测试结果数据的肺克菌株，根据assembly id信息，从ncbi的genome数据库(ftp://ftp.ncbi.nlm.nih.gov/genomes)批量下载基因组序列。
100.从patric平台数据库下载：打开网址https://patricbrc.org，在搜索窗口browse栏选择点击bacteria按钮，首先选择”amr phenotypes”，在keywords栏输入”klebsiella pneumoniae”进行筛选，同时过滤掉evidence列为”computational method”仅保留”laboratory method”条目，得到肺炎克雷伯菌菌株药敏信息并下载数据表格，然后选择”genome”,添加”assembly accession”列信息后下载数据表格。根据下载的两个表格里genome id对应找到具有药敏测试结果数据的菌株基因组patric id或assembly id，然后从patric或ncbi genome数据库(ftp://ftp.ncbi.nlm.nih.gov/genomes)批量下载基因组序列。
101.合并从ncbipathogen detetction和patric数据库搜集下载的基因组，并过滤掉冗余基因组，最后得到共3072株肺炎克雷伯菌基因组及其药敏测试结果数据，同时以随机方式，按照8：2的比例，将3072个阉割版你随机划分成两个子集，以分别用作模型训练集(n＝2410)和验证集(n＝662)。其中，针对碳青霉烯药物的菌株数目如下：
[0102][0103]
步骤2、基于下载的肺炎克雷伯菌基因组contig序列比对耐药数据库进行耐药基因检测。采用blastn(版本2.9.0 )将contig序列与card耐药基因参考序列库(v3.1.0)进行比对，过滤掉identity小于90％或参比基因覆盖度小于90％的hit，然后针对各contig上比对上的区域选择best hit(first hit)作为该contig区域的最终比对结果，并添加耐药基因的注释信息，统计各菌株耐药基因检出情况，最后汇总成0-1矩阵表格，其中0表示未检出耐药基因，1表示有检出该耐药基因。
[0104]
步骤3、基于训练样本进行基因型和抗生素耐药表型数据关联分析。以亚胺培南为例，其他药物类似。依据card库里记录的各耐药基因对应的抗生素药物分类信息，基于步骤2得到的基因检出矩阵表格，挑选出碳青霉烯类相关耐药基因的子矩阵表格，并过滤掉检出频率较低(优选的，检出频率小于3个)、ppv较低(优选的，ppv小于0.6)的基因，然后用过滤后的表格数据进行与亚胺培南药敏结果的关联分析。子矩阵表格数据(记为x)格式如下(数据多，仅呈现部分)：
[0105]
[0106][0107]
亚胺培南药敏结果数据(记为y)格式如下：
[0108]
argastgca_000401195.1rgca_000406385.1rgca_000406405.1rgca_000406425.1rgca_000406445.1rgca_000406465.1r... [0109]
以上面两个数据集(x和y)为输入，使用r语言glmnet程序包进行关联分析，并进行10重交叉验证，以筛选得到与亚胺培南耐药相关的基因。部分运行程序代码如下：
[0110]
library(glmnet)
[0111]
cv.model《-cv.glmnet(x,y,family＝"binomial",nlambda＝100,alpha＝1,standardize＝f,nfolds＝10,type.measure＝"class")
[0112]
coefficients《-coef(cv.model,s＝cv.model$lambda.min)
[0113]
library(proc)
[0114]
predict《-predict(cv.model,trainx,s＝cv.model$lambda.min,type＝"response")
[0115]
roc.predict《-roc(y,as.numeric(predict))
[0116]
...
[0117]
运行上诉程序后，分析得到与亚胺培南耐药相关的基因。鉴于初步分析得到的基因数目往往较多，所以需要对这些基因按照重要性进行排序，然后最终挑取排序靠前的重要基因。重要基因挑取过程如下：
[0118]
基于初步分析得到的排序后的基因，设置梯度选取不同数目的基因组合，按照以上程序代码构建拉索回归模型，获得模型auc值、cv错误率，并计算auc-error差值，然后绘制曲线图。按照auc-error值逐渐增加过程中第一次出现拐点下降处对应的横坐标数，或者按照auc第一次出现最大值处或error第一次出现最小值处对应的横坐标数目以选择作为最终选取的基因数目(见图3)。其次，为避免由于模型筛选可能漏掉的低频次发生的重要耐药基因，需要对与耐药表型ppv较高(》＝0.8)的基因进行人工召回，并基于以上得到的重要基因其权重系数值，同时计算得到召回基因的权重系数(计算公式：基因其权重系数值，同时计算得到召回基因的权重系数(计算公式：其中，arg_n表示检出相应耐药基因的样本数，arg_w表示相应耐药基因的权重系数，j表示与目标基因归属同家族的有检出的j个耐药基因，ppv
arg
表示目标基因与耐药表型的阳性一致率)。
[0119]
这里，针对亚胺培南最终确定基于机器学习模型筛选出7个重要基因，并通过召回得到另外9个重要基因。
[0120]
综上，最终筛选出的与肺炎克雷伯菌耐碳青霉烯类药物相关的重要基因及其权重系数如下：
[0121]
2、oxa-181、oxa-163、oxa-232、ndm-1、ndm-7、oxa-65、imp-4、adc-214、ges-19、ges-26、ges-20、vim-27、vim-1、vim-19、vim-4、oxa-24。
[0125]
步骤4、roc分析确定基于筛选出的特征基因构建的分类模型性能
[0126]
定义score指标(其中arg_w

表示检出相应基因的权重系数值)，并以此作为阴阳性判读指标。针对亚胺培南，基于以上模型筛选得到的重要基因权重系数矩阵，结合样本耐药基因实际检出情况，计算各样本的score值，然后进行roc曲线分析，得到训练集(n＝1804)模型auc值为0.959(图4左)。然后进一步用验证集(n＝472)进行roc分析，得到auc值为0.943(图4右)。训练集和验证集模型auc较高，表明本技术方法具有较好的性能，即本技术的筛选方法和模型准确有效。
[0127]
同样的，针对美罗培南，训练集(n＝1981)和验证集(n＝540)auc值分别为0.948、0.929(图5)。
[0128]
因此，基于本发明的模型，根据表1结果数据可知，通过对如下基因进行检测时，结合模型得到的基因权重系数，则可实现预测相应抗生素的药敏结果。
[0129]
针对亚胺培南，目标耐药基因包括kpc-3、kpc-14、kpc-2、oxa-181、oxa-232、ndm-1、ndm-5、ndm-7、oxa-65、oxa-83、oxa-24、adc-214、adc-56、vim-27、vim-1、vim-19、veb-1；而综合基因权重、基因发生频次以及耐药产生的所有可能机制等因素，实践中可通过对于主要介导耐药产生的高频次发生且权重较高的kpc-3、kpc-2、oxa-181、oxa-232、ndm-1基因同时进行检测，若检测结果均为阴性，则可推测为敏感，即实现药敏，尤其是敏感检测目的；而为提升检测精度等，可进一步包括上述kpc-14、ndm-5、ndm-7、oxa-65、oxa-83、oxa-24、adc-214、adc-56、vim-27、vim-1、vim-19、veb-1中的一个或多个，优选全部。
[0130]
针对美罗培南，目标基因包括：kpc-3、kpc-4、kpc-2、oxa-181、oxa-163、oxa-232、ndm-1、ndm-7、oxa-65、imp-4、adc-214、ges-19、ges-26、ges-20、vim-27、vim-1、vim-19、vim-4、oxa-24。而综合基因权重、基因发生频次以及耐药产生的所有可能机制等因素，可知，实践中通过对于主要介导耐药产生的高频次发生且权重较高的kpc-3、kpc-2、oxa-181、oxa-232、ndm-1基因同时进行检测，若检测结果均为阴性，则可推测为敏感，即可实现药敏，尤其是敏感检测目的；而为提升检测精度等，可进一步包括上述kpc-4、oxa-163、ndm-7、oxa-65、imp-4、adc-214、ges-19、ges-26、ges-20、vim-27、vim-1、vim-19、vim-4、oxa-24中的一个或多个，优选全部。
[0131]
以上对本技术具体实施方式的描述并不限制本技术，本领域技术人员可以根据本技术作出各种改变或变形，只要不脱离本技术的精神，均应属于本技术所附权利要求的范围。