由于最近遇到一个宏基因项目,第一次分析注释到的物种数目较少,因此更换数据库试试看看!
于是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'")
参考: