【R分享|实战】科白君教你定义函数批量计算微生物α多样性

“ 突破。”   --科白君

"R分享实战"专刊·第9篇

  编辑 | 科白维尼

  500字 | 10分钟阅读

本期推送内容

最近在处理一批数据量较大的数据，需要多次读入+导出，并且还需要计算微生物α多样性指标。为了减少重复代码及减轻工作量。尝试定义函数和编写循环解决当前问题。这里主要涉及内容：批量导入数据+定义计算α多样性指标的函数+批量导出。下面直接对代码进行讲解~由于目前数据还未发表，这里重新定义了函数来生成随机矩阵模拟OTU数据（这里非常感谢师兄的指点和帮助）。

直接上代码和讲解~

Code：

#加入以下R包是为了运行read.xlsx这个函数
library(xlsx)
library(openxlsx)
library(xlsxjars)
library(rJava)

##### 定义函数随机生成N个矩阵（n行，n列）######
#由于数据还未发表，这里定义函数生成随机矩阵，来模拟otu数据
set.seed(1) #只要涉及随机数都需要设定种子
random_matrix <- function(n, nrow, ncol) { #主要的函数体
  df = list() #生成一个空的df为list()结构
  otu = paste0("otu_", 1:ncol) #设置OTU id
  sample = paste0("sample_", 1:nrow) #设置samples id
  for (i in 1:n) { #第一次循环，是为了生成矩阵
    df[[i]] = matrix(NA,nrow,ncol) #循环生成N个矩阵并装进df这个list()中
    df[[i]] = data.frame(df[[i]]) #将list转成数据框结构，以便添加行名和列名
    colnames(df[[i]]) = otu #对df所有数据框赋值列名
    rownames(df[[i]]) = sample #赋值行名
    for (j in seq_along(df[[i]])) { #嵌套循环
    df[[i]][j] = sample(0:300,nrow,replace = T) #前面生成的矩阵是NA空值，这里循环生成N行、N列0~300的随机数
    } #list()中 df[[]][] 前面的双括号表示读取到哪个矩阵，后面单括号表示读取该矩阵的哪一列
  }
  return(df) #返回给df
}

#生成6个矩阵，每个矩阵10行，200列
data <- random_matrix(n=6, nrow= 10, ncol= 200) # n nrow ncol 可以根据自己的喜好，随便设定
#直接以不同sheet的形式一次性导出至一个excel中
write.xlsx(data, "mydata.xlsx", append=T, row.names=T)

###### 批量导入同excel中的多个sheet #######
library(tidyverse)
sheet.index <- c(1:6) #1至6个sheet 6可以才成n 这里根据自己的excel中的情况修改
re <- map(sheet.index, ~ read.xlsx("mydata.xlsx", sheet=. )) #一次性导入excel中1到6个sheet
str(re) #可以查看到导入的数据集是list结构，其中包含6个data.frame，也就是6个sheet

###### 定义函数计算alpha多样性 ######

alpha_diversity <- function(x){
  library(vegan) #由于α多样性有些指标需要基于vegan包 这里加载R包
  x = data.matrix(x) #转换矩阵
  shannon = diversity(x,index="shannon") #香浓指标
  richness = specnumber(x) #丰富度
  simpson = diversity(x, index = "simpson") #辛普森指数
  Chao1 = estimateR(x)[2,] #Chao1指数
  eveness = shannon / log(richness) #均匀度
  a = data.frame(shannon, richness,simpson,Chao1,eveness) #用数据框打包计算的结果
  return(a)
}

#构建一个list来装对应不同sheet中计算alpha多样性的结果
b <- list()
#循环计算sheet的过程
for(i in seq_along(re)){
  b[[i]] <- alpha_diversity(re[[i]]) #循环计算re这个list()中每个数据框并装到b列表中
}
#查看结果
head(b)

###### 批量导出多个sheet到同一个excel ######
for(i in 1:6){
  sample <- paste0("sample_", 1:nrow(re[[i]])) #这里生成对应数据集中行行数的样本名
  rownames(b[[i]]) <- sample #并赋值给b数据集中每个data.frame作为行名
  write.xlsx(b, file = "C:/Users/Zz/Desktop/a.xlsx", append=T, row.names=T) #输出所有计算结果并导出至一个excel中
}

用R markdown输出结果：

这里是循环得到虚拟的数据集。

这里是循环得到α多样性的结果。