Accurate Binning of Metagenomic Contigs Using Composition, Coverage, and Assembly Graphs

算法感觉有一个问题啊，他只在含有marker gene的contig里面选择增加新bin，但如果一个物种的genome里面没有marker基因呢，是不是没法作为一个bin出现了，可以考虑下contig graph的结构来看，距离近的且没有被添加到任何bin中的contigs 自己组成一个新bin。这一步可以在step3 结束后进行。

ABSTRACT

Metagenomic 作用是什么
Metagenomic 的流程是什么
目前软件可以改进的地方
本文提出的MetaCoAG简述，然后方法的结果和突破。

Metagenomic 可以恢复不同物种的不同基因原料，可以给微生物社区提供有效的信息。Metagenomic bin不同组织的序列，然后将短read组装成更长的contig，然后将这些contig分类的不同的组中，代表宏基因组样本中不同的分类的组。目前的软件没有用到组装图。本文提出的MetaCoAG使用了组装图和contig的覆盖信息，用single-copy marker gene 估计初始bins的数量，迭代的将contigs分配到bin中，通过binning 过程动态调整bins的数量，比其他结果都好，而且是以第一个直接使用了组装图的方法。

INTRODUCTION

Metagenomic的背景及为什么需要
Metagenomic 目前的工具的流程
为什么 reference-free的越来越流行
reference-free的主要用的特征
目前最新的工具介绍
本文使用marker的介绍
之前工具的缺陷的点
为什么需要本文的工具
本文工具介绍

二代数据的出现，Metagenomic 的背景，直接从环境样本中测序，识别样本的组成和存在的微生物，可以用来对下游分析，为了促进这种分析，Metagenomic

在组装前 bin这些reads，因为read长度限制，所以不太可靠，因此，流行的宏基因组流程的分析是先将短reads组装成更长的contig，然后将这些contigs 分到不同组中，代表不同的分类群，contigs的bins有助于构建组装的宏基因组，代表部分生物完整的基因组。

最新的contig-binning方法主要有两类

基于reference的，通过与参考序列进行比较，对带有已知分类组特征的 contig 进行分类，
不基于reference的，将contig根据它们的genomic特征，聚类到未知标签的bins中，

之前是基于reference的，但由于reference可能不完整或者质量低，而且之前未识别的微生物参考基因组可能不可靠，因此，不基于reference的越来越多，并且它们可以不依赖参考数据库识别新物种。

不基于reference的主要使用两个特征去bin（分箱）

组成：以 k-mer 的归一化频率获得
覆盖深度，比对到contig每个碱基的read的平均覆盖深度，这些工具通过结合这两个信息，提升了性能。

最新的Vamb、LRbinner、RepBin采用机器学习去获取种类序列的信号到低纬特征，促进聚类。然而准确重建相似组成和覆盖深度的仍具挑战。

估计样本中的样本数量是宏基因组中的另外的一个挑战，最近的binning 工具采用single-copy marker gene （在基因组中只出现一次，并且出现在大多数细菌基因组中）估计种类的数量。MaxBin ，MaxBin2 and SolidBin 用到了这个信息，但他们只有了一个marker，所以值得探索如何用多个marker获取一个更好的bins数量的估计结果，获取过呢更多的contig的特征提升bin

Special assemblers known as metagenomic ，大多数组装工具使用 assembly graph作为主要的结构，简化DBG图获的contig。metaSPAdes 输出的contig中带有连接信息，但现有的bins工具忽略了contig中有效的连接信息。

最近的bin-refinement工具。GraphBin，GraphBin2，METAMVGL and GraphPlas 优化bins 结果，这些工具依靠目前的binning 工具并且不能动态的调整bins的数量。此外罪行的metabinner DAS和MetaWRAP整合和优化了多种binning方法的结果。尽管这些工具都提升了binning的性能，他们仍然需要其他现有binning工具的结果，并且有些不能动态调整bins的数量。因此值得探索一个单独的contig-bin工具，使用assembly graph信息，compositing 和 coverage 信息。

MetaCoAG 不使用reference ，binning metagenomic 的方法。

比对实验中好像没和refinement后的binning结果比较？

METHODS

Step 0: Assemble reads into contigs and construct the assembly graph

用metaSPAdes 获取contigs和assembly graph，输入到MetaCoAG。也可以使用其他工具的产生的assembly graph 像MEGAHIT和metaFlye

Step 1: Identify contigs with single-copy marker genes

Single-copy marker gene 是特殊的marker gene 在细菌的基因组中只出现一次，并且在大多数的细菌染色体中都出现（conserverd）。用FragGeneScan 和 HMMER 识别包含 marker的contigs。如果marker的50%以上可以比对到一个contig，这个marker被认为被包含在这个contig。像MaxBin，MaxBin2 和 SolidBin，MetaCoAG使用marker识别属于不同物种的contigs，如果多个contig包含相同的marker，它们分别都属于不同的物种。

Step 2: Order single-copy marker genes and estimate the number of initial bins

包含相一个marker的contigs应该来自不同的物种，理想情况下，包含相同marker的contig数量就是样本中物种的数量，实际上组装的不完整性（fragmented)和错误，会降低包含一个marker的contig数量。

将marker按照包含它的contig的数量降序排序，这个列表称为SMG，每个marker g_i，有一个包含它的contig的集合C(g_i) ，初始bins的数量设置为第一个SMG中第一个marker gene的contig数量。

Step 3: Bin contigs with single-copy marker genes

Step 3a: Initialize bins

第一个marker gene中，每个contig初始化为一个bin，bin的数量在binning的过程中会改变的。

Calculating composition and coverage similarities.

最常用的表征组成信息的基因组特征是4核苷酸频率（136个典型的4-mer 也叫 tetramers 四聚体），使用总四核苷酸数量归一化每个contig的四核苷酸频率，获得它们的向量tetra(c)。

下面公式中dist_E是欧几里得距离，使用这个公式计算两个contig之间的核苷酸组成距离。(一个值)

d_{tetra}(c,c') = dist_E(tetra(c),tetra(c'))

使用下面公式计算两个contig的相似性（是否属于同一物种），N_{intra} 和 N_{inter} 分别是参数为\mu_{intra} \sigma^2_{intra} 和 \mu_{inter} \sigma^2_{inter} 的高斯分布（根据MaxBin最新的值设定），这个是通过分析来自相同genome（intra）和来自不同genome（inter）序列之间的核苷酸频率的欧几里得距离计算得到。如果距离小，他们更相似，更有可能属于同一个genome。（来自同一个genome的contig服从一个高斯分布g1（u1，delta1 为均值，方差），来自不同genome的服从另一个高斯分布g2，将d作为x带入g1,g2得到两个概率，g1/(g1+g2)是d服从g1的概率。）

S_{comp}(c‚ c') = \frac {N_{intra}(d_{tetra}(c‚ c')| \mu_{intra}, \sigma^2_{intra} )} {N_{intra}(d_{tetra}(c‚ c')| \mu_{intra}, \sigma^2_{intra} ) + N_{inter}(d {tetra}(c‚ c')| \mu_{inter}, \sigma^2_{inter} )}

使用核苷酸序列coverage泊松分布计算两个contig的相似程度。covn(c) 是样本n中 contig c 的coverage值，M是样本数量，poisson是泊松概率质量函数。（n代表碱基，M表示有n个碱基，一个泊松分布表示，一个c中碱基的coverage，在c'的这个碱基泊松分布下，c的coverage发生的概率。）

𝑆_{𝑐⁢𝑜⁢𝑣}⁡(𝑐,𝑐')=min⁡(\prod_{n=1}^𝑀 𝑃⁢𝑜⁢𝑖⁢𝑠⁢𝑠⁢𝑜⁢𝑛⁡(𝑐⁢𝑜⁢𝑣_𝑛⁡(𝑐)|𝑐⁢𝑜⁢𝑣_𝑛⁡(𝑐')),\prod_{n=1}^𝑀 𝑃⁢𝑜⁢𝑖⁢𝑠⁢𝑠⁢𝑜⁢𝑛⁡(𝑐⁢𝑜⁢𝑣_𝑛⁡(𝑐')|𝑐⁢𝑜⁢𝑣_𝑛⁡(𝑐)))

Step 3b: Construct a weighted bipartite graph and find a minimum-weight full matching

在bins中任意一个bin 和 C(g_i)的contigs 构建二分图，然后提出了一个计算contig 和现有bin之间权重的方法。

w_{c2B}(c,B) = \frac{\sum_{c' \in B}w_{c2c}(c,c')}{|B|}

其中 wc2c(c,c')是两个contig c,c' 属于一个物种的可能性。

w_{c2c}(c,c') = -(log(S_{comp}(c,c'))+log(S_{cov}(c,c')))

找一个 minimum-weight full matching ，每个contig c都会找到正好一个bin。算法的用处是找到边的集合，这个集合中的边权重和是最小的，在这个边的集合中，每个c都正好有一个边连接到一个bin？

这不就是对所有的 contig c都找一个权重最小的边及对应的bin

Step 3c: Assign contigs to existing bins and dynamically adjust bins.

之前研究表明，在连接图中互相连接的conigs更有可能属于同一个分类。

引入d_{graph}(c,B) 测量在组装图中，contig c 和bin B的连接性如何。d_{graph}(c,B) 定义为contig c 和bin B中的所有contigs的最短路径长度的平均值。w_{c2B}(c,B) 和 d_{graph}(c,B) 都会用来将contig分配给现有bin或者动态调整bins的数量（增加bin）.

定义种内阈值 w_{intra} = -(log(p_{intra})) \times M ，种间阈值w_{inter} = -(log(p_{inter})) \times M M是数据中物种的数量，每个minimum-weight full matching 种得到的候选pair(c,B) 都在下面三种情况中。

若 w_{c2B}(c,B) \le w_{intra}\ and\ d_{graph}(c,B) \le d_{limit} 将contig c 分配给B
若w_{c2B}(c,B) > w_{inter}\ and\ d_{graph}(c,B) > d_{limit} 创建一个新bin，将contig分配给这个新bin
如果不满足任意一个条件，contig c不会分配给任何bin

默认的 p_{intra}, p_{inter}, d_{limit} 根据经验设为 0.1，0.01，20。

重复3b和3c处理包含marker的所有contig，下一步处理不包含marker gene的contigs。

Step 4: Bin remaining contigs using label propagation

完成上面步骤后，买个包含marker的contig都会有一个label对应它的bin，然后从这些contig中传播label到那些unlabeled contig。

Step 4a: Propagate labels within connected components

MetaCoAG使用composition，coverage和组装图上的距离信息，从labled contig推理同一组件内的unlabelled contig 的label，对于直接连接或通过短contig连接到 labled contig c'的长contig c（at least 1000bp 长的可靠），MettaCoAG计算一个候选传播(c',c,d(c,c'),w_{c2B}(c,B')) , d(c,c') 是c和c'之间的只使用unlabeled节点的最短路径距离，B' 是 c'分配后的bin。如果有多个候选传播，d小的有优先级高，d如果相同，w越小的优先级高。

MetaCoAG迭代的选择最高优先级的候选contig传播，然后对节点打标签。如果要打标签的节点有marker，如果目标bin中没有这个marker，执行这个候选的传播。搜索深度限制为10。

Step 4b: Propagate labels across different components.

组装图的中的某些成分可能没有labeled contigs，我们需要从labeled bin 通过成分传播 label。计算任意剩余contigs的 w_{c2c}(c,c') 比较耗时。创建了一个代表contig c(B)，有一个组成和coverage 的profile，这个profile通过平均归一化所有B中的contigs计算得到。

对于每个unlabeled contigs 找一个B 有最小的w_{c2c}(c,c(B)) ,并将c分配给B。

这个传播只用长contig（1000bp），如果unlabeled contig有marker，分配给有最小w值且不包含这个marker的B中，运行step 4a。

Step 4c: Postprocessing.

如果两个bins没有公共的marker并且 w_{c2c}(c(B),c(B')) 小于 w_intra, 它们是可以合并的。MetaCoAG创建一个图，节点对应当前的bins，边对应相应的两个bins是可以合并的。查找图中的最大团，并合并最大团的中的bins。移除bins中少于三分之一marker的bins。输出相应的结果。

EXPERIMENTAL SETUP

Datasets and tools

模拟数据集 100个细菌种类，silicoseq 模拟的 paired-end read （MiSeq error model）

人体5个不同部位的模拟数据集

3个NCBI的真实数据集，

Sharon 预产婴儿肠道宏基因组 18个样本

COPD 慢性阻塞性肺疾病 (COPD) 肺微生物组的宏基因组包含 18 个样本

Deep HMP TD 舌背的人类宏基因组样本，包含 8 个样本

metaSPAdes 进行宏基因组组装，使用CoverM计算平局coverage

和MaxBin2，MetaBAT2， Vamb比较，

MaxBin2使用了概率模型和期望最大化算法，基于conttig的归一化四核苷酸频率和coverage，迭代的对contig分bin。它是使用归一化四核苷酸频率和coverage构建图，在图上聚类对contigs分bin。Vamb使用 deep variational autoencoders 对归一化四核苷酸频率和coverage编码，然后基于编码信息cluster。

使用 Assessment of Metagenome BinnERs (AMBER)，CheckM， Genome Taxonomy Database-Toolkit（GTDB-Tk）评估 binning 结果。

Evaluation metrics

用 minimap2 map contig到reference，确定simHC数据集的ground truth AMBER 评估结果

定义了一个AMBER F1-score ,precision = AMBER purity and recall = AMBER completenes

AMBER\ F1-score = 2 \times \frac{precision \times recall}{precision + recall}

对于所有的数据集，使用CheckM 去确定每个bin的完整性和纯度。

定义CheckM F1-score，precision = 1/(1+CheckM contaminatio) and recall = CheckM completenes

CheckM\ F1-score = 2 \times \frac{precision \times recall}{precision + recall}

此外，记录了不同质量bins的数量，

高质量：>80% completeness and 90% purity

中质量：>50% completeness and 80% purity

低质量：剩下的是低质量。