学习技能(点击标题跳转)
01
02
03
04
05
06
07
在前面的文章我们用Claude code辅助完成了宏基因组测序数据的质控和去宿主Claude code做宏基因组数据质控、手把手教你指挥KneadData精准“斩断”宏基因组宿主干扰,接下来就对宏基因组数据进行组装。在宏基因组学研究中,将数以百万计的短序列(Reads)拼接成更长的连续序列(Contigs),是后续物种鉴定、功能分析和基因组分箱(Binning)的关键第一步。
本教程将涵盖三大核心环节:
为了避免软件版本和依赖冲突,推荐使用 Conda 来管理生信软件。首先,创建一个独立的分析环境并安装所有需要的工具:
# 创建名为 "metagenome" 的新环境
conda create -n metagenome
# 激活环境
conda activate metagenome
# 一次性安装所有需要的软件
# -c bioconda 和 -c conda-forge 是指定的软件源
conda install -c bioconda megahit bbmap bwa samtools metabat2 -y
安装完成后,就可以在这个 metagenome
环境中开始分析了。
假设已经完成了质控,手头有高质量的双端测序数据:S1_clean_R1.fq.gz
和 S1_clean_R2.fq.gz
。
MEGAHIT 是一款专为宏基因组设计的高效组装软件。它基于简洁的德布莱英图(Succinct de Bruijn Graph)算法,具有速度快、内存占用低的优点,尤其擅长处理高复杂度的宏基因组数据。
MEGAHIT 的命令非常直观。我们将使用多线程来加速组装,并设置最小 Contig 长度为 1000 bp,这是宏基因组分析的常用标准,可以有效过滤掉过短、无意义的片段。
# 运行 MEGAHIT 组装
megahit \
-1 S1_clean_R1.fq.gz \
-2 S1_clean_R2.fq.gz \
-o S1_megahit_out \
-t 24 \
--min-contig-len 1000
# 参数解释:
# -1, -2: 分别指定双端测序的 R1 和 R2 文件。
# -o: 指定输出结果的目录名。
# -t: 指定使用的线程数,根据你的服务器配置调整。
# --min-contig-len: 设置输出 Contig 的最小长度阈值(单位:bp)。
运行结束后,最重要的输出文件是 S1_megahit_out/final.contigs.fa
,这就是我们得到的组装好的 Contigs 文件,后续所有分析都将基于它。
组装完成后,我们需要评估其质量。BBMap 套件中的 stats.sh
是一个轻量级但功能强大的工具,可以快速计算 N50、Contig 数量、总长度等关键指标。
stats.sh
默认输出人类易读的格式,但通过 format=5
参数,我们可以得到一个制表符分隔的、机器易读的格式,方便后续脚本处理。
# 运行 stats.sh 进行统计
stats.sh in=S1_megahit_out/final.contigs.fa format=5 > S1_assembly_stats.txt
# 参数解释:
# in=: 指定输入的 FASTA 文件。
# format=5: 输出为5列格式(n_scaffolds, n_contigs, scaf_bp, contig_bp, gap_pct)。
# >: 将标准输出重定向到文件 S1_assembly_stats.txt。
打开 S1_assembly_stats.txt
文件,就会看到类似下面的内容,其中 N50 是衡量组装连续性的核心指标,值越高通常意味着组装效果越好。
#n_scaffolds n_contigs scaf_bp contig_bp gap_pct scaf_N50 scaf_L50 ...
150000 150000 350000000 350000000 0.00 25000 28000 ...
宏基因组分箱(Binning)的目的是将属于同一个物种的 Contigs 聚类成一个基因组草图(MAG)。分箱算法(如 MetaBAT2)依赖两个核心信息:序列组成特征(如四核苷酸频率)和覆盖度信息。
来自同一个物种的 Contigs,其在样本中的丰度应该相似,因此它们的测序深度(覆盖度)也应该相近。此步骤就是为了计算每个 Contig 的平均覆盖深度。
这个过程是一个经典的比对流程:
这是一个三步走的流程,我们通常用管道符 |
将它们连接起来,以提高效率。
# BWA 需要先对参考序列(这里是我们的 Contigs)建立索引
bwa index S1_megahit_out/final.contigs.fa
# 定义文件名变量,方便后续使用
CONTIGS="S1_megahit_out/final.contigs.fa"
BAM_OUT="S1_mapping.sorted.bam"
DEPTH_OUT="S1_contig_depth.txt"
THREADS=24
# 核心命令:比对 -> 转 BAM -> 排序
bwa mem -t ${THREADS} ${CONTIGS} S1_clean_R1.fq.gz S1_clean_R2.fq.gz | \
samtools view -@ ${THREADS} -bS -F 4 - | \
samtools sort -@ ${THREADS} -o ${BAM_OUT} -
# 参数解释:
# bwa mem -t: 使用多线程将 reads 比对到 contigs。
# |: 管道符,将上一步的输出作为下一步的输入。
# samtools view -bS -F 4: 将 SAM 格式转为 BAM,-F 4 过滤掉未比对上的 reads。
# samtools sort -o: 对 BAM 文件按坐标排序,并输出到指定文件。
# 3. 从排序后的 BAM 文件中计算覆盖度
# jgi_summarize_bam_contig_depths 是 MetaBAT2 自带的脚本
jgi_summarize_bam_contig_depths --outputDepth ${DEPTH_OUT} ${BAM_OUT}
# 参数解释:
# --outputDepth: 指定输出的深度文件。
# 最后一个参数是输入的 sorted BAM 文件。
运行结束后,将得到一个名为 S1_contig_depth.txt
的文件,内容格式如下:
contigName contigLen totalAvgDepth S1_mapping.sorted.bam
contig_1 54321 15.78 15.78
contig_2 43210 89.12 89.12
...
这个文件包含了每个 Contig 的名称、长度和平均覆盖深度,是进行分箱分析的关键输入之一。
至此,你已经完成了宏基因组数据从原始 Reads 到结构化 Contigs 的核心处理流程。
S1_megahit_out/final.contigs.fa
(组装好的 Contigs)S1_contig_depth.txt
(Contigs 覆盖度文件)有了这两个文件,就可以进入下游分析阶段,例如使用 MetaBAT2
、MaxBin2
等工具进行基因组分箱 (Binning),从而重构出样本中微生物的基因组草图(MAGs),进而探索它们的物种身份和功能潜力。
我们深知,科研的宝贵时间不应浪费在环境配置的反复试错与计算任务的漫长等待上。为此,我们推出的高性能计算服务器,正是为解决这些痛点而生。我们不仅为您准备了拥有大内存、多核心的强劲硬件,更将通用分析工具及所需数据库进行了预装和深度优化,为您打造一个“开箱即用”的宏基因组分析平台。让您告别繁琐配置,专注数据,加速您的科研进程。
天意云服务器产品类型包括集成分析环境、共享服务器、独享服务器三种。服务器配置科学划分,按需部署。考虑到项目需求较小、周期较短的用户,我们基于共享服务器搭建了集成分析环境,通过浏览器即可登录使用Rstudio和Jupyter。服务器产品都提前预装了常用的软件、R包、Python库,无论你是做宏基因、宏病毒还是做单细胞组学、空间转录组,都有适配的软件供你使用,帮你节省很多安装软件的时间。最重要的,选择我们的服务器就相当于选择了一个专业的辅助团队,在你使用服务器期间,我们技术人员免费提供技术支持服务。