使用R语言获得16S物种丰度

Posted 2021-04-18 科技记者

 还是获得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 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 -