Commit a80deb12 authored by peguerin's avatar peguerin
Browse files

format nature comm

parent 52f1a5f9
......@@ -121,20 +121,23 @@ write.table(pgi.corF,"14-tables/figure_S1b_underlying_data.tsv",col.names=T,row.
##########################################################################
## write pdf files
pdf("10-figures/figureS1.pdf",width=14,height=10,paper='special')
pdf("10-figures/figureS1.pdf",width=12,height=6,paper='special')
#cairo_ps(filename="FigureS1.eps",width=12,height=6,fallback_resolution=600,onefile=FALSE)
layout(matrix(c(1,1,2,2),nrow=1))
par(mar=c(4,4,4,2)) # c(bottom, left, top, right)
plot(pgi.corM,xlim=c(200,10000),ylim=c(-0.3,0.3),main="", xlab="", ylab="",cex.axis=1.6,xaxt="n")
par(mar=c(4.5,4.5,4,2)) # c(bottom, left, top, right)
plot(pgi.corM,xlim=c(200,10000),ylim=c(-0.2,0.3),main="", xlab="", ylab="",cex.axis=1.6,xaxt="n")
axis(side=1, at=c(200,2000,4000,6000,8000,10000), labels=c(200,2000,4000,6000,8000,10000),mgp=c(3,1,0),cex.axis=1.6)
mtext("(a)",side=3,line=1,at=0,cex=1.6)
mtext("Moran I",side=2,line=2,at=0.05,cex=1.4)
mtext("Distance classes (km)",side=1,line=3,cex=1.4)
mtext("a",side=3,line=1,at=-400,cex=1.6,font=2)
mtext("Moran ",side=2,line=2.4,at=0.05,cex=1.4)
mtext("I",side=2,line=2.4,at=0.09,cex=1.4,font=3)
plot(pgi.corF,xlim=c(0,10000),ylim=c(-0.3,0.3),main="", xlab="", ylab="",cex.axis=1.6,xaxt="n")
mtext("Distance classes (km)",side=1,line=3.2,cex=1.4)
plot(pgi.corF,xlim=c(200,10000),ylim=c(-0.2,0.3),main="", xlab="", ylab="",cex.axis=1.6,xaxt="n")
axis(side=1, at=c(200,2000,4000,6000,8000,10000), labels=c(200,2000,4000,6000,8000,10000),mgp=c(3,1,0),cex.axis=1.6)
mtext("(b)",side=3,line=1,at=0,cex=1.6)
mtext("Distance classes (km)",side=1,line=3,cex=1.4)
mtext("b",side=3,line=1,at=-400,cex=1.6,font=2)
mtext("Distance classes (km)",side=1,line=3.2,cex=1.4)
dev.off()
\ No newline at end of file
......@@ -50,7 +50,6 @@ grid_percentiles_upper <- function(dcV, dcV.1, gridV) {
return(gridV.up)
}
## creates a grid with lower percentiles levels data for genetic diversity
grid_percentiles_lower <- function(dcV, dcV.1, gridV) {
lower10.1=quantile(dcV.1$GD_mean, 0.1)
......@@ -59,15 +58,13 @@ grid_percentiles_lower <- function(dcV, dcV.1, gridV) {
return(gridV.low)
}
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)
plot_GD_world <- function(gridW,worldcoast,biglakes,riverlakes,titleFig,cex.pt) {
plot_world(worldcoast,biglakes,riverlakes,cex.axis=0.8)
points(gridW$lon,gridW$lat,bg=as.character(gridW$col),pch=22,cex=cex.pt,col="black",lwd=0.01)
box(lwd=0.2)
mtext(titleFig,side=2,line=0.4,at=98,las=2,font=2)
}
##########################################################################
## load data
dc = load_data("09-all_descripteurs/total_data_genetic_diversity_with_all_descripteurs.tsv")
......@@ -109,13 +106,16 @@ gridF.low=grid_percentiles_lower(dcF,dcF.1,grid_Wgs84)
##########################################################################
## write pdf files
tiff(filename="10-figures/figureS2.tiff",width=22,height=16,units="cm",res=640,compression="lzw")
#tiff(filename="10-figures/figureS2.tiff",width=22,height=16,units="cm",res=640,compression="lzw")
pdf("10-figures/figureS2.pdf",width=8,height=5,paper='special')
#cairo_ps(filename="FigureS2.eps",width=8,height=5,fallback_resolution=300,onefile=FALSE)
layout(mat=rbind(c(1,2),c(1,2),c(3,4),c(3,4)))
par(mar=c(2,2,2,2))
plot_GD_world(gridM.up,worldcoast,biglakes,riverlakes,"(a)")
plot_GD_world(gridM.low,worldcoast,biglakes,riverlakes,"(b)")
plot_GD_world(gridF.up,worldcoast,biglakes,riverlakes,"(c)")
plot_GD_world(gridF.low,worldcoast,biglakes,riverlakes,"(d)")
par(mar=c(2,3,3,2))
plot_GD_world(gridM.up,worldcoast,biglakes,riverlakes,"a",cex.pt=0.55)
plot_GD_world(gridM.low,worldcoast,biglakes,riverlakes,"b",cex.pt=0.55)
par(mar=c(3,3,2,2))
plot_GD_world(gridF.up,worldcoast,biglakes,riverlakes,"c",cex.pt=0.55)
plot_GD_world(gridF.low,worldcoast,biglakes,riverlakes,"d",cex.pt=0.55)
dev.off()
......@@ -48,29 +48,47 @@ latbandrichness <- function(dcW,gridW) {
return(dcW.alb)
}
plot_latband_richness <-function(dcW.lb, firstRow, secondRow,xlabFig,titleFig){
barplot.W <- barplot(dcW.lb[firstRow:secondRow,2],horiz=T,col="grey",xlim=c(0,1600),cex.axis=1.2,tcl=-0.25,mgp=c(3,0.4,0),plot=TRUE,offset=0.000,xaxt="n")
axis(side=1,line=0,cex.axis=1.6,lwd=0.35,tcl=-0.25,bg="white",labels=as.character(seq(0,1600,by=400)), at=seq(0,1600,by=400),mgp=c(3,0.5,0))
plot_latband_richness <-function(Data=dcM.lb, firstRow, secondRow,xlab="",ylab="",titleFig,Val,xlim=c(0,3000)){
barplot.W <- barplot(Data[firstRow:secondRow,2],horiz=T,col="grey",xlim=xlim,
cex.axis=1.2,tcl=-0.25,mgp=c(3,0.4,0),plot=TRUE,offset=0.000,xaxt="n",lwd=0.4)
axis(side=1,line=0,cex.axis=1.6,lwd=0.35,tcl=-0.25,bg="white",labels=as.character(seq(0,xlim[2],by=400)),
at=seq(0,xlim[2],by=400),mgp=c(3,0.5,0))
axis(side=2,at=barplot.W[firstRow:secondRow,1],labels=seq(-90,90,10)[-10],las=2,cex.axis=1.2,tcl=-0.25)
arrows(dcW.lb$richness_mean-dcW.lb$richness_sd,
barplot.W,
dcW.lb$richness_mean+dcW.lb$richness_sd,
barplot.W,
angle=90, code=3,length=0.05)
mtext(titleFig,side=2,line=1,at=23.5,las=2,cex=1.4)
mtext(xlabFig,side=1,line=2,at=800,cex=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.2,pch=19,col="black")
arrows(Data$richness_mean,barplot.W,Data$richness_mean+Data$richness_sd,
barplot.W,angle=90,code=3,length=0.05)
mtext(titleFig,side=2,line=1,at=23.5,las=2,cex=1.4,font=2)
mtext(ylab,side=2,line=3,at=10,cex=1.2)
mtext(xlab,side=1,line=2,at=1500,cex=1)
box()
}
boxplot_richness_env <- function(gdW,titleFig,ylabFig) {
boxplot(richness~environment,data=gdW, main="", xlab="", ylab="",col="grey",ylim=c(0,max(gdW$richness)),cex.axis=1.6,cex.lab=1,yaxt = "n")
boxplot_richness_env <- function(Data,titleFig,ylabFig) {
boxplot(richness~environment,data=Data, main="",xlab="", ylab="",col="grey",
ylim=c(0,max(Data$richness)),cex.axis=1.6,cex.lab=1,yaxt = "n")
axis(side=2,line=0,cex.axis=1.6,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))
mtext(titleFig,side=2,line=1,at=3850,las=2,cex=1.4)
mtext(titleFig,side=2,line=1,at=3850,las=2,cex=1.4,font=2)
mtext(ylabFig,side=2,line=3,at=1600,cex=1.2)
}
latvalbandrichness <- function(dcW,gridW) {
gridW.1=gridW[which(gridW$IDcell %in% dcW$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,richness=dcW$richness)
Test[order(Test[,1], decreasing=T),]
}
##########################################################################
## load data
dc = load_data("09-all_descripteurs/total_data_genetic_diversity_with_all_descripteurs.tsv")
......@@ -104,44 +122,44 @@ Tr$environment <- factor(Tr$environment,levels = c('marine','freshwater'),ordere
## latclass
dat=read.table("14-tables/latclass.tsv",header=T)
##########################################################################
## load worldcoast shapefile
grid <- readOGR("01-infos/grid_equalarea200km","gridFish.b_260418")
##########################################################################
## model to compare genetic diversity bewteen marine and FW species
mod1<-lm(scale(log(GD_mean))~is_marine+latitude,data=dat)
## model to compare species diversity bewteen marine and FW species
mod2<-lm(scale(log(richness))~is_marine+latitude,data=dat)
summary(mod1)
summary(mod2)
summary(mod1); summary(mod2)
##########################################################################
## dataframe of genetic diversity by latband
dcM.lb=latbandrichness(Mr,grid)
dcF.lb=latbandrichness(Fr,grid)
dcM.lb=latbandrichness(dcW=Mr,grid=gridW)
Val <- latvalbandrichness(dcW=Mr,grid=gridW)
dcF.lb=latbandrichness(Fr,grid)
Val2 <- latvalbandrichness(dcW=Fr,grid=grid)
##########################################################################
##write tables
write.table(Tr,"14-tables/figure_S3a_underlying_data.tsv",col.names=T,row.names=F,sep="\t")
write.table(dcM.lb,"14-tables/figure_S3b_underlying_data.tsv",col.names=T,row.names=F,sep="\t")
write.table(dcF.lb,"14-tables/figure_S3c_underlying_data.tsv",col.names=T,row.names=F,sep="\t")
##########################################################################
## write pdf files
#cairo_ps(filename="10-figures/FigureS3.eps",width=17,height=8,fallback_resolution=300,onefile=FALSE)
pdf("10-figures/figureS3.pdf",width=17,height=8,paper='special')
layout(matrix(c(1,2,3),nrow=1))
par(mar=c(4,5,4,4))
par(mar=c(4,5,4,3.5))
boxplot_richness_env(Data=Tr,"a","Species diversity")
boxplot_richness_env(Tr,"(a)","Species diversity")
plot_latband_richness(dcM.lb, 1,18,"Species diversity", "(b)")
plot_latband_richness(dcF.lb, 1,18,"Species diversity", "(c)")
par(mar=c(4,2,4,3.5))
plot_latband_richness(Data=dcM.lb,firstRow= 1,secondRow=18,xlab="Species diversity",
ylab="Latitude",titleFig="b",Val=Val)
par(mar=c(4,2,4,2))
plot_latband_richness(Data=dcF.lb, 1,18,xlab="Species diversity",
ylab="Latitude",titleFig="c",Val=Val2)
dev.off()
......@@ -15,17 +15,17 @@
##
##########################################################################
## Libraries
lib_vect <-c("MASS","hier.part","countrycode","lme4","sjPlot","raster","plotrix","rgeos","rgdal","sp","maptools","shape","parallel","png","plyr")
lib_vect <-c("MASS","hier.part","countrycode","lme4","sjPlot","raster","plotrix",
"rgeos","rgdal","sp","maptools","shape","parallel","png","plyr")
sapply(lib_vect,library,character.only=TRUE)
##########################################################################
## functions
source("00-scripts/step5/figures/functions.R")
boxplot_gd_region <- function(gdW,titleFig,ylabFig) {
boxplot(gd~regions,data=gdW, main="", xlab="", ylab="",col="grey",ylim=c(0,0.1),cex.axis=1.4,cex.lab=1.4)
mtext(titleFig,side=2,line=1,at=0.11,las=2,cex=1.6)
mtext(titleFig,side=2,line=1,at=0.11,las=2,cex=1.6,font=2)
mtext(ylabFig,side=2,line=2.6,at=0.05, cex=1.4)
}
......@@ -35,7 +35,6 @@ boxplot_gd_region <- function(gdW,titleFig,ylabFig) {
dc = load_data("09-all_descripteurs/total_data_genetic_diversity_with_all_descripteurs.tsv")
dcT=dc[[1]];dcM=dc[[2]];dcF=dc[[3]];dcM.n1=dc[[4]];dcF.n1=dc[[5]];dcM.1=dc[[6]];dcF.1=dc[[7]]
##########################################################################
## MARINE
### add region indopacific/atlantic information for MARINE
......@@ -72,7 +71,6 @@ dcM.r2=data.frame(
GD_mean=dcM.n1[which(dcM.n1$cell %in% dcM.r1$cell),]$GD_mean)
dcM.r3<-na.omit(dcM.r2)
dM_GD.r=data.frame(gd=dcM.r3$GD_mean,regions=dcM.r3$regions)
......@@ -113,8 +111,8 @@ write.table(dF_GD.r,"14-tables/figure_S4b_underlying_data.tsv",col.names=T,row.n
pdf("10-figures/figureS4.pdf",width=14,height=8,paper='special')
layout(matrix(c(1,1,2,2,2),nrow=1))
par(mar=c(4,4,4,2))
boxplot_gd_region(dM_GD.r,"(a)","Genetic diversity")
boxplot_gd_region(dF_GD.r,"(b)","")
par(mar=c(4,4.5,4,2))
boxplot_gd_region(dM_GD.r,"a","Genetic diversity")
boxplot_gd_region(dF_GD.r,"b","")
dev.off()
......@@ -103,8 +103,8 @@ pdf("10-figures/figureS5.pdf",width=12,height=12,paper='special')
layout(matrix(c(1,2),nrow=2))
par(mar=c(6,4,4,4)) # c(bottom, left, top, right)
barplot(propOrder$proportion,las=2,names.arg=propOrder$order,ylab="Proportion species per order",cex.names = 0.6 )
title("(a)", adj = 0)
title("a", adj = 0)
barplot(NSeqO$numberOfElements, las=2,names.arg=NSeqO$category,ylab="Number sequences per order",cex.names = 0.6 )
title("(b)", adj = 0)
title("b", adj = 0)
dev.off()
......@@ -102,10 +102,10 @@ grid_numberofseq <- function(gridW,dcW.1,dcW) {
## and from a river and lakes shape polygon object
## create a plot map of number of sequences
plot_nseq_world <- function(gridW,colorColon,worldcoast,biglakes,riverlakes,titleFig) {
plot_world(worldcoast,biglakes,riverlakes)
points(gridW$lon,gridW$lat,bg=as.character(gridW[,colorColon]),pch=22,cex=0.65,col="black",lwd=0.2)
plot_world(worldcoast,biglakes,riverlakes,cex.axis=1.2)
points(gridW$lon,gridW$lat,bg=as.character(gridW[,colorColon]),pch=22,cex=0.85,col="black",lwd=0.05)
box()
mtext(titleFig,side=2,line=0.4,at=90,las=2)
mtext(titleFig,side=2,line=0.4,at=90,las=2,font=2)
}
......@@ -138,32 +138,32 @@ riverlakes <- readOGR("01-infos/ne_50m_rivers_lake_centerlines_scale_rank",laye
grid_Wgs84 <- spTransform(grid, proj4string(worldcoast))
biglakes <- readShapePoly("01-infos/big_lakes/GSHHS_h_L2.shp",proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"))
##########################################################################
## add number of sequences to the grid cells
gridM.1.df=grid_numberofseq(grid_Wgs84,dcM.1,dcM)
gridF.1.df=grid_numberofseq(grid_Wgs84,dcF.1,dcF)
##########################################################################
## write pdf files
tiff(filename="10-figures/figureS6.tiff",width=27,height=24,units="cm",res=640,compression="lzw")
#tiff(filename="10-figures/figureS6.tiff",width=27,height=24,units="cm",res=640,compression="lzw")
pdf("10-figures/figureS6.pdf",width=13.5,height=12,paper='special')
layout(mat=rbind(c(1,2),c(3,4),c(5,6)))
par(mar=c(1.5,2,2,1))
par(mar=c(1.5,2.7,2,1))
## number of sequences per cell
plot_nseq_world(gridM.1.df,4,worldcoast,biglakes,riverlakes,"(a)")
plot_nseq_world(gridF.1.df,4,worldcoast,biglakes,riverlakes,"(b)")
plot_nseq_world(gridM.1.df,4,worldcoast,biglakes,riverlakes,"a")
plot_nseq_world(gridF.1.df,4,worldcoast,biglakes,riverlakes,"b")
legend("bottomleft",legend=levels(gridM.1.df$numberOfSequences),col=levels(gridM.1.df$numberOfSequences_color),
pch=15,cex=0.9,title="Number of sequences",pt.cex=1)
pch=15,cex=0.9,title="Number of sequences",pt.cex=1,bg="white")
## number of species per cell
plot_nseq_world(gridM.1.df,6,worldcoast,biglakes,riverlakes,"(c)")
plot_nseq_world(gridF.1.df,6,worldcoast,biglakes,riverlakes,"(d)")
plot_nseq_world(gridM.1.df,6,worldcoast,biglakes,riverlakes,"c")
plot_nseq_world(gridF.1.df,6,worldcoast,biglakes,riverlakes,"d")
legend("bottomleft",legend=levels(gridM.1.df$numberOfSpecies),col=levels(gridM.1.df$numberOfSpecies_color),
pch=15,cex=0.9,title="Number of species",pt.cex=1)
pch=15,cex=0.9,title="Number of species",pt.cex=1,bg="white")
## number of individuals by species per cell
plot_nseq_world(gridM.1.df,8,worldcoast,biglakes,riverlakes,"(e)")
plot_nseq_world(gridF.1.df,8,worldcoast,biglakes,riverlakes,"(f)")
plot_nseq_world(gridM.1.df,8,worldcoast,biglakes,riverlakes,"e")
plot_nseq_world(gridF.1.df,8,worldcoast,biglakes,riverlakes,"f")
legg=legend("bottomleft",legend=levels(gridM.1.df$meanNumberOfIndividualsBySpecies),col=levels(gridM.1.df$meanNumberOfIndividualsBySpecies_color), pch=15,cex=0.9,title="Mean number of\nsequences by species",pt.cex=1,bty='n')
lxleft <- legg$rect[["left"]]
lytop <- legg$rect[["top"]]
......@@ -173,3 +173,6 @@ rect(lxleft, lybottom-.1, lxright, lytop+6,col="white")
legend("bottomleft",legend=levels(gridM.1.df$meanNumberOfIndividualsBySpecies),col=levels(gridM.1.df$meanNumberOfIndividualsBySpecies_color), pch=15,cex=0.9,title="Mean number of\nsequences by species",pt.cex=1,bty='n')
dev.off()
......@@ -60,14 +60,16 @@ grid_numberofseq <- function(gridW,dcW.1,dcW) {
## and from a worldcoast shape polygon object
## and from a river and lakes shape polygon object
## create a plot map of number of sequences
plot_nseq_world <- function(gridW,colorColon,worldcoast,biglakes,riverlakes,titleFig) {
plot_world(worldcoast,biglakes,riverlakes)
points(gridW$lon,gridW$lat,bg=as.character(gridW[,colorColon]),pch=22,cex=0.65,col="black",lwd=0.2)
box()
mtext(titleFig,side=2,line=0.4,at=90,las=2)
plot_world(worldcoast,biglakes,riverlakes,cex.axis = 0.8)
points(gridW$lon,gridW$lat,bg=as.character(gridW[,colorColon]),pch=22,cex=0.65,col="black",lwd=0.05)
box()
mtext(titleFig,side=2,line=0.4,at=92,las=2,font=2,cex=1.4)
}
##########################################################################
## load data
dc <- load_data("09-all_descripteurs/total_data_genetic_diversity_with_all_descripteurs.tsv")
......@@ -121,8 +123,6 @@ row.names(numberCellSpeTaxonCov)=c("marine","freshwater")
## write CSV files
write.csv(numberCellSpeTaxonCov,file="12-taxonomic_coverage/species_taxon_coverage_number_cells.csv",sep=",")
##########################################################################
## load worldcoast shapefile
grid <- readOGR("01-infos/grid_equalarea200km","gridFish.b_260418")
......@@ -141,14 +141,15 @@ gridF.1.df=grid_numberofseq(grid_Wgs84,Fr,dcF)
##########################################################################
## write pdf files
tiff(filename="10-figures/figureS7.tiff",width=26,height=22,units="cm",res=640,compression="lzw")
#tiff(filename="10-figures/figureS7.tiff",width=26,height=22,units="cm",res=640,compression="lzw")
pdf("10-figures/figureS7.pdf",width=8,height=8,paper='special')
layout(mat=rbind(1,2))
par(mar=c(1.5,2,2,1))
par(mar=c(1.5,3,2,1))
## taxonomic coverage species cell
plot_nseq_world(gridM.1.df,4,worldcoast,biglakes,riverlakes,"(a)")
plot_nseq_world(gridF.1.df,4,worldcoast,biglakes,riverlakes,"(b)")
plot_nseq_world(gridM.1.df,4,worldcoast,biglakes,riverlakes,"a")
plot_nseq_world(gridF.1.df,4,worldcoast,biglakes,riverlakes,"b")
legend("bottomleft",legend=levels(gridM.1.df$taxonCoverageSpecies),col=levels(gridM.1.df$taxonCoverageSpecies_color),
pch=15,cex=0.9,title="% taxon cover",pt.cex=1)
pch=15,cex=0.9,title="% taxon cover",pt.cex=1,bg="white")
dev.off()
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