先自我介绍一下，小编浙江大学毕业，去过华为、字节跳动等大厂，目前阿里P7
深知大多数程序员，想要提升技能，往往是自己摸索成长，但自己不成体系的自学效果低效又漫长，而且极易碰到天花板技术停滞不前！
因此收集整理了一份《2024年最新Linux运维全套学习资料》，初衷也很简单，就是希望能够帮助到想自学提升又不知道该从何学起的朋友。
既有适合小白学习的零基础资料，也有适合3年以上经验的小伙伴深入学习提升的进阶课程，涵盖了95%以上运维知识点，真正体系化！
由于文件比较多，这里只是将部分目录截图出来，全套包含大厂面经、学习笔记、源码讲义、实战项目、大纲路线、讲解视频，并且后续会持续更新
如果你需要这些资料，可以添加V获取：vip1024b （备注运维）
正文
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
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
1. 首次使用请参照`0Install.sh`脚本，安装软件和数据库(大约1-3天，仅一次)
2. 易宏基因组(EasyMetagenome)流程`1Pipeline.sh`复制到项目文件夹，如本次为meta
3. 项目文件夹准备测序数据(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=
s
o
f
t
/
b
i
n
:
soft/bin:
soft/bin:soft/condabin:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin:
d
b
/
E
a
s
y
M
i
c
r
o
b
i
o
m
e
/
l
i
n
u
x
:
db/EasyMicrobiome/linux:
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
.
f
q
.
g
z
−
I
s
e
q
/
{i}_1.fq.gz -I seq/
i1.fq.gz −Iseq/{i}_2.fq.gz
-o temp/qc/
i
1
.
f
a
s
t
q
−
O
t
e
m
p
/
q
c
/
{i}_1.fastq -O temp/qc/
i1.fastq−Otemp/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 "KaTeX parse error: Undefined control sequence: \t at position 2: i\̲t̲" >> temp/fastp…{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/KaTeX parse error: Undefined control sequence: \/ at position 2: /\̲/̲1/' temp/qc/{}_…//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
?
.
f
a
s
t
q
>
t
e
m
p
/
c
o
n
c
a
t
/
{i}_?.fastq \ > temp/concat/
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目录
( c m d ) 与 ‘ c m d ‘ 通常是等价的； ‘ c m d ‘ 写法更简单，但要注意反引号是键盘左上角 E S C 下面的按键， (cmd) 与 `cmd` 通常是等价的；`cmd`写法更简单，但要注意反引号是键盘左上角ESC下面的按键， (cmd)与‘cmd‘通常是等价的；‘cmd‘写法更简单，但要注意反引号是键盘左上角ESC下面的按键，(cmd)更通用，适合嵌套使用
for i in
(
t
a
i
l
−
n
+
2
r
e
s
u
l
t
/
m
e
t
a
d
a
t
a
.
t
x
t
∣
c
u
t
−
f
1
)
;
d
o
m
v
t
e
m
p
/
h
u
m
a
n
n
3
/
(tail -n+2 result/metadata.txt | cut -f1); do mv temp/humann3/
(tail−n+2result/metadata.txt∣cut−f1);do mvtemp/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
.
f
a
s
t
q
t
e
m
p
/
h
u
m
a
n
n
3
/
{i}_1.fastq \ temp/humann3/
i1.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[KaTeX parse error: Expected 'EOF', got '}' at position 3: 1]}̲' result/metada…/\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 KaTeX parse error: Expected group after '_' at position 96: …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
g
e
n
e
f
a
m
i
l
i
e
s
.
t
s
v
−
g
u
n
i
r
e
f
9
0
k
o
−
o
t
e
m
p
/
h
u
m
a
n
n
3
/
{i}_genefamilies.tsv \ -g uniref90_ko \ -o temp/humann3/
igenefamilies.tsv −guniref90ko −otemp/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]}’ KaTeX parse error: Undefined control sequence: \n at position 39: …p/sampleid|tr '\̲n̲' '\t'|sed 's/\…/\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
d
b
/
k
r
a
k
e
n
2
/
p
l
u
s
p
f
16
g
/
−
−
p
a
i
r
e
d
t
e
m
p
/
q
c
/
{db}/kraken2/pluspf16g/ \ --paired temp/qc/
db/kraken2/pluspf16g/ −−pairedtemp/qc/{i}_?.fastq
--threads 2 --use-names --report-zero-counts
--report temp/kraken2/
i
.
r
e
p
o
r
t
−
−
o
u
t
p
u
t
t
e
m
p
/
k
r
a
k
e
n
2
/
{i}.report \ --output temp/kraken2/
i.report −−outputtemp/kraken2/{i}.output
多样本并行生成report，1样本8线程逐个运行，内存大但速度快，不建议用多任务并行
for i in tail -n+2 result/metadata.txt | cut -f1;do
kraken2 --db
d
b
/
k
r
a
k
e
n
2
/
p
l
u
s
p
f
16
g
−
−
p
a
i
r
e
d
t
e
m
p
/
q
c
/
{db}/kraken2/pluspf16g \ --paired temp/qc/
db/kraken2/pluspf16g −−pairedtemp/qc/{i}_?.fastq
--threads 2 --use-names --report-zero-counts
--report temp/kraken2/
i
.
r
e
p
o
r
t
−
−
o
u
t
p
u
t
t
e
m
p
/
k
r
a
k
e
n
2
/
{i}.report \ --output temp/kraken2/
i.report −−outputtemp/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
.
r
e
p
o
r
t
−
−
d
i
s
p
l
a
y
−
h
e
a
d
e
r
−
o
t
e
m
p
/
k
r
a
k
e
n
2
/
{i}.report \ --display-header -o temp/kraken2/
i.report −−display−header−otemp/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
h
e
a
d
e
r
t
a
i
l
−
n
+
2
t
e
m
p
/
k
r
a
k
e
n
2
/
header tail -n+2 temp/kraken2/
headertail−n+2temp/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
d
b
/
k
r
a
k
e
n
2
/
p
l
u
s
p
f
16
g
/
−
i
t
e
m
p
/
k
r
a
k
e
n
2
/
{db}/kraken2/pluspf16g/ \ -i temp/kraken2/
db/kraken2/pluspf16g/ −itemp/kraken2/{i}.report
-r 100 -l
t
a
x
−
t
0
−
o
t
e
m
p
/
b
r
a
c
k
e
n
/
{tax} -t 0 \ -o temp/bracken/
tax−t0 −otemp/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
d
b
/
E
a
s
y
M
i
c
r
o
b
i
o
m
e
/
s
c
r
i
p
t
/
f
i
l
t
e
r
f
e
a
t
u
r
e
t
a
b
l
e
.
R
−
i
r
e
s
u
l
t
/
k
r
a
k
e
n
2
/
b
r
a
c
k
e
n
.
{db}/EasyMicrobiome/script/filter_feature_table.R \ -i result/kraken2/bracken.
db/EasyMicrobiome/script/filterfeaturetable.R −iresult/kraken2/bracken.{tax}.txt
-p 0.01
-o result/kraken2/bracken.
t
a
x
.
0.01
c
s
v
t
k
−
t
s
t
a
t
r
e
s
u
l
t
/
k
r
a
k
e
n
2
/
b
r
a
c
k
e
n
.
{tax}.0.01 csvtk -t stat result/kraken2/bracken.
tax.0.01csvtk−tstatresult/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 "KaTeX parse error: Undefined control sequence: \t at position 2: i\̲t̲" >> result/kra…{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
.
r
e
p
o
r
t
−
o
t
e
m
p
/
b
r
a
c
k
e
n
/
{i}.report -o temp/bracken/
i.report−otemp/bracken/{i}.krona --no-intermediate-ranks
ktImportText temp/bracken/
i
.
k
r
o
n
a
−
o
r
e
s
u
l
t
/
k
r
a
k
e
n
2
/
k
r
o
n
a
.
{i}.krona -o result/kraken2/krona.
i.krona−oresult/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],KaTeX parse error: Expected 'EOF', got '}' at position 2: 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
--query result/NR/protein.fa
--threads 2 -e 1e-3 --outfmt 6 --max-target-seqs 1 --quiet
--out temp/dbcan3/gene_diamond.f6
wc -l temp/dbcan3/gene_diamond.f6
提取基因与dbcan分类对应表
format_dbcan2list.pl
-i temp/dbcan3/gene_diamond.f6
-o temp/dbcan3/gene.list
按对应表累计丰度，依赖
summarizeAbundance.py
-i result/salmon/gene.TPM
-m temp/dbcan3/gene.list
-c 2 -s ‘,’ -n raw --dropkeycolumn
-o result/dbcan3/TPM
添加注释生成STAMP的spf格式，结合metadata.txt进行差异比较
awk 'BEGIN{FS=OFS=“\t”} NR==FNR{a[$1]=$2} NR>FNR{print a[$1],KaTeX parse error: Expected 'EOF', got '}' at position 2: 0}̲' \ {db}/EasyMicrobiome/dbcan2/CAZy_description.txt result/dbcan3/TPM.CAZy.raw.txt |
sed ‘s/^\t/Unannotated\t/’
> result/dbcan3/TPM.CAZy.raw.spf
head result/dbcan3/TPM.CAZy.raw.spf
检查未注释数量，有则需要检查原因
grep ‘Unannotated’ result/dbcan3/TPM.CAZy.raw.spf|wc -l
CARD耐药基因
CARD在线分析平台：https://card.mcmaster.ca/
本地软件使用教程: https://github.com/arpcard/rgi
参考文献：http://doi.org/10.1093/nar/gkz935
结果说明：protein.json，在线可视化；protein.txt，注释基因列表
mkdir -p result/card
启动rgi环境和记录版本
conda activate rgi6
rgi main -v # 6.0.3
简化蛋白ID
cut -f 1 -d ’ ’ result/NR/protein.fa > temp/protein.fa
这个错误忽略即可，不是报错，没有任何影响 grep: 写错误: 断开的管道
grep ‘>’ result/NR/protein.fa | head -n 3
grep ‘>’ temp/protein.fa | head -n 3
蛋白层面注释ARG
rgi main -i temp/protein.fa -t protein
-n 9 -a DIAMOND --include_loose --clean
-o result/card/protein
head -n3 result/card/protein.txt
基因层面注释ARG
cut -f 1 -d ’ ’ result/NR/nucleotide.fa > temp/nucleotide.fa
grep ‘>’ temp/nucleotide.fa | head -n3
rgi main -i temp/nucleotide.fa -t contig
-n 9 -a DIAMOND --include_loose --clean
-o result/card/nucleotide
head -n3 result/card/nucleotide.txt
重叠群层面注释ARG
cut -f 1 -d ’ ’ result/megahit/final.contigs.fa > temp/contigs.fa
grep ‘>’ temp/contigs.fa | head -n3
rgi main -i temp/contigs.fa -t contig
-n 9 -a DIAMOND --include_loose --clean
-o result/card/contigs
head result/card/contigs.txt
#### 3.4 基因物种注释
Generate report in default taxid output
conda activate kraken2
kraken2 --db ${db}/kraken2/pluspf16g
result/NR/nucleotide.fa
--threads 3
--report temp/NRgene.report
--output temp/NRgene.output
Genes & taxid list
grep ‘^C’ temp/NRgene.output | cut -f 2,3 | sed ‘1 i Name\ttaxid’
> temp/NRgene.taxid
Add taxonomy
awk ‘BEGIN{FS=OFS=“\t”} NR==FNR{a[$1]=$0} NR>FNR{print $1,a[$2]}’
$db/EasyMicrobiome/kraken2/taxonomy.txt
temp/NRgene.taxid > result/NR/nucleotide.tax
conda activate eggnog
summarizeAbundance.py
-i result/salmon/gene.TPM
-m result/NR/nucleotide.tax --dropkeycolumn
-c ‘2,3,4,5,6,7,8,9’ -s ‘,+,+,+,+,+,+,+,’ -n raw
-o result/NR/tax
wc -l result/NR/tax*|sort -n
### 四、分箱挖掘单菌基因组Binning
#### 4.1 MetaWRAP混合样本分箱 Samples binning
主页：https://github.com/bxlab/metaWRAP
教程: https://github.com/bxlab/metaWRAP/blob/master/Usage\_tutorial.md
挖掘单菌基因组，需要研究对象复杂度越低、测序深度越大，结果质量越好。要求单样本6GB+，复杂样本如土壤推荐数据量30GB+，至少3个样本
演示数据12个样仅140MB，无法获得单菌基因组，这里使用官方测序数据演示讲解
软件和数据库布置需1-3天，演示数据分析过程超10h，30G样也需1-30天，由服务器性能决定。
设置并进入工作目录
wd=~/meta/binning
mkdir -p ${wd} && cd ${wd}
初始化项目
mkdir -p temp/hr seq result
启动metawrap环境
conda activate metawrap
数据和环境变量 Data and enviroment
这里基于质控clean数据和拼接好的重叠群contigs，基于上游结果继续分析。由于上游测试数据过小，分箱无结果。 本次采用软件推荐的7G数据，我们进入一个新文件夹开展分析。
输入输出文件介绍：
输入：质控后序列，文件名格式为*_1.fastq和*_2.fastq，temp/qc 目录下，如C1_1.fastq、C1_2.fastq
组装的重叠群文件：result/megahit/final.contigs.fa
输出：
Binning结果：temp/binning
提纯后的Bin统计结果：temp/bin_refinement/metawrap_50_10_bins.stats
Bin定量结果文件和图：binning/temp/bin_quant/bin_abundance_table.tab 和 bin_abundance_heatmap.png
Bin物种注释：binning/temp/bin_classify/bin_taxonomy.tab
Prokka基因预测：binning/temp/bin_annotate/prokka_out/bin.*.ffn 核酸序列
Bin可视化图表：binning/temp/bloblogy/final.contigs.binned.blobplot (数据表) 和 blobplot_figures (可视化图)
准备输入文件：原始数据+组装结果
质控后数据位于temp/qc中，此处需下载并解压
方法1. 直接拷贝
/bin/cp -rf /db/metawrap/*.fastq ~/meta/binning/temp/hr/
方法2. 在线下载
cd temp/hr
for i in seq 7 9;do
wget -c ftp.sra.ebi.ac.uk/vol1/fastq/ERR011/ERR01134
i
/
E
R
R
01134
{i}/ERR01134
i/ERR01134{i}_1.fastq.gz
wget -c ftp.sra.ebi.ac.uk/vol1/fastq/ERR011/ERR01134
i
/
E
R
R
01134
{i}/ERR01134
i/ERR01134{i}_2.fastq.gz
done
gunzip -k *.gz
批量修改扩展名fq为fastq
rename .fq .fastq *.fq
megahit拼接结果
cd ${wd}
mkdir -p temp/megahit
cd temp/megahit
可从EasyMetagenome目录复制，或链接下载
wget -c http://www.imeta.science/db/metawrap/final.contigs.fa.gz
gunzip -k *.gz
cd ${wd}
分箱Binning
加载运行环境
cd ${wd}
conda activate metawrap
metawrap -v # 1.3.2
输入文件为contig和clean reads
调用maxbin2, metabat2，8p线程2h，24p耗时1h；-concoct 3h
nohup metawrap binning -o temp/binning
-t 3 -a temp/megahit/final.contigs.fa
--metabat2 --maxbin2
temp/hr/*.fastq
–concoct > /dev/null 2>&1 增加3~10倍计算量，添加/dev/null清除海量Warning信息
分箱提纯Bin refinement
8线程2h， 24p 1.3h；2方法16p 20m
metawrap bin_refinement
-o temp/bin_refinement
-A temp/binning/metabat2_bins/
-B temp/binning/maxbin2_bins/
-c 50 -x 10 -t 8
-C temp/binning/concoct_bins/ \
统计高质量Bin的数量，2方法6个，3方法9个
tail -n+2 temp/bin_refinement/metawrap_50_10_bins.stats|wc -l
分析比较图见 temp/bin_refinement/figures/
所有分箱至同一目录All bins in one directory
mkdir -p temp/drep_in
混合组装分箱链接和重命名
ln -s pwd/temp/bin_refinement/metawrap_50_10_bins/bin.* temp/drep_in/
ls -l temp/drep_in/
改名CentOS
rename ‘bin.’ ‘Mx_All_’ temp/drep_in/bin.*
改名Ubuntu
rename s/bin./Mx_All_/ temp/drep_in/bin.*
ls temp/drep_in/Mx*
(可选Opt)单样本分箱Single sample binning
多样本受硬件、计算时间限制无法完成时，需要单样本组装、分箱。多样本信息丰度，分箱结果更多，更容易降低污染。详见：- [Nature Methods | 单样本与多样本宏基因组分箱的比较揭示了广泛存在的隐藏性污染]( )
\*\*设置全局线程、并行任务数和筛选分箱的条件\*\*
p:threads线程数,job任务数,complete完整度x:contaminate污染率
conda activate metawrap
p=16
j=3
c=50
x=10
(可选)并行需要样本列表，请提前编写metadata.txt保存于result中
快速读取文件生成样本ID列表再继续编写
ls temp/hr/ | grep 1 | cut -f 1 -d '’ | sort -u | sed ‘1 i SampleID’ > result/metadata.txt
预览
cat result/metadata.txt
\*\*组装Assemble\*\*
单样本并行组装；支持中断继续运行
tail -n+2 result/metadata.txt|cut -f1|rush -j ${j}
“metawrap assembly -m 200 -t ${p} --megahit
-1 temp/hr/{}_1.fastq -2 temp/hr/{}_2.fastq
-o temp/megahit/{}”
\*\*分箱binning\*\*
单样本并行分箱，192p, 15m (concoct使用超多线程)；16p 2d/sample, >/dev/null 16p 12h/sample
time tail -n+2 result/metadata.txt|cut -f1|rush -j ${j}
“metawrap binning
-o temp/binning/{} -t ${p}
-a temp/megahit/{}/final_assembly.fasta
--metabat2 --maxbin2 --concoct
temp/hr/{}_*.fastq > /dev/null 2>&1”
\*\*分箱提纯bin refinement\*\*
time tail -n+2 result/metadata.txt|cut -f1|rush -j ${j}
"metawrap bin_refinement
-o temp/bin_refinement/{} -t ${p}
-A temp/binning/{}/metabat2_bins/
-B temp/binning/{}/maxbin2_bins/
-C temp/binning/{}/concoct_bins/
-c ${c} -x ${x} "
分别为1,2,2个
tail -n+2 result/metadata.txt|cut -f1|rush -j 1
"tail -n+2 temp/bin_refinement/{}/metawrap_50_10_bins.stats|wc -l "
单样品分箱链接和重命名
for i in tail -n+2 result/metadata.txt|cut -f1;do
ln -s pwd/temp/bin_refinement/KaTeX parse error: Expected 'EOF', got '#' at position 48: …mp/drep_in/ #̲ CentOS rena…{i}" temp/drep_in/bin.*
# Ubuntu
rename "s/bin./Sg${i}_/" temp/drep_in/bin.*
done
删除空白中无效链接
/bin/rm -f temp/drep_in/**
统计混合和单样本来源数据，10个混，5个单；不同系统结果略有差异
ls temp/drep_in/|cut -f 1 -d ‘_’|uniq -c
统计混合批次/单样本来源
ls temp/drep_in/|cut -f 2 -d ‘_’|cut -f 1 -d ‘.’ |uniq -c
(可选Opt)分组分箱 Subgroup binning
样本>30或数据量>300G在1TB内存胖结点上完成混合组装和分箱可能内存不足、且时间>1周甚至1月，需要对研究相近条件、地点进行分小组，且每组编写一个metadata??.txt。
conda activate metawrap
小组ID: A1/A2/A3
g=A1
组装Assemble：<30个或<300G样本，~12h
metawrap assembly -m 600 -t 32 --megahit
-1 tail -n+2 result/metadata${g}.txt|cut -f1|sed 's/^/temp\/hr\//;s/$/_1.fastq/'|tr '\n' ','|sed 's/,$//'
-2 tail -n+2 result/metadata${g}.txt|cut -f1|sed 's/^/temp\/hr\//;s/$/_2.fastq/'|tr '\n' ','|sed 's/,$//'
-o temp/megahit_${g}
分箱Binning，~18h
链接文件到临时位置
mkdir -p temp/
g
/
f
o
r
i
i
n
‘
t
a
i
l
−
n
+
2
r
e
s
u
l
t
/
m
e
t
a
d
a
t
a
{g}/ for i in `tail -n+2 result/metadata
g/foriin‘tail−n+2result/metadata{g}.txt|cut -f1;do ln -s pwd`/temp/hr/
i
∗
.
f
a
s
t
q
t
e
m
p
/
{i}*.fastq temp/
i∗.fastqtemp/{g}/
done
按组分箱
metawrap binning -o temp/binning_KaTeX parse error: Expected group after '_' at position 30: …-a temp/megahit_̲{g}/final_assembly.fasta
--metabat2 --maxbin2
temp/${g}/*.fastq
\*\*分箱提纯Bin refinement\*\*
metawrap bin_refinement
-o temp/bin_refinement_KaTeX parse error: Expected group after '_' at position 24: …-A temp/binning_̲{g}/metabat2_bins/
-B temp/binning_${g}/maxbin2_bins/
-c 50 -x 10 -t 32
统计高质量Bin的数量
wc -l temp/bin_refinement_${g}/metawrap_50_10_bins.stats
改名汇总 Rename & merge
mkdir -p temp/drep_in
混合组装分箱链接和重命名
ln -s pwd/temp/bin_refinement_${g}/metawrap_50_10_bins/bin.* temp/drep_in/
改名
rename “s/bin./Gp_${g}_/” temp/drep_in/bin.* # Ubuntu
rename ‘bin.’ “Gp_${g}_” temp/drep_in/bin.* # CentOS
统计
mkdir -p result/bin
echo -n KaTeX parse error: Expected group after '_' at position 47: …temp/drep_in/Gp_̲{g}_*|wc>> result/bin/groupNo.txt
cat result/bin/groupNo.txt
#### 4.2 dRep去冗余种/株基因组集
进入虚拟环境drep和工作目录
conda activate drep
cd ${wd}
按种水平去冗余：6~40min，15个为10个，8个来自混拼，2个来自单拼
mkdir -p temp/drep95
/bin/rm -rf temp/drep95/data/checkM
time dRep dereplicate temp/drep95/
-g temp/drep_in/*.fa
-sa 0.95 -nc 0.30 -comp 50 -con 10 -p 5
报错日志在temp/drep95/log/cmd_logs中查看，加-d显示更多
ls temp/drep95/dereplicated_genomes/|cut -f 1 -d ‘_’|sort|uniq -c
主要结果temp/drep95中：
\* 非冗余基因组集：temp/drep95/dereplicated\_genomes/\*.fa
\* 聚类信息表：temp/drep95/data\_tables/Cdb.csv
\* 聚类和质量图：temp/drep95/figures/\*clustering\*
(可选)按株水平99%去冗余，20-30min，本处也为10个
mkdir -p temp/drep99
time dRep dereplicate temp/drep99/
-g temp/drep_in/*.fa
-sa 0.99 -nc 0.30 -comp 50 -con 10 -p 5
ls -l temp/drep99/dereplicated_genomes/ | grep ‘.fa’ | wc -l
#### 4.3 CoverM基因组定量
启动环境
conda activate coverm
mkdir -p temp/coverm
(可选)单样本测试, 3min
i=ERR011347
time coverm genome --coupled temp/hr/
i
1
.
f
a
s
t
q
t
e
m
p
/
h
r
/
{i}_1.fastq temp/hr/
i1.fastqtemp/hr/{i}_2.fastq
--genome-fasta-directory temp/drep95/dereplicated_genomes/ -x fa
-o temp/coverm/
i
.
t
x
t
c
a
t
t
e
m
p
/
c
o
v
e
r
m
/
{i}.txt cat temp/coverm/
i.txtcattemp/coverm/{i}.txt
并行计算, 4min
tail -n+2 result/metadata.txt|cut -f1|rush -j 2
“coverm genome --coupled temp/hr/{}_1.fastq temp/hr/{}_2.fastq
--genome-fasta-directory temp/drep95/dereplicated_genomes/ -x fa
-o temp/coverm/{}.txt”
结果合并
mkdir -p result/coverm
conda activate humann3
sed -i ‘s/_1.fastq Relative Abundance (%)//’ temp/coverm/*.txt
humann_join_tables --input temp/coverm
--file_name txt
--output result/coverm/abundance.tsv
按组求均值，需要metadata中有3列且每个组有多个样本
Rscript ${db}/EasyMicrobiome/script/otu_mean.R --input result/coverm/abundance.tsv
--metadata result/metadata.txt
--group Group --thre 0
--scale TRUE --zoom 100 --all TRUE --type mean
--output result/coverm/group_mean.txt
https://www.bic.ac.cn/ImageGP/ 直接选择热图可视化
#### 4.4 GTDB-tk物种注释和进化树
启动软件所在虚拟环境
conda activate gtdbtk2.3
export GTDBTK_DATA_PATH=“${db}/gtdb”
gtdbtk -v # 2.3.2
代表性细菌基因组物种注释
mkdir -p temp/gtdb_classify
10个基因组，24p，100min 152G内存; 6p, 22基因组，1h
gtdbtk classify_wf
--genome_dir temp/drep95/dereplicated_genomes
--out_dir temp/gtdb_classify
--extension fa --skip_ani_screen
--prefix tax
--cpus 6
less -S按行查看，按q退出
less -S temp/gtdb_classify/tax.bac120.summary.tsv
less -S temp/gtdb_classify/tax.ar53.summary.tsv
代表种注释：以上面鉴定的10个种为例，注意扩展名要与输入文件一致，可使用压缩格式gz。主要结果文件描述：此9个细菌基因组在tax.bac120.summary.tsv。古菌在tax.ar53开头的文件中。
(可选)所有MAG物种注释
mkdir -p temp/gtdb_all
10000个基因组，32p，100min
time gtdbtk classify_wf
--genome_dir temp/drep_in/
--out_dir temp/gtdb_all
--extension fa --skip_ani_screen
--prefix tax
--cpus 6
less -S temp/gtdb_all/tax.bac120.summary.tsv
less -S temp/gtdb_all/tax.ar53.summary.tsv
多序列对齐结果建树
# 以9个细菌基因组的120个单拷贝基因建树，1s
mkdir -p temp/gtdb_infer
gtdbtk infer --msa_file temp/gtdb_classify/align/tax.bac120.user_msa.fasta.gz
--out_dir temp/gtdb_infer --prefix tax --cpus 3
树文件`tax.unrooted.tree`可使用iTOL在线美化，也可使用GraphLan本地美化。
制作树注释文件：以gtdb-tk物种注释(tax.bac120.summary.tsv)和drep基因组评估(Widb.csv)信息为注释信息
mkdir -p result/itol
制作分类学表
tail -n+2 temp/gtdb_classify/tax.bac120.summary.tsv|cut -f 1-2|sed ‘s/;/\t/g’|sed ‘1 s/^/ID\tDomain\tPhylum\tClass\tOrder\tFamily\tGenus\tSpecies\n/’
> result/itol/tax.txt
head result/itol/tax.txt
基因组评估信息
sed ‘s/,/\t/g;s/.fa//’ temp/drep95/data_tables/Widb.csv|cut -f 1-7,11|sed ‘1 s/genome/ID/’
> result/itol/genome.txt
head result/itol/genome.txt
整合注释文件
awk ‘BEGIN{OFS=FS=“\t”} NR==FNR{a[$1]=$0} NR>FNR{print $0,a[$1]}’ result/itol/genome.txt result/itol/tax.txt|cut -f 1-8,10- > result/itol/annotation.txt
head result/itol/annotation.txt
添加各样本相对丰度(各组替换均值)
awk ‘BEGIN{OFS=FS=“\t”} NR==FNR{a[$1]=$0} NR>FNR{print $0,a[$1]}’ <(sed ‘1 s/Genome/ID/’ result/coverm/abundance.tsv) result/itol/annotation.txt|cut -f 1-15,17- > result/itol/annotation2.txt
head result/itol/annotation2.txt
CheckM2重新评估
conda activate checkm2
mkdir -p temp/checkm2 result/checkm2
10 genomes, 2m
time checkm2 predict --input temp/drep95/dereplicated_genomes/* --output-directory temp/checkm2 --threads 8
ln temp/checkm2/quality_report.tsv result/checkm2/
查看结果
less result/checkm2//quality_report.tsv
### 五、(可选)单菌基因组
#### 5.1 Fastp质量控制
每个样本~30s，173个j2共
mkdir -p temp/qc/
time tail -n+2 result/metadata.txt | cut -f1 | rush -j 2
“time 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”
#### 5.2 metaspades组装
conda activate megahit
spades.py -v # v3.15.4，>3.14.0才支持–isolate模式
mkdir -p temp/spades result/spades
127 genoms, 1m17s
time tail -n+2 result/metadata.txt|cut -f1|rush -j 3
“spades.py --pe1-1 temp/qc/{1}_1.fastq
--pe1-2 temp/qc/{1}_2.fastq
-t 16 --isolate --cov-cutoff auto
-o temp/spades/{1}”
筛选>1k的序列并汇总、统计
time tail -n+2 result/metadata.txt|cut -f1|rush -j 3
“seqkit seq -m 1000 temp/spades/{1}/contigs.fasta
> temp/spades/{1}.fa”
seqkit stat temp/spades/*.fa | sed ‘s/temp/spades///;s/.fa//’ > result/spades/stat1k.txt
#### 5.3 checkm质量评估
checkm评估质量
conda activate drep
checkm # CheckM v1.1.2
mkdir -p temp/checkm result/checkm
127 genoms, 1m17s
time checkm lineage_wf -t 8 -x fa temp/spades/ temp/checkm
format checkm jason to tab
checkmJason2tsv.R -i temp/checkm/storage/bin_stats_ext.tsv
-o temp/checkm/bin_stats.txt
csvtk -t pretty temp/checkm/bin_stats.txt | less
(可选)checkm2评估(测试中...)
conda activate checkm2
mkdir -p temp/checkm2
time checkm2 predict --threads 8 --input temp/spades/ --output-directory temp/checkm2
筛选污染和高质量基因组 >5% contamination and high quailty
awk ‘$5<90 || $10>5’ temp/checkm/bin_stats.txt | csvtk -t cut -f 1,5,10,4,9,2 > temp/checkm/contamination5.txt
tail -n+2 temp/checkm/contamination5.txt|wc -l
筛选高质量用于下游分析 <5% high-quality for down-stream analysis
awk ‘$5>=90 && $10<=5’ temp/checkm/bin_stats.txt | csvtk -t cut -f 1,5,10,4,9,2 | sed ‘1 i ID\tCompleteness\tContamination\tGC\tN50\tsize’ > result/checkm/Comp90Cont5.txt
tail -n+2 result/checkm/Comp90Cont5.txt|wc -l
链接高质量基因组至新目录，单菌完整度通常>99%
mkdir -p temp/drep_in/
for n in tail -n+2 result/checkm/Comp90Cont5.txt|cut -f 1;do
ln temp/spades/${n}.fa temp/drep_in/
done
#### 5.4 混菌metawarp分箱
分箱和提纯binning & refinement
conda activate metawrap
mkdir -p temp/binning temp/bin
time tail -n+2 temp/checkm/contamination5.txt|cut -f1|rush -j 3
“metawrap binning
-o temp/binning/{} -t 8
-a temp/spades/{}/contigs.fasta
--metabat2 --maxbin2
temp/qc/{}_*.fastq”
time tail -n+2 temp/checkm/contamination5.txt|cut -f1|rush -j 15
“metawrap bin_refinement
-o temp/bin/{} -t 8
-A temp/binning/{}/metabat2_bins/
-B temp/binning/{}/maxbin2_bins/
-c 50 -x 10”
分箱结果汇总
echo -n -e “” > temp/bin/metawrap.stat
for m in tail -n+2 temp/checkm/contamination5.txt|cut -f1;do
echo
m
>
>
t
e
m
p
/
b
i
n
/
m
e
t
a
w
r
a
p
.
s
t
a
t
c
u
t
−
f
1
−
4
,
6
−
7
t
e
m
p
/
b
i
n
/
{m} >> temp/bin/metawrap.stat cut -f1-4,6-7 temp/bin/
m>>temp/bin/metawrap.stat cut−f1−4,6−7temp/bin/{m}/metawrap_50_10_bins.stats >> temp/bin/metawrap.stat
done
分箱后的按b1,b2,b3重命名共培养，单菌也可能减少污染
for m in tail -n+2 temp/checkm/contamination5.txt|cut -f1;do
c=1
for n in tail -n+2 temp/bin/$m/metawrap_50_10_bins.stats|cut -f 1;do
cp temp/bin/
m
/
m
e
t
a
w
r
a
p
5
0
1
0
b
i
n
s
/
m/metawrap_50_10_bins/
m/metawrap5010bins/{n}.fa temp/drep_in/
m
b
{m}b
mb{c}.fa
((c++))
done
done
分箱前后统计比较
如107个测序分箱为352个基因组，共418个基因组
tail -n+2 temp/checkm/contamination5.txt|wc -l
ls temp/drep_in/b?.fa | wc -l
ls temp/drep_in/.fa | wc -l
重建新ID列表，A代表所有，B代表Bin分箱过的单菌
ls temp/drep_in/*.fa|cut -f 3 -d ‘/’|sed ‘s/.fa//’|sed ‘1 i ID’|less -S>result/metadataA.txt
ls temp/drep_in/*b?.fa|cut -f 3 -d ‘/’|sed ‘s/.fa//’|sed ‘1 i ID’|less -S>result/metadataB.txt
可视化混菌中覆盖度分布，以第一污染菌为例
i=tail -n+2 temp/checkm/contamination5.txt|head -n1|cut -f1
grep ‘>’ temp/spades/KaTeX parse error: Expected group after '_' at position 25: …t -f 2,4,6 -d '_̲'|sed 's/^/C/;s…{i}.len
grep ‘>’ temp/drep_in/KaTeX parse error: Expected group after '_' at position 61: …-f 1,2,4,6 -d '_̲'|sed 's/_/\t/g…{i}.cov
sp_lines.sh -f temp/drep_in/${i}.cov -m TRUE -A TRUE -a Contig -H Genome
#### 5.5 drep基因组去冗余
mkdir -p temp/drep95/
conda activate drep
相似度sa 0.99995 去重复, 0.99 株水平, 0.95 种水平
dRep dereplicate
-g temp/drep_in/*.fa
-sa 0.95 -nc 0.3 -p 16 -comp 50 -con 10
temp/drep95
统计总和使用基因组数据，确保没有异常丢弃stat total and used genomes, keep no abnormal discard
ls temp/drep_in/.fa|wc -l
grep ‘passed checkM’ temp/drep95/log/logger.log|sed 's/[ ][ ]/ /g’|cut -f 4 -d ’ ’
去冗余后数量，418变为49种
ls temp/drep95/dereplicated_genomes/|wc -l
唯一和重复的基因组unique and duplicate genome
csvtk cut -f 11 temp/drep95/data_tables/Widb.csv | sort | uniq -c
整理种列表
echo “SampleID” > result/metadataS.txt
ls temp/drep95/dereplicated_genomes/|sed ‘s/.fa//’ >> result/metadataS.txt
#### 5.6 gtdb物种注释
conda activate gtdbtk2.3
所有基因组注释，400g, 1h, 1T
mkdir -p temp/gtdb_all result/gtdb_all
memusg -t gtdbtk classify_wf
--genome_dir temp/drep_in/
--out_dir temp/gtdb_all/
--extension fa --skip_ani_screen
--prefix tax
--cpus 16
基因组信息genomeInfo.csv
sed ‘s/,/\t/g;s/.fa//’ temp/drep95/data_tables/genomeInfo.csv |sed ‘1 s/genome/ID/’ > result/gtdb_all/genome.txt
95%聚类种基因组注释，40g, 1h, 500G
mkdir -p temp/gtdb_95 result/gtdb_95
Taxonomy classify
memusg -t gtdbtk classify_wf
--genome_dir temp/drep95/dereplicated_genomes/
--out_dir temp/gtdb_95
--extension fa --skip_ani_screen
--prefix tax
--cpus 8
Phylogenetic tree infer
memusg -t gtdbtk infer
--msa_file temp/gtdb_95/align/tax.bac120.user_msa.fasta.gz
--out_dir temp/gtdb_95
--cpus 8 --prefix g >> temp/gtdb_95/infer.log 2>&1
ln pwd/temp/gtdb_95/infer/intermediate_results/g.unrooted.tree result/gtdb_95/
细菌format to standard 7 levels taxonomy
tail -n+2 temp/gtdb_95/classify/tax.bac120.summary.tsv|cut -f 1-2|sed ‘s/;/\t/g’|sed ‘1 s/^/ID\tKingdom\tPhylum\tClass\tOrder\tFamily\tGenus\tSpecies\n/’ > result/gtdb_95/tax.bac.txt
古菌(可选)
tail -n+2 temp/gtdb_95/classify/tax.ar122.summary.tsv|cut -f 1-2|sed ‘s/;/\t/g’|sed ‘1 s/^/ID\tKingdom\tPhylum\tClass\tOrder\tFamily\tGenus\tSpecies\n/’ > result/gtdb_95/tax.ar.txt
cat result/gtdb_95/tax.bac.txt <(tail -n+2 result/gtdb_95/tax.ar.txt) > result/gtdb_95/tax.txt
Widb.csv 非冗余基因组信息
sed ‘s/,/\t/g;s/.fa//’ temp/drep95/data_tables/Widb.csv|cut -f 1-7,11|sed ‘1 s/genome/ID/’ > result/gtdb_95/genome.txt
整合物种注释和基因组信息 Integrated taxonomy and genomic info
awk ‘BEGIN{OFS=FS=“\t”} NR==FNR{a[$1]=$0} NR>FNR{print $0,a[$1]}’ result/gtdb_95/genome.txt result/gtdb_95/tax.txt|cut -f 1-8,10- > result/gtdb_95/annotation.txt
csvtk -t headers -v result/gtdb_95/annotation.txt
制作itol files
cd result/gtdb_95
table2itol.R -D plan1 -a -c double -i ID -l Genus -t %s -w 0.5 annotation.txt
table2itol.R -D plan2 -a -d -c none -b Phylum -i ID -l Genus -t %s -w 0.5 annotation.txt
table2itol.R -D plan3 -c keep -i ID -t %s annotation.txt
table2itol.R -D plan4 -a -c factor -i ID -l Genus -t %s -w 0 annotation.txt
Stat each level
echo -e ‘Taxonomy\tKnown\tNew’ > tax.stat
for i in seq 2 8;do
head -n1 tax.txt|cut -f ${i}|tr ‘\n’ ‘\t’ >> tax.stat
tail -n+2 tax.txt|cut -f KaTeX parse error: Expected group after '_' at position 14: {i}|grep -v '_̲_‘|sort|uniq -c|wc -l|tr ‘\n’ ‘\t’ >> tax.stat
tail -n+2 tax.txt|cut -f KaTeX parse error: Expected group after '_' at position 11: {i}|grep '_̲_’|wc -l >> tax.stat; done
cat tax.stat
tail -n+2 tax.txt|cut -f3|sort|uniq -c|awk ‘{print $2"\t"$1}’|sort -k2,2nr > count.phylum
cat count.phylum
cd …/…
#### 5.7 功能注释eggnog/dbcan/arg/antismash
基因注释
mkdir -p temp/prodigal
conda activate eggnog
prodigal -v # V2.6.3
50g, 31s, 4m
time tail -n+2 result/metadataS.txt|cut -f1|rush -j 10
“prodigal
-i temp/drep95/dereplicated_genomes/{1}.fa
-o temp/prodigal/{1}.gff
-a temp/prodigal/{1}.faa
-d temp/prodigal/{1}.ffn
-p single -f gff”
seqkit stat temp/prodigal/*.ffn | sed ‘s/temp/prodigal///;s/.ffn//;s/[[:blank:]]{1,}/\t/g’ | cut -f 1,3-
> result/prodigal.txt
碳水化合物注释
mkdir -p temp/dbcan3 result/dbcan3
time tail -n+2 result/metadataS.txt|cut -f1|rush -j 9
“diamond blastp
--db ${db}/dbcan3/CAZyDB
--query temp/prodigal/{1}.faa
--outfmt 6 --threads 8 --quiet --log
--evalue 1e-102 --max-target-seqs 1 --sensitive
--block-size 6 --index-chunks 1
--out temp/dbcan3/{1}_diamond.f6”
wc -l temp/dbcan3/*.f6|head -n-1|awk ‘{print $2"\t"$1}’|cut -f3 -d ‘/’|sed ‘s/_diamond.f6//’|sed ‘1 i ID\tCAZy’|less -S > result/dbcan3/gene.count
format blast2genelist
for i in tail -n+2 result/metadataS.txt|cut -f1;do
format_dbcan3list.pl
-i temp/dbcan3/
i
d
i
a
m
o
n
d
.
f
6
−
o
t
e
m
p
/
d
b
c
a
n
3
/
{i}_diamond.f6 \ -o temp/dbcan3/
idiamond.f6 −otemp/dbcan3/{i}.list
done
CAZy type count
for i in tail -n+2 result/metadataS.txt|cut -f1;do
tail -n+2 temp/dbcan3/${i}.list|cut -f2|sort|uniq -c|awk '{print $2"\t"KaTeX parse error: Expected 'EOF', got '}' at position 2: 1}̲'|sed "1 i CAZy…{i}"|less -S > temp/dbcan3/${i}_CAZy.tsv
done
merge2table
conda activate humann3
humann_join_tables
--input temp/dbcan3/ --file_name CAZy
--output result/dbcan3/cazy.txt
csvtk -t stat result/dbcan3/cazy.txt
merge to level1
paste <(cut -f1 result/dbcan3/cazy.txt) <(cut -f1 result/dbcan3/cazy.txt|tr ‘0-9’ ’ '|sed ‘s/ //g’) | sed ‘1 s/\tCAZy/\tLevel1/’ > result/dbcan3/cazy.L1
summarizeAbundance.py
-i result/dbcan3/cazy.txt
-m result/dbcan3/cazy.L1
-c 2 -s ‘,’ -n raw --dropkeycolumn
-o result/dbcan3/sum
基因相似度
echo -e ‘Name\tCAZy\tIdentity\tGenome’ > result/dbcan3/identity.txt
for i in tail -n+2 result/metadataS.txt|cut -f1;do
csvtk -t replace -f 2 -p “\d+” -r “” temp/dbcan3/
i
.
l
i
s
t
∣
u
n
i
q
∣
t
a
i
l
−
n
+
2
∣
s
e
d
"
s
/
{i}.list | uniq | tail -n+2 | sed "s/
i.list∣uniq∣tail−n+2∣sed"s//\t${i}/" >> result/dbcan3/identity.txt
done
csvtk -t stat result/dbcan3/identity.txt
sp_boxplot.sh -f result/dbcan3/identity.txt -m T -F CAZy -d Identity
耐药基因
mkdir -p temp/card result/card
conda activate rgi6
load database 加载数据库
rgi load -i ${db}/card/card.json
--card_annotation ${db}/card/card.fasta --local
Annotation 蛋白注释
默认为0, --include_loose 可极大增加结果，519/4657=11.14%; --exclude_nudge结果不变，但jason为空
time for i in tail -n+2 result/metadataS.txt|cut -f1;do
i=X004b2
cut -f 1 -d ’ ’ temp/prodigal/KaTeX parse error: Undefined control sequence: \* at position 18: …}.faa | sed 's/\̲*̲//' > temp/prod…{i}.fa
rgi main
--input_sequence temp/prodigal/protein_
i
.
f
a
−
−
o
u
t
p
u
t
f
i
l
e
t
e
m
p
/
c
a
r
d
/
{i}.fa \ --output_file temp/card/
i.fa −−outputfiletemp/card/{i}
--input_type protein --clean
--num_threads 8 --alignment_tool DIAMOND > temp/log 2>&1
done
# 附录：常见分析问题和补充代码
## 计算时间统计表
在60核(p)及以上服务器，单样本3个并行推荐16p，混合组装分箱推荐32p。
小数据：2个7.5M样本，在72个核、512GB服务器上测试。
大数据：20个10G样本，在192个核、2TB大内存服务器上测试。
| 步骤 | 数据小(6Gx3) | 数据(10Gx20) | 备注 |
| --- | --- | --- | --- |
| fastqc | 33s | 15m | |
| seqkit | 2s | 15m | |
| fastp | 2s | 40m | |
| kneaddata | 25s | 5h | |
| humann | 34m | 30h | |
| megahit | 39s | 15h | |
| binning | 6h | 16h | -concoct |
| binrefine | 2h | 3h | |
| coverm | 4m | 30m | |
注：m为分，h为小时，d为天
## 补充代码Supplementary scripts
\*\*for循环批量处理样本列表\*\*
# 基于样本元数据提取样本列表命令解析
# 去掉表头
tail -n+2 result/metadata.txt
# 提取第一列样本名
tail -n+2 result/metadata.txt|cut -f1
# 循环处理样本
for i in `tail -n+2 result/metadata.txt|cut -f1`;do echo "Processing "$i; done
# ` 反引号为键盘左上角Esc键下面的按键，一般在数字1的左边，代表运行命令返回结果
## 质控去宿主KneadData
\*\*双端序列质控后是否配对的检查\*\*
双端序列质控后序列数量不一致是肯定出错了。但即使序列数量一致，也可能序列不对。在运行metawrap分箱时会报错。可以kneaddata运行时添加--reorder来尝试解决。以下提供了检查双端序列ID是否配对的比较代码
**网上学习资料一大堆，但如果学到的知识不成体系，遇到问题时只是浅尝辄止，不再深入研究，那么很难做到真正的技术提升。**
**需要这份系统化的资料的朋友，可以添加V获取：vip1024b （备注运维）**
**一个人可以走的很快，但一群人才能走的更远！不论你是正从事IT行业的老鸟或是对IT行业感兴趣的新人，都欢迎加入我们的的圈子（技术交流、学习资源、职场吐槽、大厂内推、面试辅导），让我们一起学习成长！**
an3/cazy.L1
summarizeAbundance.py \
-i result/dbcan3/cazy.txt \
-m result/dbcan3/cazy.L1 \
-c 2 -s ',' -n raw --dropkeycolumn \
-o result/dbcan3/sum
# 基因相似度
echo -e 'Name\tCAZy\tIdentity\tGenome' > result/dbcan3/identity.txt
for i in `tail -n+2 result/metadataS.txt|cut -f1`;do
csvtk -t replace -f 2 -p "\d+" -r "" temp/dbcan3/${i}.list | uniq | tail -n+2 | sed "s/$/\t${i}/" >> result/dbcan3/identity.txt
done
csvtk -t stat result/dbcan3/identity.txt
sp_boxplot.sh -f result/dbcan3/identity.txt -m T -F CAZy -d Identity
耐药基因
mkdir -p temp/card result/card
conda activate rgi6
# load database 加载数据库
rgi load -i ${db}/card/card.json \
--card_annotation ${db}/card/card.fasta --local
# Annotation 蛋白注释
# 默认为0, --include_loose 可极大增加结果，519/4657=11.14%; --exclude_nudge结果不变，但jason为空
time for i in `tail -n+2 result/metadataS.txt|cut -f1`;do
# i=X004b2
cut -f 1 -d ' ' temp/prodigal/${i}.faa | sed 's/\*//' > temp/prodigal/protein_${i}.fa
rgi main \
--input_sequence temp/prodigal/protein_${i}.fa \
--output_file temp/card/${i} \
--input_type protein --clean \
--num_threads 8 --alignment_tool DIAMOND > temp/log 2>&1
done
附录：常见分析问题和补充代码
计算时间统计表
在60核§及以上服务器，单样本3个并行推荐16p，混合组装分箱推荐32p。
小数据：2个7.5M样本，在72个核、512GB服务器上测试。
大数据：20个10G样本，在192个核、2TB大内存服务器上测试。
| 步骤 | 数据小(6Gx3) | 数据(10Gx20) | 备注 |
|---|---|---|---|
| fastqc | 33s | 15m | |
| seqkit | 2s | 15m | |
| fastp | 2s | 40m | |
| kneaddata | 25s | 5h | |
| humann | 34m | 30h | |
| megahit | 39s | 15h | |
| binning | 6h | 16h | -concoct |
| binrefine | 2h | 3h | |
| coverm | 4m | 30m |
注：m为分，h为小时，d为天
补充代码Supplementary scripts
**for循环批量处理样本列表**
基于样本元数据提取样本列表命令解析
去掉表头
tail -n+2 result/metadata.txt
提取第一列样本名
tail -n+2 result/metadata.txt|cut -f1
循环处理样本
for i in tail -n+2 result/metadata.txt|cut -f1;do echo "Processing "$i; done
` 反引号为键盘左上角Esc键下面的按键，一般在数字1的左边，代表运行命令返回结果
质控去宿主KneadData
**双端序列质控后是否配对的检查**
双端序列质控后序列数量不一致是肯定出错了。但即使序列数量一致，也可能序列不对。在运行metawrap分箱时会报错。可以kneaddata运行时添加–reorder来尝试解决。以下提供了检查双端序列ID是否配对的比较代码
网上学习资料一大堆，但如果学到的知识不成体系，遇到问题时只是浅尝辄止，不再深入研究，那么很难做到真正的技术提升。
需要这份系统化的资料的朋友，可以添加V获取：vip1024b （备注运维）
[外链图片转存中…(img-yLqdMhnE-1713321937656)]
一个人可以走的很快，但一群人才能走的更远！不论你是正从事IT行业的老鸟或是对IT行业感兴趣的新人，都欢迎加入我们的的圈子（技术交流、学习资源、职场吐槽、大厂内推、面试辅导），让我们一起学习成长！
扫一扫
2万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
