R语言之简单批量处理

示例数据
提取码：p3pn

1. 循环计算获取结果

展示并设定工作路径

list.dirs()
setwd('YT')

数据准备

map = read.table('map.txt',row.names = 1, header = T)
nti = read.table('bNTI.txt',row.names = 1, header = T)
rc = read.table('RCbray.txt',row.names = 1,header = T)

设定要考察的所有组合

grps = list('T1','T2','T4','T8','T16',
            c('T1','T2'),c('T2','T4'),c('T4','T8'),c('T8','T16'),
            c('T1','T2','T4'),c('T2','T4','T8'),c('T4','T8','T16'),
            c('T1','T2','T4','T8'),c('T2','T4','T8','T16'),
            c('T1','T2','T4','T8','T16')) #there are 15 grps in total, wtf.

新建空数据框用于存储计算结果

result = matrix(0,nrow = 0, ncol = 6,dimnames = list(c(),
                                                     c("Variable selection",
                                                       "Homogeneous selection",
                                                       "Dispersal limitation",
                                                       "Homogenizing dispersal",
                                                       "Undominated process","summary")))

新建一个函数将两个对称矩阵转换为row-col-value1-value2的线性形式

mtrx2cols = function(m1,m2,val1,val2){
  lt = lower.tri(m1)
  res = data.frame(row = rownames(m1)[row(m1)[lt]],
                   col = rownames(m1)[col(m1)[lt]],
                   val1 = m1[lt], val2= m2[lt])
names(res)[3:4] = c(val1,val2)
return(res)
}

循环提取各组合对应的样品

for (i in seq(15)){
  # extract samples from the whole data
  submap = subset(map,time %in% grps[[i]])  #A %in% B，表示从A中提取出B对应的部分样品
  subnti = nti[rownames(submap),rownames(submap)]
  subrc  = rc[rownames(submap),rownames(submap)]
  #now the samples in subnti and subrc are in same order
  #将对称矩阵转换为row-col-value的线性形式
  res = matrix2cols(subnti,subrc,'betaNTI','RCbray')
  
  #循环过程中的中间变量，用于存储本循环的计算结果，后续会将这个中间变量与恒定变量合并，此中间变量的使命就完成了。
  a = matrix(0,nrow = 1, ncol = 6,dimnames = list(c(),
                                                  c("Variable selection",
                                                    "Homogeneous selection",
                                                    "Dispersal limitation",
                                                    "Homogenizing dispersal",
                                                    "Undominated process","summary")))
  
 #根据两列值的大小，分别对应至相应的过程，并计算相对比例
  num = dim(res)[1]#获取总的行数
  a[,1] = sum(res$betaNTI>2)/num#获取betaNTI > 2的数据行所占的比例
  a[,2] = sum(res$betaNTI<(-2))/num
  a[,3] = sum(abs(res$betaNTI)<2 & res$RCbray>0.95)/num
  a[,4] = sum(abs(res$betaNTI)<2 & res$RCbray<(-0.95))/num
  a[,5] = sum(abs(res$betaNTI)<2 & abs(res$RCbray)<0.95)/num
  a[,6] = num
 #每个组合均包含多个元素，这里讲该组合中的所有元素拼接在一起，作为该组的标记
  rownames(a) = paste(grps[[i]][1:length(grps[[i]])],collapse = '-')
  
  result = rbind(result,a)#将中间变量和恒定变量按行合并
  
}
#导出本次循环的计算结果
write.csv(result,'../YT.summary.csv')
#返回上一级路径，重新开始执行另一个文件
setwd('../')

---

2. 后续分析，可视化

数据准备

cs = read.csv('CS.summary.csv')
cq = read.csv('CQ.summary.csv')
yt = read.csv('YT.summary.csv')
#合并，并定义组别
ecology = rbind(cs,cq,yt)
ecology$site = rep(c('CS','CQ','YT'),each = 10)

library('reshape2')
dat = melt(ecology)#变为列表的形式，便于绘图
names(dat) = c('pairs','site','process','value')
#去除部分不需要的数据，`!=`表示不等于
dat = subset(dat, process != 'summary')

初步可视化

library('ggplot2')
library('scales')

p <- ggplot(dat,aes(pairs,value,fill = process)) +
  theme_bw() +
  coord_flip() + #坐标轴反转
  scale_y_continuous(labels = percent) + #设定y轴为百分比形式
  facet_grid(~site) +#按*site*进行分页
  labs(x = '', y = '', fill = '') +#坐标轴标题为空
  geom_bar(stat = 'identity',position = 'fill') + #绘制堆积条形图
  scale_fill_brewer(palette = 'Paired') + #个性化配色方案
  theme(panel.spacing = unit(0.5,'cm'),#设置分页面之间的间距大小
        axis.text = element_text(face = 'bold',color = 'black',size = 6),
        legend.text = element_text(face = 'bold',color = 'black',size = 8),
        strip.text = element_text(face = 'bold',color = 'black',size = 10))

ggsave(p,filename = 'ecological-process.jpg',width = 8,height = 3,dpi = 600)