Commit 6ba3bb6b authored by peguerin's avatar peguerin
Browse files

format main figures nature comm

parent a80deb12
......@@ -74,8 +74,8 @@ grid_GD <- function(gridW,dcW.1,dcW) {
plot_GD_world <- function(gridW,worldcoast,biglakes,riverlakes,titleFig) {
plot_world(worldcoast,biglakes,riverlakes)
points(gridW$lon,gridW$lat,bg=as.character(gridW$col),pch=22,cex=0.64,col="black",lwd=0.2)
box()
mtext(titleFig,side=2,line=0.4,at=96,las=2,cex=1.2)
box(lwd=0.3)
mtext(titleFig,side=2,line=0.4,at=96,las=2,cex=1.2,font=2)
}
## convert coordinates into epsg:4326 and aggregate by latband(10degree)
......@@ -111,20 +111,54 @@ latbandGD <- function(dcW.1, dcW,gridW) {
return(dcW.albc)
}
latbandGDrichness <- function(dcW.1, dcW,gridW) {
gridW.1=gridW[which(gridW$IDcell %in% dcW.1$cell),]
dcW=dcW[which(dcW$cell %in% dcW.1$cell),]
centroW <- gCentroid(gridW.1,byid=TRUE)
spW=SpatialPoints(data.frame(as.numeric(centroW@coords[,1]),as.numeric(centroW@coords[,2])))
proj4string(spW)=proj4string(gridW)
spW.g=spTransform(spW, CRS("+init=epsg:4326"))
lat_W=matrix(spW.g@coords[,2])
latband_W=apply(lat_W, 1, FUN=function(x) coords_to_latband(x))
Test=data.frame(latband=latband_W,GD=dcW$GD_mean)
Test[order(Test[,1], decreasing=T),]
}
get_pt_band <- function ()
gridW.1=gridW[which(gridW$IDcell %in% dcW.1$cell),]
dcW=dcW[which(dcW$cell %in% dcW.1$cell),]
centroW <- gCentroid(gridW.1,byid=TRUE)
spW=SpatialPoints(data.frame(as.numeric(centroW@coords[,1]),as.numeric(centroW@coords[,2])))
proj4string(spW)=proj4string(gridW)
spW.g=spTransform(spW, CRS("+init=epsg:4326"))
lat_W=matrix(spW.g@coords[,2])
latband_W=apply(lat_W, 1, FUN=function(x) coords_to_latband(x))
dcW.latband=data.frame(latband=latband_W,GD=dcW$GD_mean)
## barplot of the distribution of genetic diversity mean by cell with IC
plot_latbandGD <-function(dcW.lb, firstRow, secondRow,xlabFig,titleFig,labels.lb){
barplot.W <- barplot(dcW.lb[firstRow:secondRow,2],horiz=T,col=as.character(dcW.lb[firstRow:secondRow,5]),xlim=c(0,0.05),yaxs="i",cex.axis=1.2,tcl=-0.25,mgp=c(3,0.4,0),plot=TRUE,offset=0.000)
axis(side=2,at=barplot.W[firstRow:secondRow,1],labels=labels.lb,las=2,cex.axis=1.2,tcl=-0.25)
arrows(dcW.lb[firstRow:secondRow,2]-dcW.lb[firstRow:secondRow,3],
plot_latbandGD <-function(dcW.lb, firstRow, secondRow,xlabFig,titleFig,labels.lb,Val,xlim=c(0,0.3)){
barplot.W <- barplot(dcW.lb[firstRow:secondRow,2],horiz=T,col=as.character(dcW.lb[firstRow:secondRow,5]),xlim=xlim,yaxs="i",cex.axis=1.2,tcl=-0.25,mgp=c(3,0.4,0),plot=TRUE,offset=0.000,lwd=0.1)
axis(side=2,at=barplot.W[firstRow:secondRow,1],labels=labels.lb,las=2,cex.axis=1.2,tcl=-0.25,lwd=0.1)
t<- cbind(barplot.W[firstRow:secondRow,1],seq(-90,90,10)[-10])
b<- sapply(1:nrow(Val),function(i){t[which(t[,2]==Val[i,1]),1]})
points(x=Val[,2],as.numeric(b),cex=0.15,pch=19,col="black",lwd=0.1)
arrows(dcW.lb[firstRow:secondRow,2],
barplot.W,
dcW.lb[firstRow:secondRow,2]+dcW.lb[firstRow:secondRow,3],
barplot.W,
angle=90, code=3,length=0.05)
mtext(titleFig,side=2,line=0.4,at=21.5,las=2,cex=1.2)
mtext(xlabFig,side=1,line=1.6,at=0.02,cex=0.85)
box()
angle=90, code=3,length=0.05,lwd=0.1)
mtext(titleFig,side=2,line=0.4,at=21.5,las=2,cex=1.2,font=2)
mtext(xlabFig,side=1,line=1.6,at=0.075,cex=0.85)
box(lwd=0.3)
}
#plot_latbandGD(dcM.lb, 1,17,labels=labels.lb,"", "b",Val=Val)
##########################################################################
## load data
......@@ -166,43 +200,63 @@ gridM=grid_GD(grid_Wgs84,dcM.1,dcM)
gridF=grid_GD(grid_Wgs84,dcF.1,dcF)
## dataframe of genetic diversity by latband
dcM.lb=latbandGD(dcM.1,dcM,grid)
dcF.lb=latbandGD(dcF.1,dcF,grid)
dcM.lb=latbandGD(dcW.1 = dcM.1,dcW=dcM,gridW=grid)
Val <- latbandGDrichness(dcW.1 = dcM.1,dcW=dcM,gridW=grid)
dcF.lb=latbandGD(dcF.1,dcF,grid)
Val2 <- latbandGDrichness(dcW.1 =dcF.1,dcW=dcF,gridW=grid)
##########################################################################
##write tables
write.table(dcM.lb,"14-tables/figure_1b_underlying_data.tsv",col.names=T,row.names=F,sep="\t")
write.table(dcF.lb,"14-tables/figure_1d_underlying_data.tsv",col.names=T,row.names=F,sep="\t")
##########################################################################
## write pdf files
labels.lb=c("-80","","-60","","-40","","-20","","0","","20","","40","","60","","80")
tiff(filename="10-figures/figure1.tiff",width=20,height=16,units="cm",res=640,compression="lzw")
#cairo_ps(filename="10-figures/figure1.eps",width=7.87,height=6.3,fallback_resolution=600,onefile=FALSE)
#tiff(filename="10-figures/figure1.tiff",width=20,height=16,units="cm",res=640,compression="lzw")
cairo_ps(filename="10-figures/figure1bis.eps",width=7.87,height=6.3,fallback_resolution=600,onefile=FALSE)
layout(mat=rbind(c(1,1,2),c(1,1,2),c(1,1,2),c(3,3,4),c(3,3,4),c(3,3,4)))
par(mar=c(3,3,3,1))
### marin
plot_GD_world(gridM,worldcoast,biglakes,riverlakes,"(a)")
plot_latbandGD(dcM.lb, 1,17,labels=labels.lb,"", "(b)")
plot_GD_world(gridM,worldcoast,biglakes,riverlakes,"a")
plot_latbandGD(dcM.lb, 1,17,labels=labels.lb,"", "b",Val=Val,xlim=c(0,0.05))
rasterImage(Mraster_imag,xleft=0.038,ybottom=17.5,xright=0.049,ytop=20)
### freshwater
plot_GD_world(gridF,worldcoast,biglakes,riverlakes,"(c)")
plot_GD_world(gridF,worldcoast,biglakes,riverlakes,"c")
legend("bottomleft",legend=levels(gridF$GD),col=levels(gridF$col),
pch=15,cex=0.9,title="Genetic diversity",pt.cex=1)
plot_latbandGD(dcF.lb, 1,17,labels=labels.lb,"Genetic diversity", "(d)")
pch=15,cex=0.9,title="Genetic diversity",pt.cex=1,box.lwd=0.3 )
plot_latbandGD(dcF.lb, 1,17,labels=labels.lb,"Genetic diversity", "d",Val=Val2,xlim=c(0,0.05))
legend("bottomright",legend=levels(dcF.lb$countCells),col=levels(dcF.lb$col),
pch=15,cex=0.9,title="# cells",pt.cex=1)
pch=15,cex=0.8,title="# cells",pt.cex=1,box.lwd=0.3 )
rasterImage(Fraster_imag,xleft=0.031,ybottom=17.5,xright=0.049,ytop=20)
dev.off()
################################################################################
cairo_ps(filename="10-figures/figure1.eps",width=7.87,height=6.3,fallback_resolution=600,onefile=FALSE)
layout(mat=rbind(c(1,1,2),c(1,1,2),c(1,1,2),c(3,3,4),c(3,3,4),c(3,3,4)))
par(mar=c(3,3,3,1))
### marin
plot_GD_world(gridM,worldcoast,biglakes,riverlakes,"a")
plot_latbandGD(dcM.lb, 1,17,labels=labels.lb,"", "b",Val=Val,xlim=c(0,0.3))
rasterImage(Mraster_imag,xleft=0.24,ybottom=17.5,xright=0.295,ytop=20)
### freshwater
par(mar=c(4,3,2,1))
plot_GD_world(gridF,worldcoast,biglakes,riverlakes,"c")
legend("bottomleft",legend=levels(gridF$GD),col=levels(gridF$col),
pch=15,cex=0.9,title="Genetic diversity",pt.cex=1,box.lwd=0.3)
plot_latbandGD(dcF.lb, 1,17,labels=labels.lb,"Genetic diversity", "d",Val=Val2,xlim=c(0,0.169))
legend("bottomright",legend=levels(dcF.lb$countCells),col=levels(dcF.lb$col),
pch=15,cex=0.9,title="# cells",pt.cex=1,box.lwd=0.3)
rasterImage(Fraster_imag,xleft=0.1,ybottom=17.5,xright=0.158,ytop=20)
dev.off()
......@@ -92,14 +92,16 @@ plot_color_lm<- function(mod,df_x,df_y,x_lab,y_lab,titleFig) {
dr_y=b+a1*dr_x
dregr=data.frame(dr_x,dr_y)
plot(df_x,df_y,col=color_value(df_x,df_y),pch=19,cex=0.65,lwd=2,axes=FALSE,xlab="",ylab="",xlim=c(0,0.065),ylim=c(0,3600))
lines(dr_x,dr_y,lwd=2)
lines(dr_x,dr_y,lwd=0.9)
axis(side=1,line=0,cex.axis=0.85,lwd=0.35,tcl=-0.25,bg="white",labels=as.character(seq(0,0.6,by=0.02)), at=seq(0,0.6,by=0.02),mgp=c(3,0.5,0))
axis(side=2,line=0,cex.axis=0.85,lwd=0.35,tcl=-0.25,bg="white",labels=as.character(seq(0,3600,by=400)), at=seq(0,3600,by=400),mgp=c(3,0.5,0))
box()
box(lwd=0.1)
mtext(x_lab,side=1,line=1.5,cex=0.7)
mtext(y_lab,side=2,line=1.5,cex=0.7)
mtext(paste("r =",corrValue),side=3,line = -1.5,at=0.05,cex=0.7)
mtext(titleFig,line=0.4,at=0)
mtext("r ",side=3,line=-1.5,at=0.045,cex=0.7,font=3)
mtext(paste("=",corrValue),side=3,line = -1.5,at=0.055,cex=0.7)
mtext(titleFig,line=0.6,at=0,font=2)
}
## add GD and species diversity information and color to the grid
......@@ -113,9 +115,9 @@ grid_GD_richness <- function(gridW,dcW) {
## world map with colored squares according to their (richness,genetic div) values
plot_GD_world <- function(gridW,worldcoast,biglakes,riverlakes,titleFig) {
plot_world(worldcoast,biglakes,riverlakes)
points(gridW$lon,gridW$lat,bg=as.character(gridW$col),pch=22,cex=0.65,col="black",lwd=0.2)
box()
mtext(titleFig,side=2,line=0.4,at=90,las=2)
points(gridW$lon,gridW$lat,bg=as.character(gridW$col),pch=22,cex=0.65,col="black",lwd=0.01)
box(lwd=0.2)
mtext(titleFig,side=2,line=0.4,at=93,las=2,font=2)
}
## plot a legend color 2D
......@@ -130,7 +132,7 @@ color_map_legend <- function (colours, labels = FALSE, borders = NULL, cex_label
size <- max(dim(colours))
plot(c(0, size), c(0, -size), type = "n",axes=F)
rect(col(colours) - 1, -row(colours) + 1, col(colours), -row(colours),
col = colours, border = borders)
col = colours, border = borders,lwd=0.05)
if (labels) {
text(col(colours) - 0.5, -row(colours) + 0.5, colours,
cex = cex_label)
......@@ -225,19 +227,19 @@ write.table(data.frame(gd=Fr$GD_mean,richness=Fr$richness)
##########################################################################
## write pdf files
#tiff(filename="10-figures/figure2_colorblind.tiff",width=20,height=16,units="cm",res=640,compression="lzw")
cairo_ps(filename="10-figures/figure2_colorblind.eps",width=7.87,height=6.3,fallback_resolution=600,onefile=FALSE)
cairo_ps(filename="figure2_colorblind.eps",width=7.87,height=6.3,fallback_resolution=600,onefile=FALSE)
layout(mat=rbind(c(1,2,2),c(1,2,2),c(1,2,2),c(3,4,4),c(3,4,4),c(3,4,4)))
### marin
par(mar=c(2.5,3.5,2,2))
plot_color_lm(modMr.final,Mr$GD_mean,Mr$richness,"Genetic diversity","Species diversity","(a)")
plot_color_lm(modMr.final,Mr$GD_mean,Mr$richness,"Genetic diversity","Species diversity","a")
par(mar=c(2.5,2,2,2))
plot_GD_world(grid.M,worldcoast,biglakes,riverlakes,"(b)")
plot_GD_world(grid.M,worldcoast,biglakes,riverlakes,"b")
### freshwater
par(mar=c(2.5,3.5,2,2))
plot_color_lm(modFr.final,Fr$GD_mean,Fr$richness,"Genetic diversity","Species diversity","(c)")
plot_color_lm(modFr.final,Fr$GD_mean,Fr$richness,"Genetic diversity","Species diversity","c")
par(mar=c(2.5,2,2,2))
plot_GD_world(grid.F,worldcoast,biglakes,riverlakes,"(d)")
plot_GD_world(grid.F,worldcoast,biglakes,riverlakes,"d")
## legend
#par( fig=c(0.4,0.45,0.15,0.2),new = T)
par( fig=c(0.39,0.46,0.11,0.18),new = T)
......
......@@ -176,17 +176,18 @@ res.mod.ols # Of course non significant
##########################################################################
## confidence intervals
nom_variables=c("Mean regional\ntemperature","Average slope","Regions","Autocor","Number of species", "Number of sequences")
nom_variables=c("Mean regional\ntemperature","Average slope","Regions","Autocor",
"Number of \nspecies", "Number of \nsequences")
names(modclimF.autof$coefficients)=c("(Intercept)",nom_variables)
p_ic_F=plot_model(modclimF.autof,colors = c("blue","red"), type = "est",sort.est=T,vline.color = "grey", line.size=1.2,show.values = TRUE, value.offset = .25, value.size=8,digits=3)
p_ic_F=p_ic_F+theme(plot.title = element_text(size = 24),
p_ic_F=plot_model(modclimF.autof,colors ="black", type="est",sort.est=T,vline.color="grey", line.size=1.2,show.values = TRUE, value.offset = .25, value.size=8,digits=3)
p_ic_F=p_ic_F+theme(plot.title = element_text(size = 24,face="bold"),
axis.title=element_text(size=20),
axis.text=element_text(size=20),
plot.margin = unit(c(0.5,2,1,1), "cm"),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.border=element_rect(fill=NA,size=1),
panel.background = element_blank())+ggtitle("(c)")
panel.background = element_blank())+ggtitle("c")
## underlying data
sink("13-analysis/freshwater_summary_lm.txt")
print(summary(modclimF.autof))
......@@ -194,7 +195,7 @@ sink()
##########################################################################
## hier.part
nom_variables=c("Number of species", "Number of sequences","Mean regional\ntemperature","Regions","Average slope","Autocor")
nom_variables=c("Number of \n species", "Number of\n sequences","Mean regional\ntemperature","Regions","Average slope","Autocor")
dcF.r3.auto=cbind(env2,autocor)
hp_F=hier.part(dcF.r3.auto$GD_mean,dcF.r3.auto[,c(2,3,4,6,10,11)])
hp_df_F=data.frame(nom_var=nom_variables,hier_part=hp_F$I.perc$I)
......@@ -206,7 +207,7 @@ p_hp_F=ggplot(data=hp_df_F, aes(x=nom_var, y=hier_part)) +
ylab("% Relative variance")+
scale_y_discrete(limits=seq(0,100,20),expand = expand_scale(mult = c(0, .1)))+
coord_cartesian(ylim = c(0, 100))+
theme(plot.title = element_text(size = 24),
theme(plot.title = element_text(size = 24,face="bold"),
axis.title.y=element_text(size=20),
axis.text.x=element_text(size=20,hjust=0.75,vjust=0.8,angle=45),
axis.text.y = element_text(size=20),
......@@ -215,7 +216,7 @@ p_hp_F=ggplot(data=hp_df_F, aes(x=nom_var, y=hier_part)) +
panel.grid.major = element_blank(),
panel.border=element_rect(fill=NA,size=1),
panel.background = element_blank())+
ggtitle("(d)")
ggtitle("d")
## write underlying data tables
write.table(hp_df_F,"14-tables/figure_3d_underlying_data.tsv",col.names=T,row.names=F,sep="\t")
......@@ -352,17 +353,17 @@ res.mod.ols.f
##########################################################################
## confidence intervals
nom_variables=c("Sea surface\ntemperature","Regions","Autocor","Number of species","Number of sequences")
nom_variables=c("Sea surface\ntemperature","Regions","Autocor","Number of\n species","Number of\n sequences")
names(modclimM.autof$coefficients)=c("(Intercept)",nom_variables)
p_ic_M=plot_model(modclimM.autof,colors = c("blue","red"), type = "est",sort.est=T,vline.color = "grey", line.size=1,show.values = TRUE, value.offset = .25, value.size=8,digits=3)
p_ic_M=p_ic_M+theme(plot.title = element_text(size = 24),
p_ic_M=plot_model(modclimM.autof,colors="black", type = "est",sort.est=T,vline.color = "grey", line.size=1,show.values = TRUE, value.offset = .25, value.size=8,digits=3)
p_ic_M=p_ic_M+theme(plot.title = element_text(size = 24,face="bold"),
axis.title=element_text(size=20),
axis.text=element_text(size=20),
plot.margin = unit(c(0.5,2,1,1), "cm"),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.border=element_rect(fill=NA,size=1),
panel.background = element_blank())+ggtitle("(a)")
panel.background = element_blank())+ggtitle("a")
## underlying data
sink("13-analysis/marine_summary_lm.txt")
print(summary(modclimM.autof))
......@@ -370,7 +371,7 @@ sink()
##########################################################################
## hier.part
nom_variables=c("Sea surface\ntemperature","Regions","Autocor","Number of species","Number of sequences")
nom_variables=c("Sea surface\ntemperature","Regions","Autocor","Number of\n species","Number of\n sequences")
dcM.auto=cbind(dcM.r4,autocor)
hp_M=hier.part(dcM.auto$GD_mean,dcM.auto[,c(10,11,13,5,6)])
hp_df_M=data.frame(nom_var=nom_variables,hier_part=hp_M$I.perc$I)
......@@ -382,7 +383,7 @@ p_hp_M=ggplot(data=hp_df_M, aes(x=nom_var, y=hier_part)) +
ylab("% Relative variance")+
scale_y_discrete(limits=seq(0,100,20),expand = expand_scale(mult = c(0, .1)))+
coord_cartesian(ylim = c(0, 100))+
theme(plot.title = element_text(size = 24),
theme(plot.title = element_text(size = 24,face="bold"),
axis.title.y=element_text(size=20),
axis.text.x=element_text(size=20,hjust=0.75,vjust=0.8,angle=45),
axis.text.y = element_text(size=20),
......@@ -391,7 +392,7 @@ p_hp_M=ggplot(data=hp_df_M, aes(x=nom_var, y=hier_part)) +
panel.grid.major = element_blank(),
panel.border=element_rect(fill=NA,size=1),
panel.background = element_blank())+
ggtitle("(b)")
ggtitle("b")
## write underlying data tables
write.table(hp_df_M,"14-tables/figure_3b_underlying_data.tsv",col.names=T,row.names=F,sep="\t")
......@@ -402,6 +403,14 @@ write.table(hp_df_M,"14-tables/figure_3b_underlying_data.tsv",col.names=T,row.na
##########################################################################
## write pdf files
pdf("10-figures/figure3.pdf",width=16,height=14,paper='special')
pdf("figure3.pdf",width=16,height=14,paper='special')
grid.arrange(p_ic_M,p_hp_M,p_ic_F,p_hp_F, nrow=2, widths=c(3,4))
dev.off()
cairo_ps(filename="figure3.eps",width=16,height=14,fallback_resolution=600,onefile=FALSE)
#pdf("figure3.pdf",width=16,height=14,paper='special')
grid.arrange(p_ic_M,p_hp_M,p_ic_F,p_hp_F, nrow=2, widths=c(3,4))
dev.off()
......@@ -101,13 +101,13 @@ ggplot_lm<- function(mod,df_x,df_y,x_lab,y_lab,titre) {
}
## from shapefiles, it "plots" a world map
plot_world<- function(worldcoast,biglakes,riverlakes) {
plot_world<- function(worldcoast,biglakes,riverlakes,cex.axis=1.2) {
plot(rnorm(100),xlim=c(-167, 167),ylim=c(-84,75),xlab="",ylab="",type="n",axes=F)
plot(worldcoast,add=TRUE,col="grey",border="black",lwd=0.2)
plot(biglakes,add=TRUE,col="white",lwd=0.2)
plot(riverlakes,add=TRUE,col="white",lwd=0.2)
axis(side=1,line=0,cex.axis=1.2,lwd=0.35,tcl=-0.25,bg="white",labels=c("-150","-100","-50","0","50","100","150"), at=c(-150,-100,-50,0,50,100,150),mgp=c(3,0.5,0))
axis(side=2,line=0,las=2,cex.axis=1.2,lwd=0.35,tcl=-0.25,bg="white",labels=c("-80","-60","-40","-20","0","20","40","60","80"), at=c(-80,-60,-40,-20,0,20,40,60,80),mgp=c(3, 0.5, 0))
axis(side=1,line=0,cex.axis=cex.axis,lwd=0.2,tcl=-0.25,bg="white",labels=c("-150","-100","-50","0","50","100","150"), at=c(-150,-100,-50,0,50,100,150),mgp=c(3,0.5,0))
axis(side=2,line=0,las=2,cex.axis=cex.axis,lwd=0.2,tcl=-0.25,bg="white",labels=c("-80","-60","-40","-20","0","20","40","60","80"), at=c(-80,-60,-40,-20,0,20,40,60,80),mgp=c(3, 0.5, 0))
}
## from a latitude coordinate, return 10-based latband value
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment