CN116825191A - 一种筛选微生物细菌关键调控qtl的方法 - Google Patents
一种筛选微生物细菌关键调控qtl的方法 Download PDFInfo
- Publication number
- CN116825191A CN116825191A CN202310749301.3A CN202310749301A CN116825191A CN 116825191 A CN116825191 A CN 116825191A CN 202310749301 A CN202310749301 A CN 202310749301A CN 116825191 A CN116825191 A CN 116825191A
- Authority
- CN
- China
- Prior art keywords
- quantitative trait
- growth
- microorganism
- individual
- quantitative
- 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.)
- Pending
Links
- 241000894006 Bacteria Species 0.000 title claims abstract description 48
- 244000005700 microbiome Species 0.000 title claims abstract description 48
- 238000000034 method Methods 0.000 title claims abstract description 47
- 230000033228 biological regulation Effects 0.000 title claims abstract description 23
- 238000012216 screening Methods 0.000 title claims abstract description 15
- 230000002068 genetic effect Effects 0.000 claims abstract description 55
- 230000000813 microbial effect Effects 0.000 claims abstract description 42
- 230000001580 bacterial effect Effects 0.000 claims abstract description 34
- 230000008569 process Effects 0.000 claims abstract description 10
- 238000004458 analytical method Methods 0.000 claims abstract description 9
- 230000001105 regulatory effect Effects 0.000 claims description 37
- 238000013507 mapping Methods 0.000 claims description 14
- 239000013598 vector Substances 0.000 claims description 12
- 230000001276 controlling effect Effects 0.000 claims description 10
- 239000002773 nucleotide Substances 0.000 claims description 6
- 125000003729 nucleotide group Chemical group 0.000 claims description 6
- 238000003766 bioinformatics method Methods 0.000 abstract description 2
- 241000191967 Staphylococcus aureus Species 0.000 description 45
- 230000006870 function Effects 0.000 description 20
- 108010059993 Vancomycin Proteins 0.000 description 12
- MYPYJXKWCTUITO-LYRMYLQWSA-N vancomycin Chemical compound O([C@@H]1[C@@H](O)[C@H](O)[C@@H](CO)O[C@H]1OC1=C2C=C3C=C1OC1=CC=C(C=C1Cl)[C@@H](O)[C@H](C(N[C@@H](CC(N)=O)C(=O)N[C@H]3C(=O)N[C@H]1C(=O)N[C@H](C(N[C@@H](C3=CC(O)=CC(O)=C3C=3C(O)=CC=C1C=3)C(O)=O)=O)[C@H](O)C1=CC=C(C(=C1)Cl)O2)=O)NC(=O)[C@@H](CC(C)C)NC)[C@H]1C[C@](C)(N)[C@H](O)[C@H](C)O1 MYPYJXKWCTUITO-LYRMYLQWSA-N 0.000 description 12
- 229960003165 vancomycin Drugs 0.000 description 12
- MYPYJXKWCTUITO-UHFFFAOYSA-N vancomycin Natural products O1C(C(=C2)Cl)=CC=C2C(O)C(C(NC(C2=CC(O)=CC(O)=C2C=2C(O)=CC=C3C=2)C(O)=O)=O)NC(=O)C3NC(=O)C2NC(=O)C(CC(N)=O)NC(=O)C(NC(=O)C(CC(C)C)NC)C(O)C(C=C3Cl)=CC=C3OC3=CC2=CC1=C3OC1OC(CO)C(O)C(O)C1OC1CC(C)(N)C(O)C(C)O1 MYPYJXKWCTUITO-UHFFFAOYSA-N 0.000 description 12
- 239000011159 matrix material Substances 0.000 description 9
- 230000007423 decrease Effects 0.000 description 8
- 230000000694 effects Effects 0.000 description 7
- 108090000623 proteins and genes Proteins 0.000 description 7
- 238000009826 distribution Methods 0.000 description 6
- 230000003247 decreasing effect Effects 0.000 description 4
- 102000054765 polymorphisms of proteins Human genes 0.000 description 4
- 238000012360 testing method Methods 0.000 description 4
- 238000004364 calculation method Methods 0.000 description 3
- 230000007614 genetic variation Effects 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 2
- 238000009395 breeding Methods 0.000 description 2
- 230000001488 breeding effect Effects 0.000 description 2
- 238000012937 correction Methods 0.000 description 2
- 230000001419 dependent effect Effects 0.000 description 2
- 239000003814 drug Substances 0.000 description 2
- 230000008303 genetic mechanism Effects 0.000 description 2
- 230000035772 mutation Effects 0.000 description 2
- 238000012546 transfer Methods 0.000 description 2
- 229920001817 Agar Polymers 0.000 description 1
- 108020000946 Bacterial DNA Proteins 0.000 description 1
- 108020004414 DNA Proteins 0.000 description 1
- 238000007400 DNA extraction Methods 0.000 description 1
- 102000004163 DNA-directed RNA polymerases Human genes 0.000 description 1
- 108090000626 DNA-directed RNA polymerases Proteins 0.000 description 1
- 241000588724 Escherichia coli Species 0.000 description 1
- 208000020584 Polyploidy Diseases 0.000 description 1
- 230000006978 adaptation Effects 0.000 description 1
- 239000008272 agar Substances 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000011681 asexual reproduction Effects 0.000 description 1
- 238000013465 asexual reproduction Methods 0.000 description 1
- 238000012098 association analyses Methods 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 238000012512 characterization method Methods 0.000 description 1
- 230000005770 chromosome separation Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 201000010099 disease Diseases 0.000 description 1
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 1
- 229940079593 drug Drugs 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 230000006801 homologous recombination Effects 0.000 description 1
- 238000002744 homologous recombination Methods 0.000 description 1
- 238000000338 in vitro Methods 0.000 description 1
- 230000002401 inhibitory effect Effects 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 238000012417 linear regression Methods 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 230000002503 metabolic effect Effects 0.000 description 1
- 230000002906 microbiologic effect Effects 0.000 description 1
- 230000000877 morphologic effect Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 230000000750 progressive effect Effects 0.000 description 1
- 238000003908 quality control method Methods 0.000 description 1
- 230000033458 reproduction Effects 0.000 description 1
- 101150090202 rpoB gene Proteins 0.000 description 1
- 238000012163 sequencing technique Methods 0.000 description 1
- 238000012070 whole genome sequencing analysis Methods 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
- G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
- G16B20/40—Population genetics; Linkage disequilibrium
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/23—Clustering techniques
- G06F18/232—Non-hierarchical techniques
- G06F18/2321—Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Theoretical Computer Science (AREA)
- Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Mathematical Physics (AREA)
- Health & Medical Sciences (AREA)
- Computational Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Optimization (AREA)
- Bioinformatics & Computational Biology (AREA)
- Mathematical Analysis (AREA)
- Evolutionary Biology (AREA)
- Genetics & Genomics (AREA)
- General Health & Medical Sciences (AREA)
- Physiology (AREA)
- Analytical Chemistry (AREA)
- Probability & Statistics with Applications (AREA)
- Chemical & Material Sciences (AREA)
- Artificial Intelligence (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Evolutionary Computation (AREA)
- Medical Informatics (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Biophysics (AREA)
- Biotechnology (AREA)
- Proteomics, Peptides & Aminoacids (AREA)
- Ecology (AREA)
- Operations Research (AREA)
- Molecular Biology (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
Abstract
本发明提供了一种筛选微生物细菌关键调控QTL的方法，涉及生物信息学分析领域，该方法包括根据微生物细菌样本个体的数量性状的表型数据得到Logistic方程估计参数；根据微生物细菌样本个体的基因型和表型数据、Logistic方程估计参数和功能作图模型框架，定位调控微生物细菌样本个体数量性状生长的显著数量性状位点；根据显著数量性状位点建立不同显著数量性状位点之间的线性相关关系，得到数量性状位点调控关系网络；根据数量性状位点调控关系网络识别解释微生物细菌样本个体数量性状生长过程的关键调控数量性状位点。本发明能够对显著的QTL调控生长的遗传控制以及QTL之间的调控网络关系进行深入分析。
Description
技术领域
本发明涉及生物信息学分析领域，特别是涉及一种筛选微生物细菌关键调控QTL的方法。
背景技术
微生物细菌的数量性状在微生物学研究中具有重要的意义，它们在环境适应、代谢能力和生物合成等方面表现出差异。这些性状往往不能通过明显的形态特征来区分，只能通过数量上的差异来揭示微生物细菌个体之间的变异。微生物细菌的数量性状受到多基因联合控制的遗传过程影响，控制这些数量性状表达的基因位点被称为数量性状位点(quantitative traitlocus，简称QTL)。建立数量性状表型与基因型之间的生物学关联，筛选调控微生物细菌性状生长的QTL对从遗传层面揭示微生物生长机制具有十分重要的作用。通过研究显著的QTL对数量性状的遗传控制模式、建立QTL之间的调控关系，可以为微生物细菌的分子遗传控制育种提供可靠的依据。这有助于深入了解微生物细菌的遗传机制，优化微生物育种策略，并提高微生物在农业、环境和工业等领域的应用效果。
基因组关联研究(GWAS)是一种用于寻找遗传变异与表型特征之间关联的方法，GWAS已成为遗传学家越来越重要的研究方法。GWAS方法使得复杂疾病的表征和与特定特征相关的常见遗传变异的鉴定成为可能。然而，与人类和其他多倍体生物相比，遗传变异在微生物细菌中表现出不同的方式。人类繁殖过程中存在同源重组和染色体分离，而微生物细菌是单倍体和无性繁殖的，其群体结构为高度结构化，存在水平基因转移和反复突变的特点。因此，在微生物细菌中开发和应用GWAS的进展相对较慢，尽管已经有一些成功应用的例子，例如：Alam等人利用GWAS鉴定与金黄色葡萄球菌的耐药表型相关的RNA聚合酶rpoB基因的突变，He等人利用GWAS探索大肠杆菌和金黄色葡萄球菌基因之间的相互作用。
但是，对微生物细菌关键数量性状的遗传机制的建立并不全面，缺少针对显著的QTL调控生长的遗传控制以及QTL之间的调控网络关系的深入分析。
发明内容
本发明的目的是提供一种筛选微生物细菌关键调控QTL的方法，能够对显著的QTL调控生长的遗传控制以及QTL之间的调控网络关系的深入分析，提供了分析基因位点之间的相互联系的有效手段。
为实现上述目的，本发明提供了如下方案：
本发明提供了一种筛选微生物细菌关键调控QTL的方法，包括：
根据具有相同生活环境的多个微生物细菌样本个体的数量性状的表型数据拟合Logistic方程，得到Logistic方程的估计参数；
将每个微生物细菌样本个体的基因型和表型数据作为分析数据，将Logistic方程的估计参数应用到功能作图的模型框架中，定位调控微生物细菌样本个体数量性状生长的显著数量性状位点；
根据显著数量性状位点，计算数量性状的遗传效应值，并根据遗传效应值建立不同显著数量性状位点之间的线性相关关系，得到数量性状位点调控关系网络；
根据数量性状位点调控关系网络识别解释微生物细菌样本个体数量性状生长过程的关键调控数量性状位点。
根据本发明提供的具体实施例，本发明公开了以下技术效果：
本发明利用Logistic方程拟合微生物细菌数量性状的生长曲线，得到Logistic方程的估计参数。同时利用功能作图的模型框架定位调控微生物细菌数量性状生长的显著数量性状位点，该方法在不损失精度的同时提高了计算效率。然后根据显著数量性状位点的遗传效应建立数量性状位点之间的遗传调控网络，用以识别解释微生物细菌数量性状生长过程的关键调控数量性状位点，因此本发明提供了分析基因位点之间的相互联系的有效手段。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例提供的一种筛选微生物细菌关键调控QTL的方法的流程图；
图2为本发明示例1中金黄色葡萄球菌丰度的Logistic生长曲线图；
图3为本发明示例1中金黄色葡萄球菌全基因组显著性检验的曼哈顿图；
图4为本发明示例1中金黄色葡萄球菌丰度的遗传效应随时间变化的趋势图；
图5为本发明示例1中调节金黄色葡萄球菌丰度生长的显著数量性状位点的调控网络图；
图6为本发明示例2中金黄色葡萄球菌丰度的Logistic生长曲线图；
图7为本发明示例2中金黄色葡萄球菌全基因组显著性检验的曼哈顿图；
图8为本发明示例2中金黄色葡萄球菌丰度的遗传效应随时间变化的趋势图；
图9为本发明示例2中调节金黄色葡萄球菌丰度生长的显著数量性状位点的调控网络图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
实施例一
如图1所示，本实施例提供了一种基于Logistic方程的筛选微生物细菌关键调控QTL的方法，包括以下步骤。
步骤100：根据具有相同生活环境的多个微生物细菌样本个体的数量性状的表型数据拟合Logistic方程，得到Logistic方程的估计参数。
在本实施例中，步骤100具体包括：
采用最小二乘法和具有相同生活环境的微生物细菌样本个体的数量性状的表型数据拟合Logistic方程，并采用BFGS拟牛顿法找到使得表型值与拟合值残差平方和最小的估计参数，然后将表型值与拟合值残差平方和最小的估计参数确定为Logistic方程的估计参数。
其中，Logistic方程如公式(1)所示。
y＝A/(1+k*e-λt)。
式中，A表示最大生长量，k表示生长速率，λ表示生长的延迟时间，t表示时间，y表示数量性状的拟合值。
在本实施例中，在执行步骤100之前，需要连续测定每个微生物细菌样本个体的生长数据，得到具有相同生活环境的多个微生物细菌样本个体的数量性状的表型数据。
步骤200：将每个微生物细菌样本个体的基因型和表型数据作为分析数据，将Logistic方程的估计参数应用到功能作图的模型框架中，定位调控微生物细菌样本个体数量性状生长的显著数量性状位点。
其中，每个微生物细菌样本个体的基因型的确定过程为：对每个微生物细菌样本个体的全基因组进行单核苷酸多态性分型，得到每个微生物细菌样本个体的基因型。
在本实施例中，步骤200具体包括：
1)建立原假设和备择假设；原假设为微生物细菌样本个体的数量性状的生长在不同基因型之间没有差异；备择假设为微生物细菌样本个体的数量性状的生长在不同基因型之间存在差异。
2)建立原假设的似然函数与备择假设的似然函数；其中，似然函数包括用来构建均值向量的Logistic方程的估计参数以及用来构建构造协方差矩阵的结构参数。
3)根据原假设的似然函数与备择假设的似然函数建立似然比统计量函数。
4)根据微生物细菌样本个体的基因型和表型数据、以及似然比统计量函数定位调控微生物细菌样本个体数量性状生长的显著数量性状位点。
在本实施例中，所述功能作图的模型框架如下：
微生物细菌样本个体i的数量性状在观测时间1,…,T的生长量表示为yi＝(yi(1),…,yi(T))；数量为n的微生物细菌样本个体总体近似服从正态分布，该正态分布的均值向量μ和协方差矩阵Σ表示为：
μ＝(μ(1),…,μ(T))；
其中，协方差矩阵Σ中对角线上的元素为微生物细菌样本个体数量性状与时间相关的方差，协方差矩阵Σ中非对角上的元素为微生物细菌样本个体数量性状在不同时间的协方差。
正态分布的概率密度函数表示为f(yi；Θ,Ψ)，构成n个微生物细菌样本个体的似然函数如公式(2)所示。
式中，Θ为用来构建均值向量μ的Logistic方程的估计参数(A,k,λ)，Ψ为用来构建构造协方差矩阵Σ的结构参数(φ,γ)。
假设微生物细菌样本个体数量性状的生长由基因组中一组显著数量性状位点调控，假设某个显著数量性状位点的基因型数量为J；基因型为j的微生物细菌样本个体的数量为nj，满足
每种基因型的微生物细菌样本个体的表型数据近似服从一组正态分布，每组正态分布的概率密度函数表示为fj(yi；Θj,Ψ)，(j＝1,…,J)，服从正态分布的基因型为j的微生物细菌样本个体的均值向量(μj)如公式(3)所示。
μj＝(μj(1),…,μj(T)) (3)。
不同基因型微生物细菌样本个体的概率密度函数fj(yi；Θj,Ψ)构成似然函数，如公式(4)所示。
式中，Θj表示Logistic方程参数(Aj,kj,λj)，用来构建均值向量μj。
协方差矩阵Σ优选由一阶向前依赖结构模型的创新方差γ2和一阶前依赖参数φ构建得到；协方差矩阵Σ中对角线上表示方差的元素如公式(5)所示，协方差矩阵Σ中非对角线上表示不同时间点性状协方差的元素如公式(6)所示。
上述似然值L0(y)和L1(y)的计算，根据如下假设性检验，即公式(7)实现。
H0:Θj＝Θversus H1:Θj≠Θ，j＝1,…,J (7)。
原假设H0假设数量性状的生长在基因型之间不存在差异，因而表示生长的Logistic参数在不同基因型下是相等的。而备择假设H1假设在该基因点具有不同基因型的微生物细菌样本个体数量性状的生长存在差异，因而不同基因型的Logistic生长参数不相等。
根据似然值L0(y)和L1(y)建立似然比统计量LR，具体如公式(8)所示。
LR＝-2(logL0(y)-logL1(y)) (8)。
LR是近似服从χ2分布的统计量，其自由度为原假设和备择假设参数个数的差异。LR的拒绝域表示为W＝{LR≥c}，其中临界值c满足P(LR≥c)≤α，α为检验水平。如果p≤α，则LR属于拒绝域，即拒绝原假设H0接受备择假设H1，可以判定在该基因位点的不同基因型性状生长存在差异，该基因位点是显著的数量性状位点。
步骤300：根据显著数量性状位点，计算数量性状的遗传效应值，并根据遗传效应值建立不同显著数量性状位点之间的线性相关关系，得到数量性状位点调控关系网络。
在本实施例中，数量性状的遗传效应值的计算公式为：
其中，a(t)观测时间为t时数量性状的遗传效应值，J为个显著数量性状位点的基因型数量，nj为基因型为j的微生物细菌样本个体的数量，μj(t)为观测时间为t时基因型为j的微生物细菌样本个体的均值。
在本实施例中，不同显著数量性状位点之间的线性相关关系如公式(10)所示，表示第i个显著数量性状位点与其他i-1个数量性状位点之间的调控关系：
Ei＝β1E1+…+βi-1Ei-1+βi+1Ei+1+…+βpEp+ε (10)。
其中，Ei表示第i个显著数量性状位点的遗传效应值，i＝1,…,p；ε～N(0,σ2)，β1,…,βi-1,βi+1,…,βp和σ2是未知参数；
(E1(t),…,Ei-1(t),Ei(t),Ei+1(t),…,Ep(t))(t＝1,…,T)是遗传效应在T个时间点的计算值，则公式(10)可以表示为矩阵形式，具体如下：
其中εi～N(0,σ2)，且独立同分布；ε＝[ε1,…,εi-1,εi+1,…,εp]'，εi是p维误差向量，并且满足：E(ε)＝0，Var(ε)＝σ2In。
[β1，…，βi-1，βi+1，…，βp]′是p维参数向量，通过检验回归方程的参数向量判断显著数量性状位点之间是否存在调控关系：当参数向量中的某个元素βj,j∈(1,…,i-1,i+1,…,p)为0时，显著数量性状位点i受到数量性状位点j的调控作用，反之则存在调控关系。
数量性状位点调控关系网络的建立是由R语言中多元线性回归函数lm实现。
步骤400：根据数量性状位点调控关系网络识别解释微生物细菌样本个体数量性状生长过程的关键调控数量性状位点。
在本实施例中，微生物细菌样本个体数量性状生长过程的关键调控数量性状位点的识别方法为：将每一个数量性状位点调控其余数量性状位点的数量进行统计，选择其中调控数量较大的数量性状位点称为关键调控数量性状位点，关键调控数量性状位点数量与数量性状位点调控关系网络中全部的数量性状位点数量相关；在本实施例中，选择标准为在数量性状位点调控关系网络中全部的数量性状位点数量的前5％。
本实施例提供的方法对所有种类的微生物细菌均适用。为了具体说明本实施例提供的方法，下面示例以金黄色葡萄球菌作为微生物细菌的代表进行分析，但这不能理解为对本实施例提供的方法有限制作用。所述数量性状指相对性状之间没有明显的界限，个体间表现的差异只能用数量来区别，变异呈连续性的性状。
示例1
采用https://github.com/CCBBeijing/PPMultilayerNetwork.中公开的金黄色葡萄球菌数据为例对本实施例做进一步详细描述。通过对亲本菌株60天的体外万古霉素处理，获得了所有金黄色葡萄球菌株。所有原始亲本菌株来源于中国工业微生物菌种保藏中心、中国农业菌种保藏中心、中国林业微生物菌种保藏与管理中心、中国一般微生物菌种保藏中心、中国药品微生物菌种保藏与管理中心、中国典型培养物保藏中心、中国医学微生物菌种保藏与管理中心、中国农业大学和北京朝阳医院。这些菌株被存储在-80℃的低温冰箱中。每个亲本菌株在含有1/2最初最小抑菌浓度(MIC)的万古霉素的BHI琼脂平板上接种，并在相同浓度的培养基中进行24小时转移培养，连续进行4天。每4天重新测定菌株的MIC，并使用更新的MIC重复处理。在对照组和压力组中接种每个菌株，测量其在培养后1小时的OD600值，随后每2小时测量一次，直到12小时后每4小时测量一次，24小时后每6小时测量一次，直到48小时，即总共进行了14个时间点的测量。使用TIANamp细菌DNA提取试剂盒(天津生物，中国北京)提取菌株的基因组DNA，按照制造商的协议进行操作。使用Illumina HiSeq4000测序仪(Illumina Inc.，美国圣地亚哥，加利福尼亚州)在Allwegene(中国北京)进行了全基因组测序。为了获得初始的比对结果，使用BWA mapper v0.7.899将测序数据与参考基因组(金黄色葡萄球菌亚种NCTC 8325)进行比对。使用SAMtools v0.1.18100对比对结果进行排序，得到110678个SNP，其中25173个通过质量控制用于关联分析。
1.利用拟牛顿法对0μg/ml万古霉素环境中的金黄色葡萄球菌样本的生长数据进行拟合，得到具有相同生活环境的多个金黄色葡萄球菌样本的数量性状的表型数据；利用金黄色葡萄球菌样本的数量性状的表型数据拟合Logistic方程，搜索得到Logistic方程的估计参数。
其中，通过最小二乘法拟合Logistic方程；估计参数的搜索方法为BFGS拟牛顿法，找到使得金黄色葡萄球菌表型值与拟合值残差平方和最小的Logistic估计参数，具体通过R语言中通用优化方法optim函数实现。
图2为在0μg/ml万古霉素环境中金黄色葡萄球菌丰度的Logistic生长曲线，其中浅色曲线是99个样本的生长曲线，深色曲线是平均生长曲线。金黄色葡萄球菌丰度数据拟合Logistic方程的决定系数为0.9979，说明Logistic生长曲线对样本数据的拟合程度较好，曲线具有良好的拟合优度。平均生长曲线的Logistic参数A，k，λ分别为1.5295，29.0908，0.2898，平均曲线拟合的残差平方和为0.0105。
2.基于功能作图的模型框架对调控金黄色葡萄球菌丰度生长的显著数量性状位点进行定位。
以每个金黄色葡萄球菌样本的基因型和表型数据作为分析数据，以上述采用功能作图的模型框架定位调控金黄色葡萄球菌样本数量性状生长的显著数量性状位点。其中，功能作图的模型框架如上述实施例所述。
图3为根据样本生长曲线进行功能作图，对金黄色葡萄球菌全基因组单核苷酸多态性(SNPs)进行显著性检验的曼哈顿图像。水平线为0.05水平下FDR矫正后的阈值。筛选出69个单核苷酸多态性调控金黄色葡萄球菌丰度的生长。显著数量性状位点的SNP序号以及p值信息均在表1中。
3.根据检测得到的69个显著数量性状位点，计算数量性状的遗传效应。
图4为0μg/ml万古霉素环境中金黄色葡萄球菌丰度的遗传效应随时间变化的趋势，图4中实线表示69个显著数量性状位点的遗传效应时间变化趋势曲线。这些显著数量性状位点的遗传效应曲线主要可以分为效应值先增加后减小、效应值先增后减再增再减这两类。其中SNP 3370的遗传效应在初始的4小时内减少，在4到12小时内逐渐增大达到峰值后，在12-48小时内缓慢下降。其在第10小时的遗传效应显著高于其他数量性状位点，在30小时后的遗传效应也远高于其他数量性状位点。SNP 2772的遗传效应曲线先增加后减小，在最开始的6小时内显著高于其他数量性状位点的遗传效应。
4.建立显著数量性状位点之间基于遗传效应的调控关系网络，识别解释0μg/ml万古霉素环境下金黄色葡萄球菌生长过程的关键调控数量性状位点。
根据调节金黄色葡萄球菌丰度生长的显著数量性状位点的遗传效应，建立数量性状位点之间的线性方程组。由方程的系数表示数量性状位点之间的调控关系。图5为69个数量性状位点的遗传效应的基因调控网络，其中4个数量性状位点(1、2、7、8)：SNP 19481、SNP21485、SNP 32318、SNP 43762为遗传效应网络中的关键调控数量性状位点，对金黄色葡萄球菌丰度生长的遗传结构中起着关键作用。网络关系中数量性状位点的编号表示在表1中。
表1调控0μg/ml万古霉素环境下金黄色葡萄球菌样本生长的显著数量性状位点信息表
| No. | SNP | P-value | No. | SNP | P-value | No. | SNP | P-value |
| 1 | 19481 | 1.12×10-10 | 24 | 220967 | 5.76×10-13 | 47 | 1937320 | 3.82×10-11 |
| 2 | 21485 | 4.09×10-12 | 25 | 229384 | 1.53×10-26 | 48 | 1937407 | 1.87×10-15 |
| 3 | 25261 | 6.62×10-17 | 26 | 232510 | 5.68×10-11 | 49 | 1942424 | 5.26×10-15 |
| 4 | 28321 | 2.99×10-11 | 27 | 250761 | 7.91×10-12 | 50 | 1961496 | 1.83×10-11 |
| 5 | 31826 | 5.26×10-21 | 28 | 258855 | 1.35×10-12 | 51 | 2072412 | 8.79×10-19 |
| 6 | 31877 | 1.53×10-12 | 29 | 303477 | 4.31×10-11 | 52 | 2122872 | 1.44×10-10 |
| 7 | 32318 | 7.30×10-15 | 30 | 303479 | 9.22×10-12 | 53 | 2572670 | 5.40×10-11 |
| 8 | 43762 | 1.44×10-11 | 31 | 306103 | 4.44×10-14 | 54 | 2644158 | 3.47×10-25 |
| 9 | 56487 | 7.26×10-12 | 32 | 316281 | 1.70×10-10 | 55 | 2667695 | 5.85×10-11 |
| 10 | 73888 | 1.32×10-11 | 33 | 319914 | 4.15×10-13 | 56 | 2673960 | 2.94×10-14 |
| 11 | 89719 | 7.49×10-15 | 34 | 324105 | 9.22×10-12 | 57 | 2728755 | 1.62×10-17 |
| 12 | 92210 | 9.69×10-24 | 35 | 405662 | 3.10×10-11 | 58 | 2748895 | 1.22×10-18 |
| 13 | 94818 | 7.40×10-15 | 36 | 448470 | 2.39×10-14 | 59 | 2749269 | 7.13×10-14 |
| 14 | 96356 | 5.24×10-11 | 37 | 454194 | 1.45×10-11 | 60 | 2755050 | 7.19×10-11 |
| 15 | 119306 | 1.14×10-12 | 38 | 550323 | 6.62×10-17 | 61 | 2773236 | 2.85×10-12 |
| 16 | 155745 | 5.15×10-14 | 39 | 590287 | 1.22×10-10 | 62 | 2783126 | 1.21×10-14 |
| 17 | 164483 | 3.25×10-12 | 40 | 768337 | 2.35×10-12 | 63 | 2783127 | 2.85×10-12 |
| 18 | 183439 | 1.07×10-11 | 41 | 768340 | 1.32×10-16 | 64 | 2785600 | 1.96×10-10 |
| 19 | 183485 | 3.54×10-17 | 42 | 775710 | 7.12×10-12 | 65 | 2789423 | 2.85×10-12 |
| 20 | 183488 | 1.36×10-13 | 43 | 942441 | 5.34×10-11 | 66 | 2804537 | 3.12×10-11 |
| 21 | 193712 | 8.29×10-18 | 44 | 1076587 | 1.77×10-13 | 67 | 2804607 | 1.43×10-12 |
| 22 | 199527 | 1.48×10-10 | 45 | 1924463 | 9.43×10-11 | 68 | 2810211 | 2.54×10-12 |
| 23 | 213345 | 5.34×10-11 | 46 | 1928590 | 1.72×10-14 | 69 | 2815318 | 5.77×10-16 |
示例2
1.利用拟牛顿法对6μg/ml万古霉素环境中的金黄色葡萄球菌样本的生长数据进行拟合，得到具有相同生活环境的多个金黄色葡萄球菌样本的数量性状的表型数据；利用金黄色葡萄球菌样本的数量性状的表型数据拟合Logistic方程，搜索得到Logistic方程的估计参数。
图6为在6μg/ml万古霉素环境中金黄色葡萄球菌丰度的Logistic生长曲线，其中浅色曲线是99个样本的生长曲线，深色曲线是平均生长曲线。金黄色葡萄球菌丰度数据拟合Logistic方程的决定系数为0.9861，说明Logistic生长曲线对样本数据的拟合程度较好，曲线具有良好的拟合优度。平均生长曲线的Logistic参数A，k，λ分别为0.9567，55.114114，0.198217，平均曲线拟合的残差平方和为0.0254。
2.基于功能作图的模型框架对调控金黄色葡萄球菌丰度生长的显著数量性状位点进行定位。
以每个金黄色葡萄球菌样本的基因型和表型数据作为分析数据，以上述采用功能作图的模型框架定位调控金黄色葡萄球菌样本数量性状生长的显著数量性状位点。其中，功能作图的模型框架如上述实施例所述。
图7为根据样本生长曲线进行功能作图，对金黄色葡萄球菌全基因组单核苷酸多态性(SNPs)进行显著性检验的曼哈顿图像。水平线为0.05水平下FDR矫正后的阈值。筛选出69个单核苷酸多态性调控金黄色葡萄球菌丰度的生长。显著数量性状位点的SNP序号以及p值信息均在表1中。
3.根据检测得到的49个显著数量性状位点，计算数量性状的遗传效应。
图8为6μg/ml万古霉素环境中金黄色葡萄球菌丰度的遗传效应随时间变化的趋势，图8中实线表示49个显著数量性状位点的遗传效应时间变化趋势曲线。这些显著数量性状位点的遗传效应曲线主要可以分为效应值先增加后减小再增加、效应值先增后减再增再减、效应值先增再减这三类。其中SNP 1749325的遗传效应在最开始的8小时内显著高于其他数量性状位点的遗传效应，SNP 213585在30小时后的遗传效应远高于其他数量性状位点。这两个SNP的遗传效应的趋势都是先增加后减小再增加再减小。
4.建立显著数量性状位点之间基于遗传效应的调控关系网络，识别解释6μg/ml万古霉素环境下金黄色葡萄球菌生长过程的关键调控数量性状位点。
根据调节金黄色葡萄球菌丰度生长的显著数量性状位点的遗传效应，建立数量性状位点之间的线性方程组。由方程的系数表示数量性状位点之间的调控关系。图9为49个数量性状位点的遗传效应的基因调控网络，其中3个数量性状位点(1、3、4)：SNP 15210、SNP54383、SNP 58404为遗传效应网络中的关键调控数量性状位点，对金黄色葡萄球菌丰度生长的遗传结构中起着关键作用。网络关系中数量性状位点的编号表示在表2中。
表2调控6μg/ml万古霉素环境下金黄色葡萄球菌样本生长的显著数量性状位点信息
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。
Claims (6)
1.一种筛选微生物细菌关键调控QTL的方法，其特征在于，包括：
根据具有相同生活环境的多个微生物细菌样本个体的数量性状的表型数据拟合Logistic方程，得到Logistic方程的估计参数；
将每个微生物细菌样本个体的基因型和表型数据作为分析数据，将Logistic方程的估计参数应用到功能作图的模型框架中，定位调控微生物细菌样本个体数量性状生长的显著数量性状位点；
根据显著数量性状位点，计算数量性状的遗传效应值，并根据遗传效应值建立不同显著数量性状位点之间的线性相关关系，得到数量性状位点调控关系网络；
根据数量性状位点调控关系网络识别解释微生物细菌样本个体数量性状生长过程的关键调控数量性状位点。
2.根据权利要求1所述的一种筛选微生物细菌关键调控QTL的方法，其特征在于，根据具有相同生活环境的多个微生物细菌样本个体的数量性状的表型数据拟合Logistic方程，得到Logistic方程的估计参数，具体包括：
采用最小二乘法和具有相同生活环境的微生物细菌样本个体的数量性状的表型数据拟合Logistic方程，并采用BFGS拟牛顿法找到使得表型值与拟合值残差平方和最小的估计参数；
将表型值与拟合值残差平方和最小的估计参数确定为Logistic方程的估计参数。
3.根据权利要求1所述的一种筛选微生物细菌关键调控QTL的方法，其特征在于，所述Logistic方程为：
y＝A/(l+k*e-λt)；
式中，A表示最大生长量，k表示生长速率，λ表示生长的延迟时间，t表示时间，y表示数量性状的拟合值。
4.根据权利要求1所述的一种筛选微生物细菌关键调控QTL的方法，其特征在于，每个微生物细菌样本个体的基因型的确定过程为：
对每个微生物细菌样本个体的全基因组进行单核苷酸多态性分型，得到每个微生物细菌样本个体的基因型。
5.根据权利要求1所述的一种筛选微生物细菌关键调控QTL的方法，其特征在于，将每个微生物细菌样本个体的基因型和表型数据作为分析数据，将Logistic方程的估计参数应用到功能作图的模型框架中，定位调控微生物细菌样本个体数量性状生长的显著数量性状位点，具体包括：
建立原假设和备择假设；原假设为微生物细菌样本个体的数量性状的生长在不同基因型之间没有差异；备择假设为微生物细菌样本个体的数量性状的生长在不同基因型之间存在差异；
建立原假设的似然函数与备择假设的似然函数；其中，似然函数包括用来构建均值向量的Logistic方程的估计参数以及用来构建构造协方差矩阵的结构参数；
根据原假设的似然函数与备择假设的似然函数建立似然比统计量函数；
根据微生物细菌样本个体的基因型和表型数据、以及似然比统计量函数定位调控微生物细菌样本个体数量性状生长的显著数量性状位点。
6.根据权利要求1所述的一种筛选微生物细菌关键调控QTL的方法，其特征在于，数量性状的遗传效应值的计算公式为：
其中，a(t)观测时间为t时数量性状的遗传效应值，J为个显著数量性状位点的基因型数量，nj为基因型为j的微生物细菌样本个体的数量，μj(t)为观测时间为t时基因型为j的微生物细菌样本个体的均值。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310749301.3A CN116825191A (zh) | 2023-06-25 | 2023-06-25 | 一种筛选微生物细菌关键调控qtl的方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310749301.3A CN116825191A (zh) | 2023-06-25 | 2023-06-25 | 一种筛选微生物细菌关键调控qtl的方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN116825191A true CN116825191A (zh) | 2023-09-29 |
Family
ID=88142314
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310749301.3A Pending CN116825191A (zh) | 2023-06-25 | 2023-06-25 | 一种筛选微生物细菌关键调控qtl的方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116825191A (zh) |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103810403A (zh) * | 2014-03-05 | 2014-05-21 | 北京林业大学 | 一种基于基因互作的植物生长预测调控方法及系统 |
| CN104361263A (zh) * | 2014-10-10 | 2015-02-18 | 北京林业大学 | 基于基因与蛋白调控网络的林木育种方法和系统 |
| CN109182538A (zh) * | 2018-09-29 | 2019-01-11 | 南京农业大学 | 奶牛乳腺炎关键SNPs位点rs88640083及2b-RAD基因分型和分析方法 |
| CN113345520A (zh) * | 2021-06-07 | 2021-09-03 | 北京林业大学 | 一种基于Richards方程的树木数量性状的QTL定位框架的方法 |
-
2023
- 2023-06-25 CN CN202310749301.3A patent/CN116825191A/zh active Pending
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103810403A (zh) * | 2014-03-05 | 2014-05-21 | 北京林业大学 | 一种基于基因互作的植物生长预测调控方法及系统 |
| CN104361263A (zh) * | 2014-10-10 | 2015-02-18 | 北京林业大学 | 基于基因与蛋白调控网络的林木育种方法和系统 |
| CN109182538A (zh) * | 2018-09-29 | 2019-01-11 | 南京农业大学 | 奶牛乳腺炎关键SNPs位点rs88640083及2b-RAD基因分型和分析方法 |
| CN113345520A (zh) * | 2021-06-07 | 2021-09-03 | 北京林业大学 | 一种基于Richards方程的树木数量性状的QTL定位框架的方法 |
Non-Patent Citations (1)
| Title |
|---|
| DENGCHENG YANG ET AL.: "Inferring multilayer interactome networks shaping phenotypic plasticity and evolution", 《NATURE COMMUNICATIONS》, pages 13 - 14 * |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Liu et al. | Tetrad analysis in plants and fungi finds large differences in gene conversion rates but no GC bias | |
| Dong et al. | Accurate identification of single-nucleotide variants in whole-genome-amplified single cells | |
| Herrera et al. | Jack of all nectars, master of most: DNA methylation and the epigenetic basis of niche width in a flower‐living yeast | |
| Gerstein et al. | Genomic convergence toward diploidy in Saccharomyces cerevisiae | |
| Ness et al. | Extensive de novo mutation rate variation between individuals and across the genome of Chlamydomonas reinhardtii | |
| Chan et al. | The complex genetic architecture of the metabolome | |
| Bosshard et al. | Accumulation of deleterious mutations during bacterial range expansions | |
| Rockman et al. | Genetics of global gene expression | |
| Hottes et al. | Bacterial adaptation through loss of function | |
| Jeon et al. | Genome-wide profiling of DNA methylation provides insights into epigenetic regulation of fungal development in a plant pathogenic fungus, Magnaporthe oryzae | |
| Charron et al. | Spontaneous whole-genome duplication restores fertility in interspecific hybrids | |
| Chen et al. | Antagonistic pleiotropy conceals molecular adaptations in changing environments | |
| Saxer et al. | Whole genome sequencing of mutation accumulation lines reveals a low mutation rate in the social amoeba Dictyostelium discoideum | |
| Kumar et al. | Characterization of genetic diversity and population structure in wheat using array based SNP markers | |
| Volkers et al. | Gene-environment and protein-degradation signatures characterize genomic and phenotypic diversity in wild Caenorhabditis elegans populations | |
| Liu et al. | Genome-wide association study of 29 morphological traits in Aegilops tauschii | |
| Roncallo et al. | Genetic diversity and linkage disequilibrium using SNP (KASP) and AFLP markers in a worldwide durum wheat (Triticum turgidum L. var durum) collection | |
| Carlson et al. | MIPSTR: a method for multiplex genotyping of germline and somatic STR variation across many individuals | |
| Lyon et al. | Population structure of Peronospora effusa in the southwestern United States | |
| Diop et al. | A pseudomolecule‐scale genome assembly of the liverwort Marchantia polymorpha | |
| Jiang et al. | Impacts of mutation effects and population size on mutation rate in asexual populations: a simulation study | |
| Jeon et al. | Experimental evolution reveals genome-wide spectrum and dynamics of mutations in the rice blast fungus, Magnaporthe oryzae | |
| Amaradasa et al. | Effects of sublethal fungicides on mutation rates and genomic variation in fungal plant pathogen, Sclerotinia sclerotiorum | |
| US20020119451A1 (en) | System and method for predicting chromosomal regions that control phenotypic traits | |
| Wei et al. | Environment-dependent pleiotropic effects of mutations on the maximum growth rate r and carrying capacity K of population growth |
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 |