基于conda环境下的宏基因组学分析利器MetaWRAP 1.3.2 安装和使用,序列分析基本流程自动分析脚本
以下是基于上述内容整理的逐步指导说明:
MetaWRAP宏基因组分析自动化脚本指南
目标
完成一个基于MetaWRAP的宏基因组分析自动化脚本(MTWRP_pipeline.sh),用于对单个或多个样本进行完整的分析流程。所需工具
MetaWRAP工具库
BLAST数据库
CarbonCONDOR服务
HUMAnN数据库
MetaPhlAn4数据库数据准备
输入数据:
- 必要输入:包含不同样品对的fastq文件(如sample1forward.fastq, sample1reverse.fastq, sample2forward.fastq, sample2reverse.fastq)。
- 保证fastq文件路径正确,并且文件名格式符合预期(如{i}1.fastq和{i}2.fastq)。
输出目录:- 创建工作目录并设置为输出目录(如/path/to/working_dir/)。
- 使用管道将输出结果写入该目录下。
自动化脚本结构
`bash
#!/bin/bash
设置工作目录和输出路径
WORKDIR=WORKDIR INPUT_FASTQ PathToInputFastQ/ 创建结果文件夹 mkdir -p WORKDIR/readqc/i/proj_metawrap/i/metagenome/i/abundance/ mkdir -p WORKDIR/readqc/i/proj_metawrap/i/b bubble/
阶段性清理旧文件(可选)
find WORKDIR(readqc) -type f -name "*.fastq" | xargs rm \ && find WORKDIR/readqc -type d ! -name "*.log" | xargs rm
&& find WORKDIR/projmetawrap ! -name "*.log" | xargs rm echo "开始质量控制..." >> readqcoutput.log 质量控制模块 for i in (ls INPUT_FASTQ); do metawrap read_qc \ -t 160 \ --skip empty_p pairing \ --skip quality_file skip QC \ --skip output quality report \ --outputdir WORKDIR/read_qc/{i} \ {INPUTFASTQ}/{i}/forward*.fastq \ {INPUTFASTQ}/{i}/reverse*.fastq >> readqcoutput.log & done echo "质量控制完成" >> readqcoutput.log 组装模块 echo "开始组装..." >> assemble_output.log for i in (ls INPUT_FASTQ); do metawrap assembly \ --method Bowtie2 \ --k 500 \ --odir WORKDIR/assemble/{i} \ {INPUTFASTQ}/{i}/forward.fastq {INPUTFASTQ}/{i}/reverse.fastq >> assemble_output.log & done echo "组装完成" >> assemble_output.log Kraken2分类模块 echo "开始Kraken2分类..." >> kraken2_output.log for i in (ls $INPUT_FASTQ);
do
metawrap kraken2
--t 4000000000kMv3.5e9.txt
--k 3219999999.76e3MmzN5mGmzCNZmLwmtIv4LMZBNfsgFkIIVoD8onccrlFxmGwlVXNmd
介绍:

MetaWRAP是一种基于宏基因组学的综合性分析平台,专为解析宏基因组测序数据而设计。该平台整合了一套功能全面的数据处理系统,能够完成宏基因组数据的拼接、注释以及各类功能特性分析任务。
MetaWRAP的功能包括:
数据质量控制:涵盖剔除含有污染的reads、过滤出高质量的样本以及去除非编码区片段等。
基因组构建:该系统MetaWRAP提供了多种基因组构建方法和工具包选项。它不仅包含基于短读长的SPAdes以及基于长读长的MEGAHIT等多种算法选择,在实际应用中可以根据用户的特定需求灵活配置不同的构建方案以满足研究目标的实现
MetaWRAP具备基因组注释功能,在功能上涵盖基因预测、功能注释以及通路预测等多种服务。该工具可利用多种数据库如KEGG、COG及NR来进行注释。
基因组间的对比分析:MetaWRAP系统支持多基因组间的对比分析,并提供物种组成分析功能。该系统能够帮助研究者深入解析不同样本间的异同点。
生态位分析:MetaWRAP整合了一系列生态位分析工具,并且这些工具能够帮助用户更好地掌握样本中微生物的功能特性和代谢途径。
MetaWRAP的优势包括:
多样化的功能:该系统集成了多种分析工具,并支持从原始测序数据到生物信息学分析的完整工作流程。
多样化的操作界面:MetaWRAP提供命令行界面和Python Application Programming Interface(API)两种主要接口供用户选择;允许用户提供不同的偏好以选择最适合自身需求的解决方案。
快速实现性能:MetaWRAP通过多线程架构和并行计算技术显著提升了分析效率。
该工具库不仅具备强大的功能优势,在宏基因组数据处理方面也展现出卓越的应用价值。它通过辅助构建、注释与功能鉴定等技术, 为用户提供全面的数据分析解决方案。应用具有高度灵活性, 同时其计算效能显著提升, 完全能够满足多种研究场景的需求。
在了解相关文章之前,请先熟悉以下内容:[MetaWRAP——支持全面、细致的基因组级代谢组学数据分析的灵活分析管道 | Microbiome | 全文
github目录:https://github.com/ursky/metaWRAP
Anaconda登录页面:访问Anaconda
安装
这是教大家如何进行conda或mamba的安装过程吧?其他版本可能会有更新问题,在配置过程中有时会遇到一些挑战
conda config --add channels defaults
conda config --add channels conda-forge
conda config --add channels bioconda
conda config --add channels ursky
#
mamba config --add channels defaults
mamba config --add channels conda-forge
mamba config --add channels bioconda
mamba config --add channels ursky
安装的时候注意了,直接将所有channel都放上,不然缺包错误:
mamba create -y --name metawrap132 -c ursky -c bioconda -c conda-forge metawrap-mg=1.3.2
安装完最后的提示:

查看安装结果,主要看其中metawrap-mg的版本,这里是1.3.2,来自ursky:
mamba list

配置建议数据库
这里包括各种类型的数据库系统架构设计中涉及的数据量大小问题,并根据具体业务需求对各模块可能涉及的数据库进行相应的规划与配置;如果缺少相应的数据库配置项,则应在后续开发过程中针对具体应用场景进行补充;若当前模块未设置相关参数,则默认跳过该参数设置并继续执行后续代码逻辑

taxonomy数据库:
# 先删除原配置文件夹
rm -rf /miniconda3/envs/metawrap132/opt/krona/taxonomy
# 自己创建指定文件夹
mkdir /path/on/big/disk/taxonomy
# 创建软链接
ln -s /path/on/big/disk/taxonomy /miniconda3/envs/metawrap132/opt/krona/taxonomy
# 自动下载更新数据库,会自动下载到自己指定的文件夹
ktUpdateTaxonomy.sh
直接查看ktUpdateTaxonomy.sh文件内容,直接下载来自ncbi数据库:


下载解压完成后,在目标目录生成一个taxonomy.tab的文件:
head taxonomy.tab
1 0 1 no rank root
2 2 131567 superkingdom Bacteria
6 7 335928 genus Azorhizobium
7 8 6 species Azorhizobium caulinodans
9 8 32199 species Buchnera aphidicola
10 7 1706371 genus Cellvibrio
11 9 1707 species Cellulomonas gilvus
13 7 203488 genus Dictyoglomus
14 8 13 species Dictyoglomus thermophilum
16 7 32011 genus Methylophilus
GRIDSS\SILVA 16S rRNA\BUSCO数据库
quast-download-gridss
quast-download-silva
quast-download-busco
下载后位于目录:
/miniconda3/envs/metawrap/lib/python2.7/site-packages/quast_libs/
下载日志:
envs/metawrap/lib/python2.7/site-packages/quast_libs/silva/blastdb.log
中间几个数据库下载不下来,更换链接吧,应该是版本变了,地址不对:
busco 的数据目录存储于/v5/data/线系表中
busco的官网:BUSCO - from QC to gene prediction and phylogenomics
需要下载的文件:
https://busco-data.ezlab.org/v5/data/lineages/fungi_odb10.2021-06-28.tar.gz
https://busco-data.ezlab.org/v5/data/lineages/eukaryota_odb10.2020-09-10.tar.gz
https://busco-data.ezlab.org/v5/data/lineages/bacteria_odb10.2020-03-06.tar.gz
cd miniconda3/envs/metawrap132/lib/python2.7/site-packages/quast_libs/busco/
mv fungi_odb10.2021-06-28.tar.gz fungi.tar.gz
mv bacteria_odb10.2020-03-06.tar.gz bacteria.tar.gz
mv eukaryota_odb10.2020-09-10.tar.gz eukaryota.tar.gz
后面启动quast程序时应该就自动解压了。
主流数据库配置
喜欢原网站信息的个大家可以参考官网:
在GitHub仓库[https://github.com/bxlab/metaWRAP]中找到installation文件并下载。
将下载的文件进行解压缩操作。
按照指导文档中的说明设置环境变量。
执行包含数据库安装过程的脚本。
确认操作完成并检查相关日志文件的存在性以及是否存在任何错误输出。
CheckM
mkdir MY_CHECKM_FOLDER
# Now manually download the database:
cd MY_CHECKM_FOLDER
wget https://data.ace.uq.edu.au/public/CheckM_databases/checkm_data_2015_01_16.tar.gz
tar -xvf *.tar.gz
rm *.gz
cd ../
# Now you need to tell CheckM where to find this data before running anything:
checkm data setRoot # CheckM will prompt to to chose your storage location
# On newer versions of CheckM, you would run:
checkm data setRoot /path/to/your/dir/MY_CHECKM_FOLDER
KRAKEN1 (这个在新版本中可以忽略了,作者建议kraken2)
KRAKEN2数据库
为了深入研究这一领域的内容,请您参考以下论文:Metagenome analysis using the Kraken software suite | Nature Protocols。
该配置对系统资源要求较高,占用较多内存资源及存储空间;支持用户根据需求自行指定存储路径;完成构建后可通过config_matwrap进行参数配置。
This system requires approximately 120GB of RAM and roughly 128GB of storage. This stipulation is specifically relevant when employing the KRAKEN2 module.
kraken2-build --standard --threads 24 --db MY_KRAKEN2_DB
# 这个MY_KRAKEN2_DB和下方的路径配置建议都使用自己指定的绝对路径
# 时间较长,建议使用后台运行,后查看krakendb.log看是否完成
nohup kraken2-build --standard --threads 180 --db /refdb/KRAKEN2_DB >krakendb.log 2>&1 &
默认情况下会报错
rsync_from_ncbi.pl encountered an unforeseen file transfer path (What’s the new server?) when attempting to connect to https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/762/265/GCF_000762265.1_ASM76226v1.
由于kraken2的下载脚本路径头存在配置问题,在此之前使用的是FTP协议(ftpd),而现在NCBI已将其转换为HTTPS协议。因此必须重新配置该下载脚本的位置参数,并将路径指向NCBI提供的HTTPS资源。
# 大概第45行:
if (! ($full_path =~ s#^ftp://${qm_server}${qm_server_path}/##)) {
die "$PROG: unexpected FTP path (new server?) for $ftp_path\n";
}
# ftp 更改为https
if (! ($full_path =~ s#^https://${qm_server}${qm_server_path}/##)) {
die "$PROG: unexpected FTP path (new server?) for $ftp_path\n";
}
不过说起来也不一定非要自行搭建数据库呢?毕竟数据库规模较大确实会带来一定的麻烦且费时费力的选择压力。为了方便大家使用,在这里给大家推荐一个可以直接访问链接并下载所需版本的预建库到指定目录下进行解压即可的操作指南:通常情况下优先选择的标准版本可能会比某些特定需求版本更大一些(比如标准版可能比某个特定小而精的需求版更大),因此可以根据实际需求选择合适的预建包进行解压操作。
列表索引由BenLangmead创建
# standard bacteria, viral, plasmid, human1, UniVec_Core ,差不多了。
wget -c https://genome-idx.s3.amazonaws.com/kraken/k2_standard_20231009.tar.gz
# PlusPF 这个原核加真菌
wget -c https://genome-idx.s3.amazonaws.com/kraken/k2_pluspf_20231009.tar.gz
# PlusPFP , 这个还包含植物的。。。
wget -c https://genome-idx.s3.amazonaws.com/kraken/k2_pluspfp_20231009.tar.gz
# 还有个更大的 nt数据库的(710GB)。。。。准备好硬盘吧,也准备好计算资源,不然没时间等
wget -c https://genome-idx.s3.amazonaws.com/kraken/k2_nt_20231129.tar.gz
如果自己研究范围较小的话可以直接下载指定范围的,比如说:

建成或者下载解压好后设置数据库目录:
#查看当前metawrap环境下的配置文件在哪里
which config-metawrap
#cd 进入配置环境目录,或直接编辑config-metawrap文件
vim /miniconda3/envs/metawrap132/bin/config-metawrap
KRAKEN2_DB=/path/to/my/database/MY_KRAKEN2_DATABASE
NCBI nt数据库:
直接下载解压:
mkdir NCBI_nt
cd NCBI_nt
wget "ftp://ftp.ncbi.nlm.nih.gov/blast/db/nt.*.tar.gz"
for a in nt.*.tar.gz; do tar xzf $a; done
设置NCBI_nt数据库目录:
vim /miniconda3/envs/metawrap132/bin/config-metawrap
BLASTDB=/path/to/my/database/NCBI_nt
NCBI taxonomy
这个也是直接下载:
mkdir NCBI_tax
cd NCBI_tax
wget ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdump.tar.gz
tar -xvf taxdump.tar.gz
别忘记设置TAXDUMP目录:
vim /miniconda3/envs/metawrap132/bin/config-metawrap
TAXDUMP=/path/to/my/database/NCBI_tax
bmtagger(选配)
此方法旨在去除宿主基因序列或其他需排除的关键序列。在此处采用人类基因组数据库hg38作为基准进行比较分析。此系统允许用户下载最新版本的数据库供进一步研究使用。
如果自己研究设计到其他特定物种,可选择下载相应的全基因组序列:
UCSC Genome Viewer’s Downloads](https://hgdownload.soe UCSC Genome Viewer’sDownloads). This link provides access to the UCSC Genome Browser for downloading genomic data.

下载:
mkdir BMTAGGER_INDEX
cd BMTAGGER_INDEX
wget ftp://hgdownload.soe.ucsc.edu/goldenPath/hg38/chromosomes/*fa.gz
gunzip *fa.gz
cat *fa > hg38.fa
rm chr*.fa
建立数据库索引,耗内存资源高,大家根据能力来做
bmtool -d hg38.fa -o hg38.bitmask
srprism mkindex -i hg38.fa -o hg38.srprism -M 100000
最后设置数据库目录
BMTAGGER_DB=/path/to/your/index/BMTAGGER_INDEX
数据库设置建议
这里的几款 databases 均为 common 的工具。在某些场景下,所需的 reference databases 不需要过于庞大且完全不适合当前 needs。我们可以根据个人要求定制并构建适合的 databases 集合。
例如,在这里可以选择替换NCBI_nt数据库为nr的相关数据库之一;也可以直接使用若干个物种自己的特定数据库。
例如,在本例中所使用的宿主基因参考数据库是hg38版本。对于水稻田的研究者来说,可以直接获取并进行处理。
KRAKEN2数据库便于获取官网编译好的版本,请根据个人需求选择相应的版本进行下载。
METAWRAP主要模块介绍
先看图

再看帮助信息
metawrap --help
MetaWRAP v=1.3.2
Usage: metaWRAP [module]
Modules:
read_qc Raw read QC module (read trimming and contamination removal)
assembly Assembly module (metagenomic assembly)
kraken KRAKEN module (taxonomy annotation of reads and assemblies)
kraken2 KRAKEN2 module (taxonomy annotation of reads and assemblies)
blobology Blobology module (GC vs Abund plots of contigs and bins)
binning Binning module (metabat, maxbin, or concoct)
bin_refinement Refinement of bins from binning module
reassemble_bins Reassemble bins using metagenomic reads
quant_bins Quantify the abundance of each bin across samples
classify_bins Assign taxonomy to genomic bins
annotate_bins Functional annotation of draft genomes
--help | -h show this help message
--version | -v show metaWRAP version
--show-config show where the metawrap configuration files are stored
MetaWRAP系统整合了来自不同生物科的高通量测序数据,并通过一系列功能完善的工具模块进行解析与分析。该系统不仅提供了多样化的功能组合,而且能够有效支持宏基因组学分析流程中的各项关键步骤。其中每一个功能模块都承担着其独特的作用,在整个研究流程中发挥着不可或缺的作用。
Quality Control(质量控制) :
- 功能模块:该步骤负责校验与清理原始序列数据集中的低品质数据段、头部片段或干扰片段等异常内容,并确保后续分析工作能够基于高质量的数据执行。
- MetaWRAP组件:
metawrap read_qc。
2.
Assembly(组装) :
- 功能:通过将碎片序列整合形成较长的连续序列(contigs),从而获取更完整的微生物基因组信息。
- MetaWRAP模块:采用metawrap assembly方法,并结合MEGAHIT和metaSPAdes等组装器进行操作。
Binning(装箱) :
功能
MetaWRAP模块
Bin Refinement(装箱优化)
功能
Reassemble(重组装)
Bin Refinement(装箱优化)
Blobology(可视化) :
- 核心功能:Blobology 通过可视化展示的方式实现对样本数据集的处理与分析功能。该系统能够生成基于不同包裹类型及其基因多样性特征和化学易损性特性的图表,并提供与包裹相关的详细注释信息。从而帮助研究人员轻松实现对多种包裹类型间的对比分析。
- MetaWRAP 模块:
Metawrap Blobology。
Quant Bins(分箱丰度分析) :
* **功能** :计算不同分箱丰度。
* **MetaWRAP模块** :`metawrap quant_bins`
Classify Bins (功能注释) :
功能定位在于完成物种的分类标记。该系统模块主要负责通过机器学习算法对生物样本进行精准分类,并生成相应的分类标记数据。涵盖该功能的核心组件有多个关键子模块
Annotate Bins**(功能注释)** :
- 核心功能:在经过数据处理、打包、整合和分析之后,系统将自动完成数据可视化展示及深入的统计学分析。
- MetaWRAP模块:该模块配置项metawrap annotate_bins 通常需要依赖于外部工具或自定义脚本配置来实现完整的工作流程中的数据可视化展示及深入的数据分析。
MetaWRAP各模块之间分析流程关系
同样,先看图

查看官方提供的分析流程说明:...
主要模块使用方法介绍
1. read_qc
这个模块用于对测序数据进行质量控制和预处理。
# Usage: metaWRAP read_qc [options] -1 reads_1.fastq -2 reads_2.fastq -o output_dir
# Note: the read files have to be named in the name_1.fastq/name_2.fastq convention
# 名称必须为 *_1.fastq, *_2.fastq, 所以自己解压修改吧。
metawrap read_qc \
-1 input_R_1.fastq \
-2 input_R_2.fastq \
-o $OUTPUT_DIR \
-t 160
# -1 sample_R1.fastq 和 -2 sample_R2.fastq:输入的测序数据,分别是 paired-end 数据的两个文件,根据实际情况替换成您的文件名。
# --out_dir read_qc_output:输出目录,存储质量控制和预处理后的数据和报告。您可以将 read_qc_output 替换为您想要的输出目录名称。
2. assembly
用于将测序数据组装成宏基因组。
metawrap assembly -1 sample_R1.fastq -2 sample_R2.fastq -m 100 -o assembly_output
#metawrap assembly: 运行 MetaWRAP 工具中的组装模块。
#-1 sample_R1.fastq -2 sample_R2.fastq: 指定样本的 paired-end 测序数据。
#-1: 第一端测序文件。
#-2: 第二端测序文件。
#-m 100: 设置最小 contig 长度为 100bp。
#-o assembly_output: 指定输出目录,将组装结果和相关文件保存在此目录中。
metawrap assembly \
-1 read_qc/${i}/final_pure_reads_1.fastq \
-2 read_qc/${i}/final_pure_reads_2.fastq \
-o assembly/$i \
-m 500 \
-t 160
3. kraken2
Kraken2 用于进行序列的分类和分类物种鉴定。
metawrap kraken2 -t 4 -o kraken2_output --db_dir database_dir assembly_output/final.contigs.fa
# metawrap kraken2: 运行 MetaWRAP 中的 kraken2 模块。
# -t 4: 指定线程数为 4,根据您的计算资源进行调整。
# -o kraken2_output: 指定输出文件夹的名称。
# --db_dir /path/to/kraken2_database: 指定 kraken2 数据库的路径。
# /path/to/your_assembly.fasta: 输入文件为您的组装后的 contigs 或者序列文件(fasta 格式)的路径。
metawrap kraken2 \
-o $OUTPUT_DIR/kraken2/$i \
-t 160 \
-s 10000000 \
read_qc/${i}/final_pure_reads_*.fastq \
assembly/${i}/final_assembly.fasta
4. blobology
Blobology 用于可视化物种分类和注释结果。
请确保您的实验步骤已完成并已生成 Kraken ₂ 的分类报告。假设您提供的读物路径是 assembly\_output/final.contigs.fa ,而 Kraken ₂ 分析的结果则包含于两个文件中:一个是 kraken\_output/kraken\_report.tsv ,另一个则是 kraken\_output/kraneous\_taxonomy.tsv 。
metawrap blobology -a kraken2_output/kraken2_report.tsv \
-t kraken2_output/kraken2_taxonomy.tsv \
-c assembly_output/final.contigs.fa \
-o blobology_output
#
5. binning
用于将组装后的 contigs 进行 binning。
metawrap binning -o binning_output -t 4 --metabat2 assembly_output/final.contigs.fa
# -o binning_output: 指定输出结果的文件夹路径。
# -t 4: 指定线程数。
# --metabat2 assembly_output/final.contigs.fa: 使用 metabat2 方法进行 binning,assembly_output/final.contigs.fa 是组装后的 contigs 文件路径。
除了现有的参数外,还可以通过选择其他 binning 方法来进一步优化. 比如, 如果需要同时应用多种方法,可以在命令行中增加相应的参数设置, 如 --metabat2 assembly_output/final.contigs.fa --maxbin2 assembly_output/final.contigs.fa 等.
6. bin_refinement
对初始的 bin 进行进一步的优化和改进。
metawrap bin_refinement -o bin_refinement_output \
-t <threads> \
--bin_folder <bin_folder_path> \
--assembly <assembly_path> \
--proposals <proposals_folder_path> \
--max_edges <max_edges> \
--min_contig_length <min_contig_length> \
--min_bin_size <min_bin_size> \
--contig_edges <contig_edges_path> \
--depth <depth_file_path> \
--log <log_file_path>
# -t <threads>: 线程数。
# --bin_folder <bin_folder_path>: 包含初始 bin 的文件夹路径。
# --assembly <assembly_path>: 原始组装的 contigs 文件路径。
# --proposals <proposals_folder_path>: 候选 bin 的文件夹路径(可选参数)。
# --max_edges <max_edges>: 两个 bin 之间的最大边数 (默认为 5)。
# --min_contig_length <min_contig_length>: 最小 contig 长度 (默认为 1000)。
# --min_bin_size <min_bin_size>: 最小 bin 大小 (默认为 200000)。
# --contig_edges <contig_edges_path>: contig 之间的连接关系文件路径 (可选参数)。
# --depth <depth_file_path>: 深度信息文件路径 (可选参数)。
# --log <log_file_path>: 日志文件路径。
metawrap bin_refinement -o bin_refinement_output \
-t 8 \
--bin_folder initial_bins/ \
--assembly assembly/final.contigs.fa \
--proposals proposals/ \
--max_edges 10 \
--min_contig_length 1500 \
--min_bin_size 250000 \
--contig_edges contig_edges.txt \
--depth depth_info.txt \
--log bin_refinement.log
7. reassemble_bins
Reconstructing the contiguous sequences within bins. The reassemble_bins module is being utilized for reconstructing the contiguous sequences within bins in the MetaWRAP software. Its primary objective is to enhance the quality and integrity of the reconstructed sequences within bins, ensuring a more accurate representation of genomic data.
metawrap reassemble_bins \
-o reassemble_output \
-t 100 \
-b binning_output_bins \
-1 reads_1.fastq \
-2 reads_2.fastq \
--parallel \
-m 500
8. quant_bins
用于对 bin 进行定量。
metawrap quant_bins -o quant_bins_output -t 4 --bin_folder binning_output_bins -a sample_reads_1.fastq,sample_reads_2.fastq
# -o quant_bins_output:指定输出结果的目录。
# -t 4:指定线程数(这里设置为4,可以根据系统资源进行调整)。
# --bin_folder binning_output_bins:指定包含 bin 文件的文件夹路径。
# -a sample_reads_1.fastq,sample_reads_2.fastq:指定用于定量的原始测序数据,以逗号分隔的双端 fastq 文件。
9. classify_bins
对 bin 进行分类。
metawrap classify_bins -o classify_bins_output -t 4 \
--bin_folder binning_output_bins \
--kraken2_reports kraken2_output/kraken2_report.tsv \
--kraken2_taxonomy kraken2_output/kraken2_taxonomy.tsv
# -o classify_bins_output:指定输出文件夹名称。
# -t 4:指定线程数。
# --bin_folder binning_output_bins:指定包含 bin 文件的文件夹路径。
# --kraken2_reports kraken2_output/kraken2_report.tsv:指定 Kraken2 的报告文件路径。
# --kraken2_taxonomy kraken2_output/kraken2_taxonomy.tsv:指定 Kraken2 的分类学信息文件路径。
10. annotate_bins
用于注释 bin 中的基因。
metawrap annotate_bins -o annotate_bins_output \
-t 4 \
--bins binning_output_bins/*.fa \
--db_dir /path/to/your/database_directory \
--min_length 500 \
--min_identity 0.5 \
--min_coverage 0.5 \
--megahit_assembly assembly_output/final.contigs.fa
# -o annotate_bins_output:指定输出目录。
# -t 4:指定线程数。
# --bins binning_output_bins/*.fa:指定需要注释的 bin 文件路径。这里假设 bin 文件以 .fa 结尾。
# --db_dir /path/to/your/database_directory:指定用于注释的数据库目录,例如 NCBI 的 NR 数据库。
# --min_length 500:设定注释的基因序列最小长度。
# --min_identity 0.5:设定用于比对的最小相似性阈值。
# --min_coverage 0.5:设定比对序列的最小覆盖率。
# --megahit_assembly assembly_output/final.contigs.fa:指定用于注释的组装文件(Megahit 组装的最终 contigs 文件)。
基于MetaWRAP的序列处理自动分析脚本(主菜)
该程序具备以下功能:它能够自动生成后续阶段的数据存储结构,在默认模式下会启动后续的数据收集与组织模块。各个阶段产生的数据结果会被独立存储在相应的目录中,并按照预设的命名规范进行命名以避免混淆。
准备,创建work目录,并进入work目录,可以自己分析工作名称为目录。
进行宏基因组测序,并建立fastq文件夹结构。随后为每个样本生成相应的存储路径。解压配对测序数据后按照样本名称分别命名为‘样本1 fastq.sra’及‘样本2 fastq.sra’将其放置在与该样本相关的存储位置对于大量样本的情况建议采用循环处理方式
# fastq文件夹内容:
sample1 sample2 samle3 sample4 ... ...
# fastq文件夹下sample2里面的文件:
sample2_1.fastq sample2_2.fastq
mtwrp_pipeline.sh
#!/bin/bash
# 设置工作目录和文件路径,自行修改
echo "设置工作目录等"
WORK_DIR=/path/to/your/work_directory
INPUT_FASTQ=/path/to/your/input/fastq
OUTPUT_DIR=$WORK_DIR/metawrap_results
echo 创建结果文件夹
echo 建议先删除清空,再创建文件夹
rm -rf $OUTPUT_DIR
mkdir -p $OUTPUT_DIR
# read_qc 模块:质量控制 自己加控制参数,可以 --skip参数,自己使用帮助查看再修改脚本
rm -rf $OUTPUT_DIR/read_qc
mkdir -p $OUTPUT_DIR/read_qc
# 环境名称自己定,这里是metawrap, 一般source,新版本conda 或 mamba
# 下面使用echo命令按需添加日志提示,方便了解计算进度。
echo "激活metawrap环境"
source activate metawrap
# 建议原fastq目录以样品名为文件夹,不要有其他任何目录和文件,下面放clean_reads
for i in $(ls $INPUT_FASTQ);
do
metawrap read_qc \
-1 $INPUT_FASTQ/${i}/${i}_1.fastq \
-2 $INPUT_FASTQ/${i}/${i}_2.fastq \
-o $OUTPUT_DIR/read_qc/$i \
-t 160
done
# assembly 模块:组装, 计算资源足够或者基因组较小的建议使用metaspades,否则默认megahit
# 根据自己资源调整-t和-m参数,分别设置最大线程数和最大内存使用量。
rm -rf $OUTPUT_DIR/assembly
mkdir -p $OUTPUT_DIR/assembly
for i in $(ls $OUTPUT_DIR/read_qc);
do
metawrap assembly \
-1 $OUTPUT_DIR/read_qc/${i}/final_pure_reads_1.fastq \
-2 $OUTPUT_DIR/read_qc/${i}/final_pure_reads_2.fastq \
-o $OUTPUT_DIR/assembly/$i \
-m 500 \
-t 160
done
# kraken2 模块:分类学分析
# Run on any number of fasta assembly files and/or or paired-end reads.
# Usage: metaWRAP kraken2 [options] -o output_dir assembly.fasta reads_1.fastq reads_2.fastq ...
mkdir -p $OUTPUT_DIR/kraken2
for i in $(ls $OUTPUT_DIR/assembly);
do
metawrap kraken2 \
-o $OUTPUT_DIR/kraken2/$i \
-t 160 \
-s 10000000 \
$OUTPUT_DIR/read_qc/${i}/final_pure_reads_*.fastq \
$OUTPUT_DIR/assembly/${i}/final_assembly.fasta
done
# binning 模块:Bin 完整度和分离
# Usage: metaWRAP binning [options] -a assembly.fa -o output_dir readsA_1.fastq readsA_2.fastq ... [readsX_1.fastq readsX_2.fastq]
# Note1: Make sure to provide all your separately replicate read files, not the joined file.
# Note2: You may provide single end or interleaved reads as well with the use of the correct option
#Note3: If the output already has the .bam alignments files from previous runs, the module will skip re-aligning the reads
mkdir -p $OUTPUT_DIR/binning
for i in $(ls $OUTPUT_DIR/assembly);
do
metawrap binning \
-o $OUTPUT_DIR/binning/$i \
--metabat2 \
--concoct \
--maxbin2 \
-a $OUTPUT_DIR/assembly/${i}/final_assembly.fa \
$OUTPUT_DIR/read_qc/${i}/final_pure_reads_1.fastq \
$OUTPUT_DIR/read_qc/${i}/final_pure_reads_2.fastq \
--run-checkm \
--universal \
-t 160
done
# bin_refinement 模块:Bin 重构
# Usage: metaWRAP bin_refinement [options] -o output_dir -A bin_folderA [-B bin_folderB -C bin_folderC]
# Note: the contig names in different bin folders must be consistant (must come from the same assembly).
mkdir -p $OUTPUT_DIR/bin_refinement
for i in $(ls $OUTPUT_DIR/binning);
do
metawrap bin_refinement \
-o $OUTPUT_DIR/refined_bins/$i \
-t 160 \
-m 500 \
-A $OUTPUT_DIR/binning/$i/metabat2_bins \
-B $OUTPUT_DIR/binning/$i/maxbin2_bins \
-C $OUTPUT_DIR/binning/$i/concoct_bins
done
# reassemble_bins 模块:重新组装 Bin
# Usage: metaWRAP reassemble_bins [options] -o output_dir -b bin_folder -1 reads_1.fastq -2 reads_2.fastq
mkdir -p $OUTPUT_DIR/reassemble_bins
for i in $(ls $OUTPUT_DIR/bin_refinement);
do
metawrap reassemble_bins \
-o $OUTPUT_DIR/reassemble_bins/$i \
-b $OUTPUT_DIR/bin_refinement/$i \
-1 $OUTPUT_DIR/read_qc/${i}/final_pure_reads_1.fastq \
-2 $OUTPUT_DIR/read_qc/${i}/final_pure_reads_2.fastq \
-t 160
done
# blobology 模块:基于组装的可视化
# Usage: metaWRAP blobology [options] -a assembly.fasta -o output_dir readsA_1.fastq readsA_2.fastq [readsB_1.fastq readsB_2.fastq ... ]
mkdir -p $OUTPUT_DIR/blobology
for i in $(ls $OUTPUT_DIR/assembly);
do
metawrap blobology \
-t 160 \
-o $OUTPUT_DIR/blobology/$i \
-a $OUTPUT_DIR/assembly/$i/final_assembly.fasta \
--bins $OUTPUT_DIR/reassemble_bins/$i \
$OUTPUT_DIR/read_qc/${i}/final_pure_reads_1.fastq \
$OUTPUT_DIR/read_qc/${i}/final_pure_reads_2.fastq
done
# quant_bins 模块:Bin 定量
# Usage: metaWRAP quant_bins [options] -b bins_folder -o output_dir -a assembly.fa readsA_1.fastq readsA_2.fastq ... [readsX_1.fastq readsX_2.fastq]
mkdir -p $OUTPUT_DIR/quant_bins
for i in $(ls $OUTPUT_DIR/reassemble_bins);
do
metawrap quant_bins \
-o $OUTPUT_DIR/quant_bins \
-b $OUTPUT_DIR/reassemble_bins/$i \
-a $OUTPUT_DIR/assembly/$i/final_assembly.fasta \
$OUTPUT_DIR/read_qc/${i}/final_pure_reads_1.fastq \
$OUTPUT_DIR/read_qc/${i}/final_pure_reads_2.fastq
done
# classify_bins 模块:Bin 分类
# Usage: metaWRAP classify_bins [options] -b bin_folder -o output_dir
mkdir -p $OUTPUT_DIR/classify_bins
for i in $(ls $OUTPUT_DIR/reassemble_bins);
do
metawrap classify_bins \
-o $OUTPUT_DIR/classify_bins \
-b $OUTPUT_DIR/reassemble_bins/$i \
-t 160
done
# annotate_bins 模块:Bin 注释
# Usage: metaWRAP annotate_bins [options] -o output_dir -b bin_folder
mkdir -p $OUTPUT_DIR/annotate_bins
for i in $(ls $OUTPUT_DIR/reassemble_bins);
do
metawrap annotate_bins \
-o $OUTPUT_DIR/annotate_bins \
-b $OUTPUT_DIR/reassemble_bins/$i \
-t 160
done
echo ------------MetaWRAP流程基本模块分析完成-------------------
# 继续开始模块外必要序列处理分析
# orf预测,基因预测,减少不必要的序列比对资源消耗,当然大家可以直接使用assembly或reassembly结果直接进行基因注释
rm -rf $OUTPUT_DIR/prodigal
mkdir -p $OUTPUT_DIR/prodigal
for i in $(ls $OUTPUT_DIR/assembly);
do
prodigal \
-i $OUTPUT_DIR/assembly/${i}/final_assembly.fa \
-a $OUTPUT_DIR/assembly/${i}.protein_seq.fa \
-d $OUTPUT_DIR/assembly/${i}.nucleotide_seq.fa \
-o $OUTPUT_DIR/assembly/${i}.genes.gff
done
#
mamba deactivate
# 输出分析完成信息
echo "MetaWRAP 分析完成!结果存储在 $OUTPUT_DIR 中。"
运行:
#建议后台
nohup sh mtwrp_pipeline.sh > pp.log 2>&1 &
# 跟踪进度和查看日志:
tail -f pp.log
该流程可应用于单个计算节点处理一批样品的自动化操作。当涉及多台计算节点时,请根据集群当前状态添加节点切换与结果汇总功能。
MetaWRAP工作流注意事项:
请确认metawrap环境的安装及其相关数据库配置已正确设置,在执行流程操作前,请采用小型fastq数据集作为初始测试阶段,请确保上述配置均处于有效状态。
为了维护命名的一致性与稳定性,请遵循以下原则:首先将数据类型的文件和目录名称保持稳定以避免混淆;其次,在进行新目录或新文件的创建之前,请通过预先判断来决定是否需要清除现有数据;最后,在完成所有相关操作后,请删除现有数据后再进行新建操作以确保命名的一致性和完整性。
掌握宏基因组学分析工作流程的基本步骤:该MetaWRAP工作流程包含了一系列关键步骤,如质量控制、序列组装与基因注释等环节。建议深入学习每个步骤的功能及其输出结果,以便全面理解整个分析流程及其结果的解读。
掌握metawrap各个模块的参数配置,在需要的情况下参考系统帮助信息来个别地调节各个模块的参数设置。
在大规模的数据处理场景中需要注意事项尤其是涉及宏基因组级别的数据分析时可能会有一些必要的配置调整以及额外资源的需求。比如优化并行作业的数量通过分布式架构分散计算资源或者采用高性价比的高性能计算集群来进行复杂的分析工作以提升效率和效果。特别是在涉及宏基因组级别的数据分析时可能会有一些必要的配置调整以及额外资源的需求。
程序中断处理机制:将每个for循环视为独立的任务,在程序执行过程中若发生突然中断,请参考日志信息分析当前的工作流程状态。对于那些已经明确完成的任务,并且拥有有效结果的任务,在后续程序重启时建议对该已完成的部分及相关目录进行注释式删除,并保留相关代码标记起始点;而对于仅部分完成的模块,则应采取全部清除现有构建目录并重新初始化的方式。
同时使用echo等命令或其他指定工具标记关键起始点。
工作流的个性化配置:在对工作流程进行删减等修改时,请务必特别注意各运行单元的数据流转衔接。建议在彻底测试确认无误后才将其整合到整体流程中。
基于assembly结果或基于prodigal结果进行序列比对大家可参考:
diamond这一款高效的大基因序列比对工具的操作指南详细解析,并涵盖超级计算集群下的多节点并行计算操作流程。此外,在博客平台上的应用说明也一并提供。
bioinformatics analysis - BLAST sequence alignment and detailed results explanation on blog
碳水化合物库的分析参考,基于assembly结果或:
基于dbCAN系统的碳水化合物酶基因数据库及其相关功能实现研究与应用-博客
humann参考分析,直接基于read_qc结果:
源自生物巴克力实验室的一项重要研究成果——宏基因组学微生物群落代谢功能分析工具HUMAnN 3.0——为科研人员提供了便捷的安装指南及使用说明。该软件通过整合多组学数据,显著提升了对微生物生态系统的解析能力,特别适用于环境微生物学领域的研究与应用
该文详细介绍了基于QC校准后的读取数据进行宏基组学物种分析的方法与操作流程。文中提供了一份针对Anaconda3平台运行于centos9系统上的MetaPhlAn4安装与使用指南,并通过 streaming 管道实现了数据处理过程的具体实现细节。该文为实验研究人员提供了便捷的操作指南以完成宏基因组学物种级别的分析工作
至此为止,我们现在已经能够应用Metawrap对宏基因组序列展开全面处理工作.后续阶段中,我们将根据metawrap工具输出的结果系统性地进行统计分析以及制作详实的科学图表.
如有任何问题,欢迎大家留言。
