From 021f1ed4c32e097f0f38e0229d211bbdab2bcb78 Mon Sep 17 00:00:00 2001 From: Francesco Sabatini <francesco.sabatini@idiv.de> Date: Sat, 14 Nov 2020 13:41:11 +0100 Subject: [PATCH] Aligned Figure numbers - Output for all combinations --- 02_Mesobromion_ExamineOutput.R | 206 ++++++++++++++++++++-- 03_Figures_Simulations.R | 310 +++++++++++++++++++++++---------- 04_Additional_Figs.R | 88 ++++++---- 3 files changed, 462 insertions(+), 142 deletions(-) diff --git a/02_Mesobromion_ExamineOutput.R b/02_Mesobromion_ExamineOutput.R index 906ef04..66ee16a 100644 --- a/02_Mesobromion_ExamineOutput.R +++ b/02_Mesobromion_ExamineOutput.R @@ -348,7 +348,7 @@ mydata.best <- mydata %>% write_csv(mydata.best %>% dplyr::select(Trait.comb:sign_plus), - path = "_output/S5_BestSolutionTiers.cov.csv") + path = "_output/S9_BestSolutionTiers.cov.csv") ### Graph of all the best combinations with text legend (top.all <- ggplot(data=(mydata.best %>% @@ -397,6 +397,188 @@ ggsave(filename = "_pics/Fig5_Best_AllCombinations_CI_cov.png", dpi=400, width=6, height=2, topall.leg) + + +###### ________R1___________ ###### +###### R1.2. Import output from Cluster #### +##### R1.2.0. Trait labs for significant traits +traits.sign.cov <- read_delim(file="_data/Mesobromion/traits.v2.10perc.cov.sign.txt", delim="\t") + +traits.sign.cov <- traits.sign.cov %>% + as.data.frame() %>% + mutate_if(~is.character(.), .funs=~as.factor(.)) %>% + column_to_rownames("species0") + + +## adapt trait labs to sign traits only + trait.labs.sign.cov <- trait.labs %>% + filter(trait.name %in% colnames(traits.sign.cov)) %>% + arrange(match(trait.name, colnames(traits.sign.cov))) %>% + rename(Trait.comb.new=Trait.comb) %>% + mutate(Trait.comb=1:n()) %>% + dplyr::select(Trait.comb, everything(), -Trait.comb.new) + + + +##### R1.2.1 Cover values ###### +### sequential trait combo +myfilelist1 <- list.files(path="_derived/Mesobromion/Cover/R1_all/", pattern="HIDDENcov-nona2_[0-9]+_.RData", full.names = T) +dataFiles1 = purrr::map(myfilelist1, function(x){get(load(x))}) +#load("_derived/Mesobromion/PresAbs/HIDDEN_round_11.RData") +corXY.all = bind_rows(dataFiles1) %>% + as_tibble() %>% + distinct() +corXY.all.ci <- get.ci(corXY.all) +corXY.all.ci <- corXY.all.ci %>% + mutate(Trait.comb2=Trait.comb) %>% + separate(Trait.comb2, into=paste0("trait", 1:7)) %>% + mutate_at(.vars=vars(trait1:trait7), + .funs=~factor(., + levels=trait.labs.sign.cov$Trait.comb, + labels=trait.labs.sign.cov$trait.name)) %>% + arrange(ntraits, desc(Coef.obs)) %>% + #filter(ntraits>1) %>% + dplyr::select(Trait.comb, Test, n, ntraits, everything()) %>% + mutate(run="seq") + +rm( dataFiles1) #dataFiles0, + +### merge together +corXY.ci <- corXY.all.ci # %>% + +mydata <- corXY.ci + +######## R1.2.1.4 Best - Graph of r(XY) using best combination of traits at each level of interaction N ######## +### extract best combinations of traits +top.one.by.one <- get.best(mydata, N=1, labs=trait.labs.sign.cov) + +## Routine to extract the best combination at each level of interaction (up to max traits) +maxtraits <- 7 +for(nn in 1:maxtraits){ + if(nn==1) { + best.at.1 <- get.best(mydata, N=nn, labs=trait.labs.sign.cov) + newdata <- mydata %>% + filter_at(.vars=vars(trait1:trait7), + .vars_predicate = any_vars(. %in% best.at.1$trait.name | is.na(.))) + new.best.row <- newdata %>% + filter(Trait.comb==best.at.1$Trait.comb) + upper <- new.best.row$q975 + lower <- new.best.row$q025 + print(paste("new best at nn", nn, best.at.1$trait.name)) + best.progr <- best.at.1$Trait.comb + } + if(nn>1){ + better <- list() + better$Trait.comb <- newdata %>% + filter(ntraits==nn) %>% + filter(q025>upper) %>% + arrange(desc(Coef.obs)) %>% + slice(1) %>% + pull(Trait.comb) + + if(length(better$Trait.comb>0)){ + better$trait.name <- trait.labs.sign.cov %>% + filter(Trait.comb %in% strsplit(better$Trait.comb, split = "_")[[1]]) %>% + pull(trait.name) + + newdata <- newdata %>% + rowwise() %>% + mutate(nmatching= sum(unlist(strsplit(Trait.comb, "_")) %in% + unlist(strsplit(better$Trait.comb, "_")), + na.rm=T)) %>% + ungroup() %>% + filter(nmatching==nn) + + new.best.row <- newdata %>% + filter(Trait.comb==better$Trait.comb) + upper <- new.best.row$q975 + lower <- new.best.row$q025 + print(paste("new best at nn", nn, paste(better$trait.name, collapse=" "))) + best <- better + best.progr <- c(best.progr, better$Trait.comb) + } + } +} + +best.traits.cov <- corXY.ci %>% + filter(as.character(Trait.comb)==best.progr[length(best.progr)]) %>% + dplyr::select(trait1:trait7) %>% + mutate_all(~as.character(.)) %>% + dplyr::select(colnames(.)[which(colSums(is.na(.))==0)]) + +best.traits.cov <- as.character(best.traits.cov[1,]) +#"Leaf_Scleroph" "FP_Dur" "VP_Fragm" "Height" "SLA" + +### Create dataset with best combinations + all the one-way combinations +mydata.best <- mydata %>% + #filter_at(.vars=vars(trait1:trait7), + # # .vars_predicate = all_vars(. %in% best$trait.name | is.na(.))) %>% + # .vars_predicate = all_vars(. %in% traits.sign.alone.cov | is.na(.))) %>% + filter(ntraits>1) %>% + filter(sign_plus==T) %>% + arrange(ntraits, Coef.obs) %>% + group_by(ntraits) %>% + slice(n()) %>% + ungroup() %>% + bind_rows(corXY.ci %>% + filter(run=="seq") %>% + filter(ntraits==1)) %>% + #filter(trait1 %in% traits.sign.alone.cov)) %>% + arrange(ntraits, Coef.obs) %>% + mutate(seq=1:n()) %>% + mutate(sign_plus=factor(Trait.comb %in% best.progr)) + +write_csv(mydata.best %>% + dplyr::select(Trait.comb:sign_plus), + path = "_output/R1.S5_BestSolutionTiers.cov_allcombos.csv") + +### Graph of all the best combinations with text legend +(top.all <- ggplot(data=(mydata.best %>% + mutate(size0=.6+(as.numeric(sign_plus)-1)*.6))) + + geom_segment(aes(x=q025, xend=q975, y=seq, yend=seq, col="a", + lwd=size0)) + + geom_point(aes(x=Coef.obs, y=seq), pch=15) + + scale_y_continuous(breaks=mydata.best$seq, + labels=mydata.best$Trait.comb, name=NULL) + + scale_x_continuous(name="RD correlation") + + scale_size_identity() + + theme_bw() + + theme(panel.grid.minor = element_blank(), + axis.text = element_text(size=7), + legend.position = "none")) +# create legend of names +tt2 <- ttheme_minimal( + core = list(fg_params=list(cex = .7), + padding=unit(c(1, 1), "mm")), + colhead = list(fg_params=list(cex = .7)), + rowhead = list(fg_params=list(cex = .7))) + +ttlabs <- trait.labs.sign.cov %>% + mutate(Code=1:n()) + +tobold <- which(ttlabs$trait.name %in% best.traits.cov) +tg <- tableGrob(ttlabs %>% + dplyr::select(Code, Trait=Long_English_name) %>% + mutate(Trait=replace(x = Trait, + list = Trait=="Vegetative Propagation - Fragmentation", + values = "Veg. Propag. - Fragmentation")), + theme=tt2, rows = NULL) +## Make significant traits bold +for (i in (11 + tobold)) { + tg$grobs[[i]] <- editGrob(tg$grobs[[i]], gp=gpar(fontface="bold")) +} +#arrange into a panel +(topall.leg <- cowplot::plot_grid(top.all, tg, + nrow=1, rel_widths=c(0.60, 0.4))) +ggsave(filename = "_pics/R1/Fig5_R1_Best_AllCombinations_CI_cov.png", dpi=400, + width=6, height=2, topall.leg) + + + + + + + break() ###### ___________________ ###### @@ -755,7 +937,7 @@ PCA.fuzz1_3 <- basemap0 + ylab(paste("PC3 (", varexpl[3], "%)", sep="")) PC_fuzzy <- cowplot::plot_grid(PCA.fuzz1_2,PCA.fuzz1_3, nrow=1) -ggsave("_pics/S11_PC_Fuzzy_1-3.png", width=10, height=5, dpi=300, last_plot()) +ggsave("_pics/S13_PC_Fuzzy_1-3.png", width=10, height=5, dpi=300, last_plot()) #### 4.0.1 Alternative showing species scores #### tmp <- as.data.frame(pca.fuzz$CA$v[,1:3]*7) %>% @@ -812,8 +994,8 @@ PCAfuzz1_3.sp <- basemap0 %+% tmp + ylab(paste("PC3 (", varexpl[3], "%)", sep="")) -ggsave("_pics/S11a_PCA_Fuzzy_1-2_wSpecies.png", width=8, height=8, dpi=300, PCAfuzz1_2.sp) -ggsave("_pics/S11b_PCA_Fuzzy_1-3_wSpecies.png", width=8, height=8, dpi=300, PCAfuzz1_3.sp) +ggsave("_pics/S13a_PCA_Fuzzy_1-2_wSpecies.png", width=8, height=8, dpi=300, PCAfuzz1_2.sp) +ggsave("_pics/S13b_PCA_Fuzzy_1-3_wSpecies.png", width=8, height=8, dpi=300, PCAfuzz1_3.sp) @@ -833,7 +1015,7 @@ env.cor <- cor(env %>% scores.pca, use = "pairwise.complete.obs") #double check fuzz.cor <- cor(pca.fuzz$CA$u[,1:3], scores.pca) -sink("_output/S9_EnvFit_CWMs_env.txt") +sink("_output/S14_EnvFit_CWMs_env.txt") cwms.cor env.cor fuzz.cor @@ -953,8 +1135,8 @@ PCA3_4.sp <- basemap0 %+% tmp + ylab(paste("PC4 (", varexpl[4], "%)", sep="")) -ggsave("_pics/S10a_PCA_Beals_1-2_wSpecies.png", width=8, height=8, dpi=300, PCA1_2.sp) -ggsave("_pics/S10b_PCA_Beals_3-4_wSpecies.png", width=8, height=8, dpi=300, PCA3_4.sp) +ggsave("_pics/S14a_PCA_Beals_1-2_wSpecies.png", width=8, height=8, dpi=300, PCA1_2.sp) +ggsave("_pics/S14b_PCA_Beals_3-4_wSpecies.png", width=8, height=8, dpi=300, PCA3_4.sp) ###### _ ###### @@ -1155,7 +1337,7 @@ PCA.t2 <- baseplot + PC_traits <- cowplot::plot_grid(PCA.t1, PCA.t2, nrow=1) -ggsave("_pics/S6_PCA_Traits_1-4_only7.png", width=10, height=5, dpi=300, PC_traits) +ggsave("_pics/S10c_PCA_Traits_1-4_only7.png", width=10, height=5, dpi=300, PC_traits) ##### 4.3b Alternative version of figS6, showing the species #### tmp <- as.data.frame(pca.scores[,1:4]*.2) %>% @@ -1200,8 +1382,8 @@ PCA.t2.sp <- baseplot %+% tmp + xlab(paste("PC3 (", varexpl[3], "%)", sep="")) + ylab(paste("PC4 (", varexpl[4], "%)", sep="")) -ggsave("_pics/S6a_PCA_Traits_1-2_wSpecies.png", width=8, height=8, dpi=300, PCA.t1.sp) -ggsave("_pics/S6b_PCA_Traits_3-4_wSpecies.png", width=8, height=8, dpi=300, PCA.t2.sp) +ggsave("_pics/S10_PCA_Traits_1-2_wSpecies.png", width=8, height=8, dpi=300, PCA.t1.sp) +ggsave("_pics/S10b_PCA_Traits_3-4_wSpecies.png", width=8, height=8, dpi=300, PCA.t2.sp) #traits.dummy %>% filter(species0 %in% (tmp %>% filter(labels %in% c("Fes_pal", "Ses_alb", "Car_hum")) %>% pull(species))) %>% dplyr::select(any_of(starts_with(as.character(traits.sign.alone.cov)))) @@ -1304,7 +1486,7 @@ traits7 <- traits %>% # relocate(all_of(starts_with(as.character(best.traits.cov))), everything()) res <- cor(traits7, use = "pairwise.complete.obs") -png(file="_pics/S7_Correlations_Trait.png", width=8, height=6.5, units = "in", res=300) +png(file="_pics/S11_Correlations_Trait.png", width=8, height=6.5, units = "in", res=300) corrplot(res, type = "upper", tl.col = "black", tl.srt = 45, number.cex=0.6, addCoef.col = "black", diag=F) dev.off() @@ -1315,7 +1497,7 @@ res2 <- cor(CWM.wide %>% dplyr::select(any_of(traits.sign.alone.cov)) %>% ## caution selecting only numerical variables dplyr::select(sort(tidyselect::peek_vars())) %>% relocate(any_of(best.traits.cov), everything())) -png(file="_pics/S8_Correlations_CWMs.png", width=8, height=6.5, units = "in", res=300) +png(file="_pics/S12_Correlations_CWMs.png", width=8, height=6.5, units = "in", res=300) corrplot(res2, type = "upper", tl.col = "black", tl.srt = 45, number.cex=0.6, addCoef.col = "black", diag=F) dev.off() diff --git a/03_Figures_Simulations.R b/03_Figures_Simulations.R index fa2bdf1..aa01b8f 100644 --- a/03_Figures_Simulations.R +++ b/03_Figures_Simulations.R @@ -78,37 +78,38 @@ outp.summary <- FormatData(myfiles) mypalette <- palette(c("#e41a1c", #1 - red) "#ff7f00", #12 - orange "#984ea3", #13 - violet - "##ffed6f", #2 - yellow + "#ffed6f", #2 - yellow "#4daf4a", #23 - green "#377eb8")) #3 - blue +outp.summary2 <- outp.summary %>% + mutate(main=main/100) %>% + mutate(corr=factor(corr/10, levels=c(0, 0.4, 0.8), labels=paste0("Correlation = ", c(0, 0.4, 0.8)))) %>% + mutate(inter=factor(inter/10, levels=c(0, 0.3, 0.5), labels=paste0("Interaction = ", c(0, 0.3, 0.5)))) %>% + ungroup() %>% + dplyr::filter(stat.type=="XY") %>% + dplyr::filter(trait %in% c("1", "2", "1 2", "3", "1 3", "2 3")) %>% + mutate(trait=factor(trait, levels=c("1", "2", "3", "1 2", "1 3", "2 3"), + labels=c("t1", "t2", "tn", "t1 t2", "t1 tn", "t2 tn"))) %>% + mutate(mylinetype=as.factor(ifelse(grepl(pattern="tn", x = trait), 2, 1))) -ggplot(data=outp.summary %>% - mutate(main=main/100) %>% - mutate(corr=factor(corr/10, levels=c(0, 0.4, 0.8), labels=paste0("Correlation = ", c(0, 0.4, 0.8)))) %>% - mutate(inter=factor(inter/10, levels=c(0, 0.3, 0.5), labels=paste0("Interaction = ", c(0, 0.3, 0.5)))) %>% - ungroup() %>% - dplyr::filter(stat.type=="XY") %>% - dplyr::filter(trait %in% c("1", "2", "1 2", "3", "1 3", "2 3")) %>% - mutate(trait=factor(trait, levels=c("1", "2", "3", "1 2", "1 3", "2 3"), - labels=c("t1", "t2", "tn", "t1 t2", "t1 tn", "t2 tn"))) %>% - mutate(linetype=as.factor(ifelse(grepl(pattern="tn", x = trait), 2, 1)))) + - geom_line(aes(x=main, y=power, group=trait, col=trait, lty=linetype)) + - guides(linetype = FALSE) + - scale_color_manual("Trait\ncomb.", - values=c("#e41a1c", #1 - red) - "#e6ab02", #2 - yellow - "#377eb8", #3 - blue - "#d95f02", #12 - orange - "#984ea3", #13 - violet - "#4daf4a" #23 - green - ) - ) + +ggplot(data=outp.summary2) + + geom_line(aes(x=main, y=power, col=trait, lty=trait)) + facet_grid(corr~inter) + theme_bw() + scale_x_continuous(name="Effect of factor e1 -> trait t1") + scale_y_continuous(name="Prop. of significant r(XY)") + - theme(panel.grid = element_blank()) - + theme(panel.grid = element_blank()) + + scale_color_manual(name="Trait\ncomb.", + values = c("#e41a1c",#1 - red) + "#e6ab02",#2 - yellow + "#377eb8",#3 - blue + "#d95f02",#12 - orange + "#984ea3",#13 - violet + "#4daf4a" #23 - green + )) + + scale_linetype_manual(name="Trait\ncomb.", + values = c(1,1,2,1,2,2)) + ggsave(filename="_pics/R1/Fig2_CorrInte_02March.png", width=6, height=5, device="png", dpi = 300, last_plot()) @@ -170,23 +171,28 @@ hugepalette <- data.frame(trait=c('t1', 't2','t3', 't4', 'tn', 't1 t2 t3 t4', 't1 t2 t3 tn','t1 t2 t4 tn','t1 t3 t4 tn','t2 t3 t4 tn', 't1 t2 t3 t4 tn'), trait.col=factor(hugepalette0, levels=hugepalette0)) +hugelinetype <- c(1, 1, 1, 1, 2, + 1*1, 1*1, 1*1, 1*2, 1*1, 1*1, 1*2, 1*1, 1*2, 1*2, + 1*1*1, 1*1*1, 1*1*2, 1*1*1, 1*1*2, 1*1*2, 1*1*1, 1*1*2, 1*1*2, 1*1*2, + 1*1*1*1, 1*1*1*2, 1*1*1*2, 1*1*1*2, 1*1*1*2, + 1*1*1*1*2) outp.summary2 <- outp.summary2 %>% left_join(hugepalette, by="trait") fig3 <- ggplot(data=outp.summary2) + - geom_line(aes(x=main, y=power, group=trait, col=trait.col, lty=linetype)) + - guides(linetype = FALSE) + + geom_line(aes(x=main, y=power, col=trait, lty=trait)) + #group=trait, scale_x_continuous(name="Effect of factor e1 -> trait t1", n.breaks = 4) + scale_y_continuous(name="Prop. of significant r(XY)") + facet_grid(sel.ntraits.lab~ntraits.lab) + - scale_color_identity(guide = "legend", - labels= hugepalette$trait) + + scale_color_manual(name="trait", #guide = "legend", + values=hugepalette0) + + #labels= hugepalette$trait) + + scale_linetype_manual(name="trait", #guide="legend", + values=hugelinetype) + theme_bw() + theme(panel.grid = element_blank()) - - ncols <- c(2,4,4,3,1) @@ -196,10 +202,11 @@ for(tier in 1:5){ outp.summary.tier <- outp.summary2 %>% filter(sel.ntraits==tier) ncombinations <- length(levels(factor(outp.summary.tier$trait))) leg.list[[tier+1]] <- cowplot::get_legend(fig3 %+% outp.summary.tier + - guides(col=guide_legend(ncol=ncols[tier], byrow=TRUE))+ - scale_color_identity(name=ifelse(tier==1, "Trait combination - tier 1",paste0("tier - ", tier)), - guide = "legend", - labels= hugepalette$trait[col.used:(col.used+ncombinations-1)]) + guides(col=guide_legend(ncol=ncols[tier], byrow=TRUE))+ + scale_color_manual(name=ifelse(tier==1, "Trait combination - tier 1",paste0("tier - ", tier)), + values= hugepalette0[col.used:(col.used+ncombinations-1)]) + + scale_linetype_manual(name=ifelse(tier==1, "Trait combination - tier 1",paste0("tier - ", tier)), + values=hugelinetype[col.used:(col.used+ncombinations-1)]) ) col.used <- col.used + ncombinations } @@ -232,19 +239,17 @@ ggplot(data=outp.summary %>% mutate(trait=factor(trait, levels=c("1", "2", "3", "4"), labels=c("t1", "t2", "t3","tn"))) %>% mutate(linetype=as.factor(ifelse(grepl(pattern="tn", x = trait), 2, 1)))) + - geom_line(aes(x=main, y=power, group=trait, col=trait, lty=linetype)) + - guides(linetype = FALSE) + - #scale_colour_brewer(palette = "Dark2") + - #scale_color_manual("Trait\ncomb.", + geom_line(aes(x=main, y=power, col=trait, lty=trait)) + scale_color_manual("Trait", values=c("#e41a1c", #1 - red) "#e6ab02", #2 - yellow "#4daf4a", #23 - green - "#377eb8", #3 - blue - "#d95f02", #12 - orange - "#984ea3" #13 - violet - ) - ) + + "#377eb8" #3 - blue + #"#d95f02", #12 - orange + #"#984ea3" #13 - violet + )) + + scale_linetype_manual(name="Trait", + values = c(1,1,1,2)) + facet_grid(corr~inter) + theme_bw() + scale_x_continuous(name="Effect of factor e1 -> trait t1") + @@ -255,53 +260,119 @@ ggsave(filename="_pics/R1/FigS4_Extra_CorrInte_08Jul20.png", width=6, height=5, -### Additional Figures for JVS R1 #### -#### FIGURE SXXX - XW #### -mypath <- "_data/Experiment_30Oct2020_FactorInteraction&TraitCorr_XW" + +### Additional Figures for JVS R1 #### + +#### FIGURE R1.S5 - effect of adding neutral traits #### +#### import function +mypath <- "_data/Experiment_01Nov2020_NeutralTraitNumber/" myfiles <- list.files(path=mypath, pattern = "Summary.txt", recursive = T, full=T) outp.summary <- FormatData(myfiles) +outp.summary <- outp.summary %>% + rename(ntraits = inter) -mypalette <- palette(c("#e41a1c", #1 - red) - "#ff7f00", #12 - orange - "#984ea3", #13 - violet - "##ffed6f", #2 - yellow - "#4daf4a", #23 - green - "#377eb8")) #3 - blue +## plotting power for XY with corr +add.t.label <- function(x) { + x <- gsub(pattern=" ", replacement=" t", x=x, perl=T) + x <- paste0("t", x) + return(x) +} -ggplot(data=outp.summary %>% - mutate(main=main/100) %>% - mutate(corr=factor(corr/10, levels=c(0, 0.4, 0.8), labels=paste0("Correlation = ", c(0, 0.4, 0.8)))) %>% - mutate(inter=factor(inter/10, levels=c(0, 0.3, 0.5), labels=paste0("Interaction = ", c(0, 0.3, 0.5)))) %>% - ungroup() %>% - dplyr::filter(stat.type=="XY") %>% - dplyr::filter(trait %in% c("1", "2", "1 2", "3", "1 3", "2 3")) %>% - mutate(trait=factor(trait, levels=c("1", "2", "3", "1 2", "1 3", "2 3"), - labels=c("t1", "t2", "tn", "t1 t2", "t1 tn", "t2 tn"))) %>% - mutate(linetype=as.factor(ifelse(grepl(pattern="tn", x = trait), 2, 1)))) + - geom_line(aes(x=main, y=power, group=trait, col=trait, lty=linetype)) + - guides(linetype = FALSE) + - scale_color_manual("Trait\ncomb.", - values=c("#e41a1c", #1 - red) - "#e6ab02", #2 - yellow - "#377eb8", #3 - blue - "#d95f02", #12 - orange - "#984ea3", #13 - violet - "#4daf4a" #23 - green - )) + - facet_grid(corr~inter) + +outp.summary2 <- outp.summary %>% + ungroup() %>% + mutate(main=main/10) %>% + mutate(ntraits.lab=factor(ntraits, levels=3:5, labels=paste0("n. neutral traits = ", 1:3))) %>% + rowwise() %>% + mutate(sel.ntraits=factor(get.ntraits(trait))) %>% + mutate(sel.ntraits.lab=factor(sel.ntraits, levels=levels(sel.ntraits), labels=paste0("Comb. tier = ", levels(sel.ntraits)))) %>% + ungroup() %>% + dplyr::filter(stat.type=="XY") %>% + mutate(trait=factor(trait)) %>% + mutate(trait=factor(trait, levels=levels(trait), + labels=add.t.label(levels(trait)))) + +### rename null trait to "tn" +outp.summary2 <- outp.summary2 %>% + #mutate(tn.name=paste0("t", ntraits)) %>% + mutate(trait=str_replace(trait, pattern="t3", replacement="tn1")) %>% + mutate(trait=str_replace(trait, pattern="t4", replacement="tn2")) %>% + mutate(trait=str_replace(trait, pattern="t5", replacement="tn3")) %>% + mutate(linetype=as.factor(ifelse(grepl(pattern="tn", x = trait), 2, 1))) + +#reorder factors +outp.summary2$trait <- factor(outp.summary2$trait, levels= c('t1', 't2','tn1', 'tn2', 'tn3', + 't1 t2', 't1 tn1','t1 tn2','t1 tn3','t2 tn1','t2 tn2','t2 tn3','tn1 tn2','tn1 tn3', 'tn2 tn3', + 't1 t2 tn1', 't1 t2 tn2','t1 t2 tn3','t1 tn1 tn2', 't1 tn1 tn3','t1 tn2 tn3','t2 tn1 tn2','t2 tn1 tn3', 't2 tn2 tn3','tn1 tn2 tn3', + 't1 t2 tn1 tn2', 't1 t2 tn1 tn3','t1 t2 tn2 tn3','t1 tn1 tn2 tn3','t2 tn1 tn2 tn3', + 't1 t2 tn1 tn2 tn3')) + +hugepalette0 <- c(RColorBrewer::brewer.pal(4, "Dark2"), + gray(0.2), + RColorBrewer::brewer.pal(10, "Paired"), + RColorBrewer::brewer.pal(10, "Set3"), + RColorBrewer::brewer.pal(5, "Pastel1"), "brown") +#change tone of yellow of t1-t2-tn2 +hugepalette0[17] <- "#ccebc5" +hugepalette <- data.frame(trait=c('t1', 't2','tn1', 'tn2', 'tn3', + 't1 t2', 't1 tn1','t1 tn2','t1 tn3','t2 tn1','t2 tn2','t2 tn3','tn1 tn2','tn1 tn3', 'tn2 tn3', + 't1 t2 tn1', 't1 t2 tn2','t1 t2 tn3','t1 tn1 tn2', 't1 tn1 tn3','t1 tn2 tn3','t2 tn1 tn2','t2 tn1 tn3', 't2 tn2 tn3','tn1 tn2 tn3', + 't1 t2 tn1 tn2', 't1 t2 tn1 tn3','t1 t2 tn2 tn3','t1 tn1 tn2 tn3','t2 tn1 tn2 tn3', + 't1 t2 tn1 tn2 tn3'), + trait.col=factor(hugepalette0, levels=hugepalette0)) + +hugelinetype <- c(1, 1, 1, 1, 2, + 1*1, 1*1, 1*1, 1*2, 1*1, 1*1, 1*2, 1*1, 1*2, 1*2, + 1*1*1, 1*1*1, 1*1*2, 1*1*1, 1*1*2, 1*1*2, 1*1*1, 1*1*2, 1*1*2, 1*1*2, + 1*1*1*1, 1*1*1*2, 1*1*1*2, 1*1*1*2, 1*1*1*2, + 1*1*1*1*2) + +outp.summary2 <- outp.summary2 %>% + left_join(hugepalette, by="trait") + +figS5 <- ggplot(data=outp.summary2) + + geom_line(aes(x=main, y=power, col=trait, lty=trait)) + #group=trait, + scale_x_continuous(name="Effect of factor e1 -> trait t1", n.breaks = 4) + + scale_y_continuous(name="Prop. of significant r(XY)") + + facet_grid(sel.ntraits.lab~ntraits.lab) + + scale_color_manual(name="trait", #guide = "legend", + values=hugepalette0) + + scale_linetype_manual(name="trait", #guide="legend", + values=hugelinetype) + theme_bw() + - scale_x_continuous(name="Effect of factor e1 -> trait t1") + - scale_y_continuous(name="Prop. of significant r(XW)") + theme(panel.grid = element_blank()) -ggsave(filename="_pics/R1/FigSXXX_CorrInte_30November_XW.png", width=6, height=5, device="png", dpi = 300, last_plot()) +ncols <- c(2,4,4,3,1) + + +leg.list <- list() +col.used <- 1 +for(tier in 1:5){ + outp.summary.tier <- outp.summary2 %>% filter(sel.ntraits==tier) + ncombinations <- length(levels(factor(outp.summary.tier$trait))) + leg.list[[tier+1]] <- cowplot::get_legend(figS5 %+% outp.summary.tier + + guides(col=guide_legend(ncol=ncols[tier], byrow=TRUE))+ + scale_color_manual(name=ifelse(tier==1, "Trait combination - tier 1",paste0("tier - ", tier)), + values= hugepalette0[col.used:(col.used+ncombinations-1)]) + + scale_linetype_manual(name=ifelse(tier==1, "Trait combination - tier 1",paste0("tier - ", tier)), + values=hugelinetype[col.used:(col.used+ncombinations-1)]) + ) + col.used <- col.used + ncombinations +} +leg.list[[tier+2]] <- NULL + +figS5.panel <- cowplot::plot_grid(figS5 + theme(legend.position = "none"), + cowplot::plot_grid(plotlist = leg.list, nrow = 7, + rel_heights = c(0.05, .2,.2,.2,.2,.2, 0.15), align="hv"), + nrow=1, rel_widths = c(0.55,.45)) +ggsave(filename="_pics/R1/FigS5_Neutral_TraitNumber.png", + width=9, height=7, device="png", dpi = 300, figS5.panel) -### FIGURE SXXY - SBM #### +#### FIGURE R1.S6a - XW #### mypath <- "_data/Experiment_30Oct2020_FactorInteraction&TraitCorr_XW" myfiles <- list.files(path=mypath, pattern = "Summary.txt", recursive = T, full=T) outp.summary <- FormatData(myfiles) @@ -318,14 +389,13 @@ ggplot(data=outp.summary %>% mutate(corr=factor(corr/10, levels=c(0, 0.4, 0.8), labels=paste0("Correlation = ", c(0, 0.4, 0.8)))) %>% mutate(inter=factor(inter/10, levels=c(0, 0.3, 0.5), labels=paste0("Interaction = ", c(0, 0.3, 0.5)))) %>% ungroup() %>% - dplyr::filter(stat.type=="SbM") %>% - dplyr::filter(envir==1) %>% + dplyr::filter(stat.type=="XY") %>% dplyr::filter(trait %in% c("1", "2", "1 2", "3", "1 3", "2 3")) %>% mutate(trait=factor(trait, levels=c("1", "2", "3", "1 2", "1 3", "2 3"), - labels=c("t1", "t2", "tn", "t1 t2", "t1 tn", "t2 tn"))) %>% - mutate(linetype=as.factor(ifelse(grepl(pattern="tn", x = trait), 2, 1)))) + - geom_line(aes(x=main, y=power, group=trait, col=trait, lty=linetype)) + - guides(linetype = FALSE) + + labels=c("t1", "t2", "tn", "t1 t2", "t1 tn", "t2 tn")))# %>% + #mutate(linetype=as.factor(ifelse(grepl(pattern="tn", x = trait), 2, 1))) + ) + + geom_line(aes(x=main, y=power, col=trait, lty=trait)) + #group=trait, scale_color_manual("Trait\ncomb.", values=c("#e41a1c", #1 - red) "#e6ab02", #2 - yellow @@ -334,16 +404,19 @@ ggplot(data=outp.summary %>% "#984ea3", #13 - violet "#4daf4a" #23 - green )) + + scale_linetype_manual(name="Trait\ncomb.", + values = c(1,1,2,1,2,2)) + facet_grid(corr~inter) + theme_bw() + scale_x_continuous(name="Effect of factor e1 -> trait t1") + - scale_y_continuous(name="Prop. of significant tests") + + scale_y_continuous(name="Prop. of significant r(XW)") + theme(panel.grid = element_blank()) -ggsave(filename="_pics/R1/FigSXXY_CorrInte_30November_SbM.png", width=6, height=5, device="png", dpi = 300, last_plot()) +ggsave(filename="_pics/R1/FigS6a_CorrInte_30November_XW.png", width=6, height=5, device="png", dpi = 300, last_plot()) + -#### FIGURE SXXO - XE #### +#### FIGURE R1.S6b - XE #### mypath <- "_data/Experiment_02Mar2020_FactorInteraction&TraitCorr" myfiles <- list.files(path=mypath, pattern = "Summary.txt", recursive = T, full=T) outp.summary <- FormatData(myfiles) @@ -366,8 +439,7 @@ ggplot(data=outp.summary %>% mutate(trait=factor(trait, levels=c("1", "2", "3", "1 2", "1 3", "2 3"), labels=c("t1", "t2", "tn", "t1 t2", "t1 tn", "t2 tn"))) %>% mutate(linetype=as.factor(ifelse(grepl(pattern="tn", x = trait), 2, 1)))) + - geom_line(aes(x=main, y=power, group=trait, col=trait, lty=linetype)) + - guides(linetype = FALSE) + + geom_line(aes(x=main, y=power, col=trait, lty=trait)) + #group=trait, scale_color_manual("Trait\ncomb.", values=c("#e41a1c", #1 - red) "#e6ab02", #2 - yellow @@ -375,18 +447,22 @@ ggplot(data=outp.summary %>% "#d95f02", #12 - orange "#984ea3", #13 - violet "#4daf4a" #23 - green - ) - ) + + )) + + scale_linetype_manual(name="Trait\ncomb.", + values = c(1,1,2,1,2,2)) + facet_grid(corr~inter) + theme_bw() + scale_x_continuous(name="Effect of factor e1 -> trait t1") + scale_y_continuous(name="Prop. of significant r(XE)") + theme(panel.grid = element_blank()) -ggsave(filename="_pics/R1/FigSXXO_CorrInte_02March_XE.png", width=6, height=5, device="png", dpi = 300, last_plot()) +ggsave(filename="_pics/R1/FigS6b_CorrInte_02March_XE.png", width=6, height=5, device="png", dpi = 300, last_plot()) + + -#### FIGURE SXXP rXY & rXW vs Sample Size #### + +#### FIGURE R1.S7 rXY & rXW vs Sample Size #### mypath <- "_data/Experiment_30Oct2020_FactorInteraction&TraitCorr_XY_SampleSize_Main=040_Inter=00_Corr=00_v21169/" myfiles <- list.files(path=mypath, pattern = "Summary.txt", recursive = T, full=T) myfiles <- myfiles[!grepl("_new", x=myfiles)] ## exclude 'new' directories --> correct? @@ -419,7 +495,7 @@ ggplot(data= outp.summary %>% theme(panel.grid = element_blank(), legend.position = c(0.8, 0.2)) -ggsave(filename="_pics/R1/FigSXXP_CorrInte_02March_SampleSize.png", +ggsave(filename="_pics/R1/FigS7_CorrInte_02March_SampleSize.png", width=3, height=3, device="png", dpi = 300, last_plot()) @@ -432,3 +508,51 @@ ggsave(filename="_pics/R1/FigSXXP_CorrInte_02March_SampleSize.png", +### FIGURE SXXY - SBM - only for review #### +mypath <- "_data/Experiment_30Oct2020_FactorInteraction&TraitCorr_XW" +myfiles <- list.files(path=mypath, pattern = "Summary.txt", recursive = T, full=T) +outp.summary <- FormatData(myfiles) + +mypalette <- palette(c("#e41a1c", #1 - red) + "#ff7f00", #12 - orange + "#984ea3", #13 - violet + "##ffed6f", #2 - yellow + "#4daf4a", #23 - green + "#377eb8")) #3 - blue + +ggplot(data=outp.summary %>% + mutate(main=main/100) %>% + mutate(corr=factor(corr/10, levels=c(0, 0.4, 0.8), labels=paste0("Correlation = ", c(0, 0.4, 0.8)))) %>% + mutate(inter=factor(inter/10, levels=c(0, 0.3, 0.5), labels=paste0("Interaction = ", c(0, 0.3, 0.5)))) %>% + ungroup() %>% + dplyr::filter(stat.type=="SbM") %>% + dplyr::filter(envir==1) %>% + dplyr::filter(trait %in% c("1", "2", "1 2", "3", "1 3", "2 3")) %>% + mutate(trait=factor(trait, levels=c("1", "2", "3", "1 2", "1 3", "2 3"), + labels=c("t1", "t2", "tn", "t1 t2", "t1 tn", "t2 tn"))) %>% + mutate(linetype=as.factor(ifelse(grepl(pattern="tn", x = trait), 2, 1)))) + + geom_line(aes(x=main, y=power, col=trait, lty=trait)) + #group=trait, + scale_color_manual("Trait\ncomb.", + values=c("#e41a1c", #1 - red) + "#e6ab02", #2 - yellow + "#377eb8", #3 - blue + "#d95f02", #12 - orange + "#984ea3", #13 - violet + "#4daf4a" #23 - green + )) + + scale_linetype_manual(name="Trait\ncomb.", + values = c(1,1,2,1,2,2)) + + facet_grid(corr~inter) + + theme_bw() + + scale_x_continuous(name="Effect of factor e1 -> trait t1") + + scale_y_continuous(name="Prop. of significant tests") + + theme(panel.grid = element_blank()) + +ggsave(filename="_pics/R1/FigSXXY_CorrInte_30November_SbM.png", width=6, height=5, device="png", dpi = 300, last_plot()) + + + + + + + diff --git a/04_Additional_Figs.R b/04_Additional_Figs.R index 95b52d4..1873263 100644 --- a/04_Additional_Figs.R +++ b/04_Additional_Figs.R @@ -122,6 +122,7 @@ do.the.parse <- function(toparse) { } +### Figure R1.S8 - Comparison beta vs cor #### ## Function to create Figure SXXV create.panel <- function(x){ gg.betaW <- ggplot(data=x %>% @@ -132,7 +133,6 @@ create.panel <- function(x){ geom_density(aes(value)) + xlab("Proportional Beta Diversity (W)") + xlim(c(0,1.1)) + - ylim(c(0,60)) + theme_classic() gg.betaY <- gg.betaW %+% @@ -149,7 +149,8 @@ create.panel <- function(x){ group_by(dataset, metric, matrix) %>% summarize(value=max(abs(value))))+ xlab("Cor(WE)") + - ylab(NULL) + ylab(NULL) + + ylim(c(0,13.5)) gg.corY <- gg.betaW %+% (x %>% @@ -159,56 +160,62 @@ create.panel <- function(x){ group_by(dataset, metric, matrix) %>% summarize(value=max(abs(value)))) + xlab("Cor(YE)") + - ylab(NULL) + ylab(NULL)+ + ylim(c(0,13.5)) gg.panel <- cowplot::plot_grid(gg.betaW, gg.corW, gg.betaY, gg.corY, - nrow=2, rel_widths = c(1,1.06)) + nrow=2, rel_widths = c(1,1.06), align = "hv") return(gg.panel) } - - - -### Figure SXXV - comparison beta vs cor #### +# Set path of files to import and parse mypath <- "_data/Experiment_30Oct2020_FactorInteraction&TraitCorr_XY_SampleSize_Main=040_Inter=00_Corr=00_v21169" myfiles <- list.files(path=mypath, pattern = "FinalSimulatedData.txt", recursive = T, full.names = T) myfiles <- myfiles[grepl("_new", x=myfiles)] +# Loop over files to import and parse. Create graphs for(i in 1:length(myfiles)){ toparse <- myfiles[i] sampleN <- regmatches(toparse, gregexpr("N=[[:digit:]]+", toparse))[[1]] div.summary <- do.the.parse(toparse) gg.out <- create.panel(div.summary) - ggsave(filename = paste0("_pics/R1/FigSXXZ_Panel_BetaCor_", sampleN, ".png"), + ggsave(filename = paste0("_pics/R1/FigS8_Panel_BetaCor_", sampleN, ".png"), width=6, height=5, device="png", dpi = 300, plot = gg.out) } -#### Figure SXXK - Comparison abg across combinatios Inter X Corr (Main=0.3) #### + + +#### Figure R1.S1d - Comparison abg across combinatios Inter X Corr (Main=0.3) #### create.panel2 <- function(xx, i, tit){ require(ggpubr) - # alpha, min, mean, max - gg.alpha <- ggplot(data=xx %>% - filter(matrix=="W") %>% - filter(metric %in% c("richness", "alpha")) %>% - mutate(with=ifelse(!is.na(with), paste0("OTU Rich (", with, ")"), with)) %>% - mutate(with=ifelse(is.na(with), "Eq. OTU (mean)", with))) + - geom_density(aes(value, group=with, col=with), alpha=0.7, show.legend=FALSE)+ - stat_density(aes(x=value, colour=with), - geom="line",position="identity") + - scale_color_brewer(palette="Dark2", name=NULL) + - theme_classic() + - xlab("Alpha diversity (OTU Richness)") + - theme(legend.position = c(0.65, 0.9)#, - #legend.text = element_text(size=7) - ) + + gg.rich <- ggplot(data=xx %>% + filter(matrix=="W") %>% + filter(metric %in% c("richness"))) + + geom_density(aes(value, group=with, col=with), alpha=0.7, show.legend=FALSE)+ + stat_density(aes(x=value, colour=with), + geom="line",position="identity") + + scale_color_brewer(palette="Dark2", name=NULL) + + theme_classic() + + xlab("Species richness") + + theme(legend.position = c(0.75, 0.9)) + xlim(c(-.1,100)) + - ylim(c(0,0.25)) - + ylim(c(0,0.2)) + + # alpha, min, mean, max + gg.alpha <- gg.rich %+% + (xx %>% + filter(matrix == "W") %>% + filter(metric == "alpha")) + + geom_density(aes(value), show.legend = FALSE) + + xlab("Mean alpha diversity") + + ylim(c(0, 0.07)) + + ylab(NULL) + # beta gg.beta <- ggplot(data=xx %>% @@ -225,7 +232,7 @@ create.panel2 <- function(xx, i, tit){ (xx %>% filter(matrix=="W") %>% filter(metric %in% "propbeta")) + - xlab("Proportional Beta Diversity") + + xlab("Proportional beta Diversity") + xlim(c(0,0.7)) + ylim(c(0,4)) #gamma @@ -243,12 +250,14 @@ create.panel2 <- function(xx, i, tit){ )) if(i!=1){ - gg.alpha <- gg.alpha + + gg.rich <- gg.rich + theme(legend.position="none") - } - + } if(i!=3){ + gg.rich <- gg.rich + + xlab(NULL) + + theme(axis.text.x = element_blank()) gg.alpha <- gg.alpha + xlab(NULL) + theme(axis.text.x = element_blank()) @@ -263,27 +272,32 @@ create.panel2 <- function(xx, i, tit){ theme(axis.text.x = element_blank()) } - gg.panel <- cowplot::plot_grid(gg.title, gg.alpha, gg.beta, gg.propbeta, gg.gamma, - nrow=1, rel_widths = c(0.08, 1,.94, .94, .94)) + gg.panel <- cowplot::plot_grid(gg.title, gg.rich, gg.alpha, gg.beta, gg.propbeta, gg.gamma, + nrow=1, rel_widths = c(0.08, 1, .94, .94, .94, .94)) return(gg.panel) } - +## set path of files to import and parse mypath <- "_data/Experiment_30Oct2020_FactorInteraction&TraitCorr_XY_DataExamples/" myfiles <- list.files(path=mypath, pattern = "FinalSimulatedData.txt", recursive = T, full.names = T) +#loop of files, parse and create graphs panel.list <- list() for(i in 1:length(myfiles)){ toparse <- myfiles[i] Inter <- regmatches(toparse, gregexpr("Inter=[[:digit:]]+", toparse))[[1]] + Inter <- gsub(pattern="=0", replacement=" = 0.", x = Inter) + Inter <- gsub(pattern="0.0", replacement="0", x = Inter) Corr <- regmatches(toparse, gregexpr("Corr=[[:digit:]]+", toparse))[[1]] + Corr <- gsub(pattern="=0", replacement=" = 0.", x = Corr) + Corr <- gsub(pattern="0.0", replacement="0", x = Corr) div.summary <- do.the.parse(toparse) panel.list[[i]] <- create.panel2(div.summary, i=i,tit=paste(Inter, Corr)) } -gg.SXXK <- cowplot::plot_grid(plotlist=panel.list, nrow=3) -ggsave(filename = "_pics/R1/FigSXXK_Panel_abg.png", - width=10, height=6, device="png", dpi = 300, plot = gg.SXXK) +gg.S1d <- cowplot::plot_grid(plotlist=panel.list, nrow=3) +ggsave(filename = "_pics/R1/FigS1d_Panel_abg.png", + width=10, height=6, device="png", dpi = 300, plot = gg.S1d) -- GitLab