示例数据
提取码：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)
