专栏首页科技记者使用R语言获得16S物种丰度

使用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 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 -

本文分享自微信公众号 - 科技记者(kejijizhe),作者:zd200572

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-07-26

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • qiime2学习笔记(一)

    最近学习肠道微生物方面的知识,有一部分测序数据需要学习分析。鉴于qiime已经升级为qiime2,还有了图形版本,真是越来越人性化了,但是图形版本还处于原型阶段...

    用户1075469
  • ubiome类似数据dada2处理探索3

    我简单处理了下otu序列和表,使它们能导入qiime2,应该是一行shell代码解决的,shell水平不行,python来顶了。

    用户1075469
  • QIIME2学习笔记(二)

    这次不用测试数据了,用实际数据跑一下,所以同样重复之前的步骤,把fastq文件压缩下,然后,生成样本数据列表(ps.不知道fastq文件不压缩可不可以用,有空试...

    用户1075469
  • 聊聊skywalking的HTTPAccessLog

    skywalking-6.6.0/oap-server/server-core/src/main/java/org/apache/skywalking/oap/...

    codecraft
  • 编程小白 | 每日一练(138)

    这道理放在编程上也一并受用。在编程方面有着天赋异禀的人毕竟是少数,我们大多数人想要从编程小白进阶到高手,需要经历的是日积月累的学习,那么如何学习呢?当然是每天都...

    C语言入门到精通
  • cmake 3.5:find_package(HDF5) 指定HDF5_ROOT无效问题

    我们知道cmake提供了FindHDF5.cmake(位置:$cmake_root/Modules)模块用于搜索HDF5组件。 通过查看FindHDF5.c...

    用户1148648
  • 精品教学案例 | 金融诈骗数据分析与预测

    本案例适合作为大数据专业数据科学导引、数据清洗或机器学习实践课程的配套教学案例。通过本案例,能够达到以下教学效果:

    数据酷客
  • 互联网巨头都在争夺用户时长,为什么说百度已赢得先机?

    移动互联网已进入成熟期。当用户基数难以再大幅增长时,让存量用户花更多时间在应用上,成为众多互联网公司的新目标。最近几年,许多互联网公司都提出了内容战略,不只是百...

    罗超频道
  • 「分布式」实现分布式锁的正确姿势

    最近看到好多博主都在推分布式锁,实现方式很多,基于db、redis、zookeeper。zookeeper方式实现起来比较繁琐,这里我们就谈谈基于redis实现...

    一个程序员的成长
  • JavaScript图片隐写术 – 图片加入看不见的版权文字

    隐写术算是一种加密技术,权威的wiki说法是“隐写术是一门关于信息隐藏的技巧与科学,所谓信息隐藏指的是不让除预期的接收者之外的任何人知晓信息的传递事件或者信息的...

    Javanx

扫码关注云+社区

领取腾讯云代金券