EasyMetagenome易宏基因组——简单易用的宏基因组分析流程-来自刘永鑫团队的秘密武器

最新推荐文章于 2024-08-15 09:55:17 发布

2401_83947353

最新推荐文章于 2024-08-15 09:55:17 发布

阅读量1.1k

点赞数 7

分类专栏： 2024年程序员学习文章标签：数据库服务器网络

本文链接：https://blog.csdn.net/2401_83947353/article/details/137762303

版权

2024年程序员学习专栏收录该内容

271 篇文章 3 订阅

订阅专栏

先自我介绍一下，小编浙江大学毕业，去过华为、字节跳动等大厂，目前阿里P7

深知大多数程序员，想要提升技能，往往是自己摸索成长，但自己不成体系的自学效果低效又漫长，而且极易碰到天花板技术停滞不前！

因此收集整理了一份《2024年最新Linux运维全套学习资料》，初衷也很简单，就是希望能够帮助到想自学提升又不知道该从何学起的朋友。

既有适合小白学习的零基础资料，也有适合3年以上经验的小伙伴深入学习提升的进阶课程，涵盖了95%以上运维知识点，真正体系化！

由于文件比较多，这里只是将部分目录截图出来，全套包含大厂面经、学习笔记、源码讲义、实战项目、大纲路线、讲解视频，并且后续会持续更新

如果你需要这些资料，可以添加V获取：vip1024b （备注运维）

正文

初始化：每次开始安装必须运行下面代码 Initialization: The following code must be run every time when installation starts

安装前准备：软件和数据库位置 Before Installation: Software and Database Locations

数据库安装位置Database Locations，默认~/db目录(无需管理权限)，管理员可选/db
db=~/db
mkdir -p ${db} && cd ${db}

软件安装位置Software installation location，默认为~/miniconda3，测试服务器为/anaconda3

soft=~/miniconda3
经常使用的服务器环境，可把全文 ${db}和$ {soft}替换为绝对路径，将不再需要每次读取以上环境变量

In the frequently used server environment, you can replace the variable ${db} and ${soft} with absolute paths, and you will no longer need to run the above environment variables every time

可选：初始化环境变量，可能提高软件安装成功率

Optional: Initialize environment variables, which may improve the success rate of software installation

PATH= ${soft}/bin:$ {soft}/condabin:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin: ${db}/EasyMicrobiome/linux:$ {db}/EasyMicrobiome/script
echo $PATH

EasyMetagenome流程Pipeline

EasyMetagenome流程，包括流程安装、使用和可视化脚本，以及流程测试数据和结果正对照，网址：https://github.com/YongxinLiu/EasyMetagenome

三种下载方法：依赖尝试/任选其一至成功即可

# 方法1. 网页 https://github.com/YongxinLiu/EasyMetagenome 中Code - Download ZIP下载压缩包，上传至服务器
# 解压，Command 'unzip' not found使用sudo apt install unzip安装
unzip EasyMetagenome-master.zip 
# 改名
mv EasyMetagenome-master EasyMetagenome

# 方法2. 微生物所备用链接，可能不是最新版
wget -c ftp://download.nmdc.cn/tools/soft/EasyMetagenome.tar.gz
tar xvzf EasyMetagenome.tar.gz

# 方法3. git下载，需安装git
git clone https://github.com/YongxinLiu/EasyMetagenome
# 旧版更新
cd EasyMetagenome && git pull && cd ../

EasyMicrobiome软件和数据库合集dependencies

EasyMetagenome依赖EasyMicrobiome，其包括众多脚本、软件和数据库的集合，网址：https://github.com/YongxinLiu/EasyMicrobiome

# 方法1. 网页中下载
# https://github.com/YongxinLiu/EasyMicrobiome 中Code Download ZIP下载压缩包，上传至服务器，并解压
unzip EasyMicrobiome-master.zip
mv EasyMicrobiome-master EasyMicrobiome

# 方法2. 备用链接下载
wget -c ftp://download.nmdc.cn/tools/soft/EasyMicrobiome.tar.gz
tar -xvzf EasyMicrobiome.tar.gz

# 方法3. git下载，需安装git
git clone https://github.com/YongxinLiu/EasyMicrobiome
# 旧版更新
cd EasyMicrobiome && git pull && cd ../

软件安装

# 添加linux命令可执行权限
chmod +x `pwd`/EasyMicrobiome/linux/* `pwd`/EasyMicrobiome/script/*
# 添加环境变量
echo "export PATH=\"\$PATH:`pwd`/EasyMicrobiome/linux:`pwd`/EasyMicrobiome/script\"" >> ~/.bashrc
source ~/.bashrc
echo $PATH

软件管理器Conda

# 下载最新版miniconda3，70M
wget -c https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
# 安装，-b批量，-f无提示，-p目录，许可协议打yes
bash Miniconda3-latest-Linux-x86_64.sh -b -f -p ${soft}
# 激活，然后关闭终端重开，提示符前出现(base)即成功
${soft}/condabin/conda init
source ~/.bashrc
# 查看版本，conda 23.7.3, python 3.11.4
conda -V 
python --version
# 添加常用频道
conda config --add channels bioconda # 生物软件
conda config --add channels conda-forge # Highest priority

# conda默认配置文件为 ~/.condarc 查看配置文件位置
conda install mamba -c conda-forge -y
mamba install pandas -c conda-forge -y
mamba install conda-pack -c conda-forge -y
conda config --set channel_priority strict
conda config --show-sources
# 查看虚拟环境列表 
conda env list

质控Quality control: kneaddata/fstqc/multiqc/fastp

**注：直接安装、下载解压安装，二选一。一种方法不成功，尝试另一种。**

BioConda: https://bioconda.github.io/recipes/kneaddata/README.html

方法1.kneaddata直接安装

# 新建kneaddata环境
conda create -y -n kneaddata
conda activate kneaddata

fastqc质量评估，multiqc评估报告汇总，kneaddata质量控制流程，fastp质控工具

mamba install kneaddata fastqc multiqc fastp r-reshape2 -y

方法2.kneaddata下载解压安装

# 指定conda文件名
s=kneaddata
# 下载，可选NMDC、百度云等
wget -c ftp://download.nmdc.cn/tools/conda/${s}.tar.gz
# 指定安装目录
mkdir -p ${soft}/envs/${s}
tar -xvzf ${s}.tar.gz -C ${soft}/envs/${s}
# 启动环境
conda activate ${s}
# 初始化环境
conda unpack

kneaddata安装测试

fastqc -v # v0.12.1
kneaddata --version # 0.12.0
trimmomatic -version # 0.39
bowtie2 --version # 2.5.1
multiqc --version  # 1.15

# (可选)安装软件打包，f覆盖输出文件，ignore跳过修改检测
n=kneaddata
conda pack -f --ignore-missing-files -n ${n} -o ${n}.tar.gz

kneaddata数据库下载

# 查看可用数据库
kneaddata_database
# 包括人基因组bowtie2/bmtagger、人类转录组、核糖体RNA和小鼠基因组

# 下载人基因组bowtie2索引 3.44 GB
mkdir -p ${db}/kneaddata/human
kneaddata_database --download human_genome bowtie2 ${db}/kneaddata/human

# 备用链接下载人类基因组至上述目录，并解压
cd ${db}/kneaddata/human
wget -c ftp://download.nmdc.cn/tools/meta/kneaddata/human_genome/Homo_sapiens_hg37_and_human_contamination_Bowtie2_v0.1.tar.gz
tar xvzf Homo_sapiens_hg37_and_human_contamination_Bowtie2_v0.1.tar.gz

# 下载小鼠基因组bowtie2索引 2.83 GB
mkdir -p ${db}/kneaddata/mouse
kneaddata_database --download mouse_C57BL bowtie2 ${db}/kneaddata/mouse
# 备用手动下载
cd ${db}/kneaddata/mouse
wget -c http://huttenhower.sph.harvard.edu/kneadData_databases/mouse_C57BL_6NJ_Bowtie2_v0.1.tar.gz
tar xvzf mouse_C57BL_6NJ_Bowtie2_v0.1.tar.gz

kneaddata自定义参考基因组索引

**(任何构建好的 bowtie2 索引都可，放置位置不限；如果有多个宿主，可以把多个宿主参考基因组序列合并后构建索引)**

自定义基因组构建索引，大多数基因组可在ensembl genome下载。此处以拟南芥为例，访问 http://plants.ensembl.org/index.html ，选择Arabidopsis thaliana —— Download DNA sequence (FASTA)，选择toplevel右键复制链接，填入下面链接处

# 创建子目录
mkdir -p ${db}/kneaddata/ath
cd ${db}/kneaddata/ath
# 下载
wget -c http://ftp.ensemblgenomes.org/pub/plants/release-51/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
mv Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz tair10.fa.gz
# wget -c ftp://download.nmdc.cn/tools/meta/kneaddata/tair10.fa.gz
# 解压
gunzip tair10.fa.gz
# bowtiew建索引，输入文件，输出文件前缀，4线程2分
bowtie2-build -f tair10.fa tair10 --threads 4

二、基于读长分析 Read-based (HUMAnN3/Kraken2)

宏基因组基于读长的分析 HUMAnN3/MetaPhlAn4/GraPhlAn

HUMAnN3+MetaPhlAn4为目前最新版，目前最广泛使用的HUMAnN2安装见附录

HUMAnN3直接安装

# 安装HUMAnN3.7+MetaPhlAn4
conda create -n humann3
conda activate humann3
conda install humann=3.7 -c bioconda -c conda-forge
# 打包(可选)
conda pack -f -n humann3 -o humann3.tar.gz

HUMAnN3解包安装

# 下载
wget -c ftp://download.nmdc.cn/tools/conda/humann3.tar.gz
# 指定安装目录
mkdir -p ${soft}/envs/humann3
tar -xvzf humann3.tar.gz -C ${soft}/envs/humann3
# 启动环境
conda activate humann3
# 初始化环境
conda unpack

HUMAnN3安装测试

# 记录核心软件版本
humann --version # v3.7
metaphlan -v # 4.0.6 (1 Mar 2023)
diamond help | head -n 1 #  v2.1.8.162
# 测试
humann_test

HUMAnN3物种和功能数据库

# 显示可用分类、泛基因组和功能数据库
humann_databases

# 安装数据库
cd ${db}
mkdir -p ${db}/humann3 # 建立下载目录
# 微生物泛基因组 16 GB
humann_databases --download chocophlan full ${db}/humann3
# 功能基因diamond索引 20 GB
humann_databases --download uniref uniref90_diamond ${db}/humann3
# 输助比对数据库 2.6 GB
humann_databases --download utility_mapping full ${db}/humann3

# humann3数据库无法自动下载，备用链接下载安装
wget -c ftp://download.nmdc.cn/tools/meta/humann3/full_chocophlan.v201901_v31.tar.gz
wget -c ftp://download.nmdc.cn/tools/meta/humann3/uniref90_annotated_v201901b_full.tar.gz
wget -c ftp://download.nmdc.cn/tools/meta/humann3/full_mapping_v201901b.tar.gz
# 安装、解压
mkdir -p ${db}/humann3/chocophlan
tar xvzf full_chocophlan.v201901_v31.tar.gz -C ${db}/humann3/chocophlan
mkdir -p ${db}/humann3/uniref
tar xvzf uniref90_annotated_v201901b_full.tar.gz -C ${db}/humann3/uniref
mkdir -p ${db}/humann3/utility_mapping
tar xvzf full_mapping_v201901b.tar.gz -C ${db}/humann3/utility_mapping

# 设置数据库位置
# 显示参数
humann_config --print
# 如修改线程数，推荐3-8，根据实际情况调整
humann_config --update run_modes threads 8
# 设置核酸、蛋白和注释库位置
humann_config --update database_folders nucleotide ${db}/humann3/chocophlan
humann_config --update database_folders protein ${db}/humann3/uniref
humann_config --update database_folders utility_mapping ${db}/humann3/utility_mapping
# 核对设置结果
humann_config --print

MetaPhlAn4物种数据库

个人整理了一个详细使用流程：

202310-宏基组学物种分析工具-MetaPhlAn4安装和使用方法-Anaconda3- centos9 stream-CSDN博客

# MetaPhlAn4数据库下载2022数据和索引3G+20G
mkdir -p ${db}/metaphlan4
cd ${db}/metaphlan4

# 官网下载
wget -c http://cmprod1.cibio.unitn.it/biobakery4/metaphlan_databases/mpa_vOct22_CHOCOPhlAnSGB_202212.tar
wget -c http://cmprod1.cibio.unitn.it/biobakery4/metaphlan_databases/bowtie2_indexes/mpa_vOct22_CHOCOPhlAnSGB_202212_bt2.tar
tar xvf mpa_vOct22_CHOCOPhlAnSGB_202212.tar
tar xvf mpa_vOct22_CHOCOPhlAnSGB_202212_bt2.tar

# 官方没有压缩体积大下载慢，备用国内百度链接：https://pan.baidu.com/s/1Ikd_47HHODOqC3Rcx6eJ6Q?pwd=0315 或 微生物所FTP ftp://download.nmdc.cn/tools/meta 下载压缩包
wget -c ftp://download.nmdc.cn/tools/meta/metaphlan4/mpa_vOct22_CHOCOPhlAnSGB_202212.tar.gz
wget -c ftp://download.nmdc.cn/tools/meta/metaphlan4/mpa_vOct22_CHOCOPhlAnSGB_202212_bt2.tar.gz
tar xvzf mpa_vOct22_CHOCOPhlAnSGB_202212.tar.gz
tar xvzf mpa_vOct22_CHOCOPhlAnSGB_202212_bt2.tar.gz
# 可选(制作下载文件和md5值)
gunzip mpa_vOct22_CHOCOPhlAnSGB_202212_bt2.tar.gz
md5sum mpa_vOct22_CHOCOPhlAnSGB_202212_bt2.tar > mpa_vOct22_CHOCOPhlAnSGB_202212_bt2.md5

生物标记鉴定和可视化LEfSe

方法1. 在线ImageGP http://www.bic.ac.cn/ImageGP/ 选择左侧LEfSe

方法2. 下载conda预安装包解压

n=lefse
# 下载
wget -c ftp://download.nmdc.cn/tools/conda/${n}.tar.gz
# 指定安装目录
mkdir -p ${soft}/envs/${n}
tar -xvzf ${n}.tar.gz -C ${soft}/envs/${n}
# 启动环境
conda activate ${n}
# 初始化环境
conda unpack

方法3. 直接安装：经常无法安装，安装的代码也有可能存在问题

mamba create -n lefse lefse -c bioconda -y

物种注释Kraken2/bracken/krakentools/krona

kraken2 基于LCA算法的物种注释 https://ccb.jhu.edu/software/kraken/

Kraken2解包安装

# 下载
wget -c ftp://download.nmdc.cn/tools/conda/kraken2.tar.gz
# 指定安装目录
mkdir -p ${soft}/envs/kraken2
tar -xvzf kraken2.tar.gz -C ${soft}/envs/kraken2
# 启动环境
conda activate kraken2
# 初始化环境
conda unpack

Kraken2直接安装，居然只安了2.0.7

mamba create -n kraken2 -y -c bioconda kraken2 bracken krakentools krona r-optparse

记录软件版本

kraken2 --version # 2.1.2

Kraken2数据库安装

下载数据库(NCBI每2周更新一次)，记录下载日期和大小。需根据服务器内存、使用目的选择合适方案。–standard标准模式下只下载5种**标准数据库：古菌archaea、细菌bacteria、人类human、载体UniVec_Core、病毒viral**。也可选直接下载作者构建的索引，还包括bracken的索引。链接：https://benlangmead.github.io/aws-indexes/k2 （10/9/2023版）。注：中科院网络下载较快，家里和农科院较慢，有时新版会有错误，可以退回旧版

方案1. 下载标准+原生动物+真菌 16GB (PlusPF-16) 

v=k2_pluspf_16gb_20231009
mkdir -p ~/db/kraken2/pluspf16g
cd ~/db/kraken2
wget -c https://genome-idx.s3.amazonaws.com/kraken/${v}.tar.gz
# 备用链接
wget -c ftp://download.nmdc.cn/tools/meta/kraken2/${v}.tar.gz
tar xvzf ~/db/kraken2/${v}.tar.gz -C ~/db/kraken2/pluspf16g

方案2. 下载标准+原生动物+真菌 69GB (PlusPF) 

v=k2_pluspf_20231009
mkdir -p ~/db/kraken2/pluspf
cd ~/db/kraken2
wget -c https://genome-idx.s3.amazonaws.com/kraken/${v}.tar.gz
# 备用链接
wget -c ftp://download.nmdc.cn/tools/meta/kraken2/${v}.tar.gz
tar xvzf ~/db/kraken2/${v}.tar.gz -C ~/db/kraken2/pluspf

方案3. 下载标准+原生动物+真菌+植物完整库 144G (PlusPFP) 

指定解压目录，包括时间和类型

v=k2_pluspfp_20231009
mkdir -p ~/db/kraken2/pluspfp
cd ~/db/kraken2
wget -c https://genome-idx.s3.amazonaws.com/kraken/${v}.tar.gz
tar xvzf ${db}/kraken2/${v}.tar.gz -C pluspfp

三、组装 Assemble-based

组装、注释和定量 megahit/spades/quast/cd-hit/emboss/salmon/prodigal

megahit解包安装

# 下载
wget -c ftp://download.nmdc.cn/tools//conda/megahit.tar.gz
# 指定安装目录
mkdir -p ${soft}/envs/megahit
tar -xvzf megahit.tar.gz -C ${soft}/envs/megahit
# 启动环境
conda activate megahit
# 初始化环境
conda unpack

megahit直接安装

mamba create -y -n megahit megahit spades quast cd-hit emboss salmon prodigal
conda activate megahit

megahit安装后测试

megahit -v # MEGAHIT v1.2.9
metaspades.py -v # metaSPAdes v3.15.4
metaquast.py -v # MetaQUAST v5.0.2
cd-hit -v | grep version # CD-HIT v4.8.1
embossversion # EMBOSS v6.6
salmon -v # salmon v1.8

蛋白同源综合注释eggNOG

eggNOG http://eggnogdb.embl.de

eggNOG解包安装

# 下载
wget -c ftp://download.nmdc.cn/tools/conda/eggnog.tar.gz
# 指定安装目录
mkdir -p ${soft}/envs/eggnog
tar -xvzf eggnog.tar.gz -C ${soft}/envs/eggnog
# 启动环境
conda activate eggnog
# 初始化环境
conda unpack

eggNOG直接安装

# 新建环境并进入
conda create -n eggnog -y
conda activate eggnog
# 安装eggnog比对工具emapper
conda install eggnog-mapper -y -c bioconda -c conda-forge

eggNOG安装测试

emapper.py --version # 2.1.12
#  Expected eggNOG DB version: 5.0.2 / Installed eggNOG DB version: 5.0.2 / 
#  Diamond version found: diamond version 2.0.15 / MMseqs2 version found: 13.45111

eggNOG数据库安装

# 下载常用数据库，注意设置下载位置
mkdir -p ${db}/eggnog && cd ${db}/eggnog
# -y默认同意，-f强制下载，eggnog.db.gz 6.3G+4.9G，解压后48G
download_eggnog_data.py -y -f --data_dir ${db}/eggnog

# 百度或微生物所备用链接下载eggnog.tar.gz
wget -c ftp://download.nmdc.cn/tools/meta/eggnog/eggnog.tar.gz
tar xvzf eggnog.tar.gz
# 查看版本时会显示默认数据位置
emapper.py --version # 2.1.12
# 链接至默认目录，注意按实际情况修改
ln -sf ${db}/eggnog ${soft}/envs/eggnog/lib/python3.9/site-packages/data/
# 复制数据至内存中加速比对
# cp eggnog.* /dev/shm

碳水化合物CAZy

dbCAN3 http://bcb.unl.edu/dbCAN2

# 创建数据库存放目录并进入
mkdir -p ${db}/dbcan3 && cd ${db}/dbcan3
# 下载序列和描述
wget -c https://bcb.unl.edu/dbCAN2/download/Databases/V12/CAZyDB.07262023.fa
wget -c https://bcb.unl.edu/dbCAN2/download/Databases/V12/CAZyDB.08062022.fam-activities.txt
# 提取基因家簇对应注释
grep -v '#' CAZyDB.08062022.fam-activities.txt | sed 's/  //'| \
  sed '1 i CAZy\tDescription' > CAZy_description.txt
# 打包压缩
tar -cvzf CAZyDB.tar.gz CAZyDB.07262023.fa CAZyDB.08062022.fam-activities.txt CAZy_description.txt

# 备用数据库下载并解压(待上传)
wget -c ftp://download.nmdc.cn/tools/meta/dbcan3/CAZyDB.tar.gz
tar xvzf CAZyDB.tar.gz

# diamond建索引，1G，11s
diamond --version # 2.1.8
time diamond makedb --in CAZyDB.07262023.fa --db CAZyDB

抗生素抗性基CARD/rgi

CARD官网：https://card.mcmaster.ca
RGI Github: https://github.com/arpcard/rgi

rgi解包安装

# 下载
wget -c ftp://download.nmdc.cn/tools/conda/rgi.tar.gz
# 指定安装目录
mkdir -p ${soft}/envs/rgi
tar -xvzf rgi.tar.gz -C ${soft}/envs/rgi
# 启动环境
conda activate rgi
# 初始化环境
conda unpack

rgi直接安装

mamba create -y -n rgi6 rgi=6.0.3
conda activate rgi6
# (可选)打包，待上传
n=rgi6
conda pack -f --ignore-missing-files -n ${n} -o ${n}.tar.gz

rgi版本和数据库部署

# 查看版本 6.0.3
rgi main -v

# 数据库部署
mkdir -p ${db}/card && cd ${db}/card
# 下载最新版数据库，3.8M (2023-10-2, 3.2.8)
wget -c https://card.mcmaster.ca/latest/data
# 解压后35M
tar -xvf data ./card.json
# 加载数据库
rgi load --card_json card.json
# 宏基因组分析扩展数据库和加载
rgi card_annotation -i card.json
mv card_database_v3.2.8_all.fasta card.fasta
rgi load -i card.json --card_annotation card.fasta

四、分箱挖掘单菌基因组Binning

metawrap分箱binning

软件主页：https://github.com/bxlab/metaWRAP

metawrap下载安装

# 下载
wget -c ftp://download.nmdc.cn/tools/conda/metawrap.tar.gz
# 指定安装目录
mkdir -p ${soft}/envs/metawrap
tar -xvzf metawrap.tar.gz -C ${soft}/envs/metawrap
# 启动环境
conda activate metawrap
# 初始化环境
conda unpack

metawrap conda安装

mamba create -y --name metawrap --channel ursky metawrap-mg=1.3.2
conda activate metawrap

metawrap相关数据库

cd ${db} 

CheckM用于Bin完整和污染估计和物种注释

mkdir -p checkm && cd checkm
# 下载文件275 MB，解压后1.4 GB
wget -c https://data.ace.uq.edu.au/public/CheckM_databases/checkm_data_2015_01_16.tar.gz
tar -xvf *.tar.gz
# 设置数据库位置，直接2次回车默认为当前位置
checkm data setRoot

NCBI核酸和物种信息(可选)

# 核酸
mkdir -p ${db}/NCBI/nt
(cd ${db}/NCBI/nt; wget -c ftp://ftp.ncbi.nlm.nih.gov/blast/db/nt.*.tar.gz)
(cd ${db}/NCBI/nt; for i in *.tar.gz; do tar xzf $i; done)
# 可能会出现个别库下载不完整的情况，删了重下，不要续传
# 物种信息，压缩文件45M，解压后351M

mkdir -p ${db}/NCBI/tax
(cd ${db}/NCBI/tax; wget -c ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdump.tar.gz)
(cd ${db}/NCBI/tax; tar -xvzf taxdump.tar.gz)

数据库位置设置

which config-metawrap
# 配置文件通常为~/miniconda3/envs/metawrap/bin/config-metawrap
# 使用Rstudio/vim等文本编辑器来修改数据库的位置
# config-metawrap文件中内容如下

# Paths to metaWRAP scripts (dont have to modify)
mw_path=$(which metawrap)
bin_path=${mw_path%/*}
SOFT=${bin_path}/metawrap-scripts
PIPES=${bin_path}/metawrap-modules

# CONFIGURABLE PATHS FOR DATABASES (see 'Databases' section of metaWRAP README for details)
# path to kraken standard database
KRAKEN_DB=~/KRAKEN_DB
KRAKEN2_DB=~/db/kraken2/pluspf/

# path to indexed human (or other host) genome (see metaWRAP website for guide). This includes .bitmask and .srprism files
BMTAGGER_DB=~/BMTAGGER_DB

# paths to BLAST databases
BLASTDB=~/db/NCBI/nt
TAXDUMP=~/db/NCBI/tax

drep基因组去冗余

挑单菌测序的基因组存在大量冗余。metawrap混合分箱的结果中冗余度非常低，甚至无冗余。而单样本、分批次分箱的结果中存在大量冗余，需要采用drep获得非冗余的基因组。
GitHub: https://github.com/MrOlm/drep
Conda: https://bioconda.github.io/recipes/drep/README.html

drep 基因组去冗余解包安装

# 下载dRep v3.2.3无法安装依赖chechm，仍用旧版2.6.2(500M)，这个压缩包没有checkm且版本为3.4.2
wget -c ftp://download.nmdc.cn/tools/conda/drep.tar.gz
# 指定安装目录
mkdir -p ${soft}/envs/drep
tar -xvzf drep.tar.gz -C ${soft}/envs/drep
# 启动环境
conda activate drep
# 初始化环境
conda unpack
dRep -h

drep 基因组去冗余直接安装

# 2023/9/12尝试，仍无法安装checkm
mamba create -y -n drep drep=3.4.3
conda activate drep
# 不满足依赖关系
mamba install checkm-genome -y
dRep -h

drep 数据库构建

CheckM用于Bin完整和污染估计和物种注释，安装过metawrap已经下载完成

mkdir -p ${db}/checkm && cd ${db}/checkm
# 下载文件275 MB，解压后1.4 GB
wget -c https://data.ace.uq.edu.au/public/CheckM_databases/checkm_data_2015_01_16.tar.gz
tar -xvf *.tar.gz
# 设置数据库位置，直接2次回车默认为当前位置
checkm data setRoot `pwd`

coverm基因组定量

conda安装

conda create -n coverm -y
conda activate coverm
conda install coverm -c bioconda -y
# conda安装后打包(可选)
conda pack -f --ignore-missing-files -n coverm -o coverm.tar.gz

压缩包安装

# 指定conda文件名
s=coverm
# 下载，可选NMDC、百度云等
# wget -c ftp://download.nmdc.cn/tools/conda/${s}.tar.gz
# 指定安装目录
mkdir -p ~/miniconda3/envs/${s}
tar -xvzf ${s}.tar.gz -C ~/miniconda3/envs/${s}
# 启动环境
conda activate ${s}
# 初始化环境
conda unpack

GTDB细菌基因组注释和进化分析

Github: https://github.com/Ecogenomics/GTDBTk
GTDB-Tk是一个软件工具包，用于根据基因组数据库分类法GTDB为细菌和古细菌基因组分配客观的分类法。它旨在与最近的进展一起使用，从而可以直接对环境样本中获得数百或数千个由宏基因组组装的基因组（MAG）进行物种分类注释。它也可以用于分离和单细胞的基因组物种注释。
本次测试版本为 gtdbtk-2.2.6，Release 07-RS207v2 (11th May 2022)。
硬件要求：内存200Gb，硬盘66Gb，64核1小时可分析1000个细菌基因组

GTDB-Tk直接安装

# gtdbtk-2.3.2, 2023-7-8
n=gtdbtk2.3
mamba create -y -n ${n} -c conda-forge -c bioconda gtdbtk=2.3.2
# 检查版本
gtdbtk -v # 2.3.2

# conda pack软件打包一次
# --exclude gtdbtk-2.3.2 指定排除数据库
conda pack -n ${n} -o ${n}.tar.gz --exclude gtdbtk-2.3.2 --ignore-editable-packages --ignore-missing-files
chmod 755 *

GTDB-Tk解包安装

soft=~/miniconda3
# 下载，目前为2.1，需更新为2.3
wget -c ftp://download.nmdc.cn/tools/conda/gtdbtk.tar.gz
# 指定安装目录
mkdir -p ${soft}/envs/gtdbtk
tar -xvzf gtdbtk.tar.gz -C ${soft}/envs/gtdbtk
# 启动环境
conda activate gtdbtk
# 初始化环境
conda unpack

GTDB-Tks数据库安装

download-db.sh自动下载数据库，将下载至conda中的envs/gtdbtk/share/gtdbtk-2.3.2/db/，我们修改为~/db/gtdb中

conda activate gtdbtk2.3
# download-db.sh中，修改数据库下载位置，的 wget 建议改成wget -c 防止覆盖
sed -i 's#miniconda3/envs/gtdbtk2.3/share/gtdbtk-2.3.2/db#db/gtdb2.3#;s/wget /wget -c /' ${soft}/envs/gtdbtk2.3/bin/download-db.sh
# 下载数据,78G
download-db.sh

(备选)下面无法下载时手动下载和配置GTDB数据库

mkdir -p ${db}/gtdb2.3 && cd ${db}/gtdb2.3
# 下载解压
wget -c https://data.gtdb.ecogenomic.org/releases/release214/214.0/auxillary_files/gtdbtk_r214_data.tar.gz
# 再运行, gtdb配置数据库
download-db.sh

# 备用链接和手工解压，指定安装完整路径
wget -c ftp://download.nmdc.cn/tools/meta/gtdb/gtdbtk_r214_data.tar.gz
tar xvzf gtdbtk_r207_v2_data.tar.gz -C ./  --strip 1
conda env config vars set GTDBTK_DATA_PATH="/data/meta/db/gtdb/"

五、单菌基因组、病毒组等其他软件

CheckM2

Conda主页：https://bioconda.github.io/recipes/checkm2/README.html

软件主页：https://github.com/chklovski/CheckM2

# 软件安装
mamba create --name checkm2 checkm2
conda activate checkm2
checkm2 -h # CheckM2 v1.0.1
# 数据库安装
mkdir ~/db/checkm2
checkm2 database --download --path ~/db/checkm2
# 报错：checkm2.zenodo_backpack.ZenodoConnectionException: Connection error: HTTPSConnectionPool(host='zenodo.org', port=443): Max retries exceeded with url: /record/5571251 (Caused by NewConnectionError('<urllib3.connection.HTTPSConnection object at 0x7f335c64abe0>: Failed to establish a new connection: [Errno 111] Connection refused'))
# 在github中搜索和查找issues无解答，提新issue https://github.com/chklovski/CheckM2/issues
# 数据库 https://zenodo.org/records/5571251  下载，需要VPN

export CHECKM2DB="path/to/database"
# 测试
checkm2 testrun

# 运行，输入目录或文件列表
checkm2 predict --threads 30 --input <folder_with_bins> --output-directory <output_folder> 
checkm2 predict --threads 30 --input ../bin1.fa ../../bin2.fna /some/other/directory/bin3.fasta --output-directory <output_folder>

常见问题

软件和数据库国内备份

国家微生物科学数据中心 —— 数据下载

http://nmdc.cn/datadownload，可以使用Filezilla直接连接 ftp://download.nmdc.cn/tools

本资源由宏基因组平台发起，微生物所提供服务器，宏基因组团队负责维护的常用软件、扩增子和宏基因组数据库的国内下载链接，解决常用数据库下载慢、或无法下载的问题。同时提供定制的软件、数据库索引，节约大家下载时间，节省数据库编制索引的计算资源消耗。

# humann3为例
mkdir -p ~/db/humann3 && cd ~/db/humann3
site=ftp://download.nmdc.cn/tools
wget -c ${site}/humann3/full_chocophlan.v296_201901.tar.gz
wget -c ${site}/humann3/uniref90_annotated_v201901.tar.gz
wget -c ${site}/humann3/full_mapping_v201901.tar.gz
mkdir -p chocophlan uniref utility_mapping
tar xvzf full_chocophlan.v296_201901.tar.gz -C chocophlan/
tar xvzf uniref90_annotated_v201901.tar.gz -C uniref/
tar xvzf full_mapping_v201901.tar.gz -C utility_mapping/

百度云备份链接

https://pan.baidu.com/s/1Ikd_47HHODOqC3Rcx6eJ6Q?pwd=0315

下载的tar.gz压缩包，可放置于指定目录，使用tar -xvzf \*.tar.gz解压

# 大文件的分卷压缩和解压 以kraken2为例
cd ~/db/kraken2
# https://www.cnblogs.com/wang--lei/p/9046643.html
# 文件夹kraken2/打包压缩，1h
tar -zcvf kraken2.tar.gz kraken2/
# b分割为指定大小文件G/M/K，-d数字，a序列长度，输入和输出前缀
split -b 13G -d -a 1 kraken2.tar.gz kraken2.tar.gz.
# 一行命令打包并分割
tar -zcvf kraken2.tar.gz kraken2idx/ | split -b 19G -d -a 1 - kraken2.tar.gz.
# 分割后合并及解压缩
cat kraken2.tar.gz.* | tar -zxv

kneaddata常见问题

kneaddata运行提示java版本不支持

# 解决思路，新建虚拟环境，安装kneaddata，再安装对应java版本
# 务必指定2.7，软件依赖2.7的python，但conda会自动安装3.6，运行报错
conda create -n kneaddata python=2.7
conda activate kneaddata
conda install openjdk=8.0.152
conda install kneaddata=0.6.1

解压失败-重新下载再安装

tar -xvzf kneaddata.tar.gz -C ~/miniconda3/envs/kneaddata

解压文件提示如下错误

gzip: stdin: invalid compressed data--format violated
tar: Unexpected EOF in archive
tar: Unexpected EOF in archive
tar: Error is not recoverable: exiting now

检查md5值确认文件是否不同

md5sum kneaddata.tar.gz

当前为d26125bee1def1faa99d03a9715bf392
原文件为9fa47a364096b2c33be52a91850b2cde

删除当前文件并重新下载即可

rm kneaddata.tar.gz
wget ftp://download.nmdc.cn/tools//conda/kneaddata.tar.gz

Lefse在Rstudio中运行命令调用R版本问题的解决

# 在Rstudio中默认调用Rstudio的R，具体写在/etc/rstudio/rserver.conf
# 或在R中用Sys.getenv()["R_HOME"]，在rpy2中print(robjects.r)可以查看其调用的r版本
# 指定lefse调用的R版本，需根据conda实际目录修改
sed -i "2 i os.environ['R_HOME'] = '~/miniconda3/envs/meta/lib/R/'" \
  ~/miniconda3/envs/meta/share/lefse-1.0.8.post1-1/lefse.py

Kraken2

定制数据库

官方教程详见 https://github.com/DerrickWood/kraken2/blob/master/docs/MANUAL.markdown

本地构建最完整索引，自定义微生物数据库，如标准+真菌+原生动物+质粒+植物

mkdir -p ${db}/kraken2/kraken2_self
conda activate kraken2
# 显示帮助
kraken2-build -h
# 下载物种注释
kraken2-build --download-taxonomy --threads 24 --db ${db}/kraken2/kraken2_self
# 下载数据库，需要12-24小时
for i in archaea bacteria UniVec_Core viral human fungi plasmid protozoa plant; do
    kraken2-build --download-library $i --threads 24 --db ${db}/kraken2/kraken2_self
done
# 确定的库建索引，4p,4h
time kraken2-build --build --threads 48 --db ${db}/kraken2/kraken2_self
# bracken索引，长度推荐100/150, 24p, 1h;
time bracken-build -d ./ -t 24 -k 35 -l 100
time bracken-build -d ./ -t 24 -k 35 -l 150

Perl版本不对

常见问题：Perl版本不对，人工指定perl版本如下

PERL5LIB=~/miniconda3/envs/kraken2/lib/site_perl/5.26.2/x86_64-linux-thread-multi:~/miniconda3/envs/kraken2/lib/site_perl/5.26.2:~/miniconda3/envs/kraken2/lib/5.26.2/x86_64-linux-thread-multi:~/miniconda3/envs/kraken2/lib/5.26.2

salmon手动安装和使用

# 如不可用，尝试下载二进制和添加环境变量
wget https://github.com/COMBINE-lab/salmon/releases/download/v0.14.0/salmon-0.14.0_linux_x86_64.tar.gz
tar xvzf salmon-0.14.0_linux_x86_64.tar.gz 
cp -rf salmon-latest_linux_x86_64/ ${soft}/envs/metagenome_env/share/salmon
# 或者直接使用软件全路径
${soft}/envs/metagenome_env/share/salmon/bin/salmon -v # 0.14.0

MetaWRAP分箱

shorten_contig_names.py报错

更新 ${soft}/envs/metawrap/bin/metawrap-scripts/shorten_contig_names.py 脚本

#!/usr/bin/env python2.7
import sys
shorten=False
for line in open(sys.argv[1]):
    if line[0]!=">":
        print line.rstrip()
    else:
        if shorten==True:
            #print "_".join(line.rstrip().split("_")[:4])
            lineL = line.rstrip().split("_")
            new_line = '_'.join([lineL[0], lineL[1], lineL[3]])
            print new_line[:20]
        elif len(line)>20 and len(line.split("_"))>5:
            lineL = line.rstrip().split("_")
            new_line = '_'.join([lineL[0], lineL[1], lineL[3]])
            #print "_".join(line.rstrip().split("_")[:4])
            print new_line[:20]
            shorten=True
        else:
            print line.rstrip()

绘图plot_binning_results.py报错

更新 ${soft}/envs/metawrap/bin/metawrap-scripts/plot_binning_results.py 脚本

# 原脚本存在嵌套错误，会输出报错，修改部分如下
# Traceback (most recent call last):
#  File "/anaconda3/envs/metawrap-env/bin/metawrap-scripts/plot_binning_results.py", line 119, in <module>
#   plt.text(x_pos, y_pos, bin_set, fontsize=18, color=c)
# NameError: name 'x_pos' is not defined

    # add bin set label to plot
    for x_pos,y in enumerate(data[bin_set]):
            if y>y_pos:
                    break
            plt.text(x_pos, y_pos, bin_set, fontsize=18, color=c)
            y_pos+=y_increment
            
# add plot and axis titles and adjust the edges
plt.title("Bin contamination ranking", fontsize=26)
plt.xlabel("Acending contamination rank", fontsize=16)
plt.ylabel("Estimated bin contamination (log scale)", fontsize=16)
plt.gcf().subplots_adjust(right=0.9)

# save figure
print "Saving figures binning_results.eps and binning_results.png ..."
plt.tight_layout(w_pad=10)
plt.subplots_adjust(top=0.92, right=0.90, left=0.08)
plt.savefig("binning_results.png",format='png', dpi=300)
plt.savefig("binning_results.eps",format='eps')
#plt.show()
EOF

blast版本不兼容

更新 metawrap 中的 blast 版本，直接到https://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/，下载最新版本blastn，再到conda的metawrap环境bin目录下，替换掉旧版本的blastn

# 如果出现这个错误，BLAST Database error: Error: Not a valid version 4 database.
# 是metawrap 中 blast 版本太老了，需要更新下
wget -c https://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/ncbi-blast-2.13.0+-x64-linux.tar.gz
tar xvzf ncbi-blast-2.13.0+-x64-linux.tar.gz
mv ncbi-blast-2.13.0+/bin/* ${soft}/envs/metawrap-env/bin/

附录

Conda安装小工具

以下小工具已经整合至EasyMicrobiome项目中的linux文件夹，以下代码提供学习多种自主安装的参考方法，用于积累conda使用

并行计算管理rush/paprllel

# conda安装rush，无依赖关系更好用的并行工具
conda install rush -c bioconda
# Ubuntu下安装方法 apt install parallel
# conda安装parallel，版本有点老
conda install parallel -c bioconda
parallel --version # GNU parallel 20170422

表格统计工具csvtk和序列处理seqkit(可选中)

# 方法1. conda安装，可能有点旧
conda install csvtk -c bioconda
conda install seqkit -c bioconda

# 方法2. 直接下载最新版 https://github.com/shenwei356，如以csvtk为例手动安装
wget -c https://github.com/shenwei356/csvtk/releases/download/v0.22.0/csvtk_linux_amd64.tar.gz
tar xvzf csvtk_linux_amd64.tar.gz
cp csvtk ~/miniconda3/bin/

宿主参考基因组下载

EnsembleGenomes http://ensemblgenomes.org/
包括动物、植物、原生生物、真菌、细菌等，此外植物还 Phytozome https://phytozome-next.jgi.doe.gov/ ，以及单个物种和专用数据库

以Ensemble中拟南芥为例：Arabidopsis thaliana – Genome assembly – Download DNA sequence (无反应)，点TAIR链接跳转ENA，下载All Seq FASTA

wget https://www.ebi.ac.uk/ena/browser/api/fasta/GCA_000001735.1?download=true&gzip=true
mv GCA_000001735.1?download=true TAIR10.fa

以Ensemble中水稻为例：Oryza sativa Japonica —— IRGSP-1.0

wget https://www.ebi.ac.uk/ena/browser/api/fasta/GCA_001433935.1?download=true&gzip=true
mv GCA_001433935.1?download=true IRGSP1.0.fa

KEGG层级注释整理

己整合至EasyMicrobiome中，自己更新请访问 https://www.kegg.jp/kegg-bin/show_brite?ko00001.keg 下载htext

# 转换ABCD为列表
kegg_ko00001_htext2tsv.pl -i ko00001.keg -o ko00001.tsv
# 统计行数，2021.1月版55761行，整理后为55103个条目
wc -l ko00001.*
# 统计各级数量, /54/527/23917
for i in `seq 1 2 8`;do
    cut -f ${i} ko00001.tsv|sort|uniq|wc -l ; done
# 生成KO编号和注释列表
cut -f 7,8 ko00001.tsv|sort|uniq|sed '1 i KO\tDescription' \
  > KO_description.txt
# KO与通路(Pathway)对应表，用于合并D级为C级
awk 'BEGIN{FS=OFS="\t"} {print $7,$6}' ko00001.tsv | sed '1 i KO\tpathway' \
  > KO_path.list

毒力因子数据库VFDB

官网：http://www.mgc.ac.cn/VFs/ 数据每周更新

mkdir -p ${db}/vfdb && cd ${db}/vfdb
# 毒力因子描述文件
wget -c http://www.mgc.ac.cn/VFs/Down/VFs.xls.gz
# 核心数据库(1117K)
wget -c http://www.mgc.ac.cn/VFs/Down/VFDB_setA_pro.fas.gz
# 完整数据库(4.99M)
wget -c http://www.mgc.ac.cn/VFs/Down/VFDB_setB_pro.fas.gz
# 解压
gunzip *.gz

宏基因组基于读长的分析 HUMAnN2/Metaphlan2/graphlan

HUMAnN2解包安装

# 下载
wget -c ftp://download.nmdc.cn/tools/conda/humann2.tar.gz
# 指定安装目录
mkdir -p ${soft}/envs/humann2
tar -xvzf humann2.tar.gz -C ${soft}/envs/humann2
# 启动环境
conda activate humann2
# 初始化环境
conda unpack

HUMAnN2直接安装

# mamba 是快速版本的 conda
mamba create -n humann2 humann2 graphlan export2graphlan -c bioconda -y

HUMAnN2安装测试

conda activate humann2
# 记录核心软件版本
humann2 --version # v2.8.1
metaphlan2.py -v # 2.7.5 (6 February 2018)
diamond help | head -n 1 #  v0.8.36.98
graphlan.py --version # 1.1.3 (5 June 2018)
export2graphlan.py -h # 0.22 of 05 May

# 测试流程是否可用
humann2_test

HUMAnN2物种和功能数据库

# 显示可用分类、泛基因组和功能数据库
humann2_databases

# 安装数据库(注：数据库下载慢或失败，附录有国内备份链接)
cd ${db}
mkdir -p ${db}/humann2 # 建立下载目录
# 输助比对数据库 593MB
humann2_databases --download utility_mapping full ${db}/humann2
# 微生物泛基因组 5.37 GB
humann2_databases --download chocophlan full ${db}/humann2
# 功能基因diamond索引 10.3 GB
humann2_databases --download uniref uniref90_diamond ${db}/humann2

# humann2数据库无法下载：附录备用链接下载后手动配置
mkdir -p ${db}/humann2/chocophlan && cd ${db}/humann2/chocophlan
tar xvzf full_chocophlan_plus_viral.v0.1.1.tar.gz
mkdir -p ${db}/humann2/uniref && cd ${db}/humann2/uniref
tar xvzf uniref90_annotated_1_1.tar.gz
mkdir -p ${db}/humann2/utility_mapping && cd ${db}/humann2/utility_mapping
tar xvzf full_mapping_1_1.tar.gz

# 设置数据库位置
# 显示参数
humann2_config --print
# 如修改线程数，推荐3-8，根据实际情况调整
humann2_config --update run_modes threads 4
humann2_config --update database_folders utility_mapping ${db}/humann2/utility_mapping
humann2_config --update database_folders nucleotide ${db}/humann2/chocophlan
humann2_config --update database_folders protein ${db}/humann2/uniref
humann2_config --print

## metaphlan2数据库下载和配置
mkdir -p ${db}/humann2 && cd ${db}/humann2
wget -c ftp://download.nmdc.cn/tools/humann2/metaphlan2.tar.gz
tar xvzf metaphlan2.tar.gz
# 链接到软件安装目录
mkdir -p ${soft}/envs/humann2/bin/databases
ln -s ${db}/humann2/metaphlan2/* ${soft}/envs/humann2/bin/databases/

n=kneaddata
conda pack -f --ignore-missing-files -n ${n} -o ${n}.tar.gz

2、序列处理分析流程

易宏基因组流程 EasyMetagenome Pipeline

# 版本Version: 1.20, 2023/11/23
# 操作系统Operation System: Linux Ubuntu 22.04+ / CentOS 7.7+

一、数据预处理 Data preprocessing

1.1 准备工作 Preparing

首次使用请参照0Install.sh脚本，安装软件和数据库(大约1-3天，仅一次)
易宏基因组(EasyMetagenome)流程1Pipeline.sh复制到项目文件夹，如本次为meta
项目文件夹准备测序数据(seq/*.fq.gz)和样本元数据(result/metadata.txt)

**环境变量设置 Environment variable settings**
**分析前必须运行，设置数据库、软件和工作目录**

# Conda软件安装目录，`conda env list`查看，如/anaconda3
soft=~/miniconda3
# 数据库database(db)位置，如管理员/db，个人~/db
db=~/db
# 设置工作目录work directory(wd)，如meta
wd=~/meta
# 创建并进入工作目录
mkdir -p $wd && cd $wd
# 创建3个常用子目录：序列，临时文件和结果
mkdir -p seq temp result
# 添加分析所需的软件、脚本至环境变量，添加至~/.bashrc中自动加载
PATH=$soft/bin:$soft/condabin:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin:$db/EasyMicrobiome/linux:$db/EasyMicrobiome/script
echo $PATH

**元数据和序列文件 Metadata and Sequence Files**

元数据

# 上传元数据metadata.txt至result目录，此处下载并重命名
wget http://www.imeta.science/github/EasyMetagenome/result/metadata.txt
mv metadata.txt result/metadata.txt
# 检查文件格式，^I为制表符，$为Linux换行，^M$为Windows回车，^M为Mac换行符
cat -A result/metadata.txt

# 根据样本文件生成元数据，可筛选子集，如EB开头
ls seq/EB*|grep '_1'|cut -f1 -d '_'|cut -f 2 -d '/'|sed'1 i SampleID'>result/metadataEB.txt
cp result/metadataEB.txt result/metadata.txt

# 元数据细节优化
# 转换Windows回车为Linux换行
sed -i 's/\r//' result/metadata.txt
# 去除数据中的一个多余空格
sed -i 's/Male  /Male/' result/metadata.txt
cat -A result/metadata.txt

序列文件

# 用户使用filezilla上传测序文件至seq目录，本次从网络下载
# seq 目录下已经有测试文件，下载跳过
cd seq/
awk '{system("wget -c http://www.imeta.science/github/EasyMetagenome/seq/"$1"_1.fq.gz")}' <(tail -n+2 ../result/metadata.txt)
awk '{system("wget -c http://www.imeta.science/github/EasyMetagenome/seq/"$1"_2.fq.gz")}' <(tail -n+2 ../result/metadata.txt)
cd ..
# ls查看文件大小，-l 列出详细信息 (l: list)，-sh 显示人类可读方式文件大小 (s: size; h: human readable)
ls -lsh seq/*.fq.gz

序列文件格式检查
zless/zcat查看可压缩文件，检查序列质量格式(质量值大写字母为标准Phred33格式，小写字母为Phred64，需参考附录：质量值转换)；检查双端序列ID是否重名，如重名需要改名。参考**附录 —— 质控kneaddata，去宿主后双端不匹配；序列改名**。

# 设置某个样本名为变量i，以后再无需修改
i=C1
# zless查看压缩文件，空格翻页，q退出; head指定显示行数
zless seq/${i}_1.fq.gz | head -n4

**工作目录和文件结构总结**

# ├── pipeline.sh
# ├── result
# │   └── metadata.txt
# ├── seq
# │   ├── C1_1.fq.gz
# │   ├── ...
# │   └── N1_2.fq.gz
# └── temp

* 1pipeline.sh是分析流程代码；
* seq目录中有2个样本Illumina双端测序，4个序列文件；
* temp是临时文件夹，存储分析中间文件，结束可全部删除节约空间
* result是重要节点文件和整理化的分析结果图表，
* 实验设计metadata.txt也在此

1.2 Fastp质量控制 Quality Control

# 创建目录，记录软件版本和引文
mkdir -p temp/qc result/qc
fastp

# 单样本质控
i=C1
fastp -i seq/${i}_1.fq.gz  -I seq/${i}_2.fq.gz \
  -o temp/qc/${i}_1.fastq -O temp/qc/${i}_2.fastq

# 多样本并行
# -j 2: 表示同时处理2个样本
time tail -n+2 result/metadata.txt|cut -f1|rush -j 2 \
  "fastp -i seq/{1}_1.fq.gz -I seq/{1}_2.fq.gz \
    -j temp/qc/{1}_fastp.json -h temp/qc/{1}_fastp.html \
    -o temp/qc/{1}_1.fastq  -O temp/qc/{1}_2.fastq \
    > temp/qc/{1}.log 2>&1 "

# 质控后结果汇总
echo -e "SampleID\tRaw\tClean" > temp/fastp
for i in `tail -n+2 result/metadata.txt|cut -f1`;do
    echo -e -n "$i\t" >> temp/fastp
    grep 'total reads' temp/qc/${i}.log|uniq|cut -f2 -d ':'|tr '\n' '\t' >> temp/fastp
    echo "" >> temp/fastp
done
sed -i 's/ //g;s/\t$//' temp/fastp
# 按metadata排序
awk 'BEGIN{FS=OFS="\t"}NR==FNR{a[$1]=$0}NR>FNR{print a[$1]}' temp/fastp result/metadata.txt \
  > result/qc/fastp.txt
cat result/qc/fastp.txt

1.3 KneadData去宿主 Host removal

kneaddata是流程主要依赖bowtie2比对宿主，然后筛选非宿主序列用于下游分析。

# 创建目录、启动环境、记录版本
mkdir -p temp/hr
conda activate kneaddata
kneaddata --version # 0.12.0

多样品并行去宿主，16p 4h

time tail -n+2 result/metadata.txt|cut -f1|rush -j 2 \
  "sed '1~4 s/ 1:/.1:/;1~4 s/$/\/1/' temp/qc/{}_1.fastq > /tmp/{}_1.fastq; \
  sed '1~4 s/ 2:/.1:/;1~4 s/$/\/2/' temp/qc/{}_2.fastq > /tmp/{}_2.fastq; \
  kneaddata -i1 /tmp/{1}_1.fastq -i2 /tmp/{1}_2.fastq \
  -o temp/hr --output-prefix {1} \
  --bypass-trim --bypass-trf --reorder \
  --bowtie2-options '--very-sensitive --dovetail' \
  -db ${db}/kneaddata/human/hg37dec_v0.1 \
  --remove-intermediate-output -v -t 3; \
  rm /tmp/{}_1.fastq /tmp/{}_2.fastq"

# *  匹配任意多个字符，？ 匹配任意一个字符
ls -shtr temp/hr/*_paired_?.fastq

简化改名

# Ubuntu系统改名
rename 's/paired_//' temp/hr/*.fastq
# CentOS系统改名
rename 'paired_' '' temp/hr/*.fastq

大文件清理，高宿主含量样本可节约>90%空间

# 使用命令的绝对路径确保使用无参数的命令，管理员用alias自定义命令含参数，影响操作结果
/bin/rm -rf temp/hr/*contam* temp/hr/*unmatched* temp/hr/reformatted* temp/hr/_temp*
ls -l temp/hr/

质控结果汇总

kneaddata_read_count_table --input temp/hr \
  --output temp/kneaddata.txt
# 筛选重点结果列
cut -f 1,2,5,6 temp/kneaddata.txt | sed 's/_1_kneaddata//' > result/qc/sum.txt
# 对齐方式查看表格
csvtk -t pretty result/qc/sum.txt

二、基于读长分析 Read-based (HUMAnN3+MetaPhlAn4+Kraken2)

2.1 准备HUMAnN输入文件

HUMAnN要求双端序列合并的文件作为输入，for循环根据实验设计样本名批量双端序列合并。注意星号(\*)和问号(?)，分别代表多个和单个字符。当然大家更不能溜号，行分割的代码行末有一个\

mkdir -p temp/concat
# 双端合并为单个文件
for i in `tail -n+2 result/metadata.txt|cut -f1`;do 
  cat temp/hr/${i}_?.fastq \
  > temp/concat/${i}.fq; done
# 查看样品数量和大小
ls -shl temp/concat/*.fq
# 数据太大，计算时间长，可用head对单端分析截取20M序列，即3G，行数为80M行，详见附录：HUMAnN2减少输入文件加速

2.2 HUMAnN计算物种和功能组成

* 物种组成调用MetaPhlAn4
* 输入文件：temp/concat/*.fq 每个样品质控后双端合并后的fastq序列
* 输出文件：temp/humann3/ 目录下
* C1_pathabundance.tsv
* C1_pathcoverage.tsv
* C1_genefamilies.tsv
* 整合后的输出：
* result/metaphlan4/taxonomy.tsv 物种丰度表
* result/metaphlan4/taxonomy.spf 物种丰度表（用于stamp分析）
* result/humann3/pathabundance_relab_unstratified.tsv 通路丰度表
* result/humann3/pathabundance_relab_stratified.tsv 通路物种组成丰度表
* stratified(每个菌对此功能通路组成的贡献)和unstratified(功能组成)

启动humann3环境，检查数据库配置

conda activate humann3
# 备选source加载指定环境
# source ~/miniconda3/envs/humann3/bin/activate
mkdir -p temp/humann3
humann --version # v3.7
humann_config

单样本1.25M PE150运行测试，8p，2.5M，1~2h；0.2M, 34m；0.1M，30m；0.01M，25m；16p，18m

i=C1

3p,26m; 数据库使用ssd缩短到19m；16p,8m

time humann --input temp/concat/${i}.fq --output temp/humann3 --threads 3 --metaphlan-options '--bowtie2db /db/metaphlan4 --index mpa_vOct22_CHOCOPhlAnSGB_202212 --offline'

多样本并行计算，测试数据约30m，推荐16p，3小时/样本

# 如果服务器性能好，请设置--threads值为8/16/32
tail -n+2 result/metadata.txt | cut -f1 | rush -j 2 \
  "humann --input temp/concat/{1}.fq  \
  --output temp/humann3/ --threads 3 --metaphlan-options '--bowtie2db /db/metaphlan4 --index mpa_vOct22_CHOCOPhlAnSGB_202212 --offline'"

# 移动重要文件至humann3目录
# $(cmd) 与 `cmd` 通常是等价的；`cmd`写法更简单，但要注意反引号是键盘左上角ESC下面的按键，$(cmd)更通用，适合嵌套使用
for i in $(tail -n+2 result/metadata.txt | cut -f1); do  
   mv temp/humann3/${i}_humann_temp/${i}_metaphlan_bugs_list.tsv temp/humann3/
done
# 删除临时文件，极占用空间
/bin/rm -rf temp/concat/* temp/humann3/*_humann_temp

(可选)单独运行MetaPhlAn4

mkdir -p temp/humann3
i=C1
# 仅物种注释极快4p, 2m, 1m读取数据库
time metaphlan --input_type fastq temp/qc/${i}_1.fastq \
  temp/humann3/${i}.txt --bowtie2db /db/metaphlan4 --index mpa_vOct22_CHOCOPhlAnSGB_202212 --offline \
  --nproc 4

2.3 物种组成表

**样品结果合并**

mkdir -p result/metaphlan4
# 合并、修正样本名、预览
merge_metaphlan_tables.py temp/humann3/*_metaphlan_bugs_list.tsv | \
  sed 's/_metaphlan_bugs_list//g' | tail -n+2 | sed '1 s/clade_name/ID/' | sed '2i #metaphlan4'> result/metaphlan4/taxonomy.tsv
csvtk -t stat result/metaphlan4/taxonomy.tsv
head -n5 result/metaphlan4/taxonomy.tsv

**转换为stamp的spf格式**

# metaphlan4较2增加更多unclassified和重复结果，用sort和uniq去除
metaphlan_to_stamp.pl result/metaphlan4/taxonomy.tsv \
  |sort -r | uniq > result/metaphlan4/taxonomy.spf
head result/metaphlan4/taxonomy.spf
# STAMP不支持unclassified，需要过滤掉再使用
grep -v 'unclassified' result/metaphlan4/taxonomy.spf > result/metaphlan4/taxonomy2.spf
# 下载metadata.txt和taxonomy2.spf使用stamp分析

2.4 功能组成分析

功能组成样本合并合并

mkdir -p result/humann3
humann_join_tables --input temp/humann3 \
  --file_name pathabundance \
  --output result/humann3/pathabundance.tsv
# 样本名调整：删除列名多余信息
sed -i 's/_Abundance//g' result/humann3/pathabundance.tsv
# 统计和预览
csvtk -t stat result/humann3/pathabundance.tsv
head -n5 result/humann3/pathabundance.tsv

标准化为相对丰度relab(1)或百万比cpm(1,000,000)

humann_renorm_table \
  --input result/humann3/pathabundance.tsv \
  --units relab \
  --output result/humann3/pathabundance_relab.tsv
head -n5 result/humann3/pathabundance_relab.tsv

分层结果：功能和对应物种表(stratified)和功能组成表(unstratified)

humann_split_stratified_table \
  --input result/humann3/pathabundance_relab.tsv \
  --output result/humann3/

差异比较和柱状图

两样本无法组间比较，在pcl层面替换为HMP数据进行统计和可视化。

* 输入数据：通路丰度表格 result/humann3/pathabundance.tsv和实验设计 result/metadata.txt
* 中间数据：包含分组信息的通路丰度表格文件 result/humann3/pathabundance.pcl
* 输出结果：result/humann3/associate.txt

在通路丰度中添加分组

## 提取样品列表
head -n1 result/humann3/pathabundance.tsv | sed 's/# Pathway/SampleID/' | tr '\t' '\n' > temp/header
## 对应分组，本示例分组为第2列($2)，根据实际情况修改
awk 'BEGIN{FS=OFS="\t"}NR==FNR{a[$1]=$2}NR>FNR{print a[$1]}' result/metadata.txt temp/header | tr '\n' '\t'|sed 's/\t$/\n/' > temp/group
# 合成样本、分组+数据
cat <(head -n1 result/humann3/pathabundance.tsv) temp/group <(tail -n+2 result/humann3/pathabundance.tsv) \
  > result/humann3/pathabundance.pcl
head -n5 result/humann3/pathabundance.pcl
tail -n5 result/humann3/pathabundance.pcl

组间比较，样本量少无差异，结果为4列的文件：通路名字，通路在各个分组的丰度，差异P-value，校正后的Q-value。
演示数据2样本无法统计，此处替换为HMP的结果演示统计和绘图(上传hmp\_pathabund.pcl，替换pathabundance.pcl为hmp\_pathabund.pcl)。

wget -c http://www.imeta.science/github/EasyMetagenome/result/humann2/hmp_pathabund.pcl
/bin/cp -f hmp_pathabund.pcl result/humann3/
# 设置输入文件名
pcl=result/humann3/hmp_pathabund.pcl
# 统计表格行、列数量
csvtk -t stat ${pcl}
head -n3 ${pcl} | cut -f 1-5
# 按分组KW检验，注意第二列的分组列名
humann_associate --input ${pcl} \
    --focal-metadatum Group --focal-type categorical \
    --last-metadatum Group --fdr 0.05 \
    --output result/humann3/associate.txt
wc -l result/humann3/associate.txt
head -n5 result/humann3/associate.txt

barplot展示通路的物种组成，如：腺苷核苷酸合成

# 指定差异通路，如 P163-PWY，--sort sum metadata 按丰度和分组排序
path=P163-PWY
humann_barplot \
    --input ${pcl} --focal-feature ${path} \
    --focal-metadata Group --last-metadata Group \
    --output result/humann3/barplot_${path}.pdf --sort sum metadata

KEGG注释

支持GO、PFAM、eggNOG、level4ec、KEGG的D级KO等注释，详见humann\_regroup\_table -h。

# 转换基因家族为KO(uniref90_ko)，可选eggNOG(uniref90_eggnog)或酶(uniref90_level4ec)
for i in `tail -n+2 result/metadata.txt|cut -f1`;do
  humann_regroup_table \
    -i temp/humann3/${i}_genefamilies.tsv \
    -g uniref90_ko \
    -o temp/humann3/${i}_ko.tsv
done
# 合并，并修正样本名
humann_join_tables \
  --input temp/humann3/ \
  --file_name ko \
  --output result/humann3/ko.tsv
sed -i '1s/_Abundance-RPKs//g' result/humann3/ko.tsv
tail result/humann3/ko.tsv
# 与pathabundance类似，可进行标准化renorm、分层stratified、柱状图barplot等操作

# 分层结果：功能和对应物种表(stratified)和功能组成表(unstratified)
humann_split_stratified_table \
  --input result/humann3/ko.tsv \
  --output result/humann3/ 
wc -l result/humann3/ko*

# KO合并为高层次L2, L1通路代码KO to level 1/2/3
summarizeAbundance.py \
  -i result/humann3/ko_unstratified.tsv \
  -m ${db}/EasyMicrobiome/kegg/KO1-4.txt \
  -c 2,3,4 -s ',+,+,' -n raw \
  -o result/humann3/KEGG
wc -l result/humann3/KEGG*

2.5 GraPhlAn图

metaphlan2 to graphlan

conda activate humann2
export2graphlan.py --skip_rows 1,2 -i result/metaphlan4/taxonomy.tsv \
  --tree temp/merged_abundance.tree.txt \
  --annotation temp/merged_abundance.annot.txt \
  --most_abundant 1000 --abundance_threshold 20 --least_biomarkers 10 \
  --annotations 3,4 --external_annotations 7
# 参数说明见PPT，或运行 export2graphlan.py --help
# graphlan annotation
graphlan_annotate.py --annot temp/merged_abundance.annot.txt \
  temp/merged_abundance.tree.txt  temp/merged_abundance.xml
# output PDF figure, annoat and legend
graphlan.py temp/merged_abundance.xml result/metaphlan4/graphlan.pdf \
  --external_legends 
# GraPhlAn Plot(测试中)
graphlan_plot.r --input result/metaphlan4/taxonomy.spf \
  --design result/metadata.txt --number 100 \
  --group all --type heatmap \
  --output result/metaphlan4/heatmap

2.6 LEfSe差异分析物种

* 输入文件：物种丰度表result/metaphlan2/taxonomy.tsv
* 输入文件：样品分组信息 result/metadata.txt
* 中间文件：整合后用于LefSe分析的文件 result/metaphlan2/lefse.txt，这个文件可以提供给www.ehbio.com/ImageGP 用于在线LefSE分析
* LefSe结果输出：result/metaphlan2/目录下lefse开头和feature开头的文件

前面演示数据仅有2个样本，无法进行差异比较。下面使用result12目录中由12个样本生成的结果表进行演示

# 设置结果目录，自己的数据使用result，演示用result12
result=result12
# 如果没有，请下载演示数据
wget -c http://www.imeta.science/db/EasyMetagenome/result12.zip
unzip result12.zip

准备输入文件，修改样本品为组名(可手动修改)

# 提取样本行替换为每个样本一行，修改ID为SampleID
head -n1 $result/metaphlan2/taxonomy.tsv|tr '\t' '\n'|sed '1 s/ID/SampleID/' >temp/sampleid
head -n3 temp/sampleid
# 提取SampleID对应的分组Group(假设为metadata.txt中第二列$2)，替换换行\n为制表符\t，再把行末制表符\t替换回换行
awk 'BEGIN{OFS=FS="\t"}NR==FNR{a[$1]=$2}NR>FNR{print a[$1]}' $result/metadata.txt temp/sampleid|tr '\n' '\t'|sed 's/\t$/\n/' >groupid
cat groupid
# 合并分组和数据(替换表头)
cat groupid <(tail -n+2 $result/metaphlan2/taxonomy.tsv) > $result/metaphlan2/lefse.txt
head -n3 $result/metaphlan2/lefse.txt

方法1. 推荐在线 https://www.bic.ac.cn/ImageGP/ 中LEfSe一键分析

方法2. LEfSe命令行分析

conda activate lefse
result=result12
# 格式转换为lefse内部格式
lefse-format_input.py $result/metaphlan2/lefse.txt \
  temp/input.in -c 1 -o 1000000
# 运行lefse(样本必须有重复和分组)
run_lefse.py temp/input.in temp/input.res

# 绘制物种树注释差异
lefse-plot_cladogram.py temp/input.res \
  $result/metaphlan2/lefse_cladogram.pdf --format pdf

# 绘制所有差异features柱状图
lefse-plot_res.py temp/input.res \
  $result/metaphlan2/lefse_res.pdf --format pdf
    
# 绘制单个features柱状图
# 查看显著差异features，按丰度排序
grep -v '-' temp/input.res | sort -k3,3n 
# 手动选择指定feature绘图，如Firmicutes
lefse-plot_features.py -f one --format pdf \
  --feature_name "k__Bacteria.p__Firmicutes" \
  temp/input.in temp/input.res \
  $result/metaphlan2/lefse_Firmicutes.pdf

# 批量绘制所有差异features柱状图
lefse-plot_features.py -f diff \
  --archive none --format pdf \
  temp/input.in temp/input.res \
  $result/metaphlan2/lefse_

2.7 Kraken2+Bracken物种注释和丰度估计

Kraken2可以快速完成读长(read)层面的物种注释和定量，还可以进行重叠群(contig)、基因(gene)、宏基因组组装基因组(MAG/bin)层面的序列物种注释。

# 启动kraken2工作环境
conda activate kraken2
# 记录软件版本
kraken2 --version # 2.1.2
mkdir -p temp/kraken2

Kraken2物种注释

输入：temp/qc/{1}_?.fastq 质控后的数据，{1}代表样本名；
参考数据库：-db ${db}/kraken2/pluspfp16g/
输出结果：每个样本单独输出，temp/kraken2/中的{1}_report和{1}_output
整合物种丰度表输出结果：result/kraken2/taxonomy_count.txt

(可选) 单样本注释，5m，50G大数据库较5G库注释比例提高10~20%。以C1为例，在2023/3/14版中，8g: 31.75%; 16g: 52.35%; 150g: 71.98%；同为16g，2023/10/9版本为63.88%

# 根据电脑内存由小到大选择下面3个数据库
# pluspf16g/pluspfp(55G)/plusppfp(120G)
i=C1
time kraken2 --db ${db}/kraken2/pluspf16g/ \
  --paired temp/qc/${i}_?.fastq \
  --threads 2 --use-names --report-zero-counts \
  --report temp/kraken2/${i}.report \
  --output temp/kraken2/${i}.output

多样本并行生成report，1样本8线程逐个运行，内存大但速度快，不建议用多任务并行

for i in `tail -n+2 result/metadata.txt | cut -f1`;do
  kraken2 --db ${db}/kraken2/pluspf16g \
  --paired temp/qc/${i}_?.fastq \
  --threads 2 --use-names --report-zero-counts \
  --report temp/kraken2/${i}.report \
  --output temp/kraken2/${i}.output; done

使用krakentools转换report为mpa格式

for i in `tail -n+2 result/metadata.txt | cut -f1`;do
  kreport2mpa.py -r temp/kraken2/${i}.report \
    --display-header -o temp/kraken2/${i}.mpa; done

合并样本为表格

mkdir -p result/kraken2
# 输出结果行数相同，但不一定顺序一致，要重新排序
tail -n+2 result/metadata.txt | cut -f1 | rush -j 1 \
  'tail -n+2 temp/kraken2/{1}.mpa | LC_ALL=C sort | cut -f 2 | sed "1 s/^/{1}\n/" > temp/kraken2/{1}_count '
# 提取第一样本品行名为表行名
header=`tail -n 1 result/metadata.txt | cut -f 1`
echo $header
tail -n+2 temp/kraken2/${header}.mpa | LC_ALL=C sort | cut -f 1 | \
  sed "1 s/^/Taxonomy\n/" > temp/kraken2/0header_count
head -n3 temp/kraken2/0header_count
# paste合并样本为表格
ls temp/kraken2/*count
paste temp/kraken2/*count > result/kraken2/tax_count.mpa
# 检查表格及统计
csvtk -t stat result/kraken2/tax_count.mpa
head -n 5 result/kraken2/tax_count.mpa

Bracken丰度估计

参数简介：

* -d为数据库，-i为输入kraken2报告文件
* r是读长，此处为100，通常为150，o输出重新估计的值
* l为分类级，可选域D、门P、纲C、目O、科F、属G、种S级别丰度估计
* t是阈值，默认为0，越大越可靠，但可用数据越少

循环重新估计每个样品的丰度，请修改tax分别重新计算P和S各1次

# 设置估算的分类级别D,P,C,O,F,G,S，常用门P和种S
# 测序6G起样本-t 10过滤低丰度物种
tax=S
mkdir -p temp/bracken
for i in `tail -n+2 result/metadata.txt | cut -f1`;do
    # i=C1
    bracken -d ${db}/kraken2/pluspf16g/ \
      -i temp/kraken2/${i}.report \
      -r 100 -l ${tax} -t 0 \
      -o temp/bracken/${i}.brk \
      -w temp/bracken/${i}.report; done
# bracken结果合并成表: 需按表头排序，提取第6列reads count，并添加样本名
tail -n+2 result/metadata.txt | cut -f1 | rush -j 1 \
  'tail -n+2 temp/bracken/{1}.brk | LC_ALL=C sort | cut -f6 | sed "1 s/^/{1}\n/" \
  > temp/bracken/{1}.count'
# 提取第一样本品行名为表行名
h=`tail -n1 result/metadata.txt|cut -f1`
tail -n+2 temp/bracken/${h}.brk | sort | cut -f1 | \
  sed "1 s/^/Taxonomy\n/" > temp/bracken/0header.count
# 检查文件数，为n+1
ls temp/bracken/*count | wc
# paste合并样本为表格，并删除非零行
paste temp/bracken/*count > result/kraken2/bracken.${tax}.txt
# 统计行列，默认去除表头
csvtk -t stat result/kraken2/bracken.${tax}.txt
# 按频率过滤，-r可标准化，-e过滤(microbiome_helper)
Rscript ${db}/EasyMicrobiome/script/filter_feature_table.R \
  -i result/kraken2/bracken.${tax}.txt \
  -p 0.01 \
  -o result/kraken2/bracken.${tax}.0.01
csvtk -t stat result/kraken2/bracken.${tax}.0.01

个性化结果筛选

# 门水平去除脊索动物(人)
grep 'Chordata' result/kraken2/bracken.P.0.01
grep -v 'Chordata' result/kraken2/bracken.P.0.01 > result/kraken2/bracken.P.0.01-H

# 按物种名手动去除宿主污染，以人为例(需按种水平计算相关结果)
# 种水平去除人类P:Chordata,S:Homo sapiens
grep 'Homo sapiens' result/kraken2/bracken.S.0.01
grep -v 'Homo sapiens' result/kraken2/bracken.S.0.01 \
  > result/kraken2/bracken.S.0.01-H

分析后清理每条序列的注释大文件

/bin/rm -rf temp/kraken2/*.output

多样性和可视化alpha多样性计算：Berger Parker’s (BP), Simpson’s (Si), inverse Simpson’s (ISi), Shannon’s (Sh)

# Fisher’s (F)依赖scipy.optimize包，默认未安装

mkdir -p result/kraken2
echo -e "SampleID\tBerger Parker\tSimpson\tinverse Simpson\tShannon" > result/kraken2/alpha.txt
for i in `tail -n+2 result/metadata.txt|cut -f1`;do
    echo -e -n "$i\t" >> result/kraken2/alpha.txt
    for a in BP Si ISi Sh;do
        alpha_diversity.py -f temp/bracken/${i}.brk -a $a | cut -f 2 -d ':' | tr '\n' '\t' >> result/kraken2/alpha.txt
    done
    echo "" >> result/kraken2/alpha.txt
done
cat result/kraken2/alpha.txt

beta多样性计算

beta_diversity.py -i temp/bracken/*.brk --type bracken \
  > result/kraken2/beta.txt
cat result/kraken2/beta.txt

Krona图

for i in `tail -n+2 result/metadata.txt|cut -f1`;do
    kreport2krona.py -r temp/bracken/${i}.report -o temp/bracken/${i}.krona --no-intermediate-ranks
    ktImportText temp/bracken/${i}.krona -o result/kraken2/krona.${i}.html
done

Pavian桑基图：https://fbreitwieser.shinyapps.io/pavian/ 在线可视化:，左侧菜单，Upload sample set (temp/bracken/*.report)，支持多样本同时上传；Sample查看结果，Configure Sankey配置图样式，Save Network下载图网页

多样性分析/物种组成，详见3StatPlot.sh，Kraken2结果筛选序列见附录

三、组装分析流程 Assemble-based

组装

# 启动工作环境
conda activate megahit

### MEGAHIT组装Assembly

# 删除旧文件夹，否则megahit无法运行
# /bin/rm -rf temp/megahit
# 组装，10~30m，TB级数据需几天至几周
megahit -t 3 \
    -1 `tail -n+2 result/metadata.txt|cut -f1|sed 's/^/temp\/hr\//;s/$/_1.fastq/'|tr '\n' ','|sed 's/,$//'` \
    -2 `tail -n+2 result/metadata.txt|cut -f1|sed 's/^/temp\/hr\//;s/$/_2.fastq/'|tr '\n' ','|sed 's/,$//'` \
    -o temp/megahit 
# 统计大小通常300M~5G，如果contigs太多，可以按长度筛选，降低数据量，提高基因完整度，详见附录megahit
seqkit stat temp/megahit/final.contigs.fa
# 预览重叠群最前6行，前60列字符
head -n6 temp/megahit/final.contigs.fa | cut -c1-60

# 备份重要结果
mkdir -p result/megahit/
ln -f temp/megahit/final.contigs.fa result/megahit/
# 删除临时文件
/bin/rm -rf temp/megahit/intermediate_contigs

(可选)metaSPAdes精细组装

# 精细但使用内存和时间更多，15~65m
/usr/bin/time -v -o metaspades.py.log metaspades.py -t 3 -m 100 \
  `tail -n+2 result/metadata.txt|cut -f1|sed 's/^/temp\/qc\//;s/$/_1.fastq/'|sed 's/^/-1 /'| tr '\n' ' '` \
  `tail -n+2 result/metadata.txt|cut -f1|sed 's/^/temp\/qc\//;s/$/_2.fastq/'|sed 's/^/-2 /'| tr '\n' ' '` \
  -o temp/metaspades
# 查看软件时间User time和内存Maximum resident set size
cat metaspades.py.log
# 2.3M，contigs体积更大
ls -sh temp/metaspades/contigs.fasta
seqkit stat temp/metaspades/contigs.fasta

# 备份重要结果
mkdir -p result/metaspades/
ln -f temp/metaspades/contigs.fasta result/metaspades/
# 删除临时文件
/bin/rm -rf temp/metaspades

注：metaSPAdes支持二、三代混合组装，见附录，此外还有OPERA-MS组装二、三代方案

QUAST评估

# QUAST评估，生成report文本tsv/txt、网页html、PDF等格式报告
quast.py result/megahit/final.contigs.fa \
  -o result/megahit/quast -t 2

(可选) megahit和metaspades比较

quast.py --label "megahit,metapasdes" \
    result/megahit/final.contigs.fa \
    result/metaspades/contigs.fasta \
    -o result/quast

(可选)metaquast评估，更全面，但需下载相关数据库，受网速影响可能时间很长(经常失败)

# metaquast based on silva, and top 50 species genome, 18min
time metaquast.py result/megahit/final.contigs.fa \
  -o result/megahit/metaquast

3.2 基因预测、去冗余和定量Gene prediction, cluster & quantitfy

metaProdigal基因预测Gene prediction

输入文件：组装的序列 result/megahit/final.contigs.fa

输出文件：prodigal预测的基因序列 temp/prodigal/gene.fa

基因大，可参考附录prodigal拆分基因文件，并行计算

mkdir -p temp/prodigal
# prodigal的meta模式预测基因，>和2>&1记录分析过程至gene.log
prodigal -i result/megahit/final.contigs.fa \
    -d temp/prodigal/gene.fa \
    -o temp/prodigal/gene.gff \
    -p meta -f gff > temp/prodigal/gene.log 2>&1 
# 查看日志是否运行完成，有无错误
tail temp/prodigal/gene.log
# 统计基因数量
seqkit stat temp/prodigal/gene.fa 
# 统计完整基因数量，数据量大可只用完整基因部分
grep -c 'partial=00' temp/prodigal/gene.fa 
# 提取完整基因(完整片段获得的基因全为完整，如成环的细菌基因组)
grep 'partial=00' temp/prodigal/gene.fa | cut -f1 -d ' '| sed 's/>//' > temp/prodigal/full_length.id
seqkit grep -f temp/prodigal/full_length.id temp/prodigal/gene.fa > temp/prodigal/full_length.fa
seqkit stat temp/prodigal/full_length.fa

cd-hit基因聚类/去冗余cluster & redundancy

输入文件：prodigal预测的基因序列 temp/prodigal/gene.fa

输出文件：去冗余后的基因和蛋白序列：result/NR/nucleotide.fa, result/NR/protein.fa

mkdir -p result/NR
# aS覆盖度，c相似度，G局部比对，g最优解，T多线程，M内存0不限制
# 2万基因2m，2千万需要2000h，多线程可加速
cd-hit-est -i temp/prodigal/gene.fa \
    -o result/NR/nucleotide.fa \
    -aS 0.9 -c 0.95 -G 0 -g 0 -T 0 -M 0
# 统计非冗余基因数量，单次拼接结果数量下降不大，多批拼接冗余度高
grep -c '>' result/NR/nucleotide.fa
# 翻译核酸为对应蛋白序列, --trim去除结尾的*
seqkit translate --trim result/NR/nucleotide.fa \
    > result/NR/protein.fa 
# 两批数据去冗余使用cd-hit-est-2d加速，见附录

salmon基因定量quantitfy

输入文件：去冗余后的基因序列：result/NR/nucleotide.fa

输出文件：Salmon定量：result/salmon/gene.count, gene.TPM

mkdir -p temp/salmon
salmon -v # 1.8.0

# 建索引, -t序列, -i 索引，10s
salmon index -t result/NR/nucleotide.fa \
  -p 3 -i temp/salmon/index 

# 定量，l文库类型自动选择，p线程，--meta宏基因组模式, 2个任务并行2个样
tail -n+2 result/metadata.txt | cut -f1 | rush -j 2 \
  "salmon quant -i temp/salmon/index -l A -p 3 --meta \
    -1 temp/qc/{1}_1.fastq -2 temp/qc/{1}_2.fastq \
    -o temp/salmon/{1}.quant"

合并

mkdir -p result/salmon
salmon quantmerge --quants temp/salmon/*.quant \
    -o result/salmon/gene.TPM
salmon quantmerge --quants temp/salmon/*.quant \
    --column NumReads -o result/salmon/gene.count
sed -i '1 s/.quant//g' result/salmon/gene.*

# 预览结果表格
head -n3 result/salmon/gene.*

3.3 功能基因注释Functional gene annotation

# 输入数据：上一步预测的蛋白序列 result/NR/protein.fa
# 中间结果：temp/eggnog/protein.emapper.seed_orthologs
#           temp/eggnog/output.emapper.annotations
#           temp/eggnog/output

# COG定量表：result/eggnog/cogtab.count
#            result/eggnog/cogtab.count.spf (用于STAMP)

# KO定量表：result/eggnog/kotab.count
#           result/eggnog/kotab.count.spf  (用于STAMP)

# CAZy碳水化合物注释和定量：result/dbcan3/cazytab.count
#                           result/dbcan3/cazytab.count.spf (用于STAMP)

# 抗生素抗性：result/resfam/resfam.count
#             result/resfam/resfam.count.spf (用于STAMP)

# 这部分可以拓展到其它数据库

eggNOG基因注释gene annotation(COG/KEGG/CAZy)

软件主页：https://github.com/eggnogdb/eggnog-mapper

# 运行并记录软件版本
conda activate eggnog
emapper.py --version
# emapper-2.1.7 / Expected eggNOG DB version: 5.0.2 
# Diamond version found: diamond version 2.0.15

# 运行emapper，18m，默认diamond 1e-3
mkdir -p temp/eggnog
time emapper.py --data_dir ${db}/eggnog \
  -i result/NR/protein.fa --cpu 3 -m diamond --override \
  -o temp/eggnog/output

# 格式化结果并显示表头
grep -v '^##' temp/eggnog/output.emapper.annotations | sed '1 s/^#//' \
  > temp/eggnog/output
csvtk -t headers -v temp/eggnog/output

生成COG/KO/CAZy丰度汇总表

mkdir -p result/eggnog
# 显示帮助
summarizeAbundance.py -h
# 汇总，7列COG_category按字母分隔，12列KEGG_ko和19列CAZy按逗号分隔，原始值累加
summarizeAbundance.py \
  -i result/salmon/gene.TPM \
  -m temp/eggnog/output --dropkeycolumn \
  -c '7,12,19' -s '*+,+,' -n raw \
  -o result/eggnog/eggnog
sed -i 's#^ko:##' result/eggnog/eggnog.KEGG_ko.raw.txt
sed -i '/^-/d' result/eggnog/eggnog*
head -n3 result/eggnog/eggnog*
# eggnog.CAZy.raw.txt  eggnog.COG_category.raw.txt  eggnog.KEGG_ko.raw.txt

# 添加注释生成STAMP的spf格式
awk 'BEGIN{FS=OFS="\t"} NR==FNR{a[$1]=$2} NR>FNR{print a[$1],$0}' \
  ${db}/EasyMicrobiome/kegg/KO_description.txt \
  result/eggnog/eggnog.KEGG_ko.raw.txt | \
  sed 's/^\t/Unannotated\t/' \
  > result/eggnog/eggnog.KEGG_ko.TPM.spf
head -n 5 result/eggnog/eggnog.KEGG_ko.TPM.spf
# KO to level 1/2/3
summarizeAbundance.py \
  -i result/eggnog/eggnog.KEGG_ko.raw.txt \
  -m ${db}/EasyMicrobiome/kegg/KO1-4.txt \
  -c 2,3,4 -s ',+,+,' -n raw --dropkeycolumn \
  -o result/eggnog/KEGG
head -n3 result/eggnog/KEGG*

# CAZy
awk 'BEGIN{FS=OFS="\t"} NR==FNR{a[$1]=$2} NR>FNR{print a[$1],$0}' \
   ${db}/EasyMicrobiome/dbcan2/CAZy_description.txt result/eggnog/eggnog.CAZy.raw.txt | \
  sed 's/^\t/Unannotated\t/' > result/eggnog/eggnog.CAZy.TPM.spf
head -n 3 result/eggnog/eggnog.CAZy.TPM.spf

# COG
awk 'BEGIN{FS=OFS="\t"} NR==FNR{a[$1]=$2"\t"$3} NR>FNR{print a[$1],$0}' \
  ${db}/EasyMicrobiome/eggnog/COG.anno result/eggnog/eggnog.COG_category.raw.txt > \
  result/eggnog/eggnog.COG_category.TPM.spf
head -n 3 result/eggnog/eggnog.COG_category.TPM.spf

CAZy碳水化合物酶

# 比对CAZy数据库, 用时2~18m
mkdir -p temp/dbcan3 result/dbcan3
# --sensitive慢10倍，dbcan3e值为1e-102，此处以1e-3演示
time diamond blastp \
  --db ${db}/dbcan3/CAZyDB \


**网上学习资料一大堆，但如果学到的知识不成体系，遇到问题时只是浅尝辄止，不再深入研究，那么很难做到真正的技术提升。**

**需要这份系统化的资料的朋友，可以添加V获取：vip1024b （备注运维）**
![img](https://img-blog.csdnimg.cn/img_convert/caf2d6a6d1e87dd3d4a8ac321ac29dfa.jpeg)

**一个人可以走的很快，但一群人才能走的更远！不论你是正从事IT行业的老鸟或是对IT行业感兴趣的新人，都欢迎加入我们的的圈子（技术交流、学习资源、职场吐槽、大厂内推、面试辅导），让我们一起学习成长！**
1.7 / Expected eggNOG DB version: 5.0.2 
# Diamond version found: diamond version 2.0.15

# 运行emapper，18m，默认diamond 1e-3
mkdir -p temp/eggnog
time emapper.py --data_dir ${db}/eggnog \
  -i result/NR/protein.fa --cpu 3 -m diamond --override \
  -o temp/eggnog/output

# 格式化结果并显示表头
grep -v '^##' temp/eggnog/output.emapper.annotations | sed '1 s/^#//' \
  > temp/eggnog/output
csvtk -t headers -v temp/eggnog/output

生成COG/KO/CAZy丰度汇总表

mkdir -p result/eggnog
# 显示帮助
summarizeAbundance.py -h
# 汇总，7列COG_category按字母分隔，12列KEGG_ko和19列CAZy按逗号分隔，原始值累加
summarizeAbundance.py \
  -i result/salmon/gene.TPM \
  -m temp/eggnog/output --dropkeycolumn \
  -c '7,12,19' -s '*+,+,' -n raw \
  -o result/eggnog/eggnog
sed -i 's#^ko:##' result/eggnog/eggnog.KEGG_ko.raw.txt
sed -i '/^-/d' result/eggnog/eggnog*
head -n3 result/eggnog/eggnog*
# eggnog.CAZy.raw.txt  eggnog.COG_category.raw.txt  eggnog.KEGG_ko.raw.txt

# 添加注释生成STAMP的spf格式
awk 'BEGIN{FS=OFS="\t"} NR==FNR{a[$1]=$2} NR>FNR{print a[$1],$0}' \
  ${db}/EasyMicrobiome/kegg/KO_description.txt \
  result/eggnog/eggnog.KEGG_ko.raw.txt | \
  sed 's/^\t/Unannotated\t/' \
  > result/eggnog/eggnog.KEGG_ko.TPM.spf
head -n 5 result/eggnog/eggnog.KEGG_ko.TPM.spf
# KO to level 1/2/3
summarizeAbundance.py \
  -i result/eggnog/eggnog.KEGG_ko.raw.txt \
  -m ${db}/EasyMicrobiome/kegg/KO1-4.txt \
  -c 2,3,4 -s ',+,+,' -n raw --dropkeycolumn \
  -o result/eggnog/KEGG
head -n3 result/eggnog/KEGG*

# CAZy
awk 'BEGIN{FS=OFS="\t"} NR==FNR{a[$1]=$2} NR>FNR{print a[$1],$0}' \
   ${db}/EasyMicrobiome/dbcan2/CAZy_description.txt result/eggnog/eggnog.CAZy.raw.txt | \
  sed 's/^\t/Unannotated\t/' > result/eggnog/eggnog.CAZy.TPM.spf
head -n 3 result/eggnog/eggnog.CAZy.TPM.spf

# COG
awk 'BEGIN{FS=OFS="\t"} NR==FNR{a[$1]=$2"\t"$3} NR>FNR{print a[$1],$0}' \
  ${db}/EasyMicrobiome/eggnog/COG.anno result/eggnog/eggnog.COG_category.raw.txt > \
  result/eggnog/eggnog.COG_category.TPM.spf
head -n 3 result/eggnog/eggnog.COG_category.TPM.spf

CAZy碳水化合物酶

# 比对CAZy数据库, 用时2~18m
mkdir -p temp/dbcan3 result/dbcan3
# --sensitive慢10倍，dbcan3e值为1e-102，此处以1e-3演示
time diamond blastp \
  --db ${db}/dbcan3/CAZyDB \


**网上学习资料一大堆，但如果学到的知识不成体系，遇到问题时只是浅尝辄止，不再深入研究，那么很难做到真正的技术提升。**

**需要这份系统化的资料的朋友，可以添加V获取：vip1024b （备注运维）**
[外链图片转存中...(img-UKE3zQVo-1713129723903)]

**一个人可以走的很快，但一群人才能走的更远！不论你是正从事IT行业的老鸟或是对IT行业感兴趣的新人，都欢迎加入我们的的圈子（技术交流、学习资源、职场吐槽、大厂内推、面试辅导），让我们一起学习成长！**