CN116312798B - 一种宏基因组测序数据物种验证的方法及应用 - Google Patents
一种宏基因组测序数据物种验证的方法及应用 Download PDFInfo
- Publication number
- CN116312798B CN116312798B CN202310153201.4A CN202310153201A CN116312798B CN 116312798 B CN116312798 B CN 116312798B CN 202310153201 A CN202310153201 A CN 202310153201A CN 116312798 B CN116312798 B CN 116312798B
- Authority
- CN
- China
- Prior art keywords
- reference genome
- species
- candidate reference
- calculating
- sliding window
- 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
- 238000012163 sequencing technique Methods 0.000 title claims abstract description 53
- 238000000034 method Methods 0.000 title claims abstract description 44
- 238000012795 verification Methods 0.000 title claims abstract description 42
- 238000007476 Maximum Likelihood Methods 0.000 claims abstract description 23
- 238000012216 screening Methods 0.000 claims description 14
- 238000004364 calculation method Methods 0.000 claims description 8
- 238000009825 accumulation Methods 0.000 claims description 7
- 238000004590 computer program Methods 0.000 claims description 6
- 238000010276 construction Methods 0.000 claims description 3
- 238000010200 validation analysis Methods 0.000 claims description 3
- 238000004458 analytical method Methods 0.000 abstract description 6
- 241000894007 species Species 0.000 description 112
- 241000589517 Pseudomonas aeruginosa Species 0.000 description 23
- 238000001514 detection method Methods 0.000 description 7
- 241000711573 Coronaviridae Species 0.000 description 5
- 244000052769 pathogen Species 0.000 description 5
- 108090000623 proteins and genes Proteins 0.000 description 5
- 238000011529 RT qPCR Methods 0.000 description 4
- 238000006243 chemical reaction Methods 0.000 description 4
- 238000010606 normalization Methods 0.000 description 4
- 230000001717 pathogenic effect Effects 0.000 description 3
- 241000589516 Pseudomonas Species 0.000 description 2
- 229940072293 axid Drugs 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000011156 evaluation Methods 0.000 description 2
- SGXXNSQHWDMGGP-IZZDOVSWSA-N nizatidine Chemical compound [O-][N+](=O)\C=C(/NC)NCCSCC1=CSC(CN(C)C)=N1 SGXXNSQHWDMGGP-IZZDOVSWSA-N 0.000 description 2
- 101150054338 ref gene Proteins 0.000 description 2
- 241000894006 Bacteria Species 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000003766 bioinformatics method Methods 0.000 description 1
- 238000012512 characterization method Methods 0.000 description 1
- 239000003153 chemical reaction reagent Substances 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 239000000284 extract Substances 0.000 description 1
- 239000012530 fluid Substances 0.000 description 1
- PCHJSUWPFVWCPO-UHFFFAOYSA-N gold Chemical compound [Au] PCHJSUWPFVWCPO-UHFFFAOYSA-N 0.000 description 1
- 208000015181 infectious disease Diseases 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 230000000813 microbial effect Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000007481 next generation sequencing Methods 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 238000011002 quantification Methods 0.000 description 1
- 238000003753 real-time PCR Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000002864 sequence alignment Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000009897 systematic effect Effects 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
- G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids
- G16B30/10—Sequence alignment; Homology search
-
- 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
-
- 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)
- Physics & Mathematics (AREA)
- Health & Medical Sciences (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Engineering & Computer Science (AREA)
- Biophysics (AREA)
- General Health & Medical Sciences (AREA)
- Analytical Chemistry (AREA)
- Chemical & Material Sciences (AREA)
- Bioinformatics & Computational Biology (AREA)
- Biotechnology (AREA)
- Evolutionary Biology (AREA)
- Proteomics, Peptides & Aminoacids (AREA)
- Medical Informatics (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Theoretical Computer Science (AREA)
- Genetics & Genomics (AREA)
- Molecular Biology (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本申请属于生信分析技术领域，具体涉及一种宏基因组测序数据物种验证的方法及应用。所述方法通过精确序列比对计算待验证物种的多个参考基因组比对指标，基于最大可能性原则筛选候选参考基因组及比对指标的极大似然估计，采用固定步长滑动窗口计算基因组不同区域间序列分布的基尼系数，生成可变长度滑动窗口计算基因组全长范围内序列分布的可信度和不可信度，进而对宏基因组测序数据中待验证物种进行综合评估。此外，本申请基于待验证物种各亚型比对指标、基尼系数、唯一比对序列占比，可进一步构建物种亚型概率模型并对各亚型进行统计推断。
Description
技术领域
本申请属于生物信息学技术领域，具体涉及一种宏基因组测序数据物种验证的方法及应用。
技术背景
宏基因组测序(metagenomics next generation sequencing,mNGS)是一种新型病原检测工具，在检测不明原因，难培养及共感染病原体中具有显著优势。
由于mNGS的参考序列数据库中物种类别数以万计，mNGS检测结果中可出现大量可疑或候选物种。病原培养和或定量PCR是验证候选物种的金标准，但通量低，个别病原可能会难培养，极大限制了其应用。
目前，在生物信息学层面上对mNGS检测结果中候选物种进行验证的常用方法是比对到代表性参考基因组，计算覆盖度和平均深度，绘制深度覆盖图，最后人工判读。由于序列比对过程中所选物种的代表性参考基因组通常为NCBI官网上标识为representative或reference的单个基因组，而对于真实样本mNGS数据中待验证的物种，该“代表性”的单个基因组不一定具有“代表性”。其次，覆盖度和平均深度是衡量物种可信度的宏观指标，缺乏对候选物种序列分布的精细表征和量化，如：mNGS数据在物种参考基因组各区域范围内序列分布的集中和离散特征、序列分布不均、可信度和不可信度的精确度量。最后，人工判读环节会存在一定的主观性，而且效率较低。
因此，对mNGS的检测结果中的候选物种进行系统全面的生物信息学验证，可有效减少假阳性，提高mNGS检测的特异度和准确性。另外，对宏基因组测序数据物种验证的方法还可用于物种亚型鉴定，扩展了其应用价值。
鉴于此，提出本申请。
发明内容
为解决上述技术问题，本申请通过生物信息学分析研究，建立一套宏基因组测序数据物种验证的方法，该方法可有效减少mNGS检测结果假阳性，提高特异度和准确性，提升病原检测效能，并且可扩展应用于物种亚型的鉴定。
具体的，本申请提出如下技术方案：
本申请首先提供一种宏基因组测序数据物种验证方法，所述方法包括如下步骤：
1)比对输入序列：将宏基因组测序数据分别比对到待验证物种的多个参考基因组；
2)计算比对指标，比对指标包括：序列总数、平均深度、覆盖区域和覆盖度；
3)构建打分模型：基于比对指标的统计量，构建概率打分模型；
4)筛选候选参考基因组：基于最大可能性原则和概率打分模型，筛选候选参考基因组；
5)计算候选参考基因组的基尼系数：基于候选参考基因组的比对结果，计算候选参考基因组不同区域间序列分布的基尼系数；
6)计算候选参考基因组比对结果的可信度：基于候选参考基因组的比对结果，计算候选参考基因组比对结果的可信度和不可信度；
优选的，所述方法还包括如下步骤：
7)生成验证报告：基于候选参考基因组、比对指标、基尼系数、可信度及不可信度，对待定物种进行综合判定，生成验证报告。
进一步的，所述1)中：
所述宏基因组测序数据为去人源的宏基因组测序数据；
进一步优选的，所述比对采用bwa-mem-2精确比对到待验证物种多个参考基因组上；
更进一步优选的，所述多个参考基因组具有如下特点：
a、refgene、genebank及fda-argos官网上下载的高质量组装基因组；
b、参考基因组为近期上传的前10个或20个参考基因组；
c、包括标识为representative和reference的参考基因组；
d、包括物种常见的各个亚型。
进一步的，所述2)具体为：基于比对结果计算比对到各个参考基因组的序列总数、平均深度、覆盖区域和覆盖度；
进一步优选的，所述比对指标分别定义如下：
e、序列总数：精确比对到各个参考基因组上的序列数；
f、平均深度：各个参考基因组上平均每个碱基被覆盖的次数；
g、覆盖区域：各个参考基因组上有序列覆盖的区域长度；
h、覆盖度：各个参考基因组上有序列覆盖的区域长度占比。
进一步的，所述3)具体为：对各参考基因组的序列总数、平均深度、覆盖区域及覆盖度进行排序、加权、求和、取倒和归一化，计算统计量，构建概率打分模型；
进一步优选的，所述概率打分模型为：
其中，G：参考基因组，n：参考基因组个数，i＝1,2,3,...,n，j＝1,2,3,...,n，Rank：各参考基因组各指标排序后的名次；
更进一步优选的，所述概率打分模型的构建步骤具体包括以下：
a、提取比对到各个参考基因组的序列总数、平均深度、覆盖区域和覆盖度信息；
b、分别对序列总数、平均深度、覆盖区域和覆盖度按由大到小进行排序，得到各个参考基因上述四个比对指标上的排名，分别记为：Rank_reads(Gi)、Rank_depth(Gi)、Rank_range(Gi)和Rank_coverage(Gi)；
c、依次经过加权、求和、取倒和归一化转换，构建概率打分模型。
进一步的，所述4)具体为：基于最大可能性原则和概率打分模型，筛选候选参考基因组，提取比对到该参考基因组上的序列总数、平均深度、覆盖区域及覆盖度，作为该待验证物种比对指标的极大似然估计；
进一步优选的，所述极大似然估计的具体步骤包括：
a、确定物种的候选参考基因组Gx，其中x＝argmax(P(Gi))，即P(Gi)取最大值时参考基因组的编号；
b、待验证物种以Gx为候选参考基因组时，将输入序列精确比对到Gx得到的序列总数、平均深度、覆盖区域和覆盖度作为待验证物种比对指标的极大似然估计。
进一步的，所述5)中，所述计算是采用固定步长滑动窗口计算候选参考基因组不同区域间序列分布的基尼系数；
进一步优选的，所述步骤5)具体为：基于候选参考基因组的比对结果，创建固定步长的滑动窗口，分别计算候选参考基因组各个滑动窗口的序列分布、滑动窗口累积百分比、滑动窗口累积序列百分比，绘制洛伦兹曲线，计算基尼系数，作为基因组不同区域序列分布不均匀程度的统计指标；
更进一步优选的，所述基尼系数的计算具体步骤包括：
a、基于候选参考基因组比对结果，计算每个碱基位置的深度；
b、按固定步长，默认设置为2kb，在候选基因组范围内创建滑动窗口；
c、统计每个滑动窗口区域内碱基深度总和，除以序列长度，得到每个滑动窗口的序列数；
d、按滑动窗口序列数排序，计算滑动窗口累积百分比，滑动窗口累积序列百分比；
e、绘制洛伦兹曲线，计算A/(A+B)的值，即基尼系数，作为衡量序列分配不均的指标。
进一步的，所述6)中，所述计算是采用可变长度滑动窗口计算候选参考基因组比对结果的可信度和不可信度；
进一步优选的，所述6)具体为：基于候选参考基因组的比对结果，生成(平均深度为0.95-1.05)的可变长度滑动窗口，分别计算各滑动窗口的覆盖度和可信度，以及不可信判定，统计不可信滑动窗口的长度占比，计算候选参考基因组比对结果的可信度和不可信度；
更进一步优选的，所述可信度和不可信度的计算具体包括以下步骤：
a、基于候选参考基因组比对结果，计算每个碱基位置的深度；
b、滑动窗口初始步长设置为序列长度的20倍，从左至右开始滑动；
c、计算滑动窗口的平均深度，覆盖度和可信度，可信度的计算方法定义为：
d、如果平均深度等于0，该滑动窗口标记为NA；
e、如果平均深度小于0.95，该滑动窗口标记为可信区域；
f、如果平均深度大于等于0.95并且小于1.05：在此前提下，如果覆盖度小于0.05，该滑动窗口标记为不可信区域；如果覆盖度大于等于0.05，该滑动窗口标记为可信区域；
g、如果平均深度大于等于1.05并且滑动窗口未到达参考基因组末端：按初始步长延伸滑动窗口，重复c到g的步骤；
h、如果平均深度大于等于1.05并且滑动窗口到达参考基因组末端：在此前提下，如果覆盖度小于0.05，该滑动窗口标记为不可信区域；如果覆盖度大于等于0.05，该滑动窗口标记为可信区域；
i、统计可信区域和不可信区域的长度，计算候选参考基因组从左至右方向的不可信度：
j、统计平均深度大于0的滑动窗口长度和可信度，计算候选参考基因组从左到右方向的可信度：
其中，length(i)：第i个滑动窗口的长度；confidence_left(i)：第i个滑动窗口从左至右可信度；n：平均深度大于0的可变长度滑动窗口个数；
k、参照上述步骤，从右至左滑动窗口，计算候选参考基因组从右至左方向的不可信度：
l、进一步计算候选参考基因组从右至左的可信度：
其中，length(i)：第i个滑动窗口的长度；confidence_right(i)：第i个滑动窗口从右至左可信度。
m、取候选参考基因组两个方向不可信度中极大值作为最终候选参考基因组不可信度；
n、取候选参考基因组两个方向可信度中极大值作为最终候选参考基因组可信度。
进一步的，所述7)中，基于上述统计指标对待验证物种进行综合判定，生成验证报告的步骤包括：
a、首先假定候选参考基因组可信；
b、如果候选参考基因组的平均深度大于1，覆盖度小于0.95，并且覆盖度/平均深度小于0.05，则判定该物种不可信；
c、如果候选参考基因组的序列分布基尼系数大于0.95，则判定该物种不可信；
d、如果候选参考基因组不可信度大于0.95，则判定该物种不可信。
e、如果候选参考基因组经上述步骤不能判定该物种为不可信：则给出该待验证物种的可信度；
f、生成候选基因组序列覆盖图，固定步长滑动窗口洛伦兹曲线和可变长度滑动窗口覆盖度图。
本申请还提供一种宏基因组测序数据物种亚型鉴定的方法，包括上述任一所述方法，并进一步包括如下步骤：
步骤8)物种亚型鉴定：基于物种各亚型比对指标、基尼系数、唯一比对序列，构建物种亚型概率模型对各亚型进行推断；
优选的，所述步骤8)具体为如下：
a、首先，按上述步骤1)-7)方法对待分型物种进行物种验证，如果判定为不可信，则结束后续物种亚型鉴定步骤；
b、如果不能判定待分型物种为不可信，则继续后续物种亚型鉴定步骤；
c、分别计算待分型物种各亚型参考基因组的比对指标，包括：比对序列数、平均深度、覆盖区域和覆盖度；
d、基于最大可能性原则和打分模型筛选各亚型候选参考基因组；
e、分别计算各亚型候选参考基因组的基尼系数、不可信度和可信度；
f、分别对各亚型候选参考基因组序列总数、平均深度、覆盖区域、覆盖度和可信度按由大到小进行排序，得到各个亚型候选参考基因组在上述5个统计指标上的排名，分别记为：Rank_reads(Si)、Rank_depth(Si)、Rank_range(Si)、Rank_coverage(Si)和Rank_confidence(Si)；
g、分别对各亚型候选参考基因组基尼系数按由小到大进行排序，得到各个亚型候选参考基因在基尼指数上的排名，记为：Rank_gini(Si)；
h、对上述5个统计指标依次经过加权、求和、取倒和归一化转换，构建各亚型概率打分模型，得到：
其中，
Rank(Si)＝0.8*Rank_reads(Si)+0.8*Rank_range(Si)+1.0*Rank_depth(Si)+1.0*Rank_coverage(Si)+1.0*Rank_confidence(Si)+1.0*Rank_gini(Si)，
S：物种亚型候选参考基因组，n：物种亚型候选参考基因组个数，
i＝1,2,3,...,n，j＝1,2,3,...,n。
i、对上述P(Si)采用唯一比对序列数占比进行矫正，得到：
其中，Unique_reads_abundance(Si)指各亚型候选参考基因组中，输入序列比对到Si参考基因组上的唯一比对序列占比；
j、取x＝argmax(Pc(Si))，Sx即为待验证物种亚型候选参考基因组的极大似然估计；
k、输入序列精确比对到Sx时，计算得到的序列总数、平均深度、覆盖区域、覆盖度、基尼系数、可信度和不可信度为该物种亚型比对及统计指标的极大似然估计值。
本申请首先提供一种宏基因组测序数据物种验证系统，所述系统包括如下模块：
1)比对输入序列模块：用于将宏基因组测序数据分别比对到待验证物种的多个参考基因组；
2)计算比对指标模块，所述比对指标包括：序列总数、平均深度、覆盖区域和覆盖度；
3)构建打分模型模块：用于基于比对指标的统计量，构建概率打分模型；
4)筛选候选参考基因组模块：用于基于最大可能性原则和概率打分模型，筛选候选参考基因组；
5)计算候选参考基因组的基尼系数模块：用于基于候选参考基因组的比对结果，计算候选参考基因组不同区域间序列分布的基尼系数；
6)计算候选参考基因组比对结果的可信度和不可信度模块：用于基于候选参考基因组的比对结果，计算候选参考基因组比对结果的可信度和不可信度；
优选的，所述方法还包括如下步骤：
7)生成验证报告模块：用于基于候选参考基因组、比对指标、基尼系数、可信度及不可信度，对待定物种进行综合判定，生成验证报告。
本申请还提供一种计算机存储介质，所述计算机存储介质存储有计算机程序，所述计算机程序包括程序指令，所述程序指令当被处理器执行时，实现上述任一项所述方法。
个申请还提供一种电子设备，所述电子设备包括存储器和处理器，所述存储器和处理器相连，所述存储器用于存储计算机程序，所述处理器用于调用所述计算机程序，实现上述任一项所述方法。
本申请有益技术效果：
1)本申请基于mNGS数据在待验证物种多个参考基因组上的精确比对信息对物种进行生信信息学验证，极大提高了参考基因组的完整性，全面性以及比对信息的准确性；
2)本申请基于最大可能性原则和极大似然估计对候选物种参考基因组及比对指标进行筛选，在假定mNGS数据中物种真实存在的前提下，候选物种参考基因组更具有代表性；
3)本申请在计算常见比对指标如覆盖度和平均深度的基础上，分别采用固定步长滑动窗口计算基尼系数，以及采用可变步长滑动窗口计算可信度和不可信度，深入挖掘和提取了mNGS数据在候选参考基因组上的序列分布特征及度量指标，为物种验证提供更多评估有效维度和证据；
4)本申请所述方法可扩展应用于物种亚型的鉴定，基于各亚型候选参考基因组的精确序列比对信息、序列分布特征和概率推断模型，极大提高了mNGS物种亚型鉴定的准确性。
附图说明
图1一种宏基因组测序数据物种验证的方法及应用示意图；
图2本申请方法扩展应用于物种亚型鉴定的示意图；
图3实施例1中铜绿假单胞菌真阳性测序数据覆盖度图；
图4实施例1中铜绿假单胞菌真阳性测序数据洛伦兹曲线和基尼系数；
图5实施例2中铜绿假单胞菌假阳性测序数据覆盖度图；
图6实施例2中铜绿假单胞菌假阳性测序数据洛伦兹曲线和基尼系数。
具体实施方式
下面将结合实施例对本申请的实施方案进行详细描述，但是本领域技术人员将会理解，下列实施例仅用于说明本申请，而不应视为限制本申请的范围。实施例中未注明具体条件者，按照常规条件或制造商建议的条件进行。所用试剂或仪器未注明生产厂商者，均为可以通过市场购买获得的常规产品。
部分术语定义，除非在下文中另有定义，本申请具体实施方式中所用的所有技术术语和科学术语的含义意图与本领域技术人员通常所理解的相同。虽然相信以下术语对于本领域技术人员很好理解，但仍然阐述以下定义以更好地解释本申请。
本申请中的术语“大约”表示本领域技术人员能够理解的仍可保证论及特征的技术效果的准确度区间。该术语通常表示偏离指示数值的±10％，优选±5％。
如本申请中所使用，术语“包括”、“包含”、“具有”、“含有”或“涉及”为包含性的(inclusive)或开放式的，且不排除其它未列举的元素或方法步骤。术语“由…组成”被认为是术语“包含”的优选实施方案。如果在下文中某一组被定义为包含至少一定数目的实施方案，这也应被理解为揭示了一个优选地仅由这些实施方案组成的组。
此外，说明书和权利要求书中的术语第一、第二、第三、(a)、(b)、(c)以及诸如此类，是用于区分相似的元素，不是描述顺序或时间次序必须的。应理解，如此应用的术语在适当的环境下可互换，并且本申请描述的实施方案能以不同于本申请描述或举例说明的其它顺序实施。
下面结合具体实施例来阐述本申请。
实验例：本申请物种验证和亚型鉴定方法体系建立
本申请通过生信分析优化，初步建立一套适用于宏基因组测序数据物种验证的方法，该方法具体包括了如下步骤：
步骤1)比对输入序列：将(优选去人源)宏基因组测序数据分别精确比对到待验证物种的多个参考基因组上；
其中，所述宏基因组测序数据为去人源的宏基因组测序数据；
所述比对采用bwa-mem-2精确比对到待验证物种多个参考基因组上；
优选的，所述多个参考基因组的集合定义为：
a、refgene，genebank及fda-argos官网上下载的高质量组装基因组；
b、参考基因组的上传时间新，筛选参数设置为最近上传的前10个或20个参考基因；
c、包括标识为representative和reference的参考基因组；
d、包括物种常见的各个亚型。
步骤2)计算比对指标：基于比对结果计算比对到各个参考基因组的序列总数、平均深度、覆盖区域和覆盖度；
其中，基于比对结果计算的比对指标定义如下：
a、序列总数：精确比对到各个参考基因组上的序列数；
b、平均深度：各个参考基因组上平均每个碱基被覆盖的次数；
c、覆盖区域：各个参考基因组上有序列覆盖的区域长度；
d、覆盖度：各个参考基因组上有序列覆盖的区域长度占比。
步骤3)构建打分模型：对各参考基因组的序列总数、平均深度、覆盖区域及覆盖度进行排序、加权、求和、取倒和归一化，计算统计量，构建概率打分模型；
概率打分模型的构建包括以下步骤：
a、提取比对到各个参考基因组的序列总数、平均深度、覆盖区域和覆盖度信息；
b、分别对序列总数、平均深度、覆盖区域和覆盖度按由大到小进行排序，得到各个参考基因上述四个比对指标上的排名，分别记为：Rank_reads(Gi)、Rank_depth(Gi)、Rank_range(Gi)和Rank_coverage(Gi)；
c、依次经过加权、求和、取倒和归一化转换，构建概率打分模型，得到：
其中，G：参考基因组，n：参考基因组个数，i＝1,2,3,...,n，j＝1,2,3,...,n，Rank：参考基因组各指标排序后的名次。
步骤4)筛选候选参考基因组：基于最大可能性原则和概率打分模型，筛选候选参考基因组，提取比对到该参考基因组上的序列总数、平均深度、覆盖区域及覆盖度，作为该待验证物种比对指标的极大似然估计；
具体的：
a、确定物种的候选参考基因组Gx，其中x＝argmax(P(Gi))，即P(Gi)取最大值时参考基因组的编号；
b、待验证物种以Gx为候选参考基因组时，将输入序列精确比对到Gx得到的序列总数、平均深度、覆盖区域和覆盖度即为待验证物种比对指标的极大似然估计值。
步骤5)计算基尼系数：基于候选参考基因组的比对结果，创建固定步长的滑动窗口，分别计算候选参考基因组各个滑动窗口的序列分布，滑动窗口累积百分比，滑动窗口累积序列百分比，绘制洛伦兹曲线，计算基尼系数，作为基因组不同区域序列分布不均匀程度的统计指标；
其中，基尼系数的计算步骤包括：
a、基于候选参考基因组比对结果，计算每个碱基位置的深度；
b、按固定步长，默认设置为2kb，在候选基因组范围内创建滑动窗口；
c、统计每个滑动窗口区域内碱基深度总和，除以序列长度，得到每个滑动窗口的序列数；
d、按滑动窗口序列数排序，计算滑动窗口累积百分比，滑动窗口累积序列百分比；
e、绘制洛伦兹曲线，计算A/(A+B)的值，即基尼系数，作为衡量序列分配不均的指标。
步骤6)计算可信度和不可信度：基于候选参考基因组的比对结果，生成长度可变的平均深度0.95-1.05的滑动窗口，分别计算各滑动窗口的覆盖度和可信度，以及不可信判定，统计不可信滑动窗口的长度占比，计算候选参考基因组比对结果的可信度和不可信度；
其中，不可信度和可信度的计算包括以下步骤：
a、基于候选参考基因组比对结果，计算每个碱基位置的深度；
b、滑动窗口初始步长设置为序列长度的20倍，以75bp序列长度为例，初始步长为1.5kb，从左至右开始滑动；
c、计算滑动窗口的平均深度，覆盖度和可信度，其中，平均深度和覆盖度的计算方法参照候选参考基因组的计算方法，可信度的计算方法定义为：
d、如果平均深度等于0，该滑动窗口标记为NA；
e、如果平均深度小于0.95，该滑动窗口标记为可信区域；
f、如果平均深度大于等于0.95并且小于1.05：在此前提下，如果覆盖度小于0.05，该滑动窗口标记为不可信区域；如果覆盖度大于等于0.05，该滑动窗口标记为可信区域；
g、如果平均深度大于等于1.05并且滑动窗口未到达参考基因组末端：按初始步长延伸滑动窗口，重复c到g的步骤；
h、如果平均深度大于等于1.05并且滑动窗口到达参考基因组末端：在此前提下，如果覆盖度小于0.05，该滑动窗口标记为不可信区域；如果覆盖度大于等于0.05，该滑动窗口标记为可信区域；
i、统计可信区域和不可信区域的长度，计算候选参考基因组从左至右方向的不可信度：
j、统计平均深度大于0的滑动窗口长度和可信度，计算候选参考基因组从左到右方向的可信度：
其中，length(i)：第i个滑动窗口的长度；confidence_left(i)：第i个滑动窗口从左至右可信度；n：平均深度大于0的可变长度滑动窗口个数。
k、参照上述步骤，从右至左滑动窗口，计算候选参考基因组从右至左方向的不可信度：
l、进一步计算候选参考基因组从右至左的可信度：
其中，length(i)：第i个滑动窗口的长度；confidence_right(i)：第i个滑动窗口从右至左可信度。
m、取候选参考基因组两个方向不可信度中极大值作为最终候选参考基因组不可信度；
n、取候选参考基因组两个方向可信度中极大值作为最终候选参考基因组可信度。
步骤7)生成验证报告：基于候选参考基因组、比对指标、基尼系数、可信度及不可信度，对待定物种进行综合判定，生成验证报告；
其中，基于上述统计指标对待验证物种进行综合判定，生成验证报告的步骤包括：
a、首先假定候选参考基因组可信；
b、如果候选参考基因组的平均深度大于1，覆盖度小于0.95，并且覆盖度/平均深度小于0.05，则判定该物种不可信；
c、如果候选参考基因组的序列分布基尼系数大于0.95，则判定该物种不可信；
d、如果候选参考基因组不可信度大于0.95，则判定该物种不可信。
e、如果候选参考基因组经上述步骤不能判定该物种为不可信：则给出该待验证物种的可信度；
f、生成候选基因组序列覆盖图，固定步长滑动窗口洛伦兹曲线和可变长度滑动窗口覆盖度图。
优选的，本申请中所述物种验证方法可扩展应用于物种亚型鉴定。
步骤8)物种亚型鉴定：基于物种各亚型比对指标、基尼系数、可信度和不可信度，构建物种亚型概率模型对各亚型进行推断。
其中，物种亚型鉴定的步骤包括：
a、首先，按上述方法对待分型物种进行物种验证，如果判定为不可信，则结束后续物种亚型鉴定步骤；
b、如果不能判定待分型物种为不可信，则继续后续物种亚型鉴定步骤；
c、分别计算待分型物种各亚型参考基因组的比对指标，包括：比对序列数、平均深度、覆盖区域和覆盖度；
d、基于最大可能性原则和打分模型筛选各亚型候选参考基因组；
e、分别计算各亚型候选参考基因组的基尼系数、不可信度和可信度；
f、分别对各亚型候选参考基因组序列总数、平均深度、覆盖区域、覆盖度和可信度按由大到小进行排序，得到各个亚型候选参考基因组在上述5个统计指标上的排名，分别记为：Rank_reads(Si)、Rank_depth(Si)、Rank_range(Si)、Rank_coverage(Si)和Rank_confidence(Si)；
g、分别对各亚型候选参考基因组基尼系数按由小到大进行排序，得到各个亚型候选参考基因在基尼指数上的排名，记为：Rank_gini(Si)；
h、对上述5个统计指标依次经过加权、求和、取倒和归一化转换，构建各亚型概率打分模型，得到：
其中，
Rank(Si)＝0.8*Rank_reads(Si)+0.8*Rank_range(Si)+1.0*Rank_depth(Si)+1.0*Rank_coverage(Si)+1.0*Rank_confidence(Si)+1.0*Rank_gini(Si)，
S：物种亚型候选参考基因组，n：物种亚型候选参考基因组个数，
i＝1,2,3,...,n，j＝1,2,3,...,n。
ii、对上述P(Si)采用唯一比对序列数占比进行矫正，得到：
其中，Unique_reads_abundance(Si)指各亚型候选参考基因组中，输入序列比对到Si参考基因组上的唯一比对序列占比；
j、取x＝argmax(Pc(Si))，Sx即为待验证物种亚型候选参考基因组的极大似然估计；
k、输入序列精确比对到Sx时，计算得到的序列总数、平均深度、覆盖区域、覆盖度、基尼系数、可信度和不可信度为该物种亚型比对及统计指标的极大似然估计值。
实施例1：绿脓杆菌真阳性验证
1)铜绿假单胞菌的真阳性宏基因组测序数据介绍
铜绿假单胞菌的真阳性宏基因组测序数据SEQ_1样本来源为ZymoBIOMICS公司的微生物DNA标准品，货号为D6305，含10种不同细菌，包括铜绿假单胞菌，qPCR验证结果为阳性。
2)物种分类和丰度估计结果
常规宏基因组分析流程结果提示：铜绿假单胞菌(物种编号：287)的丰度估计结果为：47031条序列。
3)按照本申请方法步骤进行物种验证。
4)铜绿假单胞菌物种验证结果如下：
如物种验证结果所示，铜绿假单胞菌在候选参考基因组GCF_026890985.1_ASM2689098v1_genomic.fna上的比对结果打分最高，精确比对上的序列数为37589条，平均深度为0.87，候选参考基因组上的覆盖区域为3004587bp，基尼系数为0.25，不可信度为0，可信度为0.55，综合判定结果为该样本中铜绿假单胞菌的可信度为0.55，与qPCR验证结果一致。
如图3所示，真阳性铜绿杆菌测序数据的序列在铜绿杆菌候选参考基因组上均匀分布，每个固定步长滑动窗口的平均深度为10条序列。如图4所示，候选参考基因组全长范围内序列分布基尼系数为0.25，提示序列覆盖相对均匀。
由实施例1可知，本申请提供的宏基因组测序数据物种验证方法为评估待验证物种是否为真阳性提供了基于精确比对结果的多维评估指标，包括：a、候选参考基因组极大似然估计；b、基础比对指标：比对序列数、平均深度、覆盖区域和覆盖度；c、序列分布不均衡量指标：基尼系数；d、序列分布可信度和不可信度。最后基于概率模型和规则给出待验证物种的综合判定结果。
实施例2：绿脓杆菌假阳性验证
1)铜绿假单胞菌假阳性宏基因组测序数据介绍
铜绿假单胞菌假阳性宏基因组测序数据SEQ_2来自临床肺泡灌洗液样本，铜绿假单胞菌qPCR验证结果为阴性。
2)物种分类和丰度估计结果
常规宏基因组分析流程结果提示：铜绿假单胞菌(物种编号：287)的丰度估计结果为：47177条序列。
3)按照本申请方法步骤进行物种验证。
4)绿脓杆菌物种验证结果如下：
如物种验证结果所示，铜绿假单胞菌在候选参考基因组GCF_026890985.1_ASM2689098v1_genomic.fna上的比对结果打分最高，精确比对上的序列数为452667条，平均深度为5.11，候选参考基因组上的覆盖区域为549189bp，基尼系数为0.81，不可信度为0，可信度为0.02，综合判定结果为不可信，与qPCR验证结果一致。
如图5所示，假阳性铜绿杆菌测序数据的序列在铜绿杆菌候选参考基因组上分布不匀，大量固定步长滑动窗口内平均深度为0。如图4所示，候选参考基因组全长范围内序列分布基尼系数高达0.81，提示序列分布严重不均。
进一步分析该假阳性物种所在测序数据可知，与铜绿假单胞菌同属异种的物种莫氏假单胞菌的丰度最高，总序列数达到300万条，模拟测序数据提示：当模拟测序数据中莫氏假单胞菌丰度到达30万条时，铜绿假单胞菌会出现上述假阳性特征，包括：低覆盖度、高平均深度、高基尼系数以及低可信度。
由实施例2可知，本申请提供的宏基因组测序数据物种验证方法为评估待验证物种是否为假阳性提供了基于精确比对结果的多维评估指标，最后基于概率模型和规则给出综合判定结果。
实施例3：新型冠状病毒亚型XBB株的验证和鉴定
1)新型冠状病毒XBB株测序数据模拟
从网站https://ngdc.cncb.ac.cn/ncov/下载新型冠状病毒较常见的不同株序列。
本实施例中下载的参考序列如下表。
采用art_illumina模拟器模拟，以OQ202035.1(XBB.1)为参考序列，模拟生成序列丰度为1000，序列读长为75bp，测序模式为单端，测序平台为HS2500的测序数据xbb-1。
其它6种株(XBB.3，B，BQ.1，B.1.617.2，BF.7，BA.5)作为不种亚型株的候选参考基因组。
XBB.1用于模拟生成xbb-1测序数据，为避免参考序列和模拟序列在比对过程中出现的过拟合，未将其纳入参考序列集合，而采用同株的XBB.3。
2)按照本申请方法步骤进行XBB物种验证和亚型鉴定。
3)亚型分型结果如下：
物种分类和物种验证结果均为冠状病毒(taxid:694009)，候选参考基因组亚型为新型冠状病毒(taxid:2697049)。
物种亚型鉴定结果如下：
如表所示，XBB.3得分最高，覆盖区域、比对序列数、平均深度和唯一比对序列数均得分最高。其中，2条唯一比对序列，即在多个株比对信息中，有2条序列唯一比对到该株上。另外，由于BQ.1与XBB同属于Omicron系，与XBB的基因组相似度最高，位于亚型鉴定结果第2位，与期望结果一致。
由实施例3可知，本申请提供的物种亚型鉴定方法可基于精确比对和多维评估指标精确定位待分型物种的具体亚型，与物种验证模块联合使用可显著提高宏基因组测序数据分析结果的特异度和准确性。
最后应说明的是：以上各实施例仅用以说明本申请的技术方案，而非对其限制；尽管参照前述各实施例对本申请进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本申请各实施例技术方案的范围。
Claims (5)
1.一种宏基因组测序数据物种验证方法，其特征在于，包括如下步骤：
1)比对输入序列：将宏基因组测序数据分别比对到待验证物种的多个参考基因组；
2)计算比对指标，比对指标包括：序列总数、平均深度、覆盖区域和覆盖度；
3)构建打分模型：基于比对指标的统计量，构建概率打分模型；
4)筛选候选参考基因组：基于最大可能性原则和概率打分模型，筛选候选参考基因组；
5)计算候选参考基因组的基尼系数：基于候选参考基因组的比对结果，计算候选参考基因组不同区域间序列分布的基尼系数；
6)计算候选参考基因组比对结果的可信度：基于候选参考基因组的比对结果，计算候选参考基因组比对结果的可信度和不可信度；
所述1)中：所述宏基因组测序数据为去人源的宏基因组测序数据；所述比对采用bwa-mem-2精确比对到多个参考基因组上；
所述2)中，所述计算为：基于比对结果计算比对到各个参考基因组的序列总数、平均深度、覆盖区域和覆盖度；
所述3)中，所述构建为：对各参考基因组的序列总数、平均深度、覆盖区域及覆盖度进行排序、加权、求和、取倒和归一化，计算统计量，构建概率打分模型；
所述概率打分模型为：
其中，G：参考基因组，n：参考基因组个数，i＝1,2,3,...,n，j＝1,2,3,...,n，Rank：各参考基因组各指标排序后的名次；
所述4)中，所述筛选候选参考基因组为：
基于最大可能性原则和概率打分模型，筛选候选参考基因组，提取比对到该参考基因组上的序列总数、平均深度、覆盖区域及覆盖度，作为该待验证物种比对指标的极大似然估计；
所述5)中，所述计算是采用固定步长滑动窗口计算候选参考基因组不同区域间序列分布的基尼系数：基于候选参考基因组的比对结果，创建固定步长的滑动窗口，分别计算候选参考基因组各个滑动窗口的序列分布、滑动窗口累积百分比、滑动窗口累积序列百分比，绘制洛伦兹曲线，计算基尼系数，作为基因组不同区域序列分布不均匀程度的统计指标；
所述6)中，所述计算是采用可变长度滑动窗口计算候选参考基因组比对结果的可信度和不可信度：基于候选参考基因组的比对结果，生成可变长度滑动窗口，分别计算各滑动窗口的覆盖度和可信度，以及不可信判定，统计不可信滑动窗口的长度占比，计算候选参考基因组比对结果的可信度和不可信度。
2.根据权利要求1所述的方法，其特征在于，进一步包括如下步骤：
7)生成验证报告：基于候选参考基因组、比对指标、基尼系数、可信度及不可信度，对待定物种进行综合判定，生成验证报告。
3.根据权利要求1-2任一所述的方法，其特征在于，所述2)中，所述比对指标分别定义如下：
a、序列总数：精确比对到各个参考基因组上的序列数；
b、平均深度：各个参考基因组上平均每个碱基被覆盖的次数；
c、覆盖区域：各个参考基因组上有序列覆盖的区域长度；
d、覆盖度：各个参考基因组上有序列覆盖的区域长度占比。
4.一种计算机存储介质，其特征在于，所述计算机存储介质存储有计算机程序，所述计算机程序包括程序指令，所述程序指令当被处理器执行时，实现权利要求1-3中任一项所述方法。
5.一种电子设备，其特征在于，所述电子设备包括存储器和处理器，所述存储器和处理器相连，所述存储器用于存储计算机程序，所述处理器用于调用所述计算机程序，实现权利要求1-3中任一项所述方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310153201.4A CN116312798B (zh) | 2023-02-22 | 2023-02-22 | 一种宏基因组测序数据物种验证的方法及应用 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310153201.4A CN116312798B (zh) | 2023-02-22 | 2023-02-22 | 一种宏基因组测序数据物种验证的方法及应用 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116312798A CN116312798A (zh) | 2023-06-23 |
| CN116312798B true CN116312798B (zh) | 2023-11-10 |
Family
ID=86821555
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310153201.4A Active CN116312798B (zh) | 2023-02-22 | 2023-02-22 | 一种宏基因组测序数据物种验证的方法及应用 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116312798B (zh) |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112599198A (zh) * | 2020-12-29 | 2021-04-02 | 上海派森诺生物科技股份有限公司 | 一种用于宏基因组测序数据的微生物物种与功能组成分析方法 |
| CN113744807A (zh) * | 2021-11-03 | 2021-12-03 | 微岩医学科技(北京)有限公司 | 一种基于宏基因组学的病原微生物检测方法及装置 |
| CN114555835A (zh) * | 2019-07-23 | 2022-05-27 | 生物梅里埃公司 | 通过宏基因组分析检测和定量感兴趣的生物物种的方法 |
| WO2022222936A1 (en) * | 2021-04-20 | 2022-10-27 | Hangzhou Matridx Biotechnology Co., Ltd. | Methods, computer-readble media, and systems for filtering noises for dna sequencing data |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2019126824A1 (en) * | 2017-12-22 | 2019-06-27 | Trace Genomics, Inc. | Metagenomics for microbiomes |
-
2023
- 2023-02-22 CN CN202310153201.4A patent/CN116312798B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114555835A (zh) * | 2019-07-23 | 2022-05-27 | 生物梅里埃公司 | 通过宏基因组分析检测和定量感兴趣的生物物种的方法 |
| CN112599198A (zh) * | 2020-12-29 | 2021-04-02 | 上海派森诺生物科技股份有限公司 | 一种用于宏基因组测序数据的微生物物种与功能组成分析方法 |
| WO2022222936A1 (en) * | 2021-04-20 | 2022-10-27 | Hangzhou Matridx Biotechnology Co., Ltd. | Methods, computer-readble media, and systems for filtering noises for dna sequencing data |
| CN113744807A (zh) * | 2021-11-03 | 2021-12-03 | 微岩医学科技(北京)有限公司 | 一种基于宏基因组学的病原微生物检测方法及装置 |
Non-Patent Citations (2)
| Title |
|---|
| Metagenomics next-generation sequencing tests take the stage in the diagnosis of lower respiratory tract infections;Zhenli Diao.et.;《Journal of Advanced Research》;第38卷;第 201-212页 * |
| 基于序列组成的未知环境宏基因组快速序列分类系统;乔于洋;《中国优秀硕士学位论文全文数据库基础科学辑》(第6期);第A006-160页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116312798A (zh) | 2023-06-23 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Li et al. | Comparing the normalization methods for the differential analysis of Illumina high-throughput RNA-Seq data | |
| Poon et al. | Detecting signatures of selection from DNA sequences using Datamonkey | |
| CN103886168A (zh) | 基于层次分析法的多渠道分析方法及装置 | |
| AU2021212155B2 (en) | Method and apparatus for estimating the quantity of microorganisms within a taxonomic unit in a sample | |
| CN104866863B (zh) | 一种生物标志物筛选方法 | |
| EP3518974A1 (en) | Noninvasive prenatal screening using dynamic iterative depth optimization | |
| Sipos et al. | Robust computational analysis of rRNA hypervariable tag datasets | |
| CN110648720B (zh) | 宏基因组测序质控预测评估方法及模型 | |
| CN116312798B (zh) | 一种宏基因组测序数据物种验证的方法及应用 | |
| CN114388062A (zh) | 基于机器学习预测抗生素抗性表型的方法、设备及应用 | |
| WO2023124779A1 (zh) | 基于三代测序数据检测点突变的分析方法和装置 | |
| Lalande et al. | A new framework to accurately quantify soil bacterial community diversity from DGGE | |
| CN110970093B (zh) | 一种筛选引物设计模板的方法、装置及应用 | |
| CN113260710A (zh) | 用于通过多个定制掺合混合物验证微生物组序列处理和差异丰度分析的组合物、系统、设备和方法 | |
| CN116825192A (zh) | 一种ncRNA基因突变的解读方法、存储介质及终端 | |
| CN108595914A (zh) | 一种烟草线粒体rna编辑位点高精度预测方法 | |
| CN112883284B (zh) | 一种基于网络和数据分析的在线学习系统及测试题推荐方法 | |
| US20220122695A1 (en) | Methods and systems for providing sample information | |
| CN115910216B (zh) | 一种基于机器学习识别基因组序列分类错误的方法和系统 | |
| CN116030983B (zh) | 一种基于机器学习的子宫内膜癌组织学等级预测方法 | |
| CN116153411B (zh) | 多病原体探针库组合的设计方法及应用 | |
| CN116646010B (zh) | 人源性病毒检测方法及装置、设备、存储介质 | |
| Wang et al. | Evaluation of normalization methods for predicting quantitative phenotypes in metagenomic data analysis | |
| Sun et al. | Statistical calibration of qRT-PCR, microarray and RNA-Seq gene expression data with measurement error models | |
| CN109671467B (zh) | 一种病原体感染损伤机理解析方法及装置 |
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 |