-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathFigure3.R
124 lines (92 loc) · 4.81 KB
/
Figure3.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
113
114
115
116
117
118
119
120
121
122
123
124
rm(list=ls());gc()
library(MCMCglmm)
load('output_files_elife//RData/Analysis_Cemee_Pop_WI_NaCl.RData')
# Load EV_plast
manova.env=new.env()
load('output_files_elife/RData/Pi_Theta_Plasticity_A6140_MANOVA.RData',envir=manova.env)
EV_plast=manova.env$EV_plast
rm(manova.env)
# Load the null matrices
null_matrix=new.env()
load("output_files_elife/RData_Random/Random_G_Analysis_Cemee_Pop_WI_A6140_NaCl_LIGHT.RData",envir=null_matrix)
A6140_NaCl_NULL = null_matrix$df_G1
rm(null_matrix);gc()
angle_theta <- function(x, y) {
dot.prod <- x %*% y
norm.x <- norm(x, type = "2")
norm.y <- norm(y, type = "2")
theta <- 180/pi * as.numeric(acos(dot.prod/(norm.x * norm.y)))
as.numeric(theta)
}
vect_rand_piE_NaCl <- NULL
vect_rand_piE_NaCl_NULL <- NULL
vect_rand_thetaE_NaCl <- NULL
vect_rand_thetaE_NaCl_2 <- NULL
vect_rand_thetaE_NaCl_3 <- NULL
n_iter=1000
spld_idx <- sample(1:nrow(VCV_mat_NGM[[1]]$VCV_Mat),n_iter)
k=0
for(i in spld_idx){
k=k+1
# PI
temp_NaCl <- matrix(VCV_mat_NaCl[[1]]$VCV_Mat[i,1:49],7,7)
temp_NaCl_NULL <- matrix(A6140_NaCl_NULL[k,],7,7)
temp_gmax_proj <- (t(EV_plast[,1])%*%(temp_NaCl/2)%*%EV_plast[,1])/sum(EV_plast[,1]^2)
temp_gmax_proj_NULL <- (t(EV_plast[,1])%*%(temp_NaCl_NULL/2)%*%EV_plast[,1])/sum(EV_plast[,1]^2)
vect_rand_piE_NaCl <- c(vect_rand_piE_NaCl , temp_gmax_proj/eigen(temp_NaCl/2)$values[1])
vect_rand_piE_NaCl_NULL <- c(vect_rand_piE_NaCl_NULL , temp_gmax_proj_NULL/eigen(temp_NaCl_NULL/2)$values[1])
## Angles
vect_rand_thetaE_NaCl <- c(vect_rand_thetaE_NaCl,angle_theta(EV_plast[,1],eigen(temp_NaCl/2)$vector[,1]))
vect_rand_thetaE_NaCl_2 <- c(vect_rand_thetaE_NaCl_2,angle_theta(EV_plast[,1],eigen(temp_NaCl/2)$vector[,2]))
vect_rand_thetaE_NaCl_3 <- c(vect_rand_thetaE_NaCl_3,angle_theta(EV_plast[,1],eigen(temp_NaCl/2)$vector[,3]))
}
Angle_rand=NULL
for(i in 1:1000){
Angle_rand=c(Angle_rand,angle_theta(runif(7,min=(-1),max=1), runif(7,min=(-1),max=1)))
}
## All angles range [0-90]
Angle_rand[Angle_rand>90]= 180-Angle_rand[Angle_rand>90]
vect_rand_thetaE_NaCl[vect_rand_thetaE_NaCl>90] = 180 - vect_rand_thetaE_NaCl[vect_rand_thetaE_NaCl>90]
vect_rand_thetaE_NaCl_2[vect_rand_thetaE_NaCl_2>90] = 180 - vect_rand_thetaE_NaCl_2[vect_rand_thetaE_NaCl_2>90]
vect_rand_thetaE_NaCl_3[vect_rand_thetaE_NaCl_3>90] = 180 - vect_rand_thetaE_NaCl_3[vect_rand_thetaE_NaCl_3>90]
pdf(file='plots_elife/Figure3.pdf',width=6, height=4)
#par(mfrow=c(1,2))
layout(mat=matrix(c(1,2),1,2),w=c(1,1.3))
par(mar=c(5,4,4,2))
plot(1,1,type="n",ylim=c(0,1),xlim=c(2.2,2.38),bty="n",las=1,yaxt="n",xlab="",ylab=expression(paste("Projection along ",delta,"p (",Pi,")")),xaxt="n")
axis(2,at=c(0,0.5,1))
mtext(side=3,"A",at=2.1,cex=2)
axis(side=1,at=c(2.25,2.3),labels=c("Null","Observed"),las=2)
temp_95 <- HPDinterval(as.mcmc(vect_rand_piE_NaCl))
temp_80 <- HPDinterval(as.mcmc(vect_rand_piE_NaCl),prob=0.83)
arrows(2.3,temp_95[1,1],2.3,temp_95[1,2],code=3,angle=90,length=0.05)
arrows(2.3,temp_80[1,1],2.3,temp_80[1,2],code=3,angle=90,length=0,col="gray",lwd=2)
points(2.3,mean(vect_rand_piE_NaCl),pch=16)
temp_95 <- HPDinterval(as.mcmc(vect_rand_piE_NaCl_NULL))
arrows(2.25,temp_95[1,1],2.25,temp_95[1,2],code=3,angle=90,length=0.05,col="orange",lwd=2)
#legend(0,2.35,c(expression(Pi[plast]),expression(Pi[plast_Randomized])),pch=c(16,NA),lwd=1,col=c("black","orange"),cex=1)
par(mar=c(5,4,4,2))
plot(1,1,type="n",ylim=c(0,90),xlab="",xlim=c(1,2.5),bty="n",las=1,yaxt="n",ylab=expression(paste("Angle with ",delta,"p (",Theta,")")),xaxt="n")
axis(2,at=c(0,45,90))
mtext(side=3,"B",at=0,cex=2)
temp_95 <- HPDinterval(as.mcmc(vect_rand_thetaE_NaCl))
temp_80 <- HPDinterval(as.mcmc(vect_rand_thetaE_NaCl),prob=0.83)
arrows(1.5,temp_95[1,1],1.5,temp_95[1,2],code=3,angle=90,length=0.05)
arrows(1.5,temp_80[1,1],1.5,temp_80[1,2],code=3,angle=90,length=0,col="gray",lwd=2)
points(1.5,mean(vect_rand_thetaE_NaCl),pch=16)
temp_95 <- HPDinterval(as.mcmc(vect_rand_thetaE_NaCl_2))
temp_80 <- HPDinterval(as.mcmc(vect_rand_thetaE_NaCl_2),prob=0.83)
arrows(1.75,temp_95[1,1],1.75,temp_95[1,2],code=3,angle=90,length=0.05)
arrows(1.75,temp_80[1,1],1.75,temp_80[1,2],code=3,angle=90,length=0,col="gray",lwd=2)
points(1.75,mean(vect_rand_thetaE_NaCl_2),pch=16)
temp_95 <- HPDinterval(as.mcmc(vect_rand_thetaE_NaCl_3))
temp_80 <- HPDinterval(as.mcmc(vect_rand_thetaE_NaCl_3),prob=0.83)
arrows(2,temp_95[1,1],2,temp_95[1,2],code=3,angle=90,length=0.05)
arrows(2,temp_80[1,1],2,temp_80[1,2],code=3,angle=90,length=0,col="gray",lwd=2)
points(2,mean(vect_rand_thetaE_NaCl_3),pch=16)
temp_95 <- HPDinterval(as.mcmc(Angle_rand))
temp_80 <- HPDinterval(as.mcmc(Angle_rand),prob=0.83)
arrows(1.1,temp_95[1,1],1.1,temp_95[1,2],code=3,angle=90,length=0.05)
arrows(1.1,temp_80[1,1],1.1,temp_80[1,2],code=3,angle=90,length=0,col="gray",lwd=2)
axis(side=1,at=c(2,1.75,1.5,1.1),c(expression(g[3]),expression(g[2]),expression(g[max]),"Null"),las=2)
dev.off()