# PLEASE NOTE that loading the datasets is quite memory intensive. # CONSIDER running this script on a machine with at least a few # gigagbytes of free memory. library(readxl) library(tidyr) library(ggplot2) library(ggpubr) library(ggrepel) library(dplyr) library(data.table) library(ape) library(ggfortify) library(factoextra) library(viridis) library(ggConvexHull) library(pls) library(ggpubr) library(phytools) library(wesanderson) library(ggtree) data_joined = fread("FILE_S1.1_data_joined.csv.gz") data_traits = fread("FILE_S1.2_data_traits.csv.gz") phylo = read.tree("FILE_S1.3_phylo.new") trait_vars = c("Amax_mean","LMA_g.cm2_mean", "Fp_N_mm_mean","Total_tannin_mg.g") vein_vars = c("VD","MST","ER","CR") focal_codes = c("BNT-T212-BSH","DAS1-T010227-BSH","SER-T462-BSH","SLF-T53-BSH") # do SI analysis with all traits trait_vars_other = c("total_K_mg.g", "total_Ca_mg.g", "total_Mg_mg.g", "total_P_mg.g", "X15N_per_mil", "C_perc", "X13C_per_mil", "N_perc", "leaf_thickness_mm_mean", "LDMC_mg.g_mean", "chla_mg.g", "chlb_mg.g", "carot_mg.g", "hemicellulose_perc", "cellulose_perc", "lignin_recalcitrants_perc", "Total_phenol_mg.g") # calculate fit of results data_joined$FBeta2 %>% mean data_joined$FBeta2 %>% sd # count samples data_for_counting = data_joined %>% select(code,species,VNE) %>% left_join(data_traits %>% rename(code=sample_code) %>% select(code,forestplots_name)) %>% group_by(forestplots_name) %>% unique %>% summarize(num.samples=n_distinct(code),num.species=n_distinct(species), num.edges=sum(VNE)) # make the main text table write.csv(data_for_counting, 'data_for_counting.csv',row.names=F) # number of species data_joined %>% select(code,species,Family) %>% summarize(species=n_distinct(species),families=n_distinct(Family),leaves=n_distinct(code)) # ranges of input data 1000*quantile(data_joined$width_threshold[data_joined$width_threshold>0], c(0.01,0.99)) quantile(data_joined$TotA, c(0.01,0.99)) # add sun-shade info data_joined$is_shade_leaf = grepl("-BSH",data_joined$code) width_low = 0 #min(data_joined$width_threshold,na.rm=T) width_high = 0.2 # this could be as high as 0.6 #max(data_joined$width_threshold,na.rm=T) # bin the axis-wise data width_bins <- seq(width_low, width_high,length.out=50) data_joined$width_threshold_binned <- cut(data_joined$width_threshold,breaks=width_bins,labels=FALSE,include.lowest=TRUE) data_joined$width_threshold_binned <- width_bins[data_joined$width_threshold_binned] all_scales = expand.grid(width_threshold_binned=width_bins) # get bin-median values data_median = data_joined %>% group_by(code, species, width_threshold_binned) %>% dplyr::filter(Boundary != 0 & width_threshold_binned > 0) %>% summarize(VD.median=median(VD), MST.median=median(MST),ER.median=median(ER),CR.median=median(CR)) %>% ungroup # and insert an NA anywhere the bin exists but the data don't all_scales_codes = expand.grid(width_threshold_binned=width_bins,code=unique(data_median$code)) all_scales_codes$code = as.character(all_scales_codes$code) data_median = data_median %>% do(right_join(x=.,y=all_scales_codes,by = c("code", "width_threshold_binned"))) %>% rename(rmin=width_threshold_binned) %>% mutate(vein.exists=as.numeric((!is.na(VD.median)) & (!is.na(MST.median)) & (!is.na(ER.median)) & (!is.na(CR.median)))) %>% select(-species) %>% left_join(data_joined %>% select(code,Family,Genus,species) %>% unique,by="code") # look up names data_joined %>% dplyr::filter(code %in% focal_codes) %>% select(code, species, Family) %>% unique %>% arrange(code) draw_figure <- function(data, y_stat, focal_codes, y_label, do_log=TRUE) { g_all = ggplot(data, aes_string(x="rmin",y=y_stat,group="code")) + geom_line(alpha=0.2) + geom_line(data=data %>% dplyr::filter(code %in% focal_codes) %>% dplyr::filter(vein.exists==1), aes_string(x="rmin",y=y_stat,group="code",col='species'),alpha=0.5,size=1) + geom_line(data=data %>% dplyr::filter(code %in% focal_codes), aes_string(x="rmin",y=y_stat,group="code",col='species'),size=1) + geom_point(data=data %>% dplyr::filter(code %in% focal_codes), aes_string(x="rmin",y=y_stat,group="code",col='species'),size=2) + theme_bw() + theme(legend.text=element_text(face="italic")) + scale_color_brewer(palette="Set1",name="Species") + ylab(y_label) + xlab(expression(paste("r"["min"], " (mm)"))) if (do_log==TRUE) { g_all = g_all + scale_y_log10() } return(g_all) } g_er = draw_figure(data=data_median,y_stat="ER.median",y_label="Elongation ratio (ER)", focal_codes=focal_codes) g_cr = draw_figure(data=data_median,y_stat="CR.median",y_label="Circularity ratio (CR)", focal_codes=focal_codes) g_vd = draw_figure(data=data_median,y_stat="VD.median",y_label=expression(paste("Vein density (VD, mm"^{-1},")")), focal_codes=focal_codes) g_mst = draw_figure(data=data_median,y_stat="MST.median",y_label="Minimum spanning tree ratio (MST)", focal_codes=focal_codes) g_multiscale_stats = ggarrange(g_er, g_cr, g_vd, g_mst, common.legend=TRUE,legend='bottom', align='hv', nrow=2,ncol=2, labels=c("(a)","(b)","(c)","(d)")) ggsave(g_multiscale_stats,file='g_FIG3.png',width=9,height=7) ########## # figure out how many scales we have sufficient data at data_at_scale = data_joined %>% mutate(bin = cut(data_joined$width_threshold,breaks=100)) %>% group_by(bin) %>% summarize(num_codes=length(unique(code))) %>% mutate(bin_label=as.numeric(gsub("\\(","",sapply(strsplit(as.character(bin),","),head,1)))) %>% dplyr::filter(bin_label > 0) g_data_available = ggplot(data_at_scale,aes(x=bin_label,y=num_codes)) + geom_line(size=1) + theme_bw() + geom_vline(xintercept = 0,color='red') + geom_vline(xintercept = 0.2,color='red') + xlab(expression(paste("r"["min"], " (mm)"))) + ylab("Number of samples") + scale_y_log10() ggsave(g_data_available, file='g_FIGS1.png',width=8,height=5) fill_gaps_linear <- function(x,y) { df_no_na = data.frame(x,y) %>% na.omit # do linear interpolation, assume extrapolation = values at edge of data af = approxfun(df_no_na$x,df_no_na$y,rule=2) yv_new = af(x) return(yv_new) } # fill gaps with approxfun data_median_linear_fill = data_median %>% group_by(code) %>% mutate(VD.median=fill_gaps_linear(x=rmin,y=VD.median)) %>% mutate(ER.median=fill_gaps_linear(x=rmin,y=ER.median)) %>% mutate(CR.median=fill_gaps_linear(x=rmin,y=CR.median)) %>% mutate(MST.median=fill_gaps_linear(x=rmin,y=MST.median)) # transform data for more normality and to filter out outliers before PCA data_for_pca = data_median %>% dplyr::filter(vein.exists==1) %>% select(VD.median,MST.median,ER.median,CR.median,rmin,code) %>% #filter(ER.median < 10) %>% #filter(MST.median > 0.5) %>% mutate(VD.median=sqrt(VD.median), ER.median=sqrt(ER.median)) %>% na.omit # run PCA pca_all = prcomp(data_for_pca %>% select_at(vars(contains("median"))), center=TRUE,scale=TRUE) # extract scores pca_trajectories = data.frame(code=data_for_pca$code, rmin=data_for_pca$rmin, pca_all$x[,c(1,2,3)]) pca_scaleup = pca_trajectories %>% group_by(rmin) %>% select(PC1, PC2, PC3) %>% summarize_if(is.numeric, mean) %>% gather(key,value,-rmin) g_pca_scaleup = ggplot(pca_scaleup, aes(x=rmin,y=value,col=key)) + geom_line(size=1) + theme_bw() + xlab(expression(paste("r"["min"], " (mm)"))) + geom_hline(yintercept = 0,linetype='dashed',color='black') + ylab("Principal component score") + scale_color_brewer(palette="Set1",name="") ggsave(g_pca_scaleup, file='g_FIGS6.png',width=5,height=4) # variance explained varexp = 100*pca_all$sdev^2/sum(pca_all$sdev^2) # get axes pca_rotations = as.data.frame(pca_all$rotation[,1:2]) pca_rotations$var = row.names(pca_rotations) pca_rotations$x0 = 0 pca_rotations$y0 = 0 row.names(pca_rotations) = NULL scale_factor = 3 # overall PCA space g_pca_biplot = ggplot(pca_trajectories %>% mutate(rmin.factor=rmin), aes(x=PC1, y=PC2,col=rmin.factor,group=rmin.factor)) + geom_point(alpha=0.1) + geom_convexhull(aes(fill=rmin.factor),alpha=0.1,color=NA) + stat_ellipse() + geom_segment(data=pca_rotations,aes(x=x0,y=y0,xend=scale_factor*PC1,yend=scale_factor*PC2),inherit.aes = FALSE,size=1.5,color='black') + geom_text(data=pca_rotations,aes(x=scale_factor*PC1,y=scale_factor*PC2,label=gsub("\\.median","",var)),inherit.aes = FALSE,size=6,color='white') + theme_bw() + xlab(sprintf("PC1 (%.2f%%)",varexp[1])) + ylab(sprintf("PC2 (%.2f%%)",varexp[2])) + scale_fill_viridis(name=expression(paste("r"["min"], " (mm)"))) + scale_color_viridis(guide=FALSE) ggsave(g_pca_biplot, file='g_FIG4.png',width=6,height=5) # add in trait data data_median_with_traits = data_median %>% left_join(data_traits,by='code') %>% dplyr::filter(rmin > 0) # keep out the smallest class data_median_linear_fill_with_traits = data_median_linear_fill %>% left_join(data_traits,by='code') %>% dplyr::filter(rmin > 0) # keep out the smallest class do_pls <- function(data, resp_var_pls) { vein_vars_final = paste(vein_vars,"median",sep=".") # select columns of interest data_pls = data %>% ungroup %>% select(code, rmin, all_of(vein_vars_final), resp_var_pls) %>% mutate_at(resp_var_pls,as.numeric) # don't know how to scale vein columns in dplyr yet, this is a workaround data_pls[,vein_vars_final] = scale(data_pls[,vein_vars_final]) # spread the data data_pls = data_pls %>% gather(key,value,-code,-resp_var_pls,-rmin) %>% mutate(scale=paste(key, round(rmin,digits=4))) %>% select(-rmin, -key) %>% spread(scale, value) m_pls = plsr(formula(sprintf("%s~.",resp_var_pls)),data=data_pls %>% select(-code),validation='CV') cat('.') return(m_pls) } extract_pls_plot <- function(m_pls) { resp_var_pls = dimnames(m_pls$coefficients)[[2]] lvals = as.data.frame(loadings(m_pls)[,1,drop=FALSE]) # pick component one names(lvals) = make.names(names(lvals)) lvals$varscale=row.names(lvals) row.names(lvals) = NULL lvals$var = sapply(strsplit(gsub('`','',lvals$varscale),' '),function(x){x[1]}) lvals$rmin = sapply(strsplit(gsub('`','',lvals$varscale),' '),function(x){as.numeric(x[2])}) lvals = lvals %>% select(-varscale) lvals$r2 = R2(m_pls)$val[,,"1 comps"] # 2nd entry = first PLS component lvals$resp = resp_var_pls return(lvals) } # generate PLS models pls_models = lapply(trait_vars, do_pls, data=data_median_linear_fill_with_traits) # make plots of each pls_plot_data = rbindlist(lapply(pls_models, extract_pls_plot)) %>% mutate(var=gsub("\\.median","",var)) # explained variation pls_plot_data %>% select(r2, resp) %>% unique g_pls = ggplot(pls_plot_data, aes(x=rmin,y=Comp.1,col=var)) + geom_hline(yintercept = 0) + geom_line(aes(size=r2)) + facet_wrap(~resp,labeller=as_labeller(c(Amax_mean="(a) Photosynthetic capacity",LMA_g.cm2_mean="(c) Leaf mass per area",Fp_N_mm_mean="(b) Force to punch",Total_tannin_mg.g="(d) Tannins"))) + theme_bw() + xlab(expression(paste("r"["min"], " (mm)"))) + ylab("Loading coefficient") + scale_color_manual(values=wes_palette("Darjeeling1",5)[c(1,2,4,5)],name="Response variable") + scale_size_continuous(name=expression("R"^2),range=c(0.5,2)) ggsave(g_pls, file='g_FIG5.png',width=7,height=7) pls_models_other = lapply(trait_vars_other, do_pls, data=data_median_linear_fill_with_traits) pls_plot_data_other = rbindlist(lapply(pls_models_other, extract_pls_plot)) %>% mutate(var=gsub("\\.median","",var)) # r2 pls_plot_data_other %>% select(r2, resp) %>% unique # plots g_pls_other = ggplot(pls_plot_data_other, aes(x=rmin,y=Comp.1,col=var)) + geom_hline(yintercept = 0) + geom_line(aes(size=r2)) + facet_wrap(~resp, labeller=as_labeller(c(total_K_mg.g = "K fraction (mg/g)", total_Ca_mg.g = "Ca fraction (mg/g)", total_Mg_mg.g = "Mg fraction (mg/g)", total_P_mg.g = "P fraction (mg/g)", X15N_per_mil = "d15N (per mil)", C_perc = "C fraction (%)", X13C_per_mil = "d13C (per mil)", N_perc = "N fraction (%)", leaf_thickness_mm_mean = "Thickness (mm)", LDMC_mg.g_mean = "LDMC (mg/g)", chla_mg.g = "Chlorophyll A (mg/g)", chlb_mg.g = "Chlorophyll B (mg/g)", carot_mg.g = "Carotenes (mg/g)", hemicellulose_perc = "Hemicellulose (%)", cellulose_perc = "Cellulose (%)", lignin_recalcitrants_perc = "Lignins (%)", Total_phenol_mg.g = "Phenols (mg/g)"))) + theme_bw() + xlab(expression(paste("r"["min"], " (mm)"))) + ylab("Loading coefficient") + scale_color_manual(values=wes_palette("Darjeeling1",5)[c(1,2,4,5)],name="Response variable") + scale_size_continuous(name=expression("R"^2),range=c(0.5,2)) ggsave(g_pls_other, file='g_FIGS9.png',width=10,height=10) # do simulated data for the scale fusion data_simulated_small = data.frame(rmin=seq(0,15,length.out=100)) %>% mutate(Method = "A") %>% mutate(stat = sin(rmin/10)*sqrt(rmin) + rnorm(n=100,sd=0.3)) data_simulated_large = data.frame(rmin=seq(9,30,length.out=100)) %>% mutate(Method = "B") %>% mutate(stat = sin(rmin/10)*sqrt(rmin) + rnorm(n=100,sd=0.5)) data_simulated = rbind(data_simulated_small, data_simulated_large) g_simulated = ggplot(data_simulated, aes(x=rmin,y=stat,group=Method, col=Method)) + geom_line() + theme_bw() + xlab(expression(paste("r"["min"]))) + ylab("Statistic") + geom_smooth(data=data_simulated,aes(x=rmin,y=stat),col='black',inherit.aes = FALSE) ggsave(g_simulated,file='g_FIGS8.png',width=7,height=5) # do phylogenetic analysis # calculate species mean vein trait values at each scale # using the gap filled data (consider not using the gap fill) data_phylo = data_median_linear_fill %>% ungroup %>% group_by(species, rmin) %>% summarize_if(is.numeric, mean,na.rm=T) %>% mutate(taxon=gsub(" ","_",species)) phylosignal_analysis <- function(data, column_name, phylo) { phylo_pruned = keep.tip(phylo, unique(data$taxon)) # get traits data_small = data[,column_name,drop=TRUE] names(data_small) = data$taxon # reorder traits data_small = data_small[phylo_pruned$tip.label] print(table(names(data_small) %in% phylo$tip.label)) k = phylosig(tree=phylo_pruned, x=data_small,method='K',test=TRUE) df_result = data.frame(K=k$K, P=k$P, var=column_name) return(df_result) } # calculate K at each scale for each trait k_results = lapply(paste(vein_vars,"median",sep="."), function(vein_var) { k_result = data_phylo %>% ungroup %>% group_by(rmin) %>% do(phylosignal_analysis(., vein_var, phylo)) }) %>% rbindlist # plot K vs scale g_phylo = ggplot(k_results %>% mutate(var=gsub("\\.median","",var)), aes(x=rmin,y=K,group=var,col=var)) + geom_point(alpha=ifelse(k_results$P<0.05,1,0.1)) + geom_line(alpha=0.1) + #ylim(0,2) + geom_hline(yintercept = 1,linetype=2) + theme_bw() + xlab(expression(paste("r"["min"], " (mm)"))) + ylab("Blomberg's K") + scale_color_manual(values=wes_palette("Darjeeling1",5)[c(1,2,4,5)],name="Variable") ggsave(g_phylo, file='g_FIG7.png',width=5,height=4) plot_trait_matrix <- function(trait) { phylo_pruned = keep.tip(phylo, gsub(" ", "_", unique(data_joined$species))) trait_matrix = data_phylo %>% ungroup %>% select("taxon", all_of(trait), "rmin") %>% dplyr::filter(rmin > 0) %>% mutate(rmin=sprintf("%.3f",rmin)) %>% spread(rmin, trait) %>% as.data.frame row.names(trait_matrix) = trait_matrix$taxon trait_matrix = trait_matrix %>% select(-taxon) # reorder trait_matrix = trait_matrix[phylo_pruned$tip.label,] # pretty it up for plotting phylo_pruned$tip.label = gsub("_", " ", phylo_pruned$tip.label) row.names(trait_matrix) = gsub("_", " ", row.names(trait_matrix)) t = ggtree(phylo_pruned) + geom_tiplab(size=1) + geom_hilight(node=372, fill="blue", alpha=0.5) + # dipterocarps geom_hilight(node=348, fill="orange", alpha=0.5) + # euphorbs geom_hilight(node=329, fill="green", alpha=0.5) + # clusiaceae / hyperiaceae geom_hilight(node=314, fill="red", alpha=0.5) # fagaceae g = gheatmap(t, trait_matrix, offset=40, colnames_angle=0, font.size=0.5, hjust=0.5) + scale_fill_viridis(name=gsub("\\.median","",trait)) + theme(legend.position='bottom',legend.justification = c(0.9, 0)) #message('note this is sqrt-ing the data') return(g) } phylo_plot_cr = plot_trait_matrix("CR.median") phylo_plot_er = plot_trait_matrix("ER.median") phylo_plot_vd = plot_trait_matrix("VD.median") phylo_plot_mst = plot_trait_matrix("MST.median") #g_phylo_heatmap = ggarrange(phylo_plot_cr, phylo_plot_er, phylo_plot_vd, phylo_plot_mst, # nrow=2,ncol=2, # align='hv', # labels=c("(a)","(b)","(c)","(d)")) ggsave(phylo_plot_cr, file='g_FIG6a.pdf',width=5,height=9) ggsave(phylo_plot_er, file='g_FIG6b.pdf',width=5,height=9) ggsave(phylo_plot_vd, file='g_FIG6c.pdf',width=5,height=9) ggsave(phylo_plot_mst, file='g_FIG6d.pdf',width=5,height=9) report_lm <- function(xvar, yvar, data) { if (nrow(data) > 3) { m_lm = NULL try(m_lm <- lm(formula(sprintf("%s~%s",yvar,xvar)), data)) if(!is.null(m_lm)) { f <- summary(m_lm)$fstatistic if(!is.null(f)) { p.val <- pf(f[1],f[2],f[3],lower.tail=F) } else { p.val <- NA } R.squared = summary(m_lm)$r.squared slope.est = coef(m_lm)[2] slope.lo = confint(m_lm)[2,1] slope.hi = confint(m_lm)[2,2] n.points = nrow(data) } else { p.val = NA R.squared = NA slope.est = NA slope.lo = NA slope.hi = NA n.points = NA } } else { p.val = NA R.squared = NA slope.est = NA slope.lo = NA slope.hi = NA n.points = NA } return(data.frame(p.val,R.squared,slope.est,slope.lo,slope.hi,n.points)) } bin_lm_data = function(xv.name, yv.name) { data_summary = data_joined %>% group_by(code, width_threshold_binned) %>% summarize(yvar=mean((!!sym(yv.name)),na.rm=T),xvar=mean((!!sym(xv.name)),na.rm=T)) %>% ungroup %>% group_by(width_threshold_binned) %>% arrange(width_threshold_binned) %>% group_by(width_threshold_binned) %>% do(report_lm(xvar="xvar",yvar="yvar",data=.)) %>% mutate(slope.est.cut = ifelse(p.val < 0.05,slope.est,NA)) %>% mutate(xvar=xv.name,yvar=yv.name) return(data_summary) } params_scaling = expand.grid(xvar=vein_vars, yvar=trait_vars, stringsAsFactors=FALSE) results_scaling = vector(mode="list",length=nrow(params_scaling)) for (i in 1:nrow(params_scaling)) { print(params_scaling[i,]) results_scaling[[i]] = bin_lm_data(xv.name=params_scaling$xvar[i],yv.name=params_scaling$yvar[i]) } results_scaling = rbindlist(results_scaling) results_scaling_plot = results_scaling %>% mutate(xvar.pretty=paste("",xvar,sep=""),yvar.pretty=paste("",yvar,sep="")) %>% # x=, y= dplyr::filter(n.points>100) # note we only keep g_spatial_scaling = ggplot(results_scaling_plot, aes(x=width_threshold_binned,y=slope.est,size=R.squared)) + #,ymin=slope.lo,ymax=slope.hi)) + #geom_ribbon(alpha=0.2) + geom_hline(yintercept = 0) + geom_point(aes(col=(p.val<0.05 & R.squared > 0.05))) + theme_bw() + facet_wrap(yvar.pretty~xvar.pretty,scales='free',ncol=length(unique(params_scaling$xvar)), labeller=as_labeller(c(Amax_mean="Photosynthetic capacity", LMA_g.cm2_mean="Leaf mass per area", Fp_N_mm_mean="Force to punch", Total_tannin_mg.g="Tannins", VD="VD", ER="ER", CR="CR", MST="MST"))) + scale_color_manual(values=c("lightgray","red"),name='P<0.05') + scale_size(name=expression('R'^2),range=c(0.5,3)) + ylab("Slope estimate") + xlab(expression(paste("r"["min"], " (mm)"))) ggsave(g_spatial_scaling,file='g_FIGS7.png',width=2.1*length(unique(params_scaling$xvar)),height=2*length(unique(params_scaling$yvar))) # show the panels panel_plot <- function(data, y_stat, y_label, do_log=TRUE) { g_all = ggplot(data %>% mutate(fullname=paste(Family,species,sep="\n")), aes_string(x="rmin",y=y_stat,group="code",col="code")) + geom_line(alpha=0.5) + geom_point(size=0.25) + #geom_line(data=data, # aes_string(x="rmin",y=y_stat,group="code",col='species'),alpha=0.5,size=1) + #geom_line(data=data, # aes_string(x="rmin",y=y_stat,group="code",col='species'),size=1) + #geom_point(data=data, # aes_string(x="rmin",y=y_stat,group="code",col='species'),size=2) + theme_bw() + theme(legend.position='none') + #theme(legend.text=element_text(face="italic")) + #scale_color_brewer(palette="Set1",name="Species") + ylab(y_label) + xlab(expression(paste("r"["min"], " (mm)"))) + facet_wrap(~fullname,ncol=13) if (do_log==TRUE) { g_all = g_all + scale_y_log10() } return(g_all) } g_er_panels = panel_plot(data=data_median,y_stat="ER.median",y_label="Elongation ratio (ER)") g_cr_panels = panel_plot(data=data_median,y_stat="CR.median",y_label="Circularity ratio (ER)") g_vd_panels = panel_plot(data=data_median,y_stat="VD.median",y_label=expression(paste("Vein density (VD, mm"^{-1},")"))) g_mst_panels = panel_plot(data=data_median,y_stat="MST.median",y_label="Minimum spanning tree ratio (MST)") ggsave(g_er_panels, file='g_FIGS2.png',width=23,height=32) ggsave(g_cr_panels, file='g_FIGS3.png',width=23,height=32) ggsave(g_vd_panels, file='g_FIGS4.png',width=23,height=32) ggsave(g_mst_panels, file='g_FIGS5.png',width=23,height=32)