基于测序数据，使用 kraken2 进行 reads 的物种注释，使用 braken进行丰度估计，之后基于丰度进行物种多样性分析及差异物种分析。

数据库构建

运行 kraken2 和 bracken 之前需要提前进行数据库构建

kraken2 数据库地址：
https://benlangmead.github.io/aws-indexes/k2

# kraken2 database

# 方法一：直接下载官方提供的数据库（包括kraken2及bracken），解压即可
wget https://genome-idx.s3.amazonaws.com/kraken/k2_pluspf_20210517.tar.gz
wget https://genome-idx.s3.amazonaws.com/kraken/k2_pluspf_8gb_20210517.tar.gz
wget https://genome-idx.s3.amazonaws.com/kraken/k2_standard_20210517.tar.gz

# 方法二：一步构建kraken标准数据库
# standard : archaea, bacteria, viral, plasmid, human1, UniVec_Core
kraken2-build --standard    --threads 4 --db ./standardDB


# 方法3：手动构建数据库
## step1 从 NCBI 下载分类文件
kraken2-build --download-taxonomy --db  ./K2db

## step2 下载序列数据
kraken2-build  --download-library UniVec_Core  --threads 4  --db  ./K2db
kraken2-build  --download-library UniVec  --threads 4  --db  ./K2db

## step3 构建数据库
kraken2-build  --build --threads 4  --db  ./K2db


# bracken database
bracken-build  -d  K2db/ -t 4

物种组成及丰度估计

我们直接使用官方提供的数据库进行物种组成分析
数据库下载自以下网址：
https://benlangmead.github.io/aws-indexes/k2
运行 kraken2，进行 reads 物种注释:

kraken2 \
--threads 4 \
--quick \ # 快速模式
--paired \ # 输入数据为双端模式
--db ../../Database/K2/ \ # 数据库路径
--report A1.kreport \ # 输出，report文件
--output A1.kraken \ # 输出，reads注释结果
./A1_1.clean.fq.gz \ # 输入，fq1
./A1_2.clean.fq.gz # 输入，fq2

运行 bracken 进行各个分类水平物种丰度估计：

# 创建输出目录，进行species水平分析
mkdir out
# 运行bracken
bracken \
-d ../Database/K2/ \ # 数据库路径
-i A1.kreport \ # 输入, 为kraken2 report文件
-o A1.bracken.S \ # 输出，物种丰度估计表格
-w A1.bracken.S.kreport \ # 输出，kraken2 report格式文件
-l S \ # 分类水平，如 D,P,C,O,F,G,S
-t 4 # 线程数

将各样品丰度结果汇总：

# report文件合并成biom格式
kraken-biom \
./out/*.kreport \ # 输入,report文件
--max D \ # 最高物种分类
-o ./out_S/S.biom # 输出， biom文件

# biom转count表格
convert \
-i ./out/S.biom \ # 输入，biom文件
-o ./out/S.count.tsv.tmp \ # 输出，丰度表格
--to-tsv \ # 指定输出格式
--header-key taxonomy # 输出分类信息

# 输出文件格式调整，补全物种名
sed 's/; g__\([ˆ;]\+\); s__/; g__\1; s__\1 /' ./out/S.count.tsv.tmp \
> ./out/S.taxID.count.tsv
## taxonID 替换回拉丁名
sed '/ˆ#/! s/ˆ[0-9]\+\t\(.*[A-Za-z]\+__\([ˆ;]\+\)\)$/\2\t\1/' \
./out/S.taxID.count.tsv > ./out/S.taxName.count.tsv
## 保留丰度信息，用于后续绘图
sed '1d; 2s/ˆ#//' ./out/S.taxName.count.tsv | \
awk -F "\t" -v 'OFS=\t' '{$NF = ""; { print $0 }}' | \
sed 's/\t$//' > ./out/S.count.tsv

# 绘制barplot 和 heatmap, 输出相对丰度
Rscript \
./Barplot.R \
out/S.count.tsv 20 \
out/S.count.out

输出结果：
S.taxName.count.tsv 为包含物种分类信息的丰度文件
S.count.tsv 为去除物种分类信息的丰度文件
S.count.out.relative.table 相对丰度表
S.count.out.barplot.pdf 丰度barplot
S.count.out.heatmap_cluster.pdf 丰度热图
S.count.out.heatmap_original.pdf 丰度热图

欢迎关注Bioinfor 生信云