kraken2注释+braken格式转换

由于最近遇到一个宏基因项目，第一次分析注释到的物种数目较少，因此更换数据库试试看看！

于是kraken上榜，OK，let's  try 一下子！

前提是已经安装了kraken和braken，运行命令比较简单

1.注释

##基于cleanreads进行注释：
cd /kraken2/test_kraken2
mkdir sample1_3
/kraken2/kraken2-2.1.2/kraken2  \
	--db /kraken2/database/ \
	--paired  --threads 24 \
	/cleandata/sample1_3_R1.clean.fq.gz \
	/cleandata/sample1_3_R2.clean.fq.gz \
	--output sample1_3/GOE1_3.kraken --report sample1_3.report

查看一下sample1_3.report文件格式：

第一列是相似百分比，最后一列物种，倒数第二列是层级

共有8个样本，预期得到的格式为：行为物种，列为样本名，中间值为物种丰度，如图所示

上代码：

2.格式修改

for i in `cat /kraken2/test_kraken2/list`;do
	source /cqproj2/Software/Binning-software/miniconda3/bin/activate braken 
	bracken -d  /kraken2/database/ -i  /kraken2/test_kraken2/sample_out/${i}.report  -o /kraken2/test_kraken2/sample_out/3/${i}.bracken.S  \
       -l S -t 4  -r 150 \
        -w /kraken2/test_kraken2/sample_out/3/out/${i}.bracken.S.kreport 
##转换为biom格式
	kraken-biom  --max D  /kraken2/test_kraken2/sample_out/3/out/*.bracken.S.kreport  -o /kraken2/test_kraken2/sample_out/3/out/S.biom
	biom convert  -i /kraken2/test_kraken2/sample_out/3/out/S.biom  --to-tsv --header-key taxonomy -o /kraken2/test_kraken2/sample_out/3/out/S.count.tsv.tmp;done

##删除第一行
sed -i '1d' S.count.tsv.tmp

备注：list包括样本信息，每行一个样本

***计算相对丰度：

import pandas as pd
import glob
df = pd.read_csv("/kraken2/test_kraken2/sample_out/3/out/count.txt", sep='\t',index_col=0)
# 提取 taxonomy 列
taxonomy = df['taxonomy']
# 去除 taxonomy 列以计算相对丰度
df_without_taxonomy = df.drop(columns=['taxonomy'])
# 计算每行的总和
row_sums = df_without_taxonomy.sum(axis=1)

# 过滤掉总和小于70的行（这个看自己需求，也可以不过滤）
df_filtered = df_without_taxonomy.loc[row_sums > 10, :]

# 计算每列的总和
column_sums = df_filtered.sum(axis=0)


# 计算相对丰度
relative_abundance_df = df_without_taxonomy.div(column_sums, axis=1)

# 添加 taxonomy 列
relative_abundance_df.insert(0, 'taxonomy', taxonomy)

relative_abundance_df.to_csv('/output//assign_tax/tax_abund.txt', sep = "\t",index=False)
print("相对丰度表已成功计算并保存为 'relative_abundance.csv'")

参考：

宏基因组测序分析（三）物种组成及丰度估计

krakentools处理kraken2+bracken结果