宏基因组测序拿到手绝大多数人第一反应是直接上组装、跑分箱然后被一大堆bin文件和五花八门的后续分析工具搞到怀疑人生。做这个项目时我踩了不少坑最后把整条流程完整走了一遍从质控、组装、分箱到bin的质量评估、物种注释和功能分析每一步的参数、命令、判断标准都摸透了。这篇文章就直接按实操顺序往下写把每个环节为什么要这么做、有哪些容易被忽略的细节、碰到问题怎么排查都讲清楚目标是让一个刚接触宏基因组的人也能照着把bin做出来并且知道下一步该分析什么、怎么分析。1. 从测序数据到组装分箱前你到底要做哪些准备1.1 为什么质控省略不得fastp一把梭但别忽略宿主污染宏基因组数据到手第一步不是组装而是质控。原始下机数据里会有接头、低质量碱基还可能混着宿主DNA这些杂质直接影响后面的组装和分箱效果。我常用的质控工具是fastp一条命令同时完成过滤和质检报告fastp -i raw_R1.fastq.gz -I raw_R2.fastq.gz \ -o clean_R1.fastq.gz -O clean_R2.fastq.gz \ --detect_adapter_for_pe --thread 16 \ -j fastp.json -h fastp.html参数解释一下--detect_adapter_for_pe是让fastp自动检测并切除双端测序的接头比手动指定接头序列省事而且对常见的Illumina接头都有效。默认的过滤规则是滑动窗口剪裁质量值低于Q15的碱基这个阈值对宏基因组来说够用如果想更严格一点可以追加--cut_front --cut_tail把reads两端质量差的碱基直接切掉。如果是肠道、皮肤这类容易混入宿主细胞的样本强烈建议在fastp之后再做一层宿主序列过滤。常用的方案是KneadData本质是把clean reads比对到宿主参考基因组上把能比对上的reads丢掉kneaddata -i clean_R1.fastq.gz -i clean_R2.fastq.gz \ -o kneaddata_out \ --trimmomatic-options SLIDINGWINDOW:4:20 MINLEN:50 \ --reference-db human_genome_index -p 8我见过不少人跳过这一步结果组装出的contig里混着大量人类序列不仅浪费计算资源分箱时还会出现高污染的假bin。所以这里值得多花点时间。1.2 组装MEGAHIT与metaSPAdes怎么选分箱的前提是先把reads拼成contig。宏基因组组装有两个绝对绕不开的工具MEGAHIT和metaSPAdes。MEGAHIT最大的优势是内存占用小、速度快适合样本量大或服务器配置一般的场景。它的本质是基于de Bruijn图迭代k-mer从短k-mer开始逐步延长对覆盖度不均的宏基因组数据适应性很好。我通常在拿到数据后先用MEGAHIT跑一版看看组装效果megahit -1 clean_R1.fastq.gz -2 clean_R2.fastq.gz \ -o megahit_out -t 16 --min-contig-len 1000--min-contig-len建议设置为1000也就是只保留长度1000bp以上的contig。这个参数非常关键因为分箱算法几乎都不喜欢太短的contig——短contig上的覆盖度估计不稳定四核苷酸频率等特征也计算不出来留太多反而干扰后续分箱。metaSPAdes的组装质量通常比MEGAHIT更好尤其是对丰度较低的物种能把那些边缘区域的基因组拼得更完整。代价是内存占用高实测下来跑一个中等复杂度样本100G内存都不一定够。如果服务器内存充足我建议直接用metaSPAdesspades.py --meta -1 clean_R1.fastq.gz -2 clean_R2.fastq.gz \ -o metaspades_out -t 16 -m 500-m 500表示给spades最多500G内存。如果跑的过程中内存不够了程序会直接报错退出所以在参数上别抠搜给够才有意义。1.3 样品数、深度与组装效果的真实关系组装效果不只是工具决定的更多由数据本身决定。宏基因组组装依赖覆盖度某个物种的基因组如果被测序reads覆盖不到一定深度再强的工具也拼不出来。这里有一个经验值平均覆盖度低于5x的物种基本不可能从短读组装里得到完整的基因组。样本数量也是个关键因素。单个样本深度不够可以尝试共组装也就是把多个同类型样本的reads合并在一起组装这样多个样本对同一物种的覆盖度会叠加contig能拼得更长。但共组装也有代价——那些只在个别样本中出现的菌株由于平均覆盖度被稀释反而可能拼不出来。所以在共组装和单样本组装之间需要有判断先做一两个样本试探对比N50和contig总数再决定全量样本怎么处理。组装结束后先用seqkit统计一下结果seqkit stats final.contigs.fa重点看N50和总长度。如果N50只有几千bp说明样本复杂度很高或者深度不够如果N50有几十k甚至上百k说明有不少高丰度物种被拼出来了。组装这一步没有绝对的好坏标准但太差的组装结果会直接导致后面一无所有。2. 分箱主流程MetaBAT2、MaxBin2、CONCOCT怎么跑2.1 先给contig算深度与覆盖度比对这一步决定了分箱效果分箱的基本依据是同一个基因组来源的contig在多个样本间的覆盖度变化模式应当是相似的。所以跑分箱之前必须把所有样本的clean reads比对回组装得到的contig算出每个contig在样本中的覆盖度。标准做法是用bwasamtools组合bwa index final.contigs.fa bwa mem -t 16 final.contigs.fa clean_R1.fastq.gz clean_R2.fastq.gz aln.sam samtools view -bS aln.sam aln.bam samtools sort - 16 -o aln.sorted.bam aln.bam samtools index aln.sorted.bam比对完成后用jgi_summarize_bam_contig_depths统计深度信息这是MetaBAT2作者开发的配套工具可以直接生成MetaBAT2需要的输入文件jgi_summarize_bam_contig_depths --outputDepth depth.txt aln.sorted.bam生成的depth.txt里包含contig名称、总长度、平均深度、变异系数等列。如果你有多个样本就把所有样本的bam文件一起传给jgi_summarize_bam_contig_depths每个样本会单独占一列。这一步是分箱质量的分水岭——深度矩阵没构建好后面所有工具都白搭。2.2 三个分箱工具的核心原理差异MetaBAT2、MaxBin2、CONCOCT是现在最常用的三个分箱工具但它们的核心思路完全不同跑出来的结果也各有侧重。MetaBAT2用的是“序列组成覆盖度”的组合特征通过图聚类算法把具有相似丰度模式和相似四核苷酸频率的contig归到一起。优点是速度快、对高质量数据分箱完整度高目前是绝大多数流程的默认首选。metabat2 -i final.contigs.fa -a depth.txt -o bins/bin -m 1500 -t 12-m 1500是最小contig长度阈值只使用1500bp以上的contig参与分箱。这个数值来自经验1500bp以下contig的组成特征不稳定参与分箱容易造成碎片化。MaxBin2的思路则是利用单拷贝标记基因。它先用Prodigal等工具预测contig上的基因然后寻找一组在细菌基因组中通常以单拷贝形式存在的标记基因通过EM算法估计出可能的基因组数量再把contig分配到各个基因组“箱”中run_MaxBin.pl -contig final.contigs.fa -out bins/maxbin \ -abund_list sample.abund.list -thread 12sample.abund.list需要你自己写每一行是一个样本的abundance文件路径也就是depth.txt里对应那一列的数据。MaxBin2的优势是保守分出的bin不容易混入杂菌但代价是很多低丰度物种根本不会被分出来。CONCOCT则完全不同它依靠的是协方差模式。CONCOCT把contig切成固定长度的小片段通常10kb基于片段在不同样本间的覆盖度协方差用高斯混合模型聚类再把聚类结果映射回原始contig。这要求你提供多个样本样本太少时效果会明显变差。MetaWRAP可以把三个工具一键跑完并做整合命令更简单适合赶项目进度时使用metawrap binning -a final.contigs.fa -o metawrap_binning_out \ -t 12 --metabat2 --maxbin2 --concoct2.3 参数选择minContig、minBinSize、线程和内存的取舍分箱过程中的参数调整是真正区分新手和老手的地方。minContig长度我一般默认1500bp但如果你研究的是低深度样本比如病毒宏基因组或超低丰度菌群可以降到1000甚至500因为长contig本来就少再卡1500就没有东西可分箱了。相反如果样本深度很高、组装质量好提高到2500或3000反而能减少短contig造成的噪声。minBinSize同理。MetaBAT2默认的bin最小大小是200000bp200kb这个数值对应细菌基因组中位数约3-4Mb来说非常宽松主要为了防止碎片化小bin混入。MaxBin2不需要特别设置它会根据标记基因自动估计bin数量。线程和内存方面MetaBAT2和MaxBin2都是多线程友好的但CONCOCT比较吃内存尤其是拼接了大量10kb片段后内存峰值可能达到几十G。实操上我建议样本数少10个以内用单样本深度列构建矩阵MinContig设1500样本数多几十个共组装后统一比对深度矩阵包含全部样本服务器内存不充足时只跑MetaBAT2和MaxBin2不跑CONCOCT分箱跑完看看bin目录里有多少个文件夹或fasta文件。如果所有工具跑出来的bin数量差异特别大通常意味着数据里存在大量低丰度菌群或者某个样本存在异常深度分布。3. 把三个结果揉到一起bin_refinement与DAS Tool3.1 为什么一定要做结果整合按我的经验三个分箱工具的结果直接取“最优”是错误思路。MetaBAT2倾向于把亲缘关系近的菌合并到一个bin里MaxBin2则可能把一个基因组拆成好几个碎片CONCOCT容易把丰度模式相似的菌混在一起。单独依赖任何一个工具都会引入特征性的错误。解决思路就是整合把多个工具的结果放到一起取并集或做择优然后用单拷贝基因评估每个候选bin的质量找出最优组合。DAS Tool就是这方面的标杆工具MetaWRAP的bin_refinement则是整合了CheckM的简化封装。DAS Tool的命令DAS_Tool -i sample_list.txt -l sample_labels.txt \ -c final.contigs.fa -o das_out \ --search_engine diamond --write_binssample_list.txt里写三个工具输出bin的目录路径sample_labels.txt写对应的标签DAS_Tool会综合所有bin的组成特征和覆盖度输出一套去冗余的最佳结果。如果只想快速得到一个干净的结果集用MetaWRAP更省事metawrap bin_refinement -o metawrap_refine -t 12 \ -A metawrap_binning_out/metabat2_bins/ \ -B metawrap_binning_out/maxbin2_bins/ \ -C metawrap_binning_out/concoct_bins/ \ -c 50 -x 10-c 50表示bin完整度不低于50%-x 10表示污染度不高于10%。这两个阈值是MetaWRAP作者推荐的默认值我当时直接用没改后面用CheckM复核结果时发现这套标准确实能筛掉大部分低质量bin。3.2 如何用CheckM对整合后的bin做质量评判不管是直接用单个工具的结果还是整合之后的bin集合都需要用CheckM做客观的质量评估。CheckM的原理是看bin里的单拷贝标记基因如果某个bin完整度高那么一大组核心单拷贝基因都应该存在污染度高则说明有多个物种的基因组混在一起导致多个“单拷贝”基因出现多拷贝的情况。checkm lineage_wf -t 12 -x fa bins checkm_out checkm qa checkm_out/lineage.ms -o 2 -t 12 --tab_table \ -f checkm_results.tsv checkm_out第一条命令会对每个bin进行谱系特异性分析并计算完整度和污染度第二条命令把详细结果输出成表格方便后面用Excel或R筛选。如果你嫌CheckM跑得慢可以考虑用CheckM2它是基于机器学习的新版本预测速度更快而且不需要判断谱系checkm2 database --download checkm2 predict --threads 12 --input bins --output-directory checkm2_out关于bin质量的判定标准学术界有比较统一的认识来自MIMAGMinimum Information about a Metagenome-Assembled Genome建议质量等级完整度污染度用途高质量≥90%≤5%物种新分类单元描述、比较基因组分析中等质量≥50%≤10%多样性分析、功能注释的粗筛低质量50%不做硬要求尽量不用于下游分析我当时筛选bin的标准是完整度≥70%、污染度≤10%这个相对宽松的阈值既保证了数量也保证了后续物种注释和功能注释的可靠性。如果你要发文章或者做新种描述建议直接按MIMAG高质量标准来。3.3 bin数量和多度分布怎么解读跑完分箱和质控你会得到几十甚至上百个bin。bin数量的多少本身没有好坏关键看与样本的生物学背景是否吻合。比如一个复杂的环境样品土壤、海洋几十个bin很正常如果是一个简单的富集培养物却分出了几十个bin大概率是分箱参数太松或者共组装时混入了其他样本的低丰度菌群。这时候建议按bin的完整度排序用seqkit stats逐个看一下bin的总长。正常细菌bin大小应该在1-6Mb之间那些不到500kb的小bin多半是碎片可以直接舍弃或标记为“需进一步处理”。4. 后续分析第一关bin的分类学注释GTDB-Tk、Kraken24.1 GTDB-Tk是怎么工作的拿到干净的bin集合第一个目标就是搞清楚它们都是什么物种。现在最权威、我最推荐的工具是GTDB-Tk它基于基因组分类数据库Genome Taxonomy Database对bin进行物种注释。GTDB-Tk的原理不算复杂把bin送到GTDB数据库中找相似基因组通过一系列串联单拷贝蛋白基因建立进化树再依据这个树的拓扑位置确定物种分类。实际操作命令如下gtdbtk classify_wf -g bins --cpus 12 --out_dir gtdb_out -x fa --prefix gtdb注意使用GTDB-Tk之前需要先下载GTDB数据库这一步相当耗时——数据库解压后有几十上百个G而且更新频繁。下载和数据库配置都不是难点关键是别把数据库目录搞乱。跑完后会在输出目录里得到两个文件gtdb.bac123.classify.tree细菌的进化树和gtdb.bac123.summary.tsv分类结果表。summary表格里每一行对应一个bin有从界到种的完整分类层级还会给出这个bin在数据库中最相近的基因组及平均核苷酸一致性ANI。4.2 其他分类工具与GTDB-Tk的区别除了GTDB-Tk不少人会用Kraken2对bin做物种注释。Kraken2本质上是把bin的contig碎片化后与参考数据库做k-mer匹配速度快但它更适合对测序reads直接做物种组成分析而不是对已分好的bin做精细分类。如果你已经用CheckM确认了bin的基因组完整性用GTDB-Tk而不是Kraken2是更合理的选择。还有个常用工具是CAT/BAT它利用蛋白质序列比对来推断bin的分类在未知物种占比较高时表现不错。不过GTDB-Tk有天然的优势GTDB数据库本身对来自宏基因组的未培养物种做了大量整理有几个新门类的命名都来源于宏基因组组装基因组所以分类结果更可靠。如果只是快速确认bin的大致门类可以先跑一个Kraken2kraken2 --db kraken2_db --threads 12 bin.fa kraken_result.txt但是做精细层面的物种注释我还是推荐GTDB-Tk它的结果经得起推敲。4.3 分类结果怎么看GTDB-Tk输出的分类结果是一个标准化的分类学谱系类似d__Bacteria;p__Bacillota;c__Clostridia;o__Oscillospirales;f__Ruminococcaceae;g__Faecalibacterium;s__Faecalibacterium prausnitzii每一个bin如果注释到了属或者种就可以直接引用。但要注意一个很常见的坑GTDB的分类框架与NCBI传统分类基于16S有一定差异同一个物种在GTDB里的属名可能与NCBI不同。比如肠道里大名鼎鼎的Eubacterium rectale在GTDB里被归到了Agathobacter属。如果你后续要和文献里的老名字对照建议在结果里同时保留ncbi分类信息避免后期解释时混乱。5. 后续分析第二关功能基因与代谢通路挖掘5.1 Prokka到eggNOG-mapper完整功能注释路径物种注释解决的是“它是谁”功能注释解决的是“它能干什么”。功能注释的标准路径是先预测bin里的基因再把预测蛋白与数据库比对。基因预测我用Prokka它整合了Prodigal、Aragorn等多个工具输出gff、gbk、蛋白序列等文件prokka --kingdom Bacteria --cpus 8 --outdir prokka_out bin.fa注意--kingdom参数宏基因组bin一般选Bacteria如果你确定某个bin是真菌来源要改成Fungi。如果注释结果出现大量hypothetical protein可以考虑加--evalue 1e-6之类的更严格阈值但不用过分在意宏基因组bin中有大量未知基因是正常的。蛋白序列注释推荐eggNOG-mapper它会把蛋白序列比对到eggNOG数据库一次性得到COG功能分类、KEGG通路、GO注释等结果conda install -c bioconda eggnog-mapper download_eggnog_data.py -y emapper.py -i prokka_out/bin.faa -o bin_eggnog --cpu 12 -m diamond这里-m diamond是用diamond做序列比对速度比blastp快几个数量级效果基本一致。输出文件里的COG_category列可以直接拿来统计功能类别数量KEGG_Pathway列可以提取通路信息这些是后续画柱状图、热图的基础。5.2 CAZyme与antiSMASH等专项挖掘除了通用的COG和KEGG很多项目会关注特定的功能类别。比如研究肠道菌群时碳水化合物的利用能力很重要这时候用dbCAN注释碳水化合物活性酶CAZymesrun_dbcan bin.fa protein --out_dir dbcan_out --use_server--use_server表示使用在线HMMER服务器如果网络条件不稳定可以改成本地数据库运行。dbCAN的输出会给出每个蛋白对应的CAZy家族比如GH糖苷水解酶、GT糖基转移酶、PL多糖裂解酶等可以直接统计各家族数量做样本间比较。次级代谢产物合成基因簇的挖掘也很常见用的工具是antiSMASHantismash prokka_out/bin.gbk --output-dir antismash_outantiSMASH需要输入gbk格式的注释文件它会在基因组里寻找聚酮合成酶PKS、非核糖体肽合成酶NRPS等次级代谢基因簇。运行时间主要取决于基因组大小和基因簇复杂度通常一个bin需要几分钟到几十分钟。5.3 样本间差异分析怎么衔接功能注释做完之后最常做的事情是比较不同样本组之间的功能丰度差异。这里要理解一个关键点bin层面功能注释的结果是“存在/不存在”不是丰度。如果要做丰度差异分析需要把reads比对回bin计算每个bin在样本中的相对丰度samtools index bin.bam # 用coverm或者bowtie2featureCounts组合 coverm genome --bam-files aln.sorted.bam --genome-files bins/*.fa --methods relative_abundance coverm_out.tsvCoverM计算相对丰度时默认考虑基因组长度能够得到一个近似于16S rRNA基因相对丰度的结果。有了每个bin在各样本中的丰度表再结合功能注释就能做“某个物种丰度升高、其携带某功能基因的bin也对应升高”之类的结论。这是宏基因组分析中非常有说服力的证据链。6. 用dRep和系统发育树整理你的bin集合6.1 dRep去冗余如果项目涉及几十个样本、上百个bin你会发现不同样本中可能存在同一物种的近缘菌株它们的基因组高度相似直接放在一起分析会造成冗余干扰。dRep就是专门解决这个问题的dRep dereplicate dRep_out -g bins/*.fa --S_algorithm ANImf \ -sa 0.95 -nc 0.30 -p 12-sa 0.95表示平均核苷酸一致性ANI在95%以上的bin会被视为同一个物种并合并-nc 0.30表示若两者的比对覆盖度低于30%则不合并。dRep会在输出目录里生成每个rep bin以及一张包含相似度聚类关系的表格。这是后续构建“非冗余的基因组集”的标准操作。我在实际项目中遇到过这样的情况同一个样本里分出了两个ANI高达99.8%的bin检查后发现它们其实是同一个菌株的两个contig碎片dRep把它们合并之后下游注释结果反而更完整、更干净。所以去冗余这一步不是可有可无。6.2 建树和可视化有了代表性的bin集合很多文章里都需要一张系统发育树来展示这些bin与其他已知物种的关系。GTDB-Tk分类学注释过程中本身就会生成一棵进化树直接用那棵树做可视化是最省事的路径。另外也可以用IQ-TREE基于串联的单拷贝蛋白序列建一棵最大似然树iqtree -s concatenated_alignment.fasta -m MFP -bb 1000 \ -nt AUTO --prefix iqtree_result-m MFP让IQ-TREE自动选择最优替代模型-bb 1000做1000次超快自举用于评估分支支持度。建树完成后用iTOL在线工具上传tree文件可以很方便地调整外观并导出出版级图片。7. 实操中踩过的坑与排查建议7.1 分箱数量和完整度不对优先查这几项分箱结果不理想时大多数人第一反应是换分箱工具但更常见的根源在前期。我会按固定顺序排查现象优先排查方向处理建议bin数量过少测序深度不足或过滤太严检查N50和clean reads保留率适当调低minContigbin污染度高10%共组装样本复杂度太高尝试降低minContig阈值或改用MaxBin2独立跑一遍再整合bin完整度低低丰度菌群基因组碎片化尝试共组装增加覆盖度或做深度binning单个bin混入近缘菌近缘菌株覆盖度模式高度相似提高覆盖度矩阵样本维度换更严格的整合阈值深度矩阵异常比对步骤存在错误检查bam文件质量和比对比率确保reads比对上组装contig7.2 环境配置与版本管理经验最后一定要提环境管理。生物信息项目的工具依赖非常琐碎MetaBAT2、CheckM、GTDB-Tk用的conda环境经常互相冲突。我的做法是用conda建独立环境每个环境对应一个流程阶段并且把软件版本固定下来不要随手升级conda create -n binning_env -c bioconda metabat2 maxbin2 concoct conda create -n checkm_env -c bioconda checkm conda create -n gtdb_env -c bioconda gtdbtk这样做的好处是以后任何时候回看项目都能知道某个结果是在哪个版本组合下产出的。宏基因组分析的结果可复现性比想象中重要尤其是合作项目或审稿人要求提供分析流程细节的时候。我个人在实际操作中最深的体会是分箱没有“一键成功”的银弹结果质量的好坏更多取决于你前期对数据的理解包括深度、样本复杂度、组装效果。宁可花时间把组装和深度矩阵做好也别在分箱参数上反复猜测。另一个建议是保留中间文件尤其是depth.txt和组装结果如果后面有新的分箱工具发布可以直接复用这些文件做重新分箱而不用从头再跑一遍比对和组装。