使用R语言获得16S物种丰度
时间:2022-07-22
本文章向大家介绍使用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 -
- JavaScript 教程
- JavaScript 编辑工具
- JavaScript 与HTML
- JavaScript 与Java
- JavaScript 数据结构
- JavaScript 基本数据类型
- JavaScript 特殊数据类型
- JavaScript 运算符
- JavaScript typeof 运算符
- JavaScript 表达式
- JavaScript 类型转换
- JavaScript 基本语法
- JavaScript 注释
- Javascript 基本处理流程
- Javascript 选择结构
- Javascript if 语句
- Javascript if 语句的嵌套
- Javascript switch 语句
- Javascript 循环结构
- Javascript 循环结构实例
- Javascript 跳转语句
- Javascript 控制语句总结
- Javascript 函数介绍
- Javascript 函数的定义
- Javascript 函数调用
- Javascript 几种特殊的函数
- JavaScript 内置函数简介
- Javascript eval() 函数
- Javascript isFinite() 函数
- Javascript isNaN() 函数
- parseInt() 与 parseFloat()
- escape() 与 unescape()
- Javascript 字符串介绍
- Javascript length属性
- javascript 字符串函数
- Javascript 日期对象简介
- Javascript 日期对象用途
- Date 对象属性和方法
- Javascript 数组是什么
- Javascript 创建数组
- Javascript 数组赋值与取值
- Javascript 数组属性和方法
- 经典的SparkSQL/Hive-SQL/MySQL面试-练习题
- codeforces 1249E(dp)
- Redis-KV数据库Java连接以及Jedis包的使用
- codeforces 1203D1(暴力)
- codeforces 1366B(线段相交)
- 一文搞懂Python自动化测试框架
- codeforces 1005D(数学)
- JSP开发简单实例演示
- Linux笔记【003】| Linux系统目录结构与基本命令
- codeforces1322A(括号匹配)
- codeforces 1296D(贪心)
- codeforces 1399D
- JSP开发之JSTL介绍和使用
- codeforces 1283E(贪心)
- codeforces1216C (矩形面积交)