-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathecological niche theory.r
112 lines (94 loc) · 4.6 KB
/
ecological niche theory.r
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
#Relationships between key factors and motifs
#ecological niche theory
#import abundance data
otu1=read.csv("otu_net1.csv",row.names=1)
otu2=read.csv("otu_net2.csv",row.names=1)
otu3=read.csv("otu_net3.csv",row.names=1)
otu4=read.csv("otu_net4.csv",row.names=1)
otu5=read.csv("otu_net5.csv",row.names=1)
otu6=read.csv("otu_net6.csv",row.names=1)
otu7=read.csv("otu_net7.csv",row.names=1)
otu8=read.csv("otu_net8.csv",row.names=1)
otu9=read.csv("otu_net9.csv",row.names=1)
otu10=read.csv("otu_net10.csv",row.names=1)
otu11=read.csv("otu_net11.csv",row.names=1)
otu12=read.csv("otu_net12.csv",row.names=1)
otu13=read.csv("otu_net13.csv",row.names=1)
otu14=read.csv("otu_net14.csv",row.names=1)
otu15=read.csv("otu_net15.csv",row.names=1)
otu16=read.csv("otu_net16.csv",row.names=1)
otu17=read.csv("otu_net17.csv",row.names=1)
otu18=read.csv("otu_net18.csv",row.names=1)
otu19=read.csv("otu_net19.csv",row.names=1)
otu20=read.csv("otu_net20.csv",row.names=1)
#save all abundance data as otulist
otulist=list(otu1,otu2,otu3,otu4,otu5,otu6,otu7,otu8,otu9,otu10,
otu11,otu12,otu13,otu14,otu15,otu16,otu17,otu18,otu19,otu20)
#rename otulist
for (i in 1:length(otulist)){
names(otulist)[i] <- paste("net", i, sep = "")
}
envnew=read.csv("metadata.csv",row.names=1)
#calculate niche breath
#take TSEA as an example of ecological amplitude
data=c()
for(j in 1:length(otulist)){
otu=otulist[[j]]
niche=c()
for(i in 1:ncol(otu)){
envnew=env[row.names(otu),]
value = weighted.mean(envnew[,"TSEA"],otu[,i])
niche=c(niche,value)
}
niche=as.data.frame(niche)
niche$net= paste("net", j, sep = "")
data=rbind(data,niche)
print(j)
}
write.csv(data,"niche.tsea.csv")
#visualization
p1=ggplot(niche.tsea,aes(x=cycfac, y=niche))+
geom_point(alpha=.7, size=2,color="#223e83")+
ggpubr::stat_cor(label.y=10000)+
ggpubr::stat_regline_equation(label.y=12000)+
geom_smooth(aes(cycfac,niche), method=lm, se=T,color="#e71f19")+
geom_smooth(aes(cycfac,niche), method=loess, se=T,color="#33a02c")+
theme(panel.grid=element_blank(), panel.background=element_rect(fill='transparent', color='black'),panel.border = element_rect(colour = "black", fill=NA, size=1))
p1
ggsave( "tsea.cycfac.pdf",p1,width=5,height=3)
p2=ggplot(niche.tsea,aes(x=facmcom, y=niche))+
geom_point(alpha=.7, size=2,color="#223e83")+
ggpubr::stat_cor(label.y=10000)+
ggpubr::stat_regline_equation(label.y=12000)+
geom_smooth(aes(facmcom,niche), method=lm, se=T,color="#e71f19")+
geom_smooth(aes(facmcom,niche), method=loess, se=T,color="#33a02c")+
theme(panel.grid=element_blank(), panel.background=element_rect(fill='transparent', color='black'),panel.border = element_rect(colour = "black", fill=NA, size=1))
p2
ggsave( "tsea.facmcom.pdf",p2,width=5,height=3)
p3=ggplot(niche.tsea,aes(x=tranfac, y=niche))+
geom_point(alpha=.7, size=2,color="#223e83")+
ggpubr::stat_cor(label.y=10000)+
ggpubr::stat_regline_equation(label.y=12000)+
geom_smooth(aes(tranfac,niche), method=lm, se=T,color="#e71f19")+
geom_smooth(aes(tranfac,niche), method=loess, se=T,color="#33a02c")+
theme(panel.grid=element_blank(), panel.background=element_rect(fill='transparent', color='black'),panel.border = element_rect(colour = "black", fill=NA, size=1))
p3
ggsave( "tsea.tranfac.pdf",p3,width=5,height=3)
p4=ggplot(niche.tsea,aes(x=trancom, y=niche))+
geom_point(alpha=.7, size=2,color="#223e83")+
ggpubr::stat_cor(label.y=10000)+
ggpubr::stat_regline_equation(label.y=12000)+
geom_smooth(aes(trancom,niche), method=lm, se=T,color="#e71f19")+
geom_smooth(aes(trancom,niche), method=loess, se=T,color="#33a02c")+
theme(panel.grid=element_blank(), panel.background=element_rect(fill='transparent', color='black'),panel.border = element_rect(colour = "black", fill=NA, size=1))
p4
ggsave( "tsea.trancom.pdf",p4,width=5,height=3)
p5=ggplot(niche.tsea,aes(x=trancomfac, y=niche))+
geom_point(alpha=.7, size=2,color="#223e83")+
ggpubr::stat_cor(label.y=10000)+
ggpubr::stat_regline_equation(label.y=12000)+
geom_smooth(aes(trancomfac,niche), method=lm, se=T,color="#e71f19")+
geom_smooth(aes(trancomfac,niche), method=loess, se=T,color="#33a02c")+
theme(panel.grid=element_blank(), panel.background=element_rect(fill='transparent', color='black'),panel.border = element_rect(colour = "black", fill=NA, size=1))
p5
ggsave( "tsea.trancomfac.pdf",p5,width=5,height=3)