WO2018086045A1 - 一种对特定群中的亚群进行定量分析的方法 - Google Patents

一种对特定群中的亚群进行定量分析的方法 Download PDF

Info

Publication number
WO2018086045A1
WO2018086045A1 PCT/CN2016/105372 CN2016105372W WO2018086045A1 WO 2018086045 A1 WO2018086045 A1 WO 2018086045A1 CN 2016105372 W CN2016105372 W CN 2016105372W WO 2018086045 A1 WO2018086045 A1 WO 2018086045A1
Authority
WO
WIPO (PCT)
Prior art keywords
matrix
frequency
snp
specific group
vector
Prior art date
Application number
PCT/CN2016/105372
Other languages
English (en)
French (fr)
Inventor
彭也
李俊桦
唐珊媚
张慧
Original Assignee
深圳华大基因研究院
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by 深圳华大基因研究院 filed Critical 深圳华大基因研究院
Priority to CN201680090780.0A priority Critical patent/CN109997193B/zh
Priority to PCT/CN2016/105372 priority patent/WO2018086045A1/zh
Publication of WO2018086045A1 publication Critical patent/WO2018086045A1/zh

Links

Images

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B25/00—ICT specially adapted for hybridisation; ICT specially adapted for gene or protein expression
    • 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
    • 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

Definitions

  • the population is selected from the group consisting of bacteria, viruses, tumor cells, fungi, chlamydia, and mycoplasma.
  • a non-diagnostic method for quantitative analysis of a subpopulation in a particular population comprising the steps of:
  • the subpopulation is selected from the group consisting of subspecies, heterogeneous cells, strains.
  • the present inventors also provide a quantitative method for realizing the ratio of strains contained in a reference matrix in the case where the number of sequencing layers of the strain to be studied (the number of times the genome is covered by the sequencing read sequence) is small;
  • the appropriate model makes the process implementation process less complex in space and time.
  • m is the total number of bases whose frequency is not zero (ie, the total number of bases in the table 3-binarized SNP matrix).

Landscapes

  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Evolutionary Biology (AREA)
  • Theoretical Computer Science (AREA)
  • Engineering & Computer Science (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Biotechnology (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • General Health & Medical Sciences (AREA)
  • Biophysics (AREA)
  • Genetics & Genomics (AREA)
  • Molecular Biology (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)

Abstract

一种对特定群中的亚群进行定量分析的方法。具体地，包括以下步骤：（1）提供对应于所述特定群的（a）参考基因组序列数据、（b）参考SNP矩阵和（c）宏基因组测序数据；（2）将所述宏基因组测序数据比对到对应于所述特定群的参考基因组数据上，以便获得比对结果；（3）根据所述参考SNP矩阵的位点信息，构建频率矩阵；（4）根据频率矩阵对参考SNP矩阵做二值化处理，得到二值化SNP矩阵；和（5）基于所述的频率矩阵、所述的二值化SNP矩阵、理论碱基频率向量f(x)和观测碱基频率向量y，通过有约束的线性模型，得出所述特定群中各亚群的相对丰度，从而获得对所述特定群中的亚群的定量分析结果。

Description

一种对特定群中的亚群进行定量分析的方法 技术领域
本发明属于生物信息学领域，更具体地涉及一种对特定群中的亚群进行定量分析的方法。
背景技术
宏基因组学(metagenomics)，其通过从微生态样品中提取全部微生物的DNA，构建宏基因组文库，进而利用基因组学的研究策略分析样品中全部微生物的遗传组成及其群落的潜在功能。该组学技术不依赖于单细菌的分离和纯培养技术，在很大程度上解决了大部分微生物因不能分离培养而难于研究的问题，同时也能反映研究生态环境中微生物组成的真实情况。而在关于人类健康的微生态研究中，基于宏基因组学测序数据对微生物进行定量(quantification)，是研究其群落构成、物种间相互作用，探索其与疾病发生发展的关系等相关规律的基础。随着科学研究的进展，越来越多的研究显示，对特定物种较低级别的分类单元，即菌株(strains)进行准确的注释愈发重要。如果简单地在一个较高的分类级别研究细菌与疾病的关系，则很可能把与疾病的发展呈正相关、无关甚至负相关的类别都加到一起，这无论在生物学还是统计学上都有明显谬误；而现有研究成果也亟待通过提高微生物定量精度来校正，或进行更深入的机制研究。
现有宏基因组物种鉴定、定量方法中，分辨率较高，即能注释到菌种或菌株水平的有基于全基因组序列比对和基于标记序列比对的方法，代表的流程分别为WG-FAST和Constrains。其中，WG-FAST可针对特定菌种，将其菌株与参考基因组进行全基因组序列比对，构建该菌种的SNP参考矩阵，进而利用这个矩阵，对宏基因组学样品中该菌种所含菌株进行鉴定；Constrains则是通过一个菌种定量流程(MetaPhlAn)得到菌种丰度和菌种所含菌株的SNP，然后对每一个样品的SNP进行聚类得到每一个样品的SNP模式，再对这些SNP模式进行聚类，通过模型学则得到最佳的SNP模式，根据这个SNP模式对每个样品中，该菌种可能含有的菌株的个数和相对丰度做出推断。这两个流程的实现步骤分别如下：
1.WG-FAST
1)确定需要鉴定菌株的菌种；
2)通过与参考基因组的序列比对，鉴定出各菌株的SNP，构建SNP参考矩阵；
3)基于实际样品与参考基因组的比对结果，生成该样品的SNP矩阵；
4)结果实际样品的SNP矩阵和参考SNP矩阵，生成系统发生树。
2.Constrains
1)利用MetaPhlAn流程进行菌种定量；
2)通过标记序列比对，鉴定某个菌种内各菌株的SNP；
3)对上述SNP进行聚类，由聚类结果得到每个样品中的SNP模式；
4)对多个样品的SNP模式进行聚类，通过模型选择得到最佳的SNP模式；
5)根据最佳的SNP模式，推断出每个样品中该菌种各菌株的相对丰度。
现有技术的上述流程中，WG-FAST虽然能鉴定出特定菌种下的菌株，但不能给出各菌株的丰度(比例)信息，且所依赖的软件和软件库较多，包括BWA-MEM，GATK，Picard-tools，DendroPy，RAxML，BioPython，Trimmomatic和SAMtools等；而Constrains只能给出菌种下菌株的一个粗略分类和相应丰度，并不能确定具体菌株信息，并且要求比到菌种水平的读序的覆盖层数在10层(10×)以上。
综上所述，本领域尚缺乏基于宏基因组对特定群(如菌群)进行定量分析的方法，现有的方法无法准确测定具体亚群(如菌株)的丰度。此外，即使有些方法即使采用了高时空复杂度的计算，但是仍无法获得的理想的结果。
因此，本领域迫切需要开发能够对特定群中的亚群进行有效定量分析的方法。
发明内容
本发明的目的在于提供一种对特定群中的亚群进行有效定量分析的方法。
在本发明的第一方面，提供了一种非诊断性的对特定群中的亚群进行定量分析的方法，包括以下步骤：
(1)提供对应于所述特定群的(a)参考基因组序列数据、(b)参考SNP矩阵和(c)宏基因组测序数据，其中所述的特定群包括N个亚群，N为≥2的正整数；
(2)将所述宏基因组测序数据比对到对应于所述特定群的参考基因组数据上，以便获得比对结果；
(3)根据所述参考SNP矩阵的位点信息，从所述比对结果中提取出各分析位点上4种碱基A、T、C、和G各自的分布频率，构建频率矩阵；
(4)根据频率矩阵对参考SNP矩阵做二值化处理，得到二值化SNP矩阵，其中在所述二值化处理时，对于某一分析位点而言，若检出的碱基与参考矩阵上的碱基相同则记1，否则记为0；并且如果某一分析位点未被覆盖，则该分析位点不参与观测模型的构建；和
(5)基于所述的频率矩阵、所述的二值化SNP矩阵、理论碱基频率向量f(x)和观测碱基频率向量y，通过有约束的线性模型，得出所述特定群中各亚群的相对丰度，从而获得对所述特定群中的亚群的定量分析结果。
在另一优选例中，在步骤(5)中，包括以下步骤：
(5a)设所述观测碱基频率向量为y，y由所述频率矩阵得到；
以及设定理论碱基频率向量为f(x),
f(x)＝ωTx
式中，x为所述二值化SNP矩阵；
ω为列向量，即所述各亚群(1,2,3，…N)的比例值(ω1,ω2,ω3…ωN)构成的列向量ω，
上标T表示转置，行向量转置后为列向量，反之亦然；
(5b)基于下式Q1、Q2、Q3、Q4，
f(x)≈y            (Q1)
ω＝(ω1,ω2,ω3...ωN)T         (Q2)
Figure PCTCN2016105372-appb-000001
0≤ωi≤1         (Q4)
(注：Q3和Q4表示：所有亚群的比例加和为1，且各亚群比例为0到1之间的正数)；
运用回归方法，对各亚群比例进行求解，使ω*满足式Q5和Q6，从而得出所述特定群中各亚群的相对丰度，
Figure PCTCN2016105372-appb-000002
Figure PCTCN2016105372-appb-000003
各式中，
m为频率不为零的碱基总个数；
T，x，y，和ω如上定义。
在另一优选例中，N为2-50，较佳地3-20，更佳地3-10。
在另一优选例中，所述的y＝(q1，q2，….,qm)，其中，各q1，q2，….,qm为频率矩阵中不等于0的各分布频率，m为所述频率矩阵中频率不为零的碱基总个数。
在另一优选例中，所述的特定群为种，而所述的亚群为亚种，或株。
在另一优选例中，在步骤(5b)中，使理论碱基频率(f(x))与观测碱基频率(y)的残差平方和最小。
在另一优选例中，在步骤(5b)中，使用最小二乘法求解。
在另一优选例中，在步骤(5b)中，使用序贯最小二乘法求解。
在另一优选例中，所述群选自下组：细菌、病毒、肿瘤细胞、真菌、衣原体、支原体。
在另一优选例中，所述亚群选自下组：亚种、异质细胞、菌株。
在另一优选例中，所述的宏基因组测序数据的测序深度的下限为0.01×(较佳地0.05×，更佳地0.1×)，而测序深度的上限为至10×(较佳地5×，更佳地2×，更佳地1×)。
在本发明的第二方面，提供了一种非诊断性的对特定群中的亚群进行定量分析的方法，包括以下步骤：
(1)提供对应于所述特定群的(a)参考基因组序列数据或参考转录组序列数据、(b)参考SNP矩阵和(c)宏转录组测序数据，其中所述的特定群包括N个亚群，N为≥2的正整数；
(2)将所述宏转录组测序数据比对到对应于所述特定群的参考基因组数据或参考转录组序列数据上，以便获得比对结果；
(3)根据所述参考SNP矩阵的位点信息，从所述比对结果中提取出各分析位点上4种碱基A、T、C、和G各自的分布频率，构建频率矩阵；
(4)根据频率矩阵对参考SNP矩阵做二值化处理，得到二值化SNP矩阵，其中在所述二值化处理时，对于某一分析位点而言，若检出的碱基与参考矩阵上的碱基相同则记1，否则记为0；并且如果某一分析位点未被覆盖，则该分析位点不参与观测模型的构建；和
(5)基于所述的频率矩阵、所述的二值化SNP矩阵、理论碱基频率向量f(x)和观测碱基频率向量y，通过有约束的线性模型，得出所述特定群中各亚群的相对丰度，从而获得对所述特定群中的亚群的定量分析结果。
在另一优选例中，所述的参考SNP矩阵包括基于基因组序列的参考SNP矩阵、或基于转录组序列的参考SNP矩阵。
在另一优选例中，在步骤(5)中，包括以下步骤：
(5a)设所述观测碱基频率向量为y，y由所述频率矩阵得到；
以及设定理论碱基频率向量为f(x),
f(x)＝ωTx
式中，x为所述二值化SNP矩阵；
ω为列向量，即所述各亚群(1,2,3，…N)的比例值(ω1,ω2,ω3…ωN)构成的列向量ω，
上标T表示转置，行向量转置后为列向量，反之亦然；
(5b)基于下式Q1、Q2、Q3、Q4，
f(x)≈y             (Q1)
ω＝(ω1,ω2,ω3...ωN)T       (Q2)
Figure PCTCN2016105372-appb-000004
0≤ωi≤1          (Q4)
(注：Q3和Q4表示：所有亚群的比例加和为1，且各亚群比例为0到1之间的正数)；
运用回归方法，对各亚群比例进行求解，使ω*满足式Q5和Q6，从而得出所述特定群中各亚群的相对丰度，
Figure PCTCN2016105372-appb-000005
Figure PCTCN2016105372-appb-000006
各式中，
m为频率不为零的碱基总个数；
T，x，y，和ω如上定义。
在另一优选例中，N为2-50，较佳地3-20，更佳地3-10。
在另一优选例中，所述的y＝(q1，q2，….,qm)，其中，各q1，q2，….,qm为频率矩阵中不等于0的各分布频率，m为所述频率矩阵中频率不为零的碱基总个数。
在另一优选例中，所述的特定群为种，而所述的亚群为亚种，或株。
在另一优选例中，在步骤(5b)中，使理论碱基频率(f(x))与观测碱基频率(y)的残差平方和最小。
在另一优选例中，在步骤(5b)中，使用最小二乘法求解。
在另一优选例中，在步骤(5b)中，使用序贯最小二乘法求解。
在另一优选例中，所述群选自下组：细菌、病毒、肿瘤细胞、真菌、衣原体、支原体。
在另一优选例中，所述亚群选自下组：亚种、异质细胞、菌株。
在另一优选例中，所述的宏转录组测序数据的测序深度的下限为0.01×(较佳地0.05×，更佳地0.1×)，而测序深度的上限为至10×(较佳地5×，更佳地2×，更佳地1×)。
应理解，在本发明范围内中，本发明的上述各技术特征和在下文(如实施例)中具体描述的各技术特征之间都可以互相组合，从而构成新的或优选的技术方案。限于篇幅，在此不再一一累述。
附图说明
图1显示了对特定菌种内的菌株进行定量的流程图，其中虚线部分为流程构建过程特有数据或步骤。
图2显示了流程对模拟简单环境样品处理结果的残差，其中各测序深度下有三个平行实验。
图3显示了流程对模拟简单环境样品处理结果的正误比例，其中各测序深度下有三个平行实验。
图4显示了流程对模拟复杂环境样品处理结果的残差，其中各测序深度下有三个平行实验。
图5显示了流程对模拟复杂环境样品处理结果的正误比例，其中各测序深度下有三个平行实验。
具体实施方式
本发明人经过广泛而深入地研究，首次开发了一种对特定群中的亚群进行有效定量分析的方法。本发明基于参考基因组或转录组序列数据、参考SNP矩阵和宏基因组或转录组测序数据，通过对各数据(库)的比对构建频率矩阵，并通过二值化处理、通过理论碱基频率向量f(x)和观测碱基频率向量y构建有约束的线性模型，从而有效地得出特定群中各亚群的定量检测结果。在此基础上，完成了本发明。
具体地，以宏基因组学为例，本发明基于特定群的参考基因组序列数据、参考SNP矩阵和宏基因组测序数据，将宏基因组测序数据与参考基因组数据进行比对，获取比对结果；根据上述参考SNP矩阵的位点信息，从所述比对结果中提取出各分析位点上4种碱基A、T、C、和G各自的分布频率，构建频率矩阵；根据频率矩阵对参考SNP矩阵做二值化处理，得到二值化SNP矩阵；基于上述频率矩阵、二值化SNP矩阵、理论碱基频率向量f(x)和观测碱基频率向量y，通过有约束的线性模型，得出所述特定群中各亚群的相对丰度，从而获得对所述特定群中的亚群的定量分析结果。
具体地，本发明人还提供了一种在待研究菌种测序层数(基因组被测序读序覆盖的次数)较少的情况下，实现对包含于参考矩阵的菌株的比例的定量方法；构建合适的模型，使流程实现过程的时空复杂度较小。
对亚群进行定量分析的方法
以下结合一个具体例子进一步描述本发明的技术方案。在该实施例中，先提供：1)特定菌种的参考基因组，2)该菌种的参考SNP矩阵，3)宏基因组测序数据。
流程处理过程包括：1)使用序列比对软件，将宏基因组测序数据比对到参 考基因组上，得到比对结果；2)根据参考SNP矩阵的位点信息，从上述比对结果中提取出各位点上4种碱基(A、T、C、G)的频率，构建频率矩阵；3)结合参考SNP矩阵和频率矩阵，构建观测模型；4)运用回归方法，对各菌株比例进行估计。该流程如附图1所示。
假设某个菌种下有3个菌株，其中菌株Ⅰ的基因组为该菌种的参考基因组(reference,ref)，其参考SNP矩阵包括n个位点，该矩阵示例如下：
表1 SNP矩阵示例
Figure PCTCN2016105372-appb-000007
流程各步骤的详细情况介绍如下：
1.序列比对及结果提取
使用序列比对软件，将宏基因组测序数据比对到参考基因组上，经过碱基质量控制后得到包含比对信息的文件。
2.频率矩阵构建
根据参考SNP矩阵上的位点信息，从上述pileup文件中提取出各个位点上4种碱基的频率，若某位点未被覆盖，则将4种碱基的频率均设为0(如位点2的情况)。得到如下所示的频率矩阵。
表2 频率矩阵示例
Figure PCTCN2016105372-appb-000008
3.模型构建
a)根据频率矩阵对参考SNP矩阵做二值化处理，即若某一位置上检出的碱基与参考矩阵上的碱基相同则记1，否则记为0；若某位点未被覆盖，则该位 点不参与观测的构建(如位点2的情况)。二值化的SNP矩阵如下表所示。
表3 二值化SNP矩阵示例
Figure PCTCN2016105372-appb-000009
为了求解菌株在菌种中的比例，现构建有约束的线性模型。
b)设观测碱基频率向量(由上述频率矩阵得到)为y，理论碱基频率向量(经计算得到)为f(x),二值化SNP矩阵为x，待估计菌株(1,2,3)的比例值(ω1,ω2,ω3)构成的列向量为ω，其中，
f(x)＝ωTx
上标T表示转置，列向量转置后为行向量，反之亦然。
使得，
f(x)≈y
其中，
ωT＝(ω1,ω2,ω3)
Figure PCTCN2016105372-appb-000010
0≤ωi≤1
即所有菌株的比例加和为1，且各菌株比例为0到1之间的正数。
在本例中，
y＝(0.5,0.5,0.1,0.9,...,0.3,0.3,0.4)
ωT＝(ω1,ω2,ω3)
Figure PCTCN2016105372-appb-000011
4.回归求解
如前所述，估计菌株比例的问题，即解有约束的优化问题，亦即求得ω*，
使
Figure PCTCN2016105372-appb-000012
Figure PCTCN2016105372-appb-000013
其中，m为频率不为零的碱基总个数(即表3-二值化SNP矩阵中碱基的总个数)。
即，使理论碱基频率(f(x))与观测碱基频率(y)的残差平方和最小。发明人使用序贯最小二乘法求解，并将比例较小的估计结果过滤掉，即认为此类菌株并不存在于样品中。
假设上述示例中仅第1、3、n个位点有测序序列覆盖，即
y＝(0.5,0.5,0.1,0.9,0.3,0.3,0.4)，
Figure PCTCN2016105372-appb-000014
则，
Figure PCTCN2016105372-appb-000015
Figure PCTCN2016105372-appb-000016
经上述回归方法求解，可得
ω*＝(0.2427，0.3859，0.3714).
得到特定物种下菌株的比例后，结合其他流程得到的菌种定量结果，即可得到各菌株在宏基因组样品中的相对丰度。设特定菌种的菌种丰度为a,各菌株相对丰度为向量
Figure PCTCN2016105372-appb-000017
即有，
Figure PCTCN2016105372-appb-000018
基于宏转录组测序数据对亚群进行定量分析的方法，与上述基于宏基因组测序数据的方法基本类似，不同点主要在于，用宏转录组测序数据替换宏基因组测序数据。此外，由于转录组数据与基因组数据存在映射关系，因此比对时可以采用参考基因组序列数据，也可以参考转录组序列数据。在本发明方法中，参考SNP矩阵可以基于基因组序列，也可以基于转录组序列。
本发明的主要优点包括：
(a)本发明可对宏基因组数据中特定菌种的菌株间的比例进行定量。例如用于对某个粪便样品中所含大肠杆菌菌株的比例进行定量，结合现有流程的菌种定量结果得到菌株的丰度，这将有利于提高相关研究的分辨率或疾病诊断、治疗的精准性。
(b)本发明开发了一套全新的、可自定义目标的菌株定量流程，能够帮助有需求的生物研究人员快速而准确地对宏基因组样品中同一菌种下的菌株进行定量。
(c)本发明提供的分析策略能对序列相似性较高的群组进行细分，因此不 仅可用于多种复杂微生物群落的宏基因组学研究流程，也可拓展至宏转录组学等相关组学研究，甚至可用于探究肿瘤异质性方面的科学问题，这将推动精准医学的发展。
下面结合具体实施例，进一步阐述本发明。应理解，这些实施例仅用于说明本发明而不用于限制本发明的范围。下列实施例中未注明具体条件的实验方法，通常按照常规条件如Sambrook等人，分子克隆：实验室手册(New York:Cold Spring Harbor Laboratory Press,1989)中所述的条件，或按照制造厂商所建议的条件。除非另外说明，否则百分比和份数按重量计算。本发明中所涉及的实验材料如无特殊说明均可从市售渠道获得。
实施例1模拟简单环境样品的分析
1.数据模拟
为了模拟简单环境样品，即仅含同一菌种的不同菌株的样品，本发明人选取基因组资源丰富的大肠杆菌作为研究对象，从NCBI基因组数据库中选取了44个代表菌株(2015年11月19日版本)用于构建参考矩阵，其中大肠杆菌K-12的基因组为参考基因组(reference,ref)；随后从中挑选了5个菌株，用于模拟简单环境样品测序数据，各菌株信息及其在复杂环境样品中的理论丰度表4所示。
表4 模拟简单样品详细信息
菌株名称 所占比例
Escherichia coli O103:H2str.12009 0.2
Escherichia coli O104:H4str.2011C-3493 0.2
Escherichia coli O127:H6str.E2348/69 0.2
Escherichia coli O145:H28str.4865/96 0.2
Escherichia coli O157:H7str.Sakai 0.2
本发明人模拟了一系列的测序数据，各系列中，上述5个菌株的模拟的测序深度(层数，×)分别为0.01×、0.1×、和1×，各个系列均平行生成3次。模拟的单端测序片段的长度为50bp，插入片段的大小为200bp，测序错误率为0。
2.序列比对及结果提取
2.1本发明人选择大肠杆菌K-12菌株的基因组作为参考序列，使用BWA序列比对软件，将宏基因组测序数据比对到参考基因组上，得到包含原始比对结果的BAM文件；
2.2使用samtools工具，将BAM文件按位点排序后，把每个位点上比对质量值大于10，即比对正确率大于90％的碱基输出，得到包含排序、过滤处理后的比对信息的pileup文件；
3.频率矩阵构建
根据参考SNP矩阵上的位点信息，从上述pileup文件中提取出各个位点上4种碱基的频率，若某位点未被覆盖，则将4种碱基的频率均设为0，得到频率矩阵。
4.模型构建
4.1根据频率矩阵对参考SNP矩阵做二值化处理，即若某一位置上检出的碱基与参考矩阵上的碱基相同则记1，否则记为0；若某位点未被覆盖，则该位点不参与观测的构建。
4.2为了求解菌株在菌种中的比例，现构建有约束的线性模型。设碱基频率向量为f(x)，二值化SNP矩阵为x，待估计比例值构成的向量为ω，三者满足，
f(x)＝ωT×x
其中，
ω＝(ω1,ω2,ω3.....ωn)
其中，
Figure PCTCN2016105372-appb-000019
0≤ωi≤1
即所有菌株(在本例中N为44株)的比例加和为1，且各菌株比例为0到1之间的正数。
5.回归求解
如前所述，估计菌株比例的问题，即解有约束的优化问题，亦即求得ω*，使，
Figure PCTCN2016105372-appb-000020
发明人使用序贯最小二乘法求解，并将比例小于0.0001的估计结果过滤掉，结果如表5-1、表5-2、表5-3及附图2所示。
6.流程估计结果分析
对流程估计结果的分析显示，各测序深度下的估计值与理论值间的残差分布随测序深度的增加逐渐收敛，在0.01×、0.1×及1×的深度下，平均残差平方和分别为0.0578、0.0119和0.0048，模拟数据所选菌株的平均总比例分别为0.7580、0.8496和0.9235。另一方面，随着测序深度的增加，各系列平行实验的估计值更加接近，即精确度更高，结果如表5-1、表5-2、表5-3、表6及附图2、附图3所示。
表5-1 模拟简单环境样品菌株估计结果
Figure PCTCN2016105372-appb-000021
Figure PCTCN2016105372-appb-000022
表5-2 模拟简单环境样品菌株估计结果
Figure PCTCN2016105372-appb-000023
Figure PCTCN2016105372-appb-000024
表5-3 模拟简单环境样品菌株估计结果
Figure PCTCN2016105372-appb-000025
Figure PCTCN2016105372-appb-000026
表6 模拟简单环境样品菌株估计的正误比例
Figure PCTCN2016105372-appb-000027
上述结果表明，对于由5个亚群构成的特定群，当测序深度为0.01×或更高时，本发明方法就可以非常有效地对特定群中的各亚群进行定量分析，并提供可靠的检测结果。其中，当测序深度为0.01×时，检测值与理论值的误差小于约30％；当测序深度为0.1×时，检测值与理论值的误差小于约20％；当测序深度为1×时，检测值与理论值的误差小于约10％。
实施例2模拟复杂环境样品的分析
1.数据模拟
为了模拟复杂环境样品，即包含多个菌种和大肠杆菌多个菌株的样品。本发明人在常见的拟杆菌门、厚壁菌门、变形菌门和放线菌门中选取了10个物种，即长双歧杆菌(Bifidobacterium longum)、迟缓埃格特菌(Eggerthella lenta)、吉氏副拟杆菌(Parabacteroides distasonis)、活泼瘤胃球菌(Ruminococcus gnavus)、艰难梭状芽胞杆菌(Peptoclostridium difficile)、韦荣球菌HPA0037(Veillonella sp.HPA0037)、链球菌I-P16(Streptococcus sp.I-P16)、柯氏柠檬酸杆菌(Citrobacter koseri)、埃希氏菌艾伯替埃希氏菌(Escherichia albertii)和大肠杆菌(Escherichia coli)，其中大肠杆菌从前述参考矩阵中选取了5个菌株，各物种选取的菌株信息及其在复杂环境样品中的理论丰度如下表所示。
表7 复杂微生物群落详细信息
菌株名称 所占比例
Escherichia coli UMN026 1/14
Escherichia coli UMEA_3318_1 1/14
Escherichia coli NA114 1/14
Escherichia coli IAI39 1/14
Escherichia coli CFT073 1/14
Escherichia albertii KF1 1/14
Citrobacter_koseri ATCC BAA-895 1/14
Streptococcus sp.I-P16 1/14
Veillonella sp.HPA0037 1/14
Peptoclostridium difficile 630 1/14
Ruminococcus gnavus AGR2154 1/14
Parabacteroides distasonis ATCC 8503 1/14
Eggerthella lenta DSM 2243 1/14
Bifidobacterium longum NCC2705 1/14
与前述方法类似，根据表7所示的物种构成，生成了系列模拟测序数据，各系列中各菌株测序深度分别为0.01×、0.1×和1×，每个系列平行重复三次。
2.流程处理过程
本实施例的序列比对及结果提取、频率矩阵构建、模型构建和回归求解过程均与实施例1一致。
3.流程处理结果分析
运行本发明处理后，对结果的分析发现，与前述探究结果类似，随着测序深度的增加，估计值与理论值的残差平方和总体上呈减小趋势，在0.01×、0.1×和1×的深度下，平均残差平方和分别为0.0744、0.0304和0.0278，模拟数据所选菌株的平均总比例分别为0.5569、0.7797和0.7700，且对应的平行实验的结果趋于一致。由于受到大肠杆菌近缘物种测序序列的干扰，该估计效果相比简单环境样品略有下降，但在菌株测序深度达到0.1×及以上时，仍有较高的准确性和精确性，结果如表8-1、表8-2、表8-3、表9及附图4、附图5所示。
表8-1 模拟复杂环境样品菌株估计结果
Figure PCTCN2016105372-appb-000028
Figure PCTCN2016105372-appb-000029
表8-2 模拟复杂环境样品菌株估计结果
Figure PCTCN2016105372-appb-000030
Figure PCTCN2016105372-appb-000031
表8-3 模拟复杂环境样品菌株估计结果
Figure PCTCN2016105372-appb-000032
Figure PCTCN2016105372-appb-000033
表9 模拟复杂环境样品菌株估计的正误比例
Figure PCTCN2016105372-appb-000034
上述结果表明，即使是对于14个亚群构成的特定群，当测序深度为0.1 ×或更高时，本发明方法可以非常有效地对特定群中的各亚群进行定量分析，并提供可靠的检测结果。当测序深度为0.1×和1×时，检测值与理论值的误差均小于25％。
在本发明提及的所有文献都在本申请中引用作为参考，就如同每一篇文献被单独引用作为参考那样。此外应理解，在阅读了本发明的上述讲授内容之后，本领域技术人员可以对本发明作各种改动或修改，这些等价形式同样落于本申请所附权利要求书所限定的范围。

Claims (10)

  1. 一种非诊断性的对特定群中的亚群进行定量分析的方法，其特征在于，包括以下步骤：
    (1)提供对应于所述特定群的(a)参考基因组序列数据、(b)参考SNP矩阵和(c)宏基因组测序数据，其中所述的特定群包括N个亚群，N为≥2的正整数；
    (2)将所述宏基因组测序数据比对到对应于所述特定群的参考基因组数据上，以便获得比对结果；
    (3)根据所述参考SNP矩阵的位点信息，从所述比对结果中提取出各分析位点上4种碱基A、T、C、和G各自的分布频率，构建频率矩阵；
    (4)根据频率矩阵对参考SNP矩阵做二值化处理，得到二值化SNP矩阵，其中在所述二值化处理时，对于某一分析位点而言，若检出的碱基与参考矩阵上的碱基相同则记1，否则记为0；并且如果某一分析位点未被覆盖，则该分析位点不参与观测模型的构建；和
    (5)基于所述的频率矩阵、所述的二值化SNP矩阵、理论碱基频率向量f(x)和观测碱基频率向量y，通过有约束的线性模型，得出所述特定群中各亚群的相对丰度，从而获得对所述特定群中的亚群的定量分析结果。
  2. 如权利要求1所述的方法，其特征在于，在步骤(5)中，包括以下步骤：
    (5a)设所述观测碱基频率向量为y，y由所述频率矩阵得到；
    以及设定理论碱基频率向量为f(x),
    f(x)＝ωTx
    式中，x为所述二值化SNP矩阵；
    ω为列向量，即所述各亚群(1,2,3，…N)的比例值(ω1,ω2,ω3…ωN)构成的列向量ω，
    上标T表示转置，行向量转置后为列向量，反之亦然；
    (5b)基于下式Q1、Q2、Q3、Q4，
    f(x)≈y           (Q1)
    ω＝(ω1,ω2,ω3...ωN)T           (Q2)
    Figure PCTCN2016105372-appb-100001
    0≤ωi≤1           (Q4)
    运用回归方法，对各亚群比例进行求解，使ω*满足式Q5和Q6，从而得出所述特定群中各亚群的相对丰度，
    Figure PCTCN2016105372-appb-100002
    Figure PCTCN2016105372-appb-100003
    各式中，
    m为频率不为零的碱基总个数；
    T，x，y，和ω如上定义。
  3. 如权利要求1所述的方法，其特征在于，N为2-50，较佳地3-20，更佳地3-10。
  4. 如权利要求1所述的方法，其特征在于，所述的y＝(q1，q2，….,qm)，其中，各q1，q2，….,qm为频率矩阵中不等于0的各分布频率，m为所述频率矩阵中频率不为零的碱基总个数。
  5. 如权利要求1所述的方法，其特征在于，所述的特定群为种，而所述的亚群为亚种，或株。
  6. 如权利要求2所述的方法，其特征在于，在步骤(5b)中，使理论碱基频率(f(x))与观测碱基频率(y)的残差平方和最小。
  7. 如权利要求2所述的方法，其特征在于，在步骤(5b)中，使用最小二乘法求解。
  8. 如权利要求2所述的方法，其特征在于，在步骤(5b)中，使用序贯最小二乘法求解。
  9. 如权利要求1或2所述的方法，其特征在于，所述群选自下组：细菌、病毒、肿瘤细胞、真菌、衣原体、支原体；
    和/或所述亚群选自下组：亚种、异质细胞、菌株。
  10. 一种非诊断性的对特定群中的亚群进行定量分析的方法，其特征在于，包括以下步骤：
    (1)提供对应于所述特定群的(a)参考基因组序列数据或参考转录组序列数据、(b)参考SNP矩阵和(c)宏转录组测序数据，其中所述的特定群包括N个亚群，N为≥2的正整数；
    (2)将所述宏转录组测序数据比对到对应于所述特定群的参考基因组数据或参考转录组序列数据上，以便获得比对结果；
    (3)根据所述参考SNP矩阵的位点信息，从所述比对结果中提取出各分析位点上4种碱基A、T、C、和G各自的分布频率，构建频率矩阵；
    (4)根据频率矩阵对参考SNP矩阵做二值化处理，得到二值化SNP矩阵，其中在所述二值化处理时，对于某一分析位点而言，若检出的碱基与参考矩阵上的碱基相同则记1，否则记为0；并且如果某一分析位点未被覆盖，则该分析位点不参与观测模型的构建；和
    (5)基于所述的频率矩阵、所述的二值化SNP矩阵、理论碱基频率向量f(x)和观测碱基频率向量y，通过有约束的线性模型，得出所述特定群中各亚群的相对丰度，从而获得对所述特定群中的亚群的定量分析结果。
PCT/CN2016/105372 2016-11-10 2016-11-10 一种对特定群中的亚群进行定量分析的方法 WO2018086045A1 (zh)

Priority Applications (2)

Application Number Priority Date Filing Date Title
CN201680090780.0A CN109997193B (zh) 2016-11-10 2016-11-10 一种对特定群中的亚群进行定量分析的方法
PCT/CN2016/105372 WO2018086045A1 (zh) 2016-11-10 2016-11-10 一种对特定群中的亚群进行定量分析的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/CN2016/105372 WO2018086045A1 (zh) 2016-11-10 2016-11-10 一种对特定群中的亚群进行定量分析的方法

Publications (1)

Publication Number Publication Date
WO2018086045A1 true WO2018086045A1 (zh) 2018-05-17

Family

ID=62109084

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2016/105372 WO2018086045A1 (zh) 2016-11-10 2016-11-10 一种对特定群中的亚群进行定量分析的方法

Country Status (2)

Country Link
CN (1) CN109997193B (zh)
WO (1) WO2018086045A1 (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112151120A (zh) * 2020-09-23 2020-12-29 易会广 用于快速转录组表达定量的数据处理方法、装置及存储介质
CN112233726A (zh) * 2020-10-23 2021-01-15 深圳未知君生物科技有限公司 一种细菌菌株的分析方法、分析装置和存储介质
CN112786102A (zh) * 2021-01-25 2021-05-11 北京大学 一种基于宏基因组学分析精准识别水体中未知微生物群落的方法
CN114300055A (zh) * 2021-12-28 2022-04-08 江苏先声医学诊断有限公司 优化的宏基因组纳米孔测序数据定量方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103955629A (zh) * 2014-02-18 2014-07-30 吉林大学 基于模糊k均值的宏基因组片段聚类方法
CN106055924A (zh) * 2016-05-19 2016-10-26 完美（中国）有限公司 微生物操作分类单元确定和序列辅助分离

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
DK3144672T3 (en) * 2007-11-21 2018-12-03 Cosmosid Inc GENOME IDENTIFICATION SYSTEM
US8478544B2 (en) * 2007-11-21 2013-07-02 Cosmosid Inc. Direct identification and measurement of relative populations of microorganisms with direct DNA sequencing and probabilistic methods
US10127346B2 (en) * 2011-04-13 2018-11-13 The Board Of Trustees Of The Leland Stanford Junior University Systems and methods for interpreting a human genome using a synthetic reference sequence
CN102952854B (zh) * 2011-08-25 2015-01-14 深圳华大基因科技有限公司 单细胞分类和筛选方法及其装置
GB2519255B (en) * 2013-02-01 2016-01-06 Univ California Methods for genome assembly and haplotype phasing
CN105095688A (zh) * 2014-08-28 2015-11-25 吉林大学 检测人体肠道宏基因组的细菌群落及丰度的方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103955629A (zh) * 2014-02-18 2014-07-30 吉林大学 基于模糊k均值的宏基因组片段聚类方法
CN106055924A (zh) * 2016-05-19 2016-10-26 完美（中国）有限公司 微生物操作分类单元确定和序列辅助分离

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
ALEXEEV, D. ET AL.: "Bacterial Rose Garden for Metagenomic SNP-Based Phylogeny Visualization", BIODATA MINING, vol. 8, no. 10, 21 March 2015 (2015-03-21), pages 1 - 12, XP021218346 *
CHEN, BO ET AL.: "Features Extraction and Dimensions Reduction in Metagenomic Binning Problem", COMPUTER SYSTEMS & APPLICATIONS, vol. 24, no. 11, 31 December 2015 (2015-12-31), pages 31 - 37 *
LUO, C.W. ET AL.: "ConStrains Identifies Microbial Strains in Metagenomic Datasets.", NAT. BIOTECHNOL., vol. 33, no. 10, October 2015 (2015-10-01), pages 1045 - 1052, XP055503134 *
NAYFACH, S. ET AL.: "An Integrated Metagenomics Pipeline for Strain Profiling Reveals Novel Patterns of Bacterial Transmission and Biogeography", GENOME RESEARCH, vol. 26, no. 11, November 2016 (2016-11-01), pages 1612 - 1625, XP055503139 *
SAHL, J.W. ET AL.: "Phylogenetically Typing Bacterial Strains from Partial SNP Genotypes Observed from Direct Sequencing of Clinical Specimen Metagenomic Data", GENOME MEDICINE, vol. 7, no. 52, 9 June 2015 (2015-06-09), pages 1 - 13, XP055503133 *

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112151120A (zh) * 2020-09-23 2020-12-29 易会广 用于快速转录组表达定量的数据处理方法、装置及存储介质
CN112151120B (zh) * 2020-09-23 2024-03-12 易会广 用于快速转录组表达定量的数据处理方法、装置及存储介质
CN112233726A (zh) * 2020-10-23 2021-01-15 深圳未知君生物科技有限公司 一种细菌菌株的分析方法、分析装置和存储介质
CN112786102A (zh) * 2021-01-25 2021-05-11 北京大学 一种基于宏基因组学分析精准识别水体中未知微生物群落的方法
CN114300055A (zh) * 2021-12-28 2022-04-08 江苏先声医学诊断有限公司 优化的宏基因组纳米孔测序数据定量方法

Also Published As

Publication number Publication date
CN109997193A (zh) 2019-07-09
CN109997193B (zh) 2023-03-14

Similar Documents

Publication Publication Date Title
Almeida et al. A new genomic blueprint of the human gut microbiota
AU2022268283B2 (en) Phenotype/disease specific gene ranking using curated, gene library and network based data structures
Dubinkina et al. Assessment of k-mer spectrum applicability for metagenomic dissimilarity analysis
Nayfach et al. Toward accurate and quantitative comparative metagenomics
Scholz et al. Strain-level microbial epidemiology and population genomics from shotgun metagenomics
Wen et al. Evaluation of the reproducibility of amplicon sequencing with Illumina MiSeq platform
Siegwald et al. Assessment of common and emerging bioinformatics pipelines for targeted metagenomics
Girotto et al. MetaProb: accurate metagenomic reads binning based on probabilistic sequence signatures
Li Microbiome, metagenomics, and high-dimensional compositional data analysis
Mangul et al. ROP: dumpster diving in RNA-sequencing to find the source of 1 trillion reads across diverse adult human tissues
Lawson et al. A molecular epidemiological and genetic diversity study of tuberculosis in Ibadan, Nnewi and Abuja, Nigeria
WO2018086045A1 (zh) 一种对特定群中的亚群进行定量分析的方法
US20200294628A1 (en) Creation or use of anchor-based data structures for sample-derived characteristic determination
Snedecor et al. Fast and accurate kinship estimation using sparse SNPs in relatively large database searches
Baran et al. Joint analysis of multiple metagenomic samples
Nayak et al. Quality Control Pipeline for Next Generation Sequencing Data Analysis
Carrieri et al. A fast machine learning workflow for rapid phenotype prediction from whole shotgun metagenomes
CN113260710A (zh) 用于通过多个定制掺合混合物验证微生物组序列处理和差异丰度分析的组合物、系统、设备和方法
Ghaddar et al. Denoising sparse microbial signals from single-cell sequencing of mammalian host tissues
Mitra et al. Short clones or long clones? A simulation study on the use of paired reads in metagenomics
CN109360603A (zh) 确定肠道细菌亚种的方法及设备
Wickramarachchi Models and Algorithms for Metagenomics Analysis and Plasmid Classification
Sengupta et al. Classification and identification of fungal sequences using characteristic restriction endonuclease cut order
Park et al. Metagenomic association analysis of gut symbiont lactobacillus reuteri without host-specific genome isolation
Oñate et al. Abundance-based reconstitution of microbial pan-genomes from whole-metagenome shotgun sequencing data

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 16921123

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 16921123

Country of ref document: EP

Kind code of ref document: A1