示例數(shù)據(jù)
提取碼：p3pn
1. 循環(huán)計算獲取結果
展示并設定工作路徑
list.dirs()
setwd('YT')
數(shù)據(jù)準備
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.
新建空數(shù)據(jù)框用于存儲計算結果
result = matrix(0,nrow = 0, ncol = 6,dimnames = list(c(),
c("Variable selection",
"Homogeneous selection",
"Dispersal limitation",
"Homogenizing dispersal",
"Undominated process","summary")))
新建一個函數(shù)將兩個對稱矩陣轉換為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)
}
循環(huán)提取各組合對應的樣品
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')
#循環(huán)過程中的中間變量，用于存儲本循環(huán)的計算結果，后續(xù)會將這個中間變量與恒定變量合并疯潭，此中間變量的使命就完成了祠锣。
a = matrix(0,nrow = 1, ncol = 6,dimnames = list(c(),
c("Variable selection",
"Homogeneous selection",
"Dispersal limitation",
"Homogenizing dispersal",
"Undominated process","summary")))
#根據(jù)兩列值的大小谓谦，分別對應至相應的過程吨凑，并計算相對比例
num = dim(res)[1]#獲取總的行數(shù)
a[,1] = sum(res$betaNTI>2)/num#獲取betaNTI > 2的數(shù)據(jù)行所占的比例
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)#將中間變量和恒定變量按行合并
}
#導出本次循環(huán)的計算結果
write.csv(result,'../YT.summary.csv')
#返回上一級路徑，重新開始執(zhí)行另一個文件
setwd('../')
2. 后續(xù)分析疏哗，可視化
數(shù)據(jù)準備
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)#變?yōu)榱斜淼男问剑阌诶L圖
names(dat) = c('pairs','site','process','value')
#去除部分不需要的數(shù)據(jù)返奉，`!=`表示不等于
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)
