Skip to content

GitLab

Commit 20a5ec16 authored by peguerin's avatar peguerin
Browse files

permutation figure S5 S6

parent ebabf1dc
Hide whitespace changes
Inline Side-by-side
00-scripts/step5/supplementary_figures/figureS5.R
View file @ 20a5ec16
......@@ -9,167 +9,102 @@
## Submited to Nature communications, 2019
##
## functions used to generate figures
## Supplementary Figure 5. Sampling effect.
## (a,b) number of sequences per cell,
## (c,d) number of fish species per cell,
## (e,f) mean number of sequences per fish species per cell for marine species
## (a, c, e) and freshwater species (b, d,f) respectively.
## Supplementary Figure S5.
## Taxonomic coverage of the sequences used by the model
## (a) Number of species per order
## (b) Number of sequences per order
##
##
##########################################################################
## Libraries
lib_vect <-c("raster","plotrix","rgeos","rgdal","sp","maptools","shape","parallel","png","plyr")
sapply(lib_vect,library,character.only=TRUE)
## no library required
##########################################################################
## functions
source("00-scripts/step5/figures/functions.R")
## convert number of sequences value into color value
color_N <- function(gd){
gdc=gd
gdc[which(gd<5)] = "< 5"
gdc[which(gd>=5 & gd<10)] = "5 - 10"
gdc[which(gd>=10 & gd<20)] ="10 - 20"
gdc[which(gd>=20 & gd<30)] ="20 - 30"
gdc[which(gd>=30)] =">30"
return(gdc)
}
color_nspec <- function(gd){
gdc=gd
gdc[which(gd<3)] = "2 - 3"
gdc[which(gd>=3 & gd<5)] = "3 - 5"
gdc[which(gd>=5 & gd<10)] ="5 - 10"
gdc[which(gd>=10)] =">10"
return(gdc)
}
color_ind <- function(gd){
gdc=gd
gdc[which(gd<4)] = "2 - 4"
gdc[which(gd>=4 & gd<6)] = "4 - 6"
gdc[which(gd>=6 & gd<8)] ="6 - 8"
gdc[which(gd>=8 & gd<10)] ="8 - 10"
gdc[which(gd>=10)] =">10"
return(gdc)
}
## from a grid and table of descriptors
## create a grid with number of sequences information for each cell of the grid
grid_numberofseq <- function(gridW,dcW.1,dcW) {
indW=dcW[which(dcW$cell %in% dcW.1$cell),]$nb_indv_mean
specW=dcW[which(dcW$cell %in% dcW.1$cell),]$nb_species
nW=indW*specW
gridW.1=gridW[which(gridW$IDcell %in% dcW.1$cell),]
centroW <- gCentroid(gridW.1,byid=TRUE)
dfW=data.frame(lon=as.numeric(centroW@coords[,1]),
lat=as.numeric(centroW@coords[,2]),
numberOfSequences=color_N(nW),
numberOfSequences_color=color_N(nW),
numberOfSpecies=color_nspec(specW),
numberOfSpecies_color=color_nspec(specW),
meanNumberOfIndividualsBySpecies=color_ind(indW),
meanNumberOfIndividualsBySpecies_color=color_ind(indW),
nseq=nW,
nspec=specW,
nind=indW)
dfW.ordered=dfW[order(dfW$nseq),]
dfW.ordered$numberOfSequences=factor(dfW.ordered$numberOfSequences,
levels=c("< 5","5 - 10", "10 - 20", "20 - 30",">30" ))
dfW.ordered=dfW.ordered[order(dfW.ordered$nspec),]
dfW.ordered$numberOfSpecies=factor(dfW.ordered$numberOfSpecies,
levels=c("2 - 3","3 - 5", "5 - 10", ">10" ))
dfW.ordered=dfW.ordered[order(dfW.ordered$nind),]
dfW.ordered$meanNumberOfIndividualsBySpecies=factor(dfW.ordered$meanNumberOfIndividualsBySpecies,
levels=c("2 - 4","4 - 6", "6 - 8", "8 - 10", ">10" ))
dfW.ordered$numberOfSequences_color=dfW.ordered$numberOfSequences
dfW.ordered$numberOfSpecies_color=dfW.ordered$numberOfSpecies
dfW.ordered$meanNumberOfIndividualsBySpecies_color=dfW.ordered$meanNumberOfIndividualsBySpecies
dfW.ordered=dfW.ordered[order(dfW.ordered$nseq),]
levels(dfW.ordered$numberOfSequences_color) <- c("#FFFFFF","#FFF0F0","#FFB0B0","#FF7575","#FF0000")
dfW.ordered=dfW.ordered[order(dfW.ordered$nspec),]
levels(dfW.ordered$numberOfSpecies_color) <- c("#FFFFFF","#FFF0F0","#FF6F6F","#FF0000")
dfW.ordered=dfW.ordered[order(dfW.ordered$nind),]
levels(dfW.ordered$meanNumberOfIndividualsBySpecies_color) <- c("#FFFFFF","#FFF0F0","#FFB0B0","#FF7575","#FF0000")
return(dfW.ordered)
}
## from a grid with number of sequences information for each cell
## 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)
## sum of number of element by 2 columns moF moM
## according to a category defined by the third column moCategory
aggregateBy2Colons <- function(moF,moM,moCategory){
numberSpecModT=aggregate(cbind(moF,moM), by=list(Category=moCategory), FUN=sum)
numbSpecMoSumColons=data.frame(category=numberSpecModT$Category, numberOfElements=numberSpecModT$moM+numberSpecModT$moF)
return(numbSpecMoSumColons)
}
##########################################################################
## load data
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]]
## remove marine data with no SST information
dcM.1=dcM.1[which(!is.na(dcM.1$SST)),]
## remove marine data with no CloMeanVal information
dcM.1=dcM.1[which(!is.na(dcM.1$cloMeanVal)),]
## remove marine data with no bathymetry information
dcM.1=dcM.1[which(!is.na(dcM.1$bathyVal)),]
## remove freshwater data with no mean regional temperature information
dcF.1=dcF.1[which(!is.na(dcF.1$SST)),]
## remove freshwater data with no bathyval information
dcF.1=dcF.1[which(!is.na(dcF.1$bathyVal)),]
## filter wrong freshwater assignation
wrongFreshwaterCells=c(3304,4741,8453,9146,11400,11474,11548,12434,12730,11860)
dcF.1=dcF.1[which(!dcF.1$cell %in% wrongFreshwaterCells),]
tc=read.csv("11-sequences_taxonomy_habitat/watertype_all_modeles_effectives_family.csv",header=T,sep=",")
## correct wrong size of family in NCBI taxonomy
tc[which(tc$taxonFamily == "Citharinidae"),]$ncbiNumberOfSpecies = tc[which(tc$taxonFamily == "Citharinidae"),]$ncbiNumberOfSpecies + 91
tc[which(tc$taxonFamily == "Gadidae"),]$ncbiNumberOfSpecies = 53
## number of species by family in two model(marine+freshwater)
monModelMARINENumbSpec=tc$dcMNumberOfSpecies
monModelFRESHWATERNumbSpec=tc$dcF346FreshwaterNumberOfSpecies
## number of sequences by family in two model(marine+freshwater)
monModelMARINENumbSeq=tc$dcMNumberOfSequences
monModelFRESHWATERNumbSeq=tc$dcF346FreshwaterNumberOfSequences
## models with at least 4 sequences by species and at least 2 species by cell
#monModelMARINENumbSpec=tc$S2I4C29MNumberOfSpecies
#monModelFRESHWATERNumbSpec=tc$S2I5C32FNumberOfSpecies
## models with at least 3 sequences by species and at least 2 species by cell
#monModelMARINENumbSpec=tc$S2I3C66MNumberOfSpecies
#monModelFRESHWATERNumbSpec=tc$S2I3C75FNumberOfSpecies
## models with at least 2 sequences by species and at least 8 species by cell
#monModelMARINENumbSpec=tc$S8I2C36MNumberOfSpecies
#monModelFRESHWATERNumbSpec=tc$S8I2C34FNumberOfSpecies
## models with at least 2 sequences by species and at least 3 species by cell
#monModelMARINENumbSpec=tc$S2I3C66MNumberOfSpecies
#monModelFRESHWATERNumbSpec=tc$S3I2C75FNumberOfSpecies
## number of species by family in two model(marine+freshwater)
NSpeF<-aggregateBy2Colons(monModelMARINENumbSpec,monModelFRESHWATERNumbSpec,tc$taxonFamily)
TotSpeFamilyNumbSeq<-aggregate(tc$ncbiNumberOfSpecies, by=list(Category=tc$taxonFamily), FUN=sum)
propFam=data.frame(proportion=(NSpeF$numberOfElements/TotSpeFamilyNumbSeq$x)*100, family=as.vector(TotSpeFamilyNumbSeq$Category))
## number of species by ORDER in two model(marine+freshwater)
NSpeO<-aggregateBy2Colons(monModelMARINENumbSpec,monModelFRESHWATERNumbSpec,tc$taxonOrder)
## Number of species order total
TotSpeOrderNumbSeq<-aggregate(tc$ncbiNumberOfSpecies, by=list(Category=tc$taxonOrder), FUN=sum)
propOrder=data.frame(proportion=(NSpeO$numberOfElements/TotSpeOrderNumbSeq$x)*100,
order=as.vector(TotSpeOrderNumbSeq$Category))
## number of sequences by family in two model(marine+freshwater)
NSeqF<-aggregateBy2Colons(monModelMARINENumbSeq,monModelFRESHWATERNumbSeq,tc$taxonFamily)
## number of sequences by order in two model(marine+freshwater)
NSeqO<-aggregateBy2Colons(monModelMARINENumbSeq,monModelFRESHWATERNumbSeq,tc$taxonOrder)
##########################################################################
## load worldcoast shapefile
grid <- readOGR("01-infos/grid_equalarea200km","gridFish.b_260418")
worldcoast <- readOGR("01-infos/ne_50m_land",layer="ne_50m_land")
## load rivers and lake centerlines shapefile
riverlakes <- readOGR("01-infos/ne_50m_rivers_lake_centerlines_scale_rank",layer="ne_50m_rivers_lake_centerlines_scale_rank")
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"))
## Quartile of the proportion of species and number of sequences
## for the 63 orders and 480 families
summary(propFam)
summary(propOrder)
##########################################################################
## 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)
## barplots
## Number of species per family
#barplot(propFam$proportion, las=2,names.arg=propFam$family,ylab="Prop species per order",main ="full models",cex.names = 0.5 )
## Number of species per order
## Number of sequences by family
#barplot(NSeqF$numberOfElements, las=2,names.arg=NSeqF$category,ylab="Number sequences per order",cex.names = 0.5 )
## Number of sequences by order
##########################################################################
## write pdf files
tiff(filename="10-figures/figureS5.tiff",width=27,height=24,units="cm",res=640,compression="lzw")
layout(mat=rbind(c(1,2),c(3,4),c(5,6)))
par(mar=c(1.5,2,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)")
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)
## 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)")
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)
## 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)")
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"]]
lybottom <- lytop - legg$rect[["h"]]
lxright <- lxleft + legg$rect[["w"]]
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')
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)
barplot(NSeqO$numberOfElements, las=2,names.arg=NSeqO$category,ylab="Number sequences per order",cex.names = 0.6 )
title("(b)", adj = 0)
dev.off()
00-scripts/step5/supplementary_figures/figureS6.R
View file @ 20a5ec16
......@@ -9,102 +9,167 @@
## Submited to Nature communications, 2019
##
## functions used to generate figures
## Supplementary Figure S6.
## Taxonomic coverage of the sequences used by the model
## (a) Number of species per order
## (b) Number of sequences per order
## Supplementary Figure 6. Sampling effect.
## (a,b) number of sequences per cell,
## (c,d) number of fish species per cell,
## (e,f) mean number of sequences per fish species per cell for marine species
## (a, c, e) and freshwater species (b, d,f) respectively.
##
##
##########################################################################
## Libraries
## no library required
lib_vect <-c("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")
## convert number of sequences value into color value
color_N <- function(gd){
gdc=gd
gdc[which(gd<5)] = "< 5"
gdc[which(gd>=5 & gd<10)] = "5 - 10"
gdc[which(gd>=10 & gd<20)] ="10 - 20"
gdc[which(gd>=20 & gd<30)] ="20 - 30"
gdc[which(gd>=30)] =">30"
return(gdc)
}
color_nspec <- function(gd){
gdc=gd
gdc[which(gd<3)] = "2 - 3"
gdc[which(gd>=3 & gd<5)] = "3 - 5"
gdc[which(gd>=5 & gd<10)] ="5 - 10"
gdc[which(gd>=10)] =">10"
return(gdc)
}
## sum of number of element by 2 columns moF moM
## according to a category defined by the third column moCategory
aggregateBy2Colons <- function(moF,moM,moCategory){
numberSpecModT=aggregate(cbind(moF,moM), by=list(Category=moCategory), FUN=sum)
numbSpecMoSumColons=data.frame(category=numberSpecModT$Category, numberOfElements=numberSpecModT$moM+numberSpecModT$moF)
return(numbSpecMoSumColons)
color_ind <- function(gd){
gdc=gd
gdc[which(gd<4)] = "2 - 4"
gdc[which(gd>=4 & gd<6)] = "4 - 6"
gdc[which(gd>=6 & gd<8)] ="6 - 8"
gdc[which(gd>=8 & gd<10)] ="8 - 10"
gdc[which(gd>=10)] =">10"
return(gdc)
}
## from a grid and table of descriptors
## create a grid with number of sequences information for each cell of the grid
grid_numberofseq <- function(gridW,dcW.1,dcW) {
indW=dcW[which(dcW$cell %in% dcW.1$cell),]$nb_indv_mean
specW=dcW[which(dcW$cell %in% dcW.1$cell),]$nb_species
nW=indW*specW
gridW.1=gridW[which(gridW$IDcell %in% dcW.1$cell),]
centroW <- gCentroid(gridW.1,byid=TRUE)
dfW=data.frame(lon=as.numeric(centroW@coords[,1]),
lat=as.numeric(centroW@coords[,2]),
numberOfSequences=color_N(nW),
numberOfSequences_color=color_N(nW),
numberOfSpecies=color_nspec(specW),
numberOfSpecies_color=color_nspec(specW),
meanNumberOfIndividualsBySpecies=color_ind(indW),
meanNumberOfIndividualsBySpecies_color=color_ind(indW),
nseq=nW,
nspec=specW,
nind=indW)
dfW.ordered=dfW[order(dfW$nseq),]
dfW.ordered$numberOfSequences=factor(dfW.ordered$numberOfSequences,
levels=c("< 5","5 - 10", "10 - 20", "20 - 30",">30" ))
dfW.ordered=dfW.ordered[order(dfW.ordered$nspec),]
dfW.ordered$numberOfSpecies=factor(dfW.ordered$numberOfSpecies,
levels=c("2 - 3","3 - 5", "5 - 10", ">10" ))
dfW.ordered=dfW.ordered[order(dfW.ordered$nind),]
dfW.ordered$meanNumberOfIndividualsBySpecies=factor(dfW.ordered$meanNumberOfIndividualsBySpecies,
levels=c("2 - 4","4 - 6", "6 - 8", "8 - 10", ">10" ))
dfW.ordered$numberOfSequences_color=dfW.ordered$numberOfSequences
dfW.ordered$numberOfSpecies_color=dfW.ordered$numberOfSpecies
dfW.ordered$meanNumberOfIndividualsBySpecies_color=dfW.ordered$meanNumberOfIndividualsBySpecies
dfW.ordered=dfW.ordered[order(dfW.ordered$nseq),]
levels(dfW.ordered$numberOfSequences_color) <- c("#FFFFFF","#FFF0F0","#FFB0B0","#FF7575","#FF0000")
dfW.ordered=dfW.ordered[order(dfW.ordered$nspec),]
levels(dfW.ordered$numberOfSpecies_color) <- c("#FFFFFF","#FFF0F0","#FF6F6F","#FF0000")
dfW.ordered=dfW.ordered[order(dfW.ordered$nind),]
levels(dfW.ordered$meanNumberOfIndividualsBySpecies_color) <- c("#FFFFFF","#FFF0F0","#FFB0B0","#FF7575","#FF0000")
return(dfW.ordered)
}
## from a grid with number of sequences information for each cell
## 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)
}
##########################################################################
## load data
tc=read.csv("11-sequences_taxonomy_habitat/watertype_all_modeles_effectives_family.csv",header=T,sep=",")
## correct wrong size of family in NCBI taxonomy
tc[which(tc$taxonFamily == "Citharinidae"),]$ncbiNumberOfSpecies = tc[which(tc$taxonFamily == "Citharinidae"),]$ncbiNumberOfSpecies + 91
tc[which(tc$taxonFamily == "Gadidae"),]$ncbiNumberOfSpecies = 53
## number of species by family in two model(marine+freshwater)
monModelMARINENumbSpec=tc$dcMNumberOfSpecies
monModelFRESHWATERNumbSpec=tc$dcF346FreshwaterNumberOfSpecies
## number of sequences by family in two model(marine+freshwater)
monModelMARINENumbSeq=tc$dcMNumberOfSequences
monModelFRESHWATERNumbSeq=tc$dcF346FreshwaterNumberOfSequences
## models with at least 4 sequences by species and at least 2 species by cell
#monModelMARINENumbSpec=tc$S2I4C29MNumberOfSpecies
#monModelFRESHWATERNumbSpec=tc$S2I5C32FNumberOfSpecies
## models with at least 3 sequences by species and at least 2 species by cell
#monModelMARINENumbSpec=tc$S2I3C66MNumberOfSpecies
#monModelFRESHWATERNumbSpec=tc$S2I3C75FNumberOfSpecies
## models with at least 2 sequences by species and at least 8 species by cell
#monModelMARINENumbSpec=tc$S8I2C36MNumberOfSpecies
#monModelFRESHWATERNumbSpec=tc$S8I2C34FNumberOfSpecies
## models with at least 2 sequences by species and at least 3 species by cell
#monModelMARINENumbSpec=tc$S2I3C66MNumberOfSpecies
#monModelFRESHWATERNumbSpec=tc$S3I2C75FNumberOfSpecies
## number of species by family in two model(marine+freshwater)
NSpeF<-aggregateBy2Colons(monModelMARINENumbSpec,monModelFRESHWATERNumbSpec,tc$taxonFamily)
TotSpeFamilyNumbSeq<-aggregate(tc$ncbiNumberOfSpecies, by=list(Category=tc$taxonFamily), FUN=sum)
propFam=data.frame(proportion=(NSpeF$numberOfElements/TotSpeFamilyNumbSeq$x)*100, family=as.vector(TotSpeFamilyNumbSeq$Category))
## number of species by ORDER in two model(marine+freshwater)
NSpeO<-aggregateBy2Colons(monModelMARINENumbSpec,monModelFRESHWATERNumbSpec,tc$taxonOrder)
## Number of species order total
TotSpeOrderNumbSeq<-aggregate(tc$ncbiNumberOfSpecies, by=list(Category=tc$taxonOrder), FUN=sum)
propOrder=data.frame(proportion=(NSpeO$numberOfElements/TotSpeOrderNumbSeq$x)*100,
order=as.vector(TotSpeOrderNumbSeq$Category))
## number of sequences by family in two model(marine+freshwater)
NSeqF<-aggregateBy2Colons(monModelMARINENumbSeq,monModelFRESHWATERNumbSeq,tc$taxonFamily)
## number of sequences by order in two model(marine+freshwater)
NSeqO<-aggregateBy2Colons(monModelMARINENumbSeq,monModelFRESHWATERNumbSeq,tc$taxonOrder)
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]]
## remove marine data with no SST information
dcM.1=dcM.1[which(!is.na(dcM.1$SST)),]
## remove marine data with no CloMeanVal information
dcM.1=dcM.1[which(!is.na(dcM.1$cloMeanVal)),]
## remove marine data with no bathymetry information
dcM.1=dcM.1[which(!is.na(dcM.1$bathyVal)),]
## remove freshwater data with no mean regional temperature information
dcF.1=dcF.1[which(!is.na(dcF.1$SST)),]
## remove freshwater data with no bathyval information
dcF.1=dcF.1[which(!is.na(dcF.1$bathyVal)),]
## filter wrong freshwater assignation
wrongFreshwaterCells=c(3304,4741,8453,9146,11400,11474,11548,12434,12730,11860)
dcF.1=dcF.1[which(!dcF.1$cell %in% wrongFreshwaterCells),]
##########################################################################
## Quartile of the proportion of species and number of sequences
## for the 63 orders and 480 families
summary(propFam)
summary(propOrder)
## load worldcoast shapefile
grid <- readOGR("01-infos/grid_equalarea200km","gridFish.b_260418")
worldcoast <- readOGR("01-infos/ne_50m_land",layer="ne_50m_land")
## load rivers and lake centerlines shapefile
riverlakes <- readOGR("01-infos/ne_50m_rivers_lake_centerlines_scale_rank",layer="ne_50m_rivers_lake_centerlines_scale_rank")
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"))
##########################################################################
## barplots
## Number of species per family
#barplot(propFam$proportion, las=2,names.arg=propFam$family,ylab="Prop species per order",main ="full models",cex.names = 0.5 )
## Number of species per order
## Number of sequences by family
#barplot(NSeqF$numberOfElements, las=2,names.arg=NSeqF$category,ylab="Number sequences per order",cex.names = 0.5 )
## Number of sequences by order
##########################################################################
## 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
pdf("10-figures/figureS6.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)
barplot(NSeqO$numberOfElements, las=2,names.arg=NSeqO$category,ylab="Number sequences per order",cex.names = 0.6 )
title("(b)", adj = 0)
tiff(filename="10-figures/figureS6.tiff",width=27,height=24,units="cm",res=640,compression="lzw")
layout(mat=rbind(c(1,2),c(3,4),c(5,6)))
par(mar=c(1.5,2,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)")
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)
## 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)")
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)
## 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)")
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"]]
lybottom <- lytop - legg$rect[["h"]]
lxright <- lxleft + legg$rect[["w"]]
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()
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