☰
宏基因组Binning实战:从组装到高质量MAG的完整流程
2026/9/29 18:54:20 网站建设 项目流程

1. 项目概述与全流程思路

宏基因组分析这几年几乎成了微生物生态学、环境科学、临床感染诊断领域的标配技能。一个典型的问题是:拿到一批测序数据后,我们想知道“这里面有哪些微生物?它们各自在做什么?有没有新的物种或者功能基因?”单纯靠扩增子测序(比如16S)只能回答“有哪些物种大致类群”,分辨率不够,也看不见功能。要真正把复杂微生物群落里每个成员的基因组草图和代谢潜力挖出来,就得走完整的宏基因组流程:质控、组装、Binning、Bin优化、注释。

我最初接触这个方向是被一个环境样本项目逼的——样本来自污水处理系统,里面混杂了几百种微生物,优势菌和低丰度菌共生存,用16S只能看到丰度前几十的属,而那些丰度虽低但功能关键的菌完全被淹没。后来切到宏基因组Binning路线后,才真正拿到每个物种的基因组拼接结果,能做物种界定、功能注释、代谢通路重构。这套流程现在已经成为环境微生物、肠道菌群、土壤微生物、海洋微生物方向的高频操作,也是许多生物信息学岗位的面试技能点。

这篇文章要讲的,是一条从测序数据开始,经过组装到分箱(Binning),再通过Bin精炼与去冗余得到高质量MAG(Metagenome-Assembled Genome)的完整实战链路。我默认读者已经会用Illumina测序平台的基础数据格式,会跑一些简单的命令行,但未必系统跑过Binning。

标题里的“binning ic”这个说法,我理解是在强调Binning过程中真正决定结果质量的核心信息维度——比如每条contig的覆盖深度、四核苷酸频率、GC含量、系统发育标记基因信息,这些维度合在一起才构成一个能分好箱的信号组合。网上有人把它叫“binning信息含量”(binning information content),也有人直接当成评估分箱效果的信息指标。后面我会结合参数和实际操作展开说。

适合看这篇文章的人:刚入门宏基因组分析的学生、想从扩增子测序转向全基因组水平分析的研究人员、以及做环境/临床样本需要标准化MAG产出流程的从业者。读完你会掌握一条可以直接复制运行的流程,并且知道每个环节为什么这么做、参数怎么根据实际数据调整、遇到问题怎么排查。

2. Binning方案选型与关键逻辑

2.1 为什么不能跳过Binning直接做物种注释

很多人刚接触宏基因组时有个误解:把测序数据拼接成长的contig,然后BLAST一下,不就完事了吗?实际操作中你会发现,宏基因组组装结果是一堆混合碎片的集合——不同物种、不同菌株、甚至不同的高变区序列全部搅在一起。你对着“contig_123”做物种注释,可能前半段来自A菌,中间混着B菌的质粒,后面又回到了A菌的基因组区域。这就是典型的嵌合(chimeric)片段。

Binning要解决的正是这个问题:把scaffold/contig按“来自同一个基因组”的原则重新分组。分组依据主要是两条:

  • 序列组成特征:四核苷酸频率(tetranucleotide frequency)、GC含量、密码子使用偏性。
  • 样本间的丰度共变:同一基因组的不同区域,在不同样本里的覆盖深度变化趋势应该一致。

这两条信号单独用都有局限。比如四核苷酸频率对高GC和低GC的菌群区分度好,但亲缘物种之间差异不大;覆盖度共变则依赖于样本数量和环境梯度,单样本数据几乎没法用。所以主流分箱工具都是把二者结合成特征向量,再用聚类算法划分。这也解释了为什么多组样本(multi-sample)的分箱效果通常优于单样本——coverage profile多了一个强有力的约束维度。

2.2 主流分箱工具怎么选

目前最常用的分箱工具有MetaBAT2、MaxBin2、CONCOCT、SEMIBIN、VAMB,以及后来整合它们的自动化流程MetaWRAP、Snakemake流程、还有比较新的GetMAG等。不同工具的建模方法完全不同,结果差距往往很大,不能迷信某一个工具。

MetaBAT2是现在跑批处理最稳的选项。它对不同深度覆盖、不同基因组大小的适应性强,速度也快,尤其适用几十个G数据的中等复杂度样本。底层用的是一种基于图模型的聚类方式,结合序列组成与丰度分布,内存消耗可控。

MaxBin2的核心是EM算法(期望最大化),利用标记基因的覆盖率估算每个bin的丰度和召回率,迭代优化分组。它在低丰度物种的回收上经常比MetaBAT2更敏感,但速度偏慢,而且对contig长度有要求,建议只喂长度大于1500bp(甚至2500bp)的contig。

CONCOCT用高斯混合模型做聚类,特征工程里做了一堆PCA降维和归一化,设计上有数学美感,但实际项目里稳定性一般,对数据量大的样本,很容易把丰度近似的两个菌分到同一个cluster里造成污染。

VAMB是一个基于变分自编码器的方法,可以用contig的k-mer特征和覆盖度特征做深度聚类。效果在某些场景(比如人类肠道、海洋样本)很强,但调参相对复杂,不适合新手第一版流程直接用。

我的建议:项目初期用MetaBAT2 + MaxBin2两条线跑,再在精炼步骤用MetaWRAP的bin_refinement模块整合两个工具的结果,保留完整性高、污染率低的bin集合。不推荐只用一个工具出结果的原因很简单:单个工具的分箱偏差是系统性的,不整合你会丢失一部分真实基因组。

2.3 从工具结果到MAG:精炼和去冗余才是关键

分箱工具输出的是“候选bins”,这些bin能不能升级为MAG,取决于质量评估和后续精炼。MAG的完整定义在MIMAG标准里写得很清楚:基于基因组的完整性(completeness)和污染率(contamination)等指标划分等级。如果完整性>90%、污染率<5%,可以算high-quality MAG;>50%完整且<10%污染算medium-quality。另外最好还检测到23S/16S/5S rRNA基因和至少18个tRNA基因,才算接近完成图级别。

实际输出里,initial binning出来的bin往往有大量冗余——同一个物种在两个工具里都被捕获,或者一个bin其实是另外两个bin的合并体。所以下一步必须做bin refinement和dereplication。MetaWRAP的refinement模块和dRep都是这一阶段的核心工具。dRep用ANIm或MASH算法计算bin之间的成对相似度,按95% ANI阈值去冗余,还能用完整性/污染率优先排序选代表基因组。

很多新手忽略的一点是:dRep不只是为了“去重”,它也是对MAG集合的一次质量控制。在dRep输出里,你能看到每个基因组簇的size(估计基因组大小)、ANI、完整度、污染率,这些信息比直接从分箱工具得到的列表更可信。因为两个本身很相似但各自只有60%完整度的bin,可能其实来自同一个基因组的不同部分,去冗余后merge(通过coverage和组成信息修正),你能拿到一个完整度更高的基因组。

3. 组装与分箱实操步骤

3.1 上游数据准备:质控、宿主过滤、组装三连

Binning质量高度依赖于上游数据质量。第一步是质控。推荐用fastp做adapter trimming和质量过滤,参数一般这样:

fastp -i sample_R1.fastq.gz -I sample_R2.fastq.gz \ -o clean_R1.fastq.gz -O clean_R2.fastq.gz \ --cut_front --cut_tail \ --length_required 75 \ --thread 16

这一步会把低质量碱基和短片段滤掉,避免后续拼接时产生大量碎contig。对于临床样本或宿主污染明显的样本(比如组织、粪便),还要做宿主序列过滤。人源样本可以用bowtie2比对到GRCh38参考基因组,把比对上的reads扔掉;植物、动物样本同理,取决于研究对象的宿主参考基因组是否能拿到。

接下来是组装。宏基因组组装的主流选择是MEGAHIT和metaSPAdes。这里有一个经验规律:数据量小(比如单样本5-10G)、想快速看结果时用MEGAHIT;想追求更高质量的contigs、尤其关注低丰度物种时用metaSPAdes,但内存消耗和运行时间要大得多。

我自己的标准流程是:

  • 单样本、快速探索:MEGAHIT,k-mer范围21-141,默认参数,内存控制在64G内。
  • 多样本、高深度:metaSPAdes,运行参数中指定每个样本的reads路径,让SPAdes利用跨样本信息构建coverage profile。

metaSPAdes常用命令形如:

metaspades.py -1 sample1_R1.fq -2 sample1_R2.fq ... \ -o spades_output \ -t 32 -m 200

注意:metaSPAdes不是简单的拼接器,它在组装过程中会构建一个多色de Bruijn图,并且对同一菌株的不同变体做局部组装合并。所以它对多样本共组的场景支持很友好,输出里包含scaffolds.fasta,这也就是我们分箱的输入之一。

组装结束后,必须检查组装质量的几个指标:

  • N50(权重中位长度)
  • 总组装长度
  • 最大contig长度
  • ≥1000bp的contig数量
  • 完整度检查:用BUSCO或CheckM assess

如果N50只有几百bp,别急着继续Binning,先去回溯质控流程或考虑换组装工具。N50过低说明组装碎片化严重,后续分箱会损失大量信号。

3.2 Binning前处理:从组装结果到输入特征

直接拿所有contig去分箱是常见的错误做法。短contig(比如小于1000bp)里四核苷酸频率的统计噪声太大,覆盖度估计也不准确,分箱时会产生大量错误归属。我一般用seqkit或awk过滤掉长度小于2500bp的contig,有些团队喜欢用1500bp阈值,但经过对比,2500bp在多数场景下能显著降低污染率。

seqkit seq -m 2500 scaffolds.fasta > scaffolds_min2500.fasta

接着是需要跑一个mapping来获得每个contig在各样本里的覆盖深度。这里有个关键点:如果你已经用多样本做了共组装,那mapping要分别对每个样本跑,不要混在一起,因为分箱工具依赖的就是每个样本独立的coverage profile。

Mapping我用bowtie2或minimap2,前者对Illumina短读长更合适,后者适合长读长或速度快。以bowt2为例,大致流程:

bowtie2-build scaffolds_min2500.fasta scaffolds_index for sample in sample1 sample2 sample3; do bowtie2 -x scaffolds_index \ -1 ${sample}_R1.fq.gz -2 ${sample}_R2.fq.gz \ -S ${sample}.sam --threads 16 samtools view -bS ${sample}.sam > ${sample}.bam samtools sort -o ${sample}.sorted.bam ${sample}.bam samtools index ${sample}.sorted.bam jgi_summarize_bam_contig_depths --outputDepth ${sample}.depth.txt ${sample}.sorted.bam done

注意jgi_summarize_bam_contig_depths是MetaBAT2自带脚本,输出文件里每行对应一条contig,包含length、totalAvgDepth、bam文件各自的覆盖深度等列。之后MetaBAT2直接用这个depth文件作为输入。

如果你用的是多组样本的共组装结果,别把不同样本的bam合并成一个,而是分别生成各自的depth列,最后在jgi_summarize_bam_contig_depths输入时一次性提供所有bam文件:

jgi_summarize_bam_contig_depths --outputDepth merged.depth.txt \ sample1.sorted.bam sample2.sorted.bam sample3.sorted.bam

这样输出里每个contig会有s1、s2、s3各自的coverage列,MetaBAT2就能利用跨样本丰度共变信息做聚类了。

3.3 三种主流分箱工具的参数实操

MetaBAT2命令行比较简洁,但参数理解很重要:

runMetaBAT.sh scaffolds_min2500.fasta merged.depth.txt \ -o bins_dir/bin \ -t 16

核心参数:

  • --minContig:默认2000,这里我手动过滤到了2500bp,传入时可以直接用原始文件,也可以在runMetaBAT里把minContig设为2500。若depth文件和fasta不一致会warning,需要保持统一。
  • --maxP:这个和概率阈值有关,调低会让分箱更保守,会产生更多小bin,调高则聚类更激进,bin数量少但可能污染高。默认是95%,一般情况下不需要改。
  • --minS:最小样本数量,用于考虑coverage profile时要至少几个样本有覆盖信号。默认1,即单样本也可。多样本时建议设成样本数的一半左右。

实际测试中,MetaBAT2对深度信息敏感,coverage跨度大的样本,它倾向于把相同深度的contig聚在一起,哪怕物种不同。所以务必要在上游质控阶段就把宿主污染、adapter污染清干净,否则分箱结果里很容易出现大量非目标序列。

MaxBin2运行要准备一个reads文件列表,还需要自己把contig长度过滤后的fasta传进去。先写一个配置:

run_MaxBin.pl -contig scaffolds_min2500.fasta \ -reads_list read_list.txt \ -out maxbin_out/bin \ -thread 16

read_list.txt每一行是一个样本的reads路径,支持fastq.gz。MaxBin2内部会自己做mapping计算coverage,不需要额外输入depth。但要留意:它要求输入reads是未经过human过滤的原始or clean reads,主要因为它用reads来估计丰度和补漏。

CONCOCT现在用的人有所减少,操作环节多、速度慢。它的核心输入除contig和coverage外,还要一个坐标文件。一个快速示例:

cut -f1,2,3 scaffolds_min2500.fasta > contig_coverage.bed

然后调用concoct命令。一般不建议新手再单独用CONCOCT,除非你有大量样本(>30组)且需要观察coverage共变结构时才考虑。

3.4 Bin优化与MAG质量评估

MetaBAT2和MaxBin2跑完,每个工具会输出几十上百个bin的fasta文件。接下来用MetaWRAP的bin_refinement来整合、评估。

MetaWRAP的安装有一点繁琐,建议用conda单独建环境:

conda create -n metawrap -c bioconda metawrap-env conda activate metawrap

bin_refinement用法:

metawrap bin_refinement -o refine_out \ -t 16 \ -A metabat2_bins/ \ -B maxbin2_bins/ \ -c 50 -x 10

其中-c和-x是“至少要求完整性<50,污染<10”之类的阈值设置参数,意思是最终保留的那些bin,其完整性不得低于50%,污染不得超过10%。实际处理时我通常把c设为70,因为MetaWRAP默认比较保守。

这个模块的核心动作是用CheckM对两个工具的bin做质量评估,然后做一种“bin splitting/merging”的优化:如果某个MaxBin2 bin内部覆盖度和四核苷酸组成明显分成两群,它会被拆分成两个;如果两个bin的组成特征高度相似、且合并后完整度提升而污染不显著增加,它会被合并。最后输出refined bins,并带有一个stats表格。

随后要用CheckM单独再跑一遍看指标:

checkm lineage_wf -t 16 -x fasta refine_out/ checkm_output/

CheckM会下载/使用内置的lineage-specific marker gene set,通过标记基因检测来估计完整性和污染率。完成后检查每个bin的Completeness、Contamination、Strain heterogeneity。

之后到dRep去冗余阶段。dRep能用更全面的比较算法(主要是ANIm),把多个条件下得到的MAG集合合并去重:

dRep dereplicate drep_output/ \ -g refined_bins/*.fasta \ -sa 0.95 \ -nc 0.30 \ -comp 50 \ -con 10 \ -p 16

参数说明:

  • -sa(secondary ANI threshold)设0.95,符合物种级别ANI界限。
  • -nc(minimum genome completeness)让dRep自动忽略完整度低于30%的bin,避免它们干扰去冗余。
  • -comp和-con配合CheckM,用来优先选择代表MAG。

最终产生的drep_output/dereplicated_genomes/里,就是你可以继续下游分析的高质量MAG集合。再配合GTDB-Tk做物种注释:

gtdbtk classify_wf --genome_dir drep_output/dereplicated_genomes/ \ --out_dir gtdb_out \ --cpus 16 \ --pplacer_cpus 16

GTDB-Tk如今比NCBI RefSeq的参考基因组更全面,尤其对未培养微生物的分类学注释效果更好。它会输出gtdbtk.bac120.summary.tsv和gtdbtk.ar53.summary.tsv,包含每个MAG的分类学路径、相对进化距离等。

4. 质量评估与MAG判定的关键指标详解

4.1 完整性、污染率、N50和菌株异质性

判定一个bin能不能算MAG,核心指标其实就那几个。完整性(Completeness)是估计这个bin覆盖了目标基因组多大比例。CheckM通过标记基因检测来估算:假设某个谱系里普遍存在n个单拷贝标记基因,你的bin里有m个,完整性就是m/n。

污染率(Contamination)则更微妙。它检测的是标记基因的多拷贝率——如果某个本该单拷贝的标记基因在你的bin里出现两个以上拷贝,说明bin里可能混入了两个不同物种的序列。这里有个细节:不是所有多拷贝都是污染。如果两个拷贝的等位基因差异很小,可能来自同一个物种的两个菌株,这叫菌株异质性(strain heterogeneity)。CheckM输出里的Strain heterogeneity会标记这个比例,处理时可以参考,不能一刀切。

N50(或contig N50)对MAG来说反映了bin内部组装的连续性。一个完整度95%但N50只有5kb的MAG,和N50有100kb的同完整度MAG相比,后者在基因簇结构、基因组排列分析上要可靠得多。很多下游工具(比如Prokka基因注释、antiSMASH次生代谢物基因簇预测)对长contig的结果输出更友好,短contig会导致基因簇被人为打断。

我自己的判断标准通常这样:high-quality MAG至少要完整性>90、污染<5,同时N50>10kb,最好有rRNA基因证据。Medium-quality只要完整性>50、污染<10,N50>5kb。如果连这个都达不到,我不建议把这类bin放进后续比较分析里,最多算探索性结果。

4.2 这些指标怎么影响后续分析

很多人跑完dRep,看着MAG列表和统计表就以为流程结束了,其实MAG质量直接决定下游结论的可信度。

举个例子:做宏基因组功能注释时,如果你用Prokka或DRAM对每个MAG注释代谢基因,一个污染率高的bin会同时包含两个物种的基因,KEGG代谢通路完整性会被严重高估——两个物种各有半个通路,拼在一起变成了一条完整的通路,这种假阳性在差异分析里非常致命。

同样,做泛基因组分析时,bin的完整度不够,会导致gene presence/absence矩阵里大量“缺失基因”,这些缺失其实不是生物学的缺失,而是技术性的假阴性。如果要发文章,审稿人大概率会要求MIMAG级别的质量标准表,这时候你每个MAG背后是90%完整还是70%完整,影响很大。

因此我强烈建议在流程终点的报告里,为每个MAG保留四项记录:

  • 分箱来源(MetaBAT2、MaxBin2、refinement合并产物)
  • CheckM完整性与污染率
  • dRep聚类后所属ANI cluster信息
  • GTDB-Tk分类学注释以及16S/tRNA检测结果

这样后续不管做统计分析还是投稿,都有完整溯源。

5. 常见问题与排查技巧实录

5.1 组装质量差导致分箱崩溃

表现:MetaBAT2只跑出几个大bin,或者大部分contig都被丢到unbinned里。

原因大概率是组装N50过低。尤其高复杂度样本(比如土壤),如果总contig数几十万,但N50只有几百bp,说明de Bruijn图断得太碎,coverage信号被截断成碎片,分箱聚类时根本没有足够的特征长度来稳定判断。

处理办法:

  • 先检查原始reads是否有大量duplicates或低质量碱基:用FastQC和multiQC确认。
  • 提高组装深度阈值。MEGAHIT有--min-contig-len参数,比如设500,丢弃短的初始输出。
  • 改用metaSPAdes重跑,内存充足时优先。
  • 如果样本复杂度实在太高,考虑先做co-assembly能不能提升N50——把同类型样本合并共组装,可以提升低丰度物种的覆盖深度。

5.2 MetaBAT2和MaxBin2结果差异很大

正常。两个工具的聚类算法完全不同,一个对coverage敏感,一个对序列组成敏感。所以它们的差异恰恰说明数据集里有模式冲突。这时候bin_refinement的合并/拆分机制非常有用。但有个前提:如果两个工具各自跑出来的bin数量差别超过3倍,你要回头检查是否有一个工具的参数设置不合理(比如MaxBin2用了过短的contig,或者MetaBAT2的depth文件里有大量0覆盖contig)。

经验法则:两个工具结果做intersection的部分通常是可信的;A有B没有的部分需要人工查看CheckM报告再做决断。

5.3 CheckM报“No marker genes found”或完整性极低

这种情况常见于非常新的谱系或者高度分化的基因组。CheckM的marker set依赖数据库,如果你的样本来自极端环境(比如嗜酸热菌)或者较偏门的生态位,数据库里的lineage-specific marker可能没有覆盖到位。

替代方案:

  • 用BUSCO的genome模式重新评估,其实思路类似,也基于单拷贝直系同源基因。
  • 用GTDB-Tk的分类结果反向定位参考基因组近缘种,手动比较保守基因存在情况。
  • 如果之后还要发文章,可以考虑补充长读长测序(Nanopore/PacBio),做hybrid assembly,这样完整度会显著改善。

5.4 污染率下不来,但bin看起来确实像单一物种

只有当bin里存在两个不同基因组的序列片段但它们的GC和覆盖度高度相似时,才会出现这种情况。比如一对共生菌,或者一对亲缘关系非常近的姐妹种。CheckM里的strain heterogeneity会把这种情况显示出来。

一个折中的处理:先把污染率高的bin拿出来,做一次更严格的二次分箱(rebinning),MetaWRAP也有reassemble_bins模块,思路是利用reads重新组装,往往能改善边界。或者用长读长技术做bin polishing,但成本高,前期尝试不建议。

5.5 GTDB-Tk的数据库下载和运行速度

GTDB-Tk需要下载庞大的数据库(几十GB),而且局域网环境容易中断。建议使用官方提供的wget方式,并设置好GTDBTK_DATA_PATH环境变量。如果只是快速看分类,可以先使用基于MASH的快速模式fastani,但最终发表时还是建议完整classify_wf。

运行慢时优先加--pplacer_cpus和--cpus,如果机器内存不够,考虑增加swap或用高内存节点。

6. 实战经验与个人心得

这一路踩下来,我最想对刚接触宏基因组Binning的人说一句:不要把流程跑通当成终点,要把每个bin的质量指标当成分数线。你最后投稿时的MAG统计表,每一个数字都必须能追溯回参数和版本号。最好从一开始就给每个过程产物记录参数:fastp版本、MEGAHIT k-mer范围、MetaBAT2的minContig阈值、CheckM数据库版本,GTDB-Tk数据库版本。这个习惯能让你在修改流程、复现实验时省下大量时间。

另外,我现在的标准流程里会保留一条“单样本组装 + 多样本共组装”并行的路线。单样本组装结果虽然N50可能低一些,但物种组成更真实,不易制造出样本间不存在的嵌合体。多样本共组装则对低丰度物种的回收有天然优势。两条路线同时跑,最后合并dRep,MAG数量和质量都比只走一条线高不少。

“binning ic”这个关键词,我之前也查过,表面上是一个浓缩术语,本质上说的是同一件事——你喂给分箱工具的信号够不够、有没有混杂;组装质量、覆盖度、序列组成特征、样本数量,这些信息维度组合起来,决定了你从数据里能回收多少高完整性、低污染率的MAG。所以不要机械地照搬某个大神的固定参数,理解每个参数影响的是哪个信号维度,才能真正把这个流程用活。

最后分享一个我常用的小技巧:跑完dRep后,我会用seqkit统计一下每个MAG的染色体保守基因(比如DNA旋转酶、RNA聚合酶亚基)是否存在,如果关键基因缺失,即使CheckM完整性看着可以,我也会对这个MAG打一个问号。这种“关键证据链检查”虽然冷门,但能拦住不少后续分析的返工。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询