使用R语言获得16S物种丰度
还是获得16S物种丰度得老问题,最近在一台新机器上安装qiime1,发现有报错,对于这种停止维护的软件,也是正常现象吧,于是想别的办法解决,恰巧最近读R几本R语言的入门书,发现prop.table()这个函数是可以实现相关功能的,于是学习使用下。可能你早已会做这个啦,还是分享一下,看看有没有人需要。
从qiime2的结果开始
qiime2可以按级别导出物种的丰富计数,通过view.qiime2.org或者view.qiime2.cn(后者是我自己镜像的前面一个网站) 查看taxa-bar-plots.qzv导出csv文件获得。或者通过下面的命令导出,结果应该是相同的。
# 按门和属水平合并,并统计
#门
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语言水平不高,一定有优化的余地,也欢迎指出。
#读取文件
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 p<- paste0(strsplit(p, ' ')[[1]][1], "属")
}
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)
就这样啦!