还是获得16S物种丰度得老问题，最近在一台新机器上安装qiime1，发现有报错，对于这种停止维护的软件，也是正常现象吧，于是想别的办法解决，恰巧最近读R几本R语言的入门书，发现prop.table()这个函数是可以实现相关功能的，于是学习使用下。可能你早已会做这个啦，还是分享一下，看看有没有人需要。
从qiime2的结果开始
qiime2可以按级别导出物种的计数，通过view.qiime2.org或者view.qiime2.cn（后者是我自己镜像的前面一个网站） 查看taxa-bar-plots.qzv导出csv文件获得。或者通过下面的命令导出，结果应该是相同的。
代码语言：javascript复制# 按门和属水平合并，并统计
#门
qiime taxa collapse
--i-table table_final1.qza
--i-taxonomy taxonomy.qza
--p-level 2
--o-collapsed-table table-l2.qza
#科
qiime taxa collapse
--i-table table_final1.qza
--i-taxonomy taxonomy.qza
--p-level 5
--o-collapsed-table table-l5.qza
#属
qiime taxa collapse
--i-table table_final1.qza
--i-taxonomy taxonomy.qza
--p-level 6
--o-collapsed-table table-l6.qza
#species
qiime taxa collapse
--i-table table_final1.qza
--i-taxonomy taxonomy.qza
--p-level 7
--o-collapsed-table table-l7.qza
#导出 tsv
for file in ./table-l*.qza
do
base=$(basename $file .qza)
echo $base
qiime tools export
--input-path $file
--output-path $base
biom convert --to-tsv -i $base/feature-table.biom -o $base.tsv
done上R语言
粗略看了下结果，主要是有一两个属竟然分属于两个高级别分类的情况，比如梭菌属等，需要先用代码合并下，并格式化命名，还有就是有的分类没有到达这个级别的，也需要稍作处理。我的代码如下，可能存在错误，欢迎交流指正。由于R语言水平不高，一定有优化的余地，也欢迎指出。
代码语言：javascript复制#读取文件
sample <- paste('table-l6','.tsv', sep = "")
df <- read.table(sample, header = TRUE, sep = 't',comment.char="",skip=1, stringsAsFactors = FALSE)
#整理属名，去除多余信息
genus <- c()
for (j in 1:length(df[, 1])) {
p <- paste(strsplit(df[, 1][j], "__")[[1]][7], j)
if(strsplit(p, ' ')[[1]][1]=="NA"){
if(strsplit(df[, 1][j], "__")[[1]][6]==';'){
p <- df[, 1][j]
}
else if(strsplit(df[, 1][j], "__")[[1]][6]==';g'){
p <- df[, 1][j]
}
else if(strsplit(p, ' ')[[1]][1]=="NA") p <- paste(strsplit(strsplit(df[, 1][j], "__")[[1]][6], ';g')[[1]][1], j)
else {p <- paste(strsplit(strsplit(df[, 1][j], "__")[[1]][6], ';')[[1]][1], j) }
}
genus <- c(genus, p)
}
#属重命名
df <- df[-1]
row.names(df) <- genus
#合并相同属
get_genus_summary <- function(df, bacterium){
Bac_name <- df[grepl(bacterium, rownames(df)),]
Bac_name <- sapply(Bac_name, as.numeric)
bac <- colSums(Bac_name)
}
df_new <- data.frame()
for (bact in 1:length(row.names(df))) {
if(grepl("\[",row.names(df)[bact])) {
bac_name <- strsplit(strsplit(strsplit(row.names(df)[bact], ' ')[[1]][1], '\[')[[1]][2], "\]")[[1]][1]
}else {
bac_name <- strsplit(row.names(df)[bact], ' ')[[1]][1]
}
if (length(row.names(df[grepl(bac_name, rownames(df)),])) > 1) {
if(!(bac_name %in% row.names(df_new))) {
bac <- get_genus_summary(df, bac_name)
df_new <- rbind(df_new, bac)
row.names(df_new)[length(row.names(df_new))] <- bac_name
}else next
}
else {
df_new <- rbind(df_new, df[bact,])
row.names(df_new)[length(row.names(df_new))] <- bac_name
}
}
#获得比例数据，先转化成矩阵，2代表以列求和
df_new <- prop.table(data.matrix(df_new), 2)就这样啦！
- END -
