Genetic correlations
Anoob Prakash
02/02/2021
Methods
lmer model
- Mixed effect model run to get the family level blups for use in the MCMC run.
- Garden is set as the fixed effect.
- Family and mBed (5 beds x 3 sites) set as random effect.
trait1 <- trait1 %>% unite(mBed, Garden, Bed, sep = "_", remove = FALSE)
trait1 <- lmer(Growth ~ Garden + (1|mBed) + (1|Family), trait1)
trait1 <- ranef(trait1)
trait1 <- trait1$Family
trait1 <- cbind(Family=rownames(trait1),trait1)
rownames(trait1) <- NULL
names(trait1)[2] <- "Growth"
trait1$Population <- gsub("\\_.*", "", trait1$Family)MCMC model - The scaled family blups of the traits are used as input MCMC run.
- Random effect of Population accounted for in this model
trait1_trait2_fam_blup_pop_corr <- MCMCglmm(cbind(scale(trait1),scale(trait2)) ~ trait - 1,
random=~us(trait):Population,
rcov=~us(trait):units, # or idh
family=c("gaussian","gaussian"),
prior=prior_corr,
# pedigree=Ped,
data=bs_19_gw_19,
pr=TRUE,
verbose=TRUE,
nitt=10000000,
burnin=10000,
thin=1000)notes
gw - growth
bb - budbreak
bs - budset
B - lmer family blups (means)
P - plasticity value for the trait
- notes from: https://stat.ethz.ch/pipermail/r-sig-mixed-models/2017q1/025532.html
- multiresponse mcmc models, page 92: https://cran.r-project.org/web/packages/MCMCglmm/vignettes/CourseNotes.pdf
packages
require(MCMCglmm)## Loading required package: MCMCglmm## Loading required package: Matrix## Loading required package: coda## Loading required package: aperequire(lmerTest)## Loading required package: lmerTest## Loading required package: lme4##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## steprequire(dplyr)## Loading required package: dplyr## Warning: package 'dplyr' was built under R version 3.6.3##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionrequire(tidyr)## Loading required package: tidyr##
## Attaching package: 'tidyr'## The following objects are masked from 'package:Matrix':
##
## expand, pack, unpackrequire(plotly)## Loading required package: plotly## Loading required package: ggplot2##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutrequire(plyr)## Loading required package: plyr## Warning: package 'plyr' was built under R version 3.6.3## ------------------------------------------------------------------------------## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)## ------------------------------------------------------------------------------##
## Attaching package: 'plyr'## The following objects are masked from 'package:plotly':
##
## arrange, mutate, rename, summarise## The following objects are masked from 'package:dplyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarizerequire(viridis)## Loading required package: viridis## Loading required package: viridisLite# require(grid)
# require(gridExtra)
# setwd("Z:/Anoob/MCMCglmm")data
# exome data
exome_meta <- read.table("./Exome/RS_Exome_metadata.txt", sep="\t", header = T)
# climateNA
climateNA <- read.csv("./ClimateNA/Family/PCA_sprucePops.txt",
header=T, sep="\t")
colnames(climateNA)[colnames(climateNA)=="PC1_sel"] <- "PC1"
Population_climateNA <- read.csv("./ClimateNA/Family/PCA_sprucePops.txt",
header=T, sep="\t")
colnames(Population_climateNA)[colnames(Population_climateNA)=="PC1_sel"] <- "PC1"
Family_climateNA <- read.csv("./ClimateNA/Family/PCA_spruceFams.txt",
header=T, sep="\t")
colnames(Family_climateNA)[colnames(Family_climateNA)=="PC1_sel"] <- "PC1"
# data
budset_19 <- read.csv("./trait_data/BudSet_2019.csv")
budset_20 <- read.csv("./trait_data/BudSet_2020.csv")
budbreak_20 <- read.csv("./trait_data/BudBreak_2020_cGDD.csv")
heightGrowth_19 <- read.csv("./trait_data/Growth_2019.csv")
heightGrowth_20 <- read.csv("./trait_data/Growth_2020.csv")Family blups
## family blups
# blups for budset 2019
budset_19 <- budset_19 %>% unite(mBed, Garden, Bed, sep = "_", remove = FALSE)
budset_19 <- lmer(BudSet ~ Garden + (1|mBed) + (1|Family), budset_19)
budset_19 <- ranef(budset_19)
budset_19 <- budset_19$Family
budset_19 <- cbind(Family=rownames(budset_19),budset_19)
rownames(budset_19) <- NULL
names(budset_19)[2] <- "BudSet"
budset_19$Population <- gsub("\\_.*", "", budset_19$Family)
# blups for budset 2020
budset_20 <- budset_20 %>% unite(mBed, Garden, Bed, sep = "_", remove = FALSE)
budset_20 <- lmer(BudSet ~ Garden + (1|mBed) + (1|Family), budset_20)
budset_20 <- ranef(budset_20)
budset_20 <- budset_20$Family
budset_20 <- cbind(Family=rownames(budset_20),budset_20)
rownames(budset_20) <- NULL
names(budset_20)[2] <- "BudSet"
budset_20$Population <- gsub("\\_.*", "", budset_20$Family)
# blups for budbreak 2020
budbreak_20 <- budbreak_20 %>% unite(mBed, Garden, Bed, sep = "_", remove = FALSE)
budbreak_20 <- lmer(cGDD ~ Garden + (1|mBed) + (1|Family), budbreak_20, na.action = na.exclude)
budbreak_20 <- ranef(budbreak_20)
budbreak_20 <- budbreak_20$Family
budbreak_20 <- cbind(Family=rownames(budbreak_20),budbreak_20)
rownames(budbreak_20) <- NULL
names(budbreak_20)[2] <- "cGDD"
budbreak_20$Population <- gsub("\\_.*", "", budbreak_20$Family)
# blups for height growth 2019
heightGrowth_19 <- heightGrowth_19 %>% unite(mBed, Garden, Bed, sep = "_", remove = FALSE)
heightGrowth_19 <- lmer(Growth ~ Garden + (1|mBed) + (1|Family), heightGrowth_19)
heightGrowth_19 <- ranef(heightGrowth_19)
heightGrowth_19 <- heightGrowth_19$Family
heightGrowth_19 <- cbind(Family=rownames(heightGrowth_19),heightGrowth_19)
rownames(heightGrowth_19) <- NULL
names(heightGrowth_19)[2] <- "Growth"
heightGrowth_19$Population <- gsub("\\_.*", "", heightGrowth_19$Family)
# blups for height growth 2020
heightGrowth_20 <- heightGrowth_20 %>% unite(mBed, Garden, Bed, sep = "_", remove = FALSE)
heightGrowth_20 <- lmer(Growth ~ Garden + (1|mBed) + (1|Family), heightGrowth_20)
heightGrowth_20 <- ranef(heightGrowth_20)
heightGrowth_20 <- heightGrowth_20$Family
heightGrowth_20 <- cbind(Family=rownames(heightGrowth_20),heightGrowth_20)
rownames(heightGrowth_20) <- NULL
names(heightGrowth_20)[2] <- "Growth"
heightGrowth_20$Population <- gsub("\\_.*", "", heightGrowth_20$Family)Phenology vs growth
Budset 2019-Growth 2019
Summary
# against current year
bs_19_gw_19_fam_blup_pop_corr <- readRDS("./genetic_correlation_outputs/bs_19_gw_19_fam_blup_pop_corr")
summary(bs_19_gw_19_fam_blup_pop_corr)##
## Iterations = 10001:9999001
## Thinning interval = 1000
## Sample size = 9990
##
## DIC: 1544.499
##
## G-structure: ~us(trait):Population
##
## post.mean l-95% CI u-95% CI eff.samp
## traitBudSet:traitBudSet.Population 0.4038 0.2241 0.5926 9990
## traitGrowth:traitBudSet.Population 0.2674 0.1277 0.4291 9990
## traitBudSet:traitGrowth.Population 0.2674 0.1277 0.4291 9990
## traitGrowth:traitGrowth.Population 0.4231 0.2491 0.6158 9935
##
## R-structure: ~us(trait):units
##
## post.mean l-95% CI u-95% CI eff.samp
## traitBudSet:traitBudSet.units 0.6212 0.5210 0.7286 9990
## traitGrowth:traitBudSet.units 0.3212 0.2451 0.4087 9990
## traitBudSet:traitGrowth.units 0.3212 0.2451 0.4087 9990
## traitGrowth:traitGrowth.units 0.5878 0.4897 0.6868 9990
##
## Location effects: cbind(scale(BudSet), scale(Growth)) ~ trait - 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## traitBudSet -0.001759 -0.183235 0.169447 9990 0.976
## traitGrowth -0.010651 -0.196784 0.163363 9990 0.912effectiveSize(bs_19_gw_19_fam_blup_pop_corr$VCV)## traitBudSet:traitBudSet.Population traitGrowth:traitBudSet.Population
## 9990.000 9990.000
## traitBudSet:traitGrowth.Population traitGrowth:traitGrowth.Population
## 9990.000 9934.643
## traitBudSet:traitBudSet.units traitGrowth:traitBudSet.units
## 9990.000 9990.000
## traitBudSet:traitGrowth.units traitGrowth:traitGrowth.units
## 9990.000 9990.000heidel.diag(bs_19_gw_19_fam_blup_pop_corr$VCV)##
## Stationarity start p-value
## test iteration
## traitBudSet:traitBudSet.Population passed 1 0.598
## traitGrowth:traitBudSet.Population passed 1 0.292
## traitBudSet:traitGrowth.Population passed 1 0.292
## traitGrowth:traitGrowth.Population passed 1 0.349
## traitBudSet:traitBudSet.units passed 1 0.790
## traitGrowth:traitBudSet.units passed 1 0.455
## traitBudSet:traitGrowth.units passed 1 0.455
## traitGrowth:traitGrowth.units passed 1 0.360
##
## Halfwidth Mean Halfwidth
## test
## traitBudSet:traitBudSet.Population passed 0.404 0.001892
## traitGrowth:traitBudSet.Population passed 0.267 0.001552
## traitBudSet:traitGrowth.Population passed 0.267 0.001552
## traitGrowth:traitGrowth.Population passed 0.423 0.001916
## traitBudSet:traitBudSet.units passed 0.621 0.001051
## traitGrowth:traitBudSet.units passed 0.321 0.000819
## traitBudSet:traitGrowth.units passed 0.321 0.000819
## traitGrowth:traitGrowth.units passed 0.588 0.000990head(bs_19_gw_19_fam_blup_pop_corr$VCV)## Markov Chain Monte Carlo (MCMC) output:
## Start = 10001
## End = 16001
## Thinning interval = 1000
## traitBudSet:traitBudSet.Population traitGrowth:traitBudSet.Population
## [1,] 0.3818943 0.2691710
## [2,] 0.2535308 0.1383171
## [3,] 0.6081488 0.3111277
## [4,] 0.3086442 0.2216499
## [5,] 0.3018752 0.2448904
## [6,] 0.2829332 0.1732858
## [7,] 0.4019098 0.2670035
## traitBudSet:traitGrowth.Population traitGrowth:traitGrowth.Population
## [1,] 0.2691710 0.4377213
## [2,] 0.1383171 0.2967979
## [3,] 0.3111277 0.3483605
## [4,] 0.2216499 0.3232103
## [5,] 0.2448904 0.4342290
## [6,] 0.1732858 0.4075192
## [7,] 0.2670035 0.4171464
## traitBudSet:traitBudSet.units traitGrowth:traitBudSet.units
## [1,] 0.5828179 0.3045284
## [2,] 0.6180179 0.3037584
## [3,] 0.5549477 0.3140133
## [4,] 0.5926414 0.2614456
## [5,] 0.6009730 0.3302068
## [6,] 0.6872240 0.4168780
## [7,] 0.6050895 0.2973151
## traitBudSet:traitGrowth.units traitGrowth:traitGrowth.units
## [1,] 0.3045284 0.5716309
## [2,] 0.3037584 0.5754469
## [3,] 0.3140133 0.6270778
## [4,] 0.2614456 0.4938555
## [5,] 0.3302068 0.5821787
## [6,] 0.4168780 0.6181461
## [7,] 0.2973151 0.5936725# genetic correlation
bs_19_gw_19_pop_gencorr <- bs_19_gw_19_fam_blup_pop_corr$VCV[,"traitBudSet:traitGrowth.Population"]/sqrt(bs_19_gw_19_fam_blup_pop_corr$VCV[,"traitBudSet:traitBudSet.Population"] * bs_19_gw_19_fam_blup_pop_corr$VCV[,"traitGrowth:traitGrowth.Population"])
posterior.mode(bs_19_gw_19_pop_gencorr)## var1
## 0.6704044HPDinterval(bs_19_gw_19_pop_gencorr)## lower upper
## var1 0.450486 0.8210264
## attr(,"Probability")
## [1] 0.9499499Data
MCMC
## setting up dataframe for correlation based on the model blups
bs_19_gw_19_coefs <- data_frame(Trait = attr(colMeans(bs_19_gw_19_fam_blup_pop_corr$Sol), "names"),
Value = colMeans(bs_19_gw_19_fam_blup_pop_corr$Sol)) %>%
separate(Trait, c("Trait","Type","Population"), sep = "\\.", fill = "right") %>%
filter(Type == "Population") %>%
filter(Trait %in% c("traitGrowth", "traitBudSet")) %>%
select(-Type) %>%
spread(Trait, Value)## Warning: `data_frame()` is deprecated as of tibble 1.1.0.
## Please use `tibble()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.bs_19_gw_19_coefs <- merge(x=bs_19_gw_19_coefs,
y=exome_meta[,c("Pop","Region")],
all.x=T,
by.x="Population",
by.y="Pop")
# rename regions
bs_19_gw_19_coefs$Region <- revalue(bs_19_gw_19_coefs$Region, c("C"="Core", "M"="Margin","E"="Edge"))
# reorder regions
bs_19_gw_19_coefs$Region <- factor(bs_19_gw_19_coefs$Region,levels = c("Core", "Margin", "Edge"))
# merge climateNA data
bs_19_gw_19_coefs <- merge(x=bs_19_gw_19_coefs,
y=Population_climateNA[,c("Population","PC1","MAT")],
by="Population")
cor.test(bs_19_gw_19_coefs$traitBudSet,bs_19_gw_19_coefs$traitGrowth) #0.6867328##
## Pearson's product-moment correlation
##
## data: bs_19_gw_19_coefs$traitBudSet and bs_19_gw_19_coefs$traitGrowth
## t = 17.369, df = 338, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6261026 0.7391087
## sample estimates:
## cor
## 0.6867328Family
bs_19_gw_19 <- merge(budset_19,heightGrowth_19[,c("Family","Growth")])
bs_19_gw_19 <- bs_19_gw_19[!is.na(bs_19_gw_19$BudSet),]
bs_19_gw_19 <- bs_19_gw_19[!is.na(bs_19_gw_19$Growth),]
# merge Region data
bs_19_gw_19 <- merge(x=bs_19_gw_19,
y=exome_meta[,c("Pop","Region")],
by.x="Population",
by.y="Pop")
# rename and reorder regions
bs_19_gw_19$Region <- plyr::revalue(bs_19_gw_19$Region,
c("C" = "Core",
"M" = "Margin",
"E" = "Edge"))
bs_19_gw_19$Region <- factor(bs_19_gw_19$Region,levels = c("Core", "Margin", "Edge"))
# merge climateNA data
bs_19_gw_19 <- merge(x=bs_19_gw_19,
y=Family_climateNA[,c("Population","PC1","MAT")],
by="Population")
bs_19_gw_19 <- distinct(bs_19_gw_19)Population
bs_19_gw_19$Population <- as.factor(bs_19_gw_19$Population)
bs_19_gw_19_pop <- bs_19_gw_19 %>%
select(Population, BudSet, Growth, Region) %>%
group_by(Population) %>%
dplyr::summarise(BudSet = mean(BudSet), Growth = mean(Growth))
bs_19_gw_19_pop <- merge(x=bs_19_gw_19_pop,
y=Population_climateNA[,c("Population","PC1","MAT")],
by="Population")
bs_19_gw_19_pop <- merge(x=bs_19_gw_19_pop,
y=exome_meta[,c("Pop","Region")],
by.x="Population",
by.y="Pop")
# rename and reorder regions
bs_19_gw_19_pop$Region <- plyr::revalue(bs_19_gw_19_pop$Region,
c("C" = "Core",
"M" = "Margin",
"E" = "Edge"))
bs_19_gw_19_pop$Region <- factor(bs_19_gw_19_pop$Region,levels = c("Core", "Margin", "Edge"))
bs_19_gw_19_pop <- distinct(bs_19_gw_19_pop)MCMC plot
# plot
ggplot(bs_19_gw_19_coefs, aes(x = traitBudSet, y = traitGrowth, group = Population, color = PC1)) +
geom_point(aes(shape = Region), size=3, alpha = 0.7) +
scale_color_viridis(option = "A", direction = -1) +
# geom_abline(intercept = 0, slope = mean(gw_19_B_gw_19_P_fit_slope)) +
labs(x = "Budset in 2019", y = "Height Growth in 2019") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) #+
# annotate("text", label = "Correlation between traits: 0.688", x = 0.7, y = -1.2)Family plot
# plot
ggplot(bs_19_gw_19, aes(x = BudSet, y = Growth, group = Population, color = PC1)) +
geom_point(aes(shape = Region), size=3, alpha = 0.7) +
scale_color_viridis(option = "A", direction = -1) +
# geom_abline(intercept = 0, slope = mean(gw_19_B_gw_19_P_fit_slope)) +
labs(x = "Budset in 2019", y = "Height Growth in 2019") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) #+
# annotate("text", label = "Correlation between traits: 0.688", x = 0.7, y = -1.2)Population plot
# plot
ggplot(bs_19_gw_19_pop, aes(x = BudSet, y = Growth, group = Population, color = PC1)) +
geom_point(aes(shape = Region), size=3, alpha = 0.7) +
scale_color_viridis(option = "A", direction = -1) +
# geom_abline(intercept = 0, slope = mean(gw_19_B_gw_19_P_fit_slope)) +
labs(x = "Budset in 2019", y = "Height Growth in 2019") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) #+
# annotate("text", label = "Correlation between traits: 0.688", x = 0.7, y = -1.2)Budset 2020-Growth 2020
Summary
# against current year
bs_20_gw_20_fam_blup_pop_corr <- readRDS("./genetic_correlation_outputs/bs_20_gw_20_fam_blup_pop_corr")
summary(bs_20_gw_20_fam_blup_pop_corr)##
## Iterations = 10001:9999001
## Thinning interval = 1000
## Sample size = 9990
##
## DIC: 1584.176
##
## G-structure: ~us(trait):Population
##
## post.mean l-95% CI u-95% CI eff.samp
## traitBudSet:traitBudSet.Population 0.1708 0.06742 0.2904 9990
## traitGrowth:traitBudSet.Population 0.2404 0.11839 0.3782 10389
## traitBudSet:traitGrowth.Population 0.2404 0.11839 0.3782 10389
## traitGrowth:traitGrowth.Population 0.5416 0.33543 0.7642 9990
##
## R-structure: ~us(trait):units
##
## post.mean l-95% CI u-95% CI eff.samp
## traitBudSet:traitBudSet.units 0.8370 0.6977 0.9785 9818
## traitGrowth:traitBudSet.units 0.2249 0.1477 0.3042 10359
## traitBudSet:traitGrowth.units 0.2249 0.1477 0.3042 10359
## traitGrowth:traitGrowth.units 0.4362 0.3627 0.5097 10440
##
## Location effects: cbind(scale(BudSet), scale(Growth)) ~ trait - 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## traitBudSet -0.002690 -0.139495 0.139808 10575 0.968
## traitGrowth -0.003049 -0.194897 0.186240 10198 0.987effectiveSize(bs_20_gw_20_fam_blup_pop_corr$VCV)## traitBudSet:traitBudSet.Population traitGrowth:traitBudSet.Population
## 9990.000 10388.596
## traitBudSet:traitGrowth.Population traitGrowth:traitGrowth.Population
## 10388.596 9990.000
## traitBudSet:traitBudSet.units traitGrowth:traitBudSet.units
## 9818.434 10359.171
## traitBudSet:traitGrowth.units traitGrowth:traitGrowth.units
## 10359.171 10439.823heidel.diag(bs_20_gw_20_fam_blup_pop_corr$VCV)##
## Stationarity start p-value
## test iteration
## traitBudSet:traitBudSet.Population passed 1 0.266
## traitGrowth:traitBudSet.Population passed 1 0.660
## traitBudSet:traitGrowth.Population passed 1 0.660
## traitGrowth:traitGrowth.Population passed 1 0.827
## traitBudSet:traitBudSet.units passed 1 0.152
## traitGrowth:traitBudSet.units passed 1 0.154
## traitBudSet:traitGrowth.units passed 1 0.154
## traitGrowth:traitGrowth.units passed 1 0.561
##
## Halfwidth Mean Halfwidth
## test
## traitBudSet:traitBudSet.Population passed 0.171 0.001172
## traitGrowth:traitBudSet.Population passed 0.240 0.001300
## traitBudSet:traitGrowth.Population passed 0.240 0.001300
## traitGrowth:traitGrowth.Population passed 0.542 0.002215
## traitBudSet:traitBudSet.units passed 0.837 0.001418
## traitGrowth:traitBudSet.units passed 0.225 0.000765
## traitBudSet:traitGrowth.units passed 0.225 0.000765
## traitGrowth:traitGrowth.units passed 0.436 0.000728head(bs_20_gw_20_fam_blup_pop_corr$VCV)## Markov Chain Monte Carlo (MCMC) output:
## Start = 10001
## End = 16001
## Thinning interval = 1000
## traitBudSet:traitBudSet.Population traitGrowth:traitBudSet.Population
## [1,] 0.1910824 0.2425390
## [2,] 0.1539166 0.1959441
## [3,] 0.1917259 0.1686148
## [4,] 0.2171373 0.2341403
## [5,] 0.1484663 0.2817071
## [6,] 0.1548874 0.1986629
## [7,] 0.3049565 0.3676478
## traitBudSet:traitGrowth.Population traitGrowth:traitGrowth.Population
## [1,] 0.2425390 0.4199850
## [2,] 0.1959441 0.5634770
## [3,] 0.1686148 0.5054699
## [4,] 0.2341403 0.4639595
## [5,] 0.2817071 0.5953026
## [6,] 0.1986629 0.5326382
## [7,] 0.3676478 0.6754438
## traitBudSet:traitBudSet.units traitGrowth:traitBudSet.units
## [1,] 0.7266947 0.1779280
## [2,] 0.9039916 0.2812020
## [3,] 0.9054705 0.2723551
## [4,] 0.8367770 0.2314403
## [5,] 0.8332138 0.2163979
## [6,] 0.8099028 0.2076021
## [7,] 0.8070284 0.2417213
## traitBudSet:traitGrowth.units traitGrowth:traitGrowth.units
## [1,] 0.1779280 0.4792074
## [2,] 0.2812020 0.4973274
## [3,] 0.2723551 0.4547770
## [4,] 0.2314403 0.4690094
## [5,] 0.2163979 0.4098387
## [6,] 0.2076021 0.4040608
## [7,] 0.2417213 0.4336318# genetic correlation
bs_20_gw_20_pop_gencorr <- bs_20_gw_20_fam_blup_pop_corr$VCV[,"traitBudSet:traitGrowth.Population"]/sqrt(bs_20_gw_20_fam_blup_pop_corr$VCV[,"traitBudSet:traitBudSet.Population"] * bs_20_gw_20_fam_blup_pop_corr$VCV[,"traitGrowth:traitGrowth.Population"])
posterior.mode(bs_20_gw_20_pop_gencorr)## var1
## 0.8518504HPDinterval(bs_20_gw_20_pop_gencorr)## lower upper
## var1 0.6039526 0.9867893
## attr(,"Probability")
## [1] 0.9499499Data
MCMC
## setting up dataframe for correlation based on the model blups
bs_20_gw_20_coefs <- data_frame(Trait = attr(colMeans(bs_20_gw_20_fam_blup_pop_corr$Sol), "names"),
Value = colMeans(bs_20_gw_20_fam_blup_pop_corr$Sol)) %>%
separate(Trait, c("Trait","Type","Population"), sep = "\\.", fill = "right") %>%
filter(Type == "Population") %>%
filter(Trait %in% c("traitGrowth", "traitBudSet")) %>%
select(-Type) %>%
spread(Trait, Value)
bs_20_gw_20_coefs <- merge(x=bs_20_gw_20_coefs,
y=exome_meta[,c("Pop","Region")],
all.x=T,
by.x="Population",
by.y="Pop")
# rename regions
bs_20_gw_20_coefs$Region <- revalue(bs_20_gw_20_coefs$Region, c("C"="Core", "M"="Margin","E"="Edge"))
# reorder regions
bs_20_gw_20_coefs$Region <- factor(bs_20_gw_20_coefs$Region,levels = c("Core", "Margin", "Edge"))
# merge climateNA data
bs_20_gw_20_coefs <- merge(x=bs_20_gw_20_coefs,
y=Population_climateNA[,c("Population","PC1","MAT")],
by="Population")
cor.test(bs_20_gw_20_coefs$traitBudSet,bs_20_gw_20_coefs$traitGrowth) #0.9357684##
## Pearson's product-moment correlation
##
## data: bs_20_gw_20_coefs$traitBudSet and bs_20_gw_20_coefs$traitGrowth
## t = 49.016, df = 338, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.9217366 0.9482352
## sample estimates:
## cor
## 0.9363054Family
# against current year
bs_20_gw_20 <- merge(budset_20,heightGrowth_20[,c("Family","Growth")])
bs_20_gw_20 <- bs_20_gw_20[!is.na(bs_20_gw_20$BudSet),]
bs_20_gw_20 <- bs_20_gw_20[!is.na(bs_20_gw_20$Growth),]
# merge Region data
bs_20_gw_20 <- merge(x=bs_20_gw_20,
y=exome_meta[,c("Pop","Region")],
by.x="Population",
by.y="Pop")
# rename and reorder regions
bs_20_gw_20$Region <- plyr::revalue(bs_20_gw_20$Region,
c("C" = "Core",
"M" = "Margin",
"E" = "Edge"))
bs_20_gw_20$Region <- factor(bs_20_gw_20$Region,levels = c("Core", "Margin", "Edge"))
# merge climateNA data
bs_20_gw_20 <- merge(x=bs_20_gw_20,
y=Family_climateNA[,c("Population","PC1","MAT")],
by="Population")
bs_20_gw_20 <- distinct(bs_20_gw_20)Population
bs_20_gw_20$Population <- as.factor(bs_20_gw_20$Population)
bs_20_gw_20_pop <- bs_20_gw_20 %>%
select(Population, BudSet, Growth, Region) %>%
group_by(Population) %>%
dplyr::summarise(BudSet = mean(BudSet), Growth = mean(Growth))
bs_20_gw_20_pop <- merge(x=bs_20_gw_20_pop,
y=Population_climateNA[,c("Population","PC1","MAT")],
by="Population")
bs_20_gw_20_pop <- merge(x=bs_20_gw_20_pop,
y=exome_meta[,c("Pop","Region")],
by.x="Population",
by.y="Pop")
# rename and reorder regions
bs_20_gw_20_pop$Region <- plyr::revalue(bs_20_gw_20_pop$Region,
c("C" = "Core",
"M" = "Margin",
"E" = "Edge"))
bs_20_gw_20_pop$Region <- factor(bs_20_gw_20_pop$Region,levels = c("Core", "Margin", "Edge"))
bs_20_gw_20_pop <- distinct(bs_20_gw_20_pop)MCMC plot
# plot
bs_20_gw_20_PC1 <- ggplot(bs_20_gw_20_coefs, aes(x = traitBudSet, y = traitGrowth, group = Population, color = PC1)) +
geom_point(aes(shape = Region), size=3, alpha = 0.7) +
scale_color_viridis(option = "A", direction = -1) +
# geom_abline(intercept = 0, slope = mean(gw_19_B_gw_19_P_fit_slope)) +
labs(x = "Budset in 2020", y = "Height Growth in 2020") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) #+
# annotate("text", label = "Correlation between traits: 0.936", x = 1.5, y = -0.9)
bs_20_gw_20_PC1
# ggplotly(bs_20_gw_20_PC1)Family plot
# plot
ggplot(bs_20_gw_20, aes(x = BudSet, y = Growth, group = Population, color = PC1)) +
geom_point(aes(shape = Region), size=3, alpha = 0.7) +
scale_color_viridis(option = "A", direction = -1) +
# geom_abline(intercept = 0, slope = mean(gw_19_B_gw_19_P_fit_slope)) +
labs(x = "Budset in 2020", y = "Height Growth in 2020") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) #+
annotate("text", label = "Correlation between traits: 0.936", x = 1.5, y = -0.9) ## mapping: x = ~x, y = ~y
## geom_text: na.rm = FALSE
## stat_identity: na.rm = FALSE
## position_identity# ggplotly(bs_20_gw_20_PC1)Population plot
# plot
ggplot(bs_20_gw_20_pop, aes(x = BudSet, y = Growth, group = Population, color = PC1)) +
geom_point(aes(shape = Region), size=3, alpha = 0.7) +
scale_color_viridis(option = "A", direction = -1) +
# geom_abline(intercept = 0, slope = mean(gw_19_B_gw_19_P_fit_slope)) +
labs(x = "Budset in 2020", y = "Height Growth in 2020") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13))
Budbreak 2020-Growth 2019
Summary
# against previous year
bb_20_gw_19_fam_blup_pop_corr <- readRDS("./genetic_correlation_outputs/bb_20_gw_19_fam_blup_pop_corr")
summary(bb_20_gw_19_fam_blup_pop_corr)##
## Iterations = 10001:9999001
## Thinning interval = 1000
## Sample size = 9990
##
## DIC: 1644.253
##
## G-structure: ~us(trait):Population
##
## post.mean l-95% CI u-95% CI eff.samp
## traitcGDD:traitcGDD.Population 0.4171 0.23479 0.6151 10865
## traitGrowth:traitcGDD.Population 0.2344 0.09358 0.3875 10553
## traitcGDD:traitGrowth.Population 0.2344 0.09358 0.3875 10553
## traitGrowth:traitGrowth.Population 0.4356 0.25605 0.6363 9990
##
## R-structure: ~us(trait):units
##
## post.mean l-95% CI u-95% CI eff.samp
## traitcGDD:traitcGDD.units 0.62425 0.52030 0.7318 9990
## traitGrowth:traitcGDD.units 0.09891 0.02431 0.1717 9693
## traitcGDD:traitGrowth.units 0.09891 0.02431 0.1717 9693
## traitGrowth:traitGrowth.units 0.58448 0.49360 0.6871 9990
##
## Location effects: cbind(scale(cGDD), scale(Growth)) ~ trait - 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## traitcGDD 0.014841 -0.173056 0.182626 9635 0.867
## traitGrowth -0.007359 -0.182190 0.176073 9621 0.938effectiveSize(bb_20_gw_19_fam_blup_pop_corr$VCV)## traitcGDD:traitcGDD.Population traitGrowth:traitcGDD.Population
## 10865.199 10552.613
## traitcGDD:traitGrowth.Population traitGrowth:traitGrowth.Population
## 10552.613 9990.000
## traitcGDD:traitcGDD.units traitGrowth:traitcGDD.units
## 9990.000 9692.877
## traitcGDD:traitGrowth.units traitGrowth:traitGrowth.units
## 9692.877 9990.000heidel.diag(bb_20_gw_19_fam_blup_pop_corr$VCV)##
## Stationarity start p-value
## test iteration
## traitcGDD:traitcGDD.Population passed 1 0.0955
## traitGrowth:traitcGDD.Population passed 1 0.7712
## traitcGDD:traitGrowth.Population passed 1 0.7712
## traitGrowth:traitGrowth.Population passed 1 0.0748
## traitcGDD:traitcGDD.units passed 1 0.4847
## traitGrowth:traitcGDD.units passed 1 0.4341
## traitcGDD:traitGrowth.units passed 1 0.4341
## traitGrowth:traitGrowth.units passed 1 0.6435
##
## Halfwidth Mean Halfwidth
## test
## traitcGDD:traitcGDD.Population passed 0.4171 0.001914
## traitGrowth:traitcGDD.Population passed 0.2344 0.001463
## traitcGDD:traitGrowth.Population passed 0.2344 0.001463
## traitGrowth:traitGrowth.Population passed 0.4356 0.001979
## traitcGDD:traitcGDD.units passed 0.6242 0.001076
## traitGrowth:traitcGDD.units passed 0.0989 0.000746
## traitcGDD:traitGrowth.units passed 0.0989 0.000746
## traitGrowth:traitGrowth.units passed 0.5845 0.000978head(bb_20_gw_19_fam_blup_pop_corr$VCV)## Markov Chain Monte Carlo (MCMC) output:
## Start = 10001
## End = 16001
## Thinning interval = 1000
## traitcGDD:traitcGDD.Population traitGrowth:traitcGDD.Population
## [1,] 0.4322287 0.2154402
## [2,] 0.3829411 0.1902089
## [3,] 0.4131572 0.1269036
## [4,] 0.4175106 0.2836771
## [5,] 0.5422818 0.3685903
## [6,] 0.3395365 0.1565171
## [7,] 0.4544345 0.2975085
## traitcGDD:traitGrowth.Population traitGrowth:traitGrowth.Population
## [1,] 0.2154402 0.3500787
## [2,] 0.1902089 0.2278045
## [3,] 0.1269036 0.4079173
## [4,] 0.2836771 0.5245097
## [5,] 0.3685903 0.4223203
## [6,] 0.1565171 0.4293298
## [7,] 0.2975085 0.5129945
## traitcGDD:traitcGDD.units traitGrowth:traitcGDD.units
## [1,] 0.6979880 0.13582270
## [2,] 0.5383521 0.07096139
## [3,] 0.6669764 0.11684992
## [4,] 0.6286008 0.13565059
## [5,] 0.5901310 0.07846901
## [6,] 0.6118179 0.11322611
## [7,] 0.6221738 0.11483953
## traitcGDD:traitGrowth.units traitGrowth:traitGrowth.units
## [1,] 0.13582270 0.5433132
## [2,] 0.07096139 0.6001118
## [3,] 0.11684992 0.6387003
## [4,] 0.13565059 0.5347718
## [5,] 0.07846901 0.5762160
## [6,] 0.11322611 0.6187398
## [7,] 0.11483953 0.5438970# genetic correlation
bb_20_gw_19_pop_gencorr <- bb_20_gw_19_fam_blup_pop_corr$VCV[,"traitcGDD:traitGrowth.Population"]/sqrt(bb_20_gw_19_fam_blup_pop_corr$VCV[,"traitcGDD:traitcGDD.Population"] * bb_20_gw_19_fam_blup_pop_corr$VCV[,"traitGrowth:traitGrowth.Population"])
posterior.mode(bb_20_gw_19_pop_gencorr)## var1
## 0.5671157HPDinterval(bb_20_gw_19_pop_gencorr)## lower upper
## var1 0.3167645 0.7710248
## attr(,"Probability")
## [1] 0.9499499Data
MCMC
## setting up dataframe for correlation based on the model blups
bb_20_gw_19_coefs <- data_frame(Trait = attr(colMeans(bb_20_gw_19_fam_blup_pop_corr$Sol), "names"),
Value = colMeans(bb_20_gw_19_fam_blup_pop_corr$Sol)) %>%
separate(Trait, c("Trait","Type","Population"), sep = "\\.", fill = "right") %>%
filter(Type == "Population") %>%
filter(Trait %in% c("traitGrowth", "traitcGDD")) %>%
select(-Type) %>%
spread(Trait, Value)
bb_20_gw_19_coefs <- merge(x=bb_20_gw_19_coefs,
y=exome_meta[,c("Pop","Region")],
all.x=T,
by.x="Population",
by.y="Pop")
# rename regions
bb_20_gw_19_coefs$Region <- revalue(bb_20_gw_19_coefs$Region, c("C"="Core", "M"="Margin","E"="Edge"))
# reorder regions
bb_20_gw_19_coefs$Region <- factor(bb_20_gw_19_coefs$Region,levels = c("Core", "Margin", "Edge"))
# merge climateNA data
bb_20_gw_19_coefs <- merge(x=bb_20_gw_19_coefs,
y=Population_climateNA[,c("Population","PC1","MAT")],
by="Population")
cor.test(bb_20_gw_19_coefs$traitcGDD,bb_20_gw_19_coefs$traitGrowth) #0.6323891##
## Pearson's product-moment correlation
##
## data: bb_20_gw_19_coefs$traitcGDD and bb_20_gw_19_coefs$traitGrowth
## t = 15.008, df = 338, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.5639602 0.6921929
## sample estimates:
## cor
## 0.6323891Family
# against previous year
bb_20_gw_19 <- merge(budbreak_20,heightGrowth_19[,c("Family","Growth")])
bb_20_gw_19 <- bb_20_gw_19[!is.na(bb_20_gw_19$cGDD),]
bb_20_gw_19 <- bb_20_gw_19[!is.na(bb_20_gw_19$Growth),]
# merge Region data
bb_20_gw_19 <- merge(x=bb_20_gw_19,
y=exome_meta[,c("Pop","Region")],
by.x="Population",
by.y="Pop")
# rename and reorder regions
bb_20_gw_19$Region <- plyr::revalue(bb_20_gw_19$Region,
c("C" = "Core",
"M" = "Margin",
"E" = "Edge"))
bb_20_gw_19$Region <- factor(bb_20_gw_19$Region,levels = c("Core", "Margin", "Edge"))
# merge climateNA data
bb_20_gw_19 <- merge(x=bb_20_gw_19,
y=Family_climateNA[,c("Population","PC1","MAT")],
by="Population")
bb_20_gw_19 <- distinct(bb_20_gw_19)Population
bb_20_gw_19$Population <- as.factor(bb_20_gw_19$Population)
bb_20_gw_19_pop <- bb_20_gw_19 %>%
select(Population, cGDD, Growth, Region) %>%
group_by(Population) %>%
dplyr::summarise(cGDD = mean(cGDD), Growth = mean(Growth))
bb_20_gw_19_pop <- merge(x=bb_20_gw_19_pop,
y=Population_climateNA[,c("Population","PC1","MAT")],
by="Population")
bb_20_gw_19_pop <- merge(x=bb_20_gw_19_pop,
y=exome_meta[,c("Pop","Region")],
by.x="Population",
by.y="Pop")
# rename and reorder regions
bb_20_gw_19_pop$Region <- plyr::revalue(bb_20_gw_19_pop$Region,
c("C" = "Core",
"M" = "Margin",
"E" = "Edge"))
bb_20_gw_19_pop$Region <- factor(bb_20_gw_19_pop$Region,levels = c("Core", "Margin", "Edge"))
bb_20_gw_19_pop <- distinct(bb_20_gw_19_pop)MCMC plot
# plot
ggplot(bb_20_gw_19_coefs, aes(x = traitcGDD, y = traitGrowth, group = Population, color = PC1)) +
geom_point(aes(shape = Region), size=3, alpha = 0.7) +
scale_color_viridis(option = "A", direction = -1) +
labs(x = "Budbreak in 2020", y = "Height Growth in 2019") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) #+
# annotate("text", label = "Correlation between traits: 0.63", x = 0.632, y = -1.4) Family plot
# plot
ggplot(bb_20_gw_19, aes(x = Growth, y = cGDD, group = Population, color = PC1)) +
geom_point(aes(shape = Region), size=3, alpha = 0.7) +
scale_color_viridis(option = "A", direction = -1) +
labs(x = "Height Growth in 2019", y = "Budbreak in 2020") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) #+
# annotate("text", label = "Correlation between traits: 0.63", x = 0.632, y = -1.4) Population plot
# plot
ggplot(bb_20_gw_19_pop, aes(x = Growth, y = cGDD, group = Population, color = PC1)) +
geom_point(aes(shape = Region), size=3, alpha = 0.7) +
scale_color_viridis(option = "A", direction = -1) +
labs(x = "Height Growth in 2019", y = "Budbreak in 2020") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) 
Budbreak 2020-Growth 2020
Summary
# against current year
bb_20_gw_20_fam_blup_pop_corr <- readRDS("./genetic_correlation_outputs/bb_20_gw_20_fam_blup_pop_corr")
summary(bb_20_gw_20_fam_blup_pop_corr)##
## Iterations = 10001:9999001
## Thinning interval = 1000
## Sample size = 9990
##
## DIC: 1538.56
##
## G-structure: ~us(trait):Population
##
## post.mean l-95% CI u-95% CI eff.samp
## traitcGDD:traitcGDD.Population 0.41857 0.2326 0.61840 9990
## traitGrowth:traitcGDD.Population -0.08789 -0.2336 0.05614 10461
## traitcGDD:traitGrowth.Population -0.08789 -0.2336 0.05614 10461
## traitGrowth:traitGrowth.Population 0.56880 0.3558 0.80790 9990
##
## R-structure: ~us(trait):units
##
## post.mean l-95% CI u-95% CI eff.samp
## traitcGDD:traitcGDD.units 0.6248 0.52349 0.7333 9990
## traitGrowth:traitcGDD.units 0.1435 0.07832 0.2100 9990
## traitcGDD:traitGrowth.units 0.1435 0.07832 0.2100 9990
## traitGrowth:traitGrowth.units 0.4356 0.36386 0.5127 9990
##
## Location effects: cbind(scale(cGDD), scale(Growth)) ~ trait - 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## traitcGDD 0.01398 -0.16261 0.19615 10337 0.874
## traitGrowth -0.00632 -0.20359 0.19154 9985 0.939effectiveSize(bb_20_gw_20_fam_blup_pop_corr$VCV)## traitcGDD:traitcGDD.Population traitGrowth:traitcGDD.Population
## 9990.00 10460.61
## traitcGDD:traitGrowth.Population traitGrowth:traitGrowth.Population
## 10460.61 9990.00
## traitcGDD:traitcGDD.units traitGrowth:traitcGDD.units
## 9990.00 9990.00
## traitcGDD:traitGrowth.units traitGrowth:traitGrowth.units
## 9990.00 9990.00heidel.diag(bb_20_gw_20_fam_blup_pop_corr$VCV)##
## Stationarity start p-value
## test iteration
## traitcGDD:traitcGDD.Population passed 1 0.7857
## traitGrowth:traitcGDD.Population passed 1 0.8362
## traitcGDD:traitGrowth.Population passed 1 0.8362
## traitGrowth:traitGrowth.Population passed 1 0.6697
## traitcGDD:traitcGDD.units passed 1 0.8574
## traitGrowth:traitcGDD.units passed 1 0.0621
## traitcGDD:traitGrowth.units passed 1 0.0621
## traitGrowth:traitGrowth.units passed 1 0.5786
##
## Halfwidth Mean Halfwidth
## test
## traitcGDD:traitcGDD.Population passed 0.4186 0.002004
## traitGrowth:traitcGDD.Population passed -0.0879 0.001410
## traitcGDD:traitGrowth.Population passed -0.0879 0.001410
## traitGrowth:traitGrowth.Population passed 0.5688 0.002355
## traitcGDD:traitcGDD.units passed 0.6248 0.001066
## traitGrowth:traitcGDD.units passed 0.1435 0.000657
## traitcGDD:traitGrowth.units passed 0.1435 0.000657
## traitGrowth:traitGrowth.units passed 0.4356 0.000744head(bb_20_gw_20_fam_blup_pop_corr$VCV)## Markov Chain Monte Carlo (MCMC) output:
## Start = 10001
## End = 16001
## Thinning interval = 1000
## traitcGDD:traitcGDD.Population traitGrowth:traitcGDD.Population
## [1,] 0.7332963 -0.19062354
## [2,] 0.3624696 -0.11983786
## [3,] 0.3975573 -0.15806714
## [4,] 0.4882200 0.04925123
## [5,] 0.4856402 -0.02837657
## [6,] 0.3152074 -0.08132773
## [7,] 0.3332120 0.01372202
## traitcGDD:traitGrowth.Population traitGrowth:traitGrowth.Population
## [1,] -0.19062354 0.5440355
## [2,] -0.11983786 0.4497206
## [3,] -0.15806714 0.6325305
## [4,] 0.04925123 0.7300534
## [5,] -0.02837657 0.5837543
## [6,] -0.08132773 0.4991884
## [7,] 0.01372202 0.6434916
## traitcGDD:traitcGDD.units traitGrowth:traitcGDD.units
## [1,] 0.6823315 0.18386770
## [2,] 0.6783677 0.18096362
## [3,] 0.5566154 0.13404684
## [4,] 0.6425861 0.18587703
## [5,] 0.6195159 0.16434974
## [6,] 0.5339721 0.08307037
## [7,] 0.6016257 0.13589455
## traitcGDD:traitGrowth.units traitGrowth:traitGrowth.units
## [1,] 0.18386770 0.4330680
## [2,] 0.18096362 0.4320629
## [3,] 0.13404684 0.4001016
## [4,] 0.18587703 0.4303359
## [5,] 0.16434974 0.4041085
## [6,] 0.08307037 0.4042373
## [7,] 0.13589455 0.4029554# genetic correlation
bb_20_gw_20_pop_gencorr <- bb_20_gw_20_fam_blup_pop_corr$VCV[,"traitcGDD:traitGrowth.Population"]/sqrt(bb_20_gw_20_fam_blup_pop_corr$VCV[,"traitcGDD:traitcGDD.Population"] * bb_20_gw_20_fam_blup_pop_corr$VCV[,"traitGrowth:traitGrowth.Population"])
posterior.mode(bb_20_gw_20_pop_gencorr)## var1
## -0.1885222HPDinterval(bb_20_gw_20_pop_gencorr)## lower upper
## var1 -0.4708898 0.09056808
## attr(,"Probability")
## [1] 0.9499499Data
MCMC
## setting up dataframe for correlation based on the model blups
bb_20_gw_20_coefs <- data_frame(Trait = attr(colMeans(bb_20_gw_20_fam_blup_pop_corr$Sol), "names"),
Value = colMeans(bb_20_gw_20_fam_blup_pop_corr$Sol)) %>%
separate(Trait, c("Trait","Type","Population"), sep = "\\.", fill = "right") %>%
filter(Type == "Population") %>%
filter(Trait %in% c("traitGrowth", "traitcGDD")) %>%
select(-Type) %>%
spread(Trait, Value)
bb_20_gw_20_coefs <- merge(x=bb_20_gw_20_coefs,
y=exome_meta[,c("Pop","Region")],
all.x=T,
by.x="Population",
by.y="Pop")
# rename regions
bb_20_gw_20_coefs$Region <- revalue(bb_20_gw_20_coefs$Region, c("C"="Core", "M"="Margin","E"="Edge"))
# reorder regions
bb_20_gw_20_coefs$Region <- factor(bb_20_gw_20_coefs$Region,levels = c("Core", "Margin", "Edge"))
# merge climateNA data
bb_20_gw_20_coefs <- merge(x=bb_20_gw_20_coefs,
y=Population_climateNA[,c("Population","PC1","MAT")],
by="Population")
cor.test(bb_20_gw_20_coefs$traitcGDD,bb_20_gw_20_coefs$traitGrowth) #-0.2778647 ##
## Pearson's product-moment correlation
##
## data: bb_20_gw_20_coefs$traitcGDD and bb_20_gw_20_coefs$traitGrowth
## t = -5.3179, df = 338, p-value = 1.913e-07
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.3731973 -0.1767255
## sample estimates:
## cor
## -0.2778647Family
# against current year
bb_20_gw_20 <- merge(budbreak_20,heightGrowth_20[,c("Family","Growth")])
bb_20_gw_20 <- bb_20_gw_20[!is.na(bb_20_gw_20$cGDD),]
bb_20_gw_20 <- bb_20_gw_20[!is.na(bb_20_gw_20$Growth),]
# merge Region data
bb_20_gw_20 <- merge(x=bb_20_gw_20,
y=exome_meta[,c("Pop","Region")],
by.x="Population",
by.y="Pop")
# rename and reorder regions
bb_20_gw_20$Region <- plyr::revalue(bb_20_gw_20$Region,
c("C" = "Core",
"M" = "Margin",
"E" = "Edge"))
bb_20_gw_20$Region <- factor(bb_20_gw_20$Region,levels = c("Core", "Margin", "Edge"))
# merge climateNA data
bb_20_gw_20 <- merge(x=bb_20_gw_20,
y=Family_climateNA[,c("Population","PC1","MAT")],
by="Population")
bb_20_gw_20 <- distinct(bb_20_gw_20)Population
bb_20_gw_20$Population <- as.factor(bb_20_gw_20$Population)
bb_20_gw_20_pop <- bb_20_gw_20 %>%
select(Population, cGDD, Growth, Region) %>%
group_by(Population) %>%
dplyr::summarise(cGDD = mean(cGDD), Growth = mean(Growth))
bb_20_gw_20_pop <- merge(x=bb_20_gw_20_pop,
y=Population_climateNA[,c("Population","PC1","MAT")],
by="Population")
bb_20_gw_20_pop <- merge(x=bb_20_gw_20_pop,
y=exome_meta[,c("Pop","Region")],
by.x="Population",
by.y="Pop")
# rename and reorder regions
bb_20_gw_20_pop$Region <- plyr::revalue(bb_20_gw_20_pop$Region,
c("C" = "Core",
"M" = "Margin",
"E" = "Edge"))
bb_20_gw_20_pop$Region <- factor(bb_20_gw_20_pop$Region,levels = c("Core", "Margin", "Edge"))
bb_20_gw_20_pop <- distinct(bb_20_gw_20_pop)MCMC plot
# plot
ggplot(bb_20_gw_20_coefs, aes(x = traitcGDD, y = traitGrowth, group = Population, color = PC1)) +
geom_point(aes(shape = Region), size=3, alpha = 0.7) +
scale_color_viridis(option = "A", direction = -1) +
labs(x = "Budbreak in 2020", y = "Height Growth in 2020") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) #+
# annotate("text", label = "Correlation between traits: -0.278", x = 0.7, y = -0.9)Family plot
# plot
ggplot(bb_20_gw_20, aes(x = cGDD, y = Growth, group = Population, color = PC1)) +
geom_point(aes(shape = Region), size=3, alpha = 0.7) +
scale_color_viridis(option = "A", direction = -1) +
labs(x = "Budbreak in 2020", y = "Height Growth in 2020") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) #+
# annotate("text", label = "Correlation between traits: -0.278", x = 0.7, y = -0.9)Population plot
# plot
ggplot(bb_20_gw_20_pop, aes(x = cGDD, y = Growth, group = Population, color = PC1)) +
geom_point(aes(shape = Region), size=3, alpha = 0.7) +
scale_color_viridis(option = "A", direction = -1) +
labs(x = "Budbreak in 2020", y = "Height Growth in 2020") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) 
Phenology vs Phenology
Budbreak 2020-Budset 2019
Summary
# against previous year
bb_20_bs_19_fam_blup_pop_corr <- readRDS("./genetic_correlation_outputs/bb_20_bs_19_fam_blup_pop_corr")
summary(bb_20_bs_19_fam_blup_pop_corr)##
## Iterations = 10001:9999001
## Thinning interval = 1000
## Sample size = 9990
##
## DIC: 1638.933
##
## G-structure: ~us(trait):Population
##
## post.mean l-95% CI u-95% CI eff.samp
## traitcGDD:traitcGDD.Population 0.4197 0.2435 0.6124 9990
## traitBudSet:traitcGDD.Population 0.3551 0.2026 0.5117 9990
## traitcGDD:traitBudSet.Population 0.3551 0.2026 0.5117 9990
## traitBudSet:traitBudSet.Population 0.4205 0.2389 0.6162 9990
##
## R-structure: ~us(trait):units
##
## post.mean l-95% CI u-95% CI eff.samp
## traitcGDD:traitcGDD.units 0.61736 0.5140 0.72133 9990
## traitBudSet:traitcGDD.units -0.09953 -0.1716 -0.02526 10398
## traitcGDD:traitBudSet.units -0.09953 -0.1716 -0.02526 10398
## traitBudSet:traitBudSet.units 0.61288 0.5143 0.71680 9990
##
## Location effects: cbind(scale(cGDD), scale(BudSet)) ~ trait - 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## traitcGDD 0.012704 -0.162639 0.194799 9990 0.884
## traitBudSet 0.006879 -0.168223 0.188259 9990 0.946effectiveSize(bb_20_bs_19_fam_blup_pop_corr$VCV)## traitcGDD:traitcGDD.Population traitBudSet:traitcGDD.Population
## 9990.00 9990.00
## traitcGDD:traitBudSet.Population traitBudSet:traitBudSet.Population
## 9990.00 9990.00
## traitcGDD:traitcGDD.units traitBudSet:traitcGDD.units
## 9990.00 10397.68
## traitcGDD:traitBudSet.units traitBudSet:traitBudSet.units
## 10397.68 9990.00heidel.diag(bb_20_bs_19_fam_blup_pop_corr$VCV)##
## Stationarity start p-value
## test iteration
## traitcGDD:traitcGDD.Population passed 1 0.326
## traitBudSet:traitcGDD.Population passed 1 0.201
## traitcGDD:traitBudSet.Population passed 1 0.201
## traitBudSet:traitBudSet.Population passed 1 0.702
## traitcGDD:traitcGDD.units passed 1 0.465
## traitBudSet:traitcGDD.units passed 1 0.536
## traitcGDD:traitBudSet.units passed 1 0.536
## traitBudSet:traitBudSet.units passed 1 0.487
##
## Halfwidth Mean Halfwidth
## test
## traitcGDD:traitcGDD.Population passed 0.4197 0.001946
## traitBudSet:traitcGDD.Population passed 0.3551 0.001601
## traitcGDD:traitBudSet.Population passed 0.3551 0.001601
## traitBudSet:traitBudSet.Population passed 0.4205 0.001962
## traitcGDD:traitcGDD.units passed 0.6174 0.001036
## traitBudSet:traitcGDD.units passed -0.0995 0.000724
## traitcGDD:traitBudSet.units passed -0.0995 0.000724
## traitBudSet:traitBudSet.units passed 0.6129 0.001018head(bb_20_bs_19_fam_blup_pop_corr$VCV)## Markov Chain Monte Carlo (MCMC) output:
## Start = 10001
## End = 16001
## Thinning interval = 1000
## traitcGDD:traitcGDD.Population traitBudSet:traitcGDD.Population
## [1,] 0.3252018 0.2629739
## [2,] 0.4190668 0.3745077
## [3,] 0.5334350 0.4563674
## [4,] 0.2391750 0.2448858
## [5,] 0.4862936 0.5060902
## [6,] 0.4595838 0.2931318
## [7,] 0.8025672 0.6471861
## traitcGDD:traitBudSet.Population traitBudSet:traitBudSet.Population
## [1,] 0.2629739 0.3297154
## [2,] 0.3745077 0.4040576
## [3,] 0.4563674 0.4945162
## [4,] 0.2448858 0.3137955
## [5,] 0.5060902 0.6161313
## [6,] 0.2931318 0.3040156
## [7,] 0.6471861 0.6590743
## traitcGDD:traitcGDD.units traitBudSet:traitcGDD.units
## [1,] 0.7455586 -0.06242558
## [2,] 0.6392679 -0.06050381
## [3,] 0.5723862 -0.11720235
## [4,] 0.7118270 -0.07310031
## [5,] 0.7397971 -0.11217070
## [6,] 0.6070517 -0.08031806
## [7,] 0.5928800 -0.01726702
## traitcGDD:traitBudSet.units traitBudSet:traitBudSet.units
## [1,] -0.06242558 0.5313842
## [2,] -0.06050381 0.5772358
## [3,] -0.11720235 0.5768493
## [4,] -0.07310031 0.6250420
## [5,] -0.11217070 0.5109948
## [6,] -0.08031806 0.5883263
## [7,] -0.01726702 0.6264233# genetic correlation
bb_20_bs_19_pop_gencorr <- bb_20_bs_19_fam_blup_pop_corr$VCV[,"traitcGDD:traitBudSet.Population"]/sqrt(bb_20_bs_19_fam_blup_pop_corr$VCV[,"traitcGDD:traitcGDD.Population"] * bb_20_bs_19_fam_blup_pop_corr$VCV[,"traitBudSet:traitBudSet.Population"])
posterior.mode(bb_20_bs_19_pop_gencorr)## var1
## 0.8844458HPDinterval(bb_20_bs_19_pop_gencorr)## lower upper
## var1 0.6961866 0.9899278
## attr(,"Probability")
## [1] 0.9499499Data
MCMC
## setting up dataframe for correlation based on the model blups
bb_20_bs_19_coefs <- data_frame(Trait = attr(colMeans(bb_20_bs_19_fam_blup_pop_corr$Sol), "names"),
Value = colMeans(bb_20_bs_19_fam_blup_pop_corr$Sol)) %>%
separate(Trait, c("Trait","Type","Population"), sep = "\\.", fill = "right") %>%
filter(Type == "Population") %>%
filter(Trait %in% c("traitBudSet", "traitcGDD")) %>%
select(-Type) %>%
spread(Trait, Value)
bb_20_bs_19_coefs <- merge(x=bb_20_bs_19_coefs,
y=exome_meta[,c("Pop","Region")],
all.x=T,
by.x="Population",
by.y="Pop")
# rename regions
bb_20_bs_19_coefs$Region <- revalue(bb_20_bs_19_coefs$Region, c("C"="Core", "M"="Margin","E"="Edge"))
# reorder regions
bb_20_bs_19_coefs$Region <- factor(bb_20_bs_19_coefs$Region,levels = c("Core", "Margin", "Edge"))
# merge climateNA data
bb_20_bs_19_coefs <- merge(x=bb_20_bs_19_coefs,
y=Population_climateNA[,c("Population","PC1","MAT")],
by="Population")
cor.test(bb_20_bs_19_coefs$traitcGDD,bb_20_bs_19_coefs$traitBudSet) #0.952888 ##
## Pearson's product-moment correlation
##
## data: bb_20_bs_19_coefs$traitcGDD and bb_20_bs_19_coefs$traitBudSet
## t = 57.756, df = 338, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.9419985 0.9617732
## sample estimates:
## cor
## 0.952888Family
# against previous year
bb_20_bs_19 <- merge(budbreak_20,budset_19[,c("Family","BudSet")])
bb_20_bs_19 <- bb_20_bs_19[!is.na(bb_20_bs_19$BudSet),]
bb_20_bs_19 <- bb_20_bs_19[!is.na(bb_20_bs_19$cGDD),]
# merge Region data
bb_20_bs_19 <- merge(x=bb_20_bs_19,
y=exome_meta[,c("Pop","Region")],
by.x="Population",
by.y="Pop")
# rename and reorder regions
bb_20_bs_19$Region <- plyr::revalue(bb_20_bs_19$Region,
c("C" = "Core",
"M" = "Margin",
"E" = "Edge"))
bb_20_bs_19$Region <- factor(bb_20_bs_19$Region,levels = c("Core", "Margin", "Edge"))
# merge climateNA data
bb_20_bs_19 <- merge(x=bb_20_bs_19,
y=Family_climateNA[,c("Population","PC1","MAT")],
by="Population")
bb_20_bs_19 <- distinct(bb_20_bs_19)Population
bb_20_bs_19$Population <- as.factor(bb_20_bs_19$Population)
bb_20_bs_19_pop <- bb_20_bs_19 %>%
select(Population, cGDD, BudSet, Region) %>%
group_by(Population) %>%
dplyr::summarise(cGDD = mean(cGDD), BudSet = mean(BudSet))
bb_20_bs_19_pop <- merge(x=bb_20_bs_19_pop,
y=Population_climateNA[,c("Population","PC1","MAT")],
by="Population")
bb_20_bs_19_pop <- merge(x=bb_20_bs_19_pop,
y=exome_meta[,c("Pop","Region")],
by.x="Population",
by.y="Pop")
# rename and reorder regions
bb_20_bs_19_pop$Region <- plyr::revalue(bb_20_bs_19_pop$Region,
c("C" = "Core",
"M" = "Margin",
"E" = "Edge"))
bb_20_bs_19_pop$Region <- factor(bb_20_bs_19_pop$Region,levels = c("Core", "Margin", "Edge"))
bb_20_bs_19_pop <- distinct(bb_20_bs_19_pop)MCMC plot
# plot
ggplot(bb_20_bs_19_coefs, aes(x = traitcGDD, y = traitBudSet, group = Population, color = PC1)) +
geom_point(aes(shape = Region), size=3, alpha = 0.7) +
scale_color_viridis(option = "A", direction = -1) +
labs(x = "Budbreak in 2020", y = "Budset in 2019") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) #+
# annotate("text", label = "Correlation between traits: 0.953", x = 0.8, y = -1.5)Family plot
# plot
ggplot(bb_20_bs_19, aes(x = BudSet, y = cGDD, group = Population, color = PC1)) +
geom_point(aes(shape = Region), size=3, alpha = 0.7) +
scale_color_viridis(option = "A", direction = -1) +
labs(x = "Budset in 2019", y = "Budbreak in 2020") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) #+
# annotate("text", label = "Correlation between traits: 0.953", x = 0.8, y = -1.5)Population plot
# plot
ggplot(bb_20_bs_19_pop, aes(x = BudSet, y = cGDD, group = Population, color = PC1)) +
geom_point(aes(shape = Region), size=3, alpha = 0.7) +
scale_color_viridis(option = "A", direction = -1) +
labs(x = "Budset in 2019", y = "Budbreak in 2020") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) 
Budbreak 2020-Budset 2020
Summary
# against current year
bb_20_bs_20_fam_blup_pop_corr <- readRDS("./genetic_correlation_outputs/bb_20_bs_20_fam_blup_pop_corr")
summary(bb_20_bs_20_fam_blup_pop_corr)##
## Iterations = 10001:9999001
## Thinning interval = 1000
## Sample size = 9990
##
## DIC: 1762.214
##
## G-structure: ~us(trait):Population
##
## post.mean l-95% CI u-95% CI eff.samp
## traitcGDD:traitcGDD.Population 0.41906 0.24138 0.6228 11121
## traitBudSet:traitcGDD.Population 0.01342 -0.09562 0.1231 9990
## traitcGDD:traitBudSet.Population 0.01342 -0.09562 0.1231 9990
## traitBudSet:traitBudSet.Population 0.18575 0.06284 0.3178 9990
##
## R-structure: ~us(trait):units
##
## post.mean l-95% CI u-95% CI eff.samp
## traitcGDD:traitcGDD.units 0.623804 0.52373 0.73344 10377
## traitBudSet:traitcGDD.units -0.009515 -0.09257 0.08476 9990
## traitcGDD:traitBudSet.units -0.009515 -0.09257 0.08476 9990
## traitBudSet:traitBudSet.units 0.837642 0.69875 0.98718 9990
##
## Location effects: cbind(scale(cGDD), scale(BudSet)) ~ trait - 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## traitcGDD 0.01451 -0.16675 0.19276 9443 0.875
## traitBudSet -0.00266 -0.14163 0.14990 9990 0.970effectiveSize(bb_20_bs_20_fam_blup_pop_corr$VCV)## traitcGDD:traitcGDD.Population traitBudSet:traitcGDD.Population
## 11121.04 9990.00
## traitcGDD:traitBudSet.Population traitBudSet:traitBudSet.Population
## 9990.00 9990.00
## traitcGDD:traitcGDD.units traitBudSet:traitcGDD.units
## 10376.86 9990.00
## traitcGDD:traitBudSet.units traitBudSet:traitBudSet.units
## 9990.00 9990.00heidel.diag(bb_20_bs_20_fam_blup_pop_corr$VCV)##
## Stationarity start p-value
## test iteration
## traitcGDD:traitcGDD.Population passed 1 0.924
## traitBudSet:traitcGDD.Population passed 1 0.102
## traitcGDD:traitBudSet.Population passed 1 0.102
## traitBudSet:traitBudSet.Population passed 1 0.596
## traitcGDD:traitcGDD.units passed 1 0.224
## traitBudSet:traitcGDD.units passed 1 0.468
## traitcGDD:traitBudSet.units passed 1 0.468
## traitBudSet:traitBudSet.units passed 1 0.686
##
## Halfwidth Mean Halfwidth
## test
## traitcGDD:traitcGDD.Population passed 0.41906 0.001870
## traitBudSet:traitcGDD.Population passed 0.01342 0.001077
## traitcGDD:traitBudSet.Population passed 0.01342 0.001077
## traitBudSet:traitBudSet.Population passed 0.18575 0.001313
## traitcGDD:traitcGDD.units passed 0.62380 0.001041
## traitBudSet:traitcGDD.units passed -0.00952 0.000875
## traitcGDD:traitBudSet.units passed -0.00952 0.000875
## traitBudSet:traitBudSet.units passed 0.83764 0.001461head(bb_20_bs_20_fam_blup_pop_corr$VCV)## Markov Chain Monte Carlo (MCMC) output:
## Start = 10001
## End = 16001
## Thinning interval = 1000
## traitcGDD:traitcGDD.Population traitBudSet:traitcGDD.Population
## [1,] 0.5938413 0.019237060
## [2,] 0.3386735 0.051503667
## [3,] 0.3909280 0.037601014
## [4,] 0.3035156 -0.050115282
## [5,] 0.3526406 0.036870643
## [6,] 0.3016790 -0.006825787
## [7,] 0.2809946 0.021874884
## traitcGDD:traitBudSet.Population traitBudSet:traitBudSet.Population
## [1,] 0.019237060 0.1493862
## [2,] 0.051503667 0.1505557
## [3,] 0.037601014 0.2430505
## [4,] -0.050115282 0.2024429
## [5,] 0.036870643 0.1424567
## [6,] -0.006825787 0.1574127
## [7,] 0.021874884 0.1531409
## traitcGDD:traitcGDD.units traitBudSet:traitcGDD.units
## [1,] 0.5751446 -0.023363690
## [2,] 0.5750166 0.033277313
## [3,] 0.6271948 0.047683487
## [4,] 0.6795201 -0.015439101
## [5,] 0.6361227 -0.002506603
## [6,] 0.6501550 -0.041714029
## [7,] 0.6464858 0.009585230
## traitcGDD:traitBudSet.units traitBudSet:traitBudSet.units
## [1,] -0.023363690 0.7736927
## [2,] 0.033277313 0.9402188
## [3,] 0.047683487 0.8286687
## [4,] -0.015439101 0.8461274
## [5,] -0.002506603 0.8373066
## [6,] -0.041714029 0.7640329
## [7,] 0.009585230 0.8102757# genetic correlation
bb_20_bs_20_pop_gencorr <- bb_20_bs_20_fam_blup_pop_corr$VCV[,"traitcGDD:traitBudSet.Population"]/sqrt(bb_20_bs_20_fam_blup_pop_corr$VCV[,"traitcGDD:traitcGDD.Population"] * bb_20_bs_20_fam_blup_pop_corr$VCV[,"traitBudSet:traitBudSet.Population"])
posterior.mode(bb_20_bs_20_pop_gencorr)## var1
## 0.08470821HPDinterval(bb_20_bs_20_pop_gencorr)## lower upper
## var1 -0.3355211 0.4265705
## attr(,"Probability")
## [1] 0.9499499Data
MCMC
## setting up dataframe for correlation based on the model blups
bb_20_bs_20_coefs <- data_frame(Trait = attr(colMeans(bb_20_bs_20_fam_blup_pop_corr$Sol), "names"),
Value = colMeans(bb_20_bs_20_fam_blup_pop_corr$Sol)) %>%
separate(Trait, c("Trait","Type","Population"), sep = "\\.", fill = "right") %>%
filter(Type == "Population") %>%
filter(Trait %in% c("traitBudSet", "traitcGDD")) %>%
select(-Type) %>%
spread(Trait, Value)
bb_20_bs_20_coefs <- merge(x=bb_20_bs_20_coefs,
y=exome_meta[,c("Pop","Region")],
all.x=T,
by.x="Population",
by.y="Pop")
# rename regions
bb_20_bs_20_coefs$Region <- revalue(bb_20_bs_20_coefs$Region, c("C"="Core", "M"="Margin","E"="Edge"))
# reorder regions
bb_20_bs_20_coefs$Region <- factor(bb_20_bs_20_coefs$Region,levels = c("Core", "Margin", "Edge"))
# merge climateNA data
bb_20_bs_20_coefs <- merge(x=bb_20_bs_20_coefs,
y=Population_climateNA[,c("Population","PC1","MAT")],
by="Population")
cor.test(bb_20_bs_20_coefs$traitcGDD,bb_20_bs_20_coefs$traitBudSet) #0.08288131##
## Pearson's product-moment correlation
##
## data: bb_20_bs_20_coefs$traitcGDD and bb_20_bs_20_coefs$traitBudSet
## t = 1.529, df = 338, p-value = 0.1272
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.02368974 0.18758983
## sample estimates:
## cor
## 0.08288131Family
# against current year
bb_20_bs_20 <- merge(budbreak_20,budset_20[,c("Family","BudSet")])
bb_20_bs_20 <- bb_20_bs_20[!is.na(bb_20_bs_20$BudSet),]
bb_20_bs_20 <- bb_20_bs_20[!is.na(bb_20_bs_20$cGDD),]
# merge Region data
bb_20_bs_20 <- merge(x=bb_20_bs_20,
y=exome_meta[,c("Pop","Region")],
by.x="Population",
by.y="Pop")
# rename and reorder regions
bb_20_bs_20$Region <- plyr::revalue(bb_20_bs_20$Region,
c("C" = "Core",
"M" = "Margin",
"E" = "Edge"))
bb_20_bs_20$Region <- factor(bb_20_bs_20$Region,levels = c("Core", "Margin", "Edge"))
# merge climateNA data
bb_20_bs_20 <- merge(x=bb_20_bs_20,
y=Family_climateNA[,c("Population","PC1","MAT")],
by="Population")
bb_20_bs_20 <- distinct(bb_20_bs_20)Population
bb_20_bs_20$Population <- as.factor(bb_20_bs_20$Population)
bb_20_bs_20_pop <- bb_20_bs_20 %>%
select(Population, cGDD, BudSet, Region) %>%
group_by(Population) %>%
dplyr::summarise(cGDD = mean(cGDD), BudSet = mean(BudSet))
bb_20_bs_20_pop <- merge(x=bb_20_bs_20_pop,
y=Population_climateNA[,c("Population","PC1","MAT")],
by="Population")
bb_20_bs_20_pop <- merge(x=bb_20_bs_20_pop,
y=exome_meta[,c("Pop","Region")],
by.x="Population",
by.y="Pop")
# rename and reorder regions
bb_20_bs_20_pop$Region <- plyr::revalue(bb_20_bs_20_pop$Region,
c("C" = "Core",
"M" = "Margin",
"E" = "Edge"))
bb_20_bs_20_pop$Region <- factor(bb_20_bs_20_pop$Region,levels = c("Core", "Margin", "Edge"))
bb_20_bs_20_pop <- distinct(bb_20_bs_20_pop)MCMC plot
# plot
ggplot(bb_20_bs_20_coefs, aes(x = traitcGDD, y = traitBudSet, group = Population, color = PC1)) +
geom_point(aes(shape = Region), size=3, alpha = 0.7) +
scale_color_viridis(option = "A", direction = -1) +
labs(x = "Budbreak in 2020", y = "Budset in 2020") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) #+
# annotate("text", label = "Correlation between traits: 0.083", x = 0.8, y = -0.9)Family plot
# plot
ggplot(bb_20_bs_20, aes(x = cGDD, y = BudSet, group = Population, color = PC1)) +
geom_point(aes(shape = Region), size=3, alpha = 0.7) +
scale_color_viridis(option = "A", direction = -1) +
labs(x = "Budbreak in 2020", y = "Budset in 2020") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) #+
# annotate("text", label = "Correlation between traits: 0.083", x = 0.8, y = -0.9)Population plot
# plot
ggplot(bb_20_bs_20_pop, aes(x = cGDD, y = BudSet, group = Population, color = PC1)) +
geom_point(aes(shape = Region), size=3, alpha = 0.7) +
scale_color_viridis(option = "A", direction = -1) +
labs(x = "Budbreak in 2020", y = "Budset in 2020") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) 
Within trait correlations
Budset 2019-Budset 2020
Summary
bs_19_bs_20_fam_blup_pop_corr <- readRDS("./genetic_correlation_outputs/bs_19_bs_20_fam_blup_pop_corr")
summary(bs_19_bs_20_fam_blup_pop_corr)##
## Iterations = 10001:9999001
## Thinning interval = 1000
## Sample size = 9990
##
## DIC: 1748.191
##
## G-structure: ~us(trait):Population
##
## post.mean l-95% CI u-95% CI
## traitBudSet2019:traitBudSet2019.Population 0.41490 0.238730 0.6233
## traitBudSet2020:traitBudSet2019.Population 0.09842 -0.007746 0.2137
## traitBudSet2019:traitBudSet2020.Population 0.09842 -0.007746 0.2137
## traitBudSet2020:traitBudSet2020.Population 0.18456 0.059982 0.3124
## eff.samp
## traitBudSet2019:traitBudSet2019.Population 10233
## traitBudSet2020:traitBudSet2019.Population 9990
## traitBudSet2019:traitBudSet2020.Population 9990
## traitBudSet2020:traitBudSet2020.Population 9990
##
## R-structure: ~us(trait):units
##
## post.mean l-95% CI u-95% CI eff.samp
## traitBudSet2019:traitBudSet2019.units 0.6213 0.51628 0.7289 10470
## traitBudSet2020:traitBudSet2019.units 0.1265 0.04174 0.2162 9990
## traitBudSet2019:traitBudSet2020.units 0.1265 0.04174 0.2162 9990
## traitBudSet2020:traitBudSet2020.units 0.8389 0.70079 0.9868 9990
##
## Location effects: cbind(scale(BudSet2019), scale(BudSet2020)) ~ trait - 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## traitBudSet2019 -0.0009861 -0.1877059 0.1749038 8934 0.988
## traitBudSet2020 -0.0021814 -0.1534925 0.1407445 9990 0.971effectiveSize(bs_19_bs_20_fam_blup_pop_corr$VCV)## traitBudSet2019:traitBudSet2019.Population
## 10233.03
## traitBudSet2020:traitBudSet2019.Population
## 9990.00
## traitBudSet2019:traitBudSet2020.Population
## 9990.00
## traitBudSet2020:traitBudSet2020.Population
## 9990.00
## traitBudSet2019:traitBudSet2019.units
## 10469.85
## traitBudSet2020:traitBudSet2019.units
## 9990.00
## traitBudSet2019:traitBudSet2020.units
## 9990.00
## traitBudSet2020:traitBudSet2020.units
## 9990.00heidel.diag(bs_19_bs_20_fam_blup_pop_corr$VCV)##
## Stationarity start p-value
## test iteration
## traitBudSet2019:traitBudSet2019.Population passed 1 0.0502
## traitBudSet2020:traitBudSet2019.Population passed 1 0.1820
## traitBudSet2019:traitBudSet2020.Population passed 1 0.1820
## traitBudSet2020:traitBudSet2020.Population passed 1 0.7749
## traitBudSet2019:traitBudSet2019.units passed 1 0.4283
## traitBudSet2020:traitBudSet2019.units passed 1 0.9023
## traitBudSet2019:traitBudSet2020.units passed 1 0.9023
## traitBudSet2020:traitBudSet2020.units passed 1000 0.1349
##
## Halfwidth Mean Halfwidth
## test
## traitBudSet2019:traitBudSet2019.Population passed 0.4149 0.001971
## traitBudSet2020:traitBudSet2019.Population passed 0.0984 0.001116
## traitBudSet2019:traitBudSet2020.Population passed 0.0984 0.001116
## traitBudSet2020:traitBudSet2020.Population passed 0.1846 0.001301
## traitBudSet2019:traitBudSet2019.units passed 0.6213 0.001038
## traitBudSet2020:traitBudSet2019.units passed 0.1265 0.000883
## traitBudSet2019:traitBudSet2020.units passed 0.1265 0.000883
## traitBudSet2020:traitBudSet2020.units passed 0.8384 0.001525head(bs_19_bs_20_fam_blup_pop_corr$VCV)## Markov Chain Monte Carlo (MCMC) output:
## Start = 10001
## End = 16001
## Thinning interval = 1000
## traitBudSet2019:traitBudSet2019.Population
## [1,] 0.4669275
## [2,] 0.3834779
## [3,] 0.4145036
## [4,] 0.3197958
## [5,] 0.3850419
## [6,] 0.3329645
## [7,] 0.4242091
## traitBudSet2020:traitBudSet2019.Population
## [1,] 0.09420649
## [2,] 0.05773232
## [3,] 0.19339156
## [4,] 0.03599759
## [5,] 0.10061068
## [6,] 0.12324162
## [7,] 0.11999759
## traitBudSet2019:traitBudSet2020.Population
## [1,] 0.09420649
## [2,] 0.05773232
## [3,] 0.19339156
## [4,] 0.03599759
## [5,] 0.10061068
## [6,] 0.12324162
## [7,] 0.11999759
## traitBudSet2020:traitBudSet2020.Population
## [1,] 0.17757150
## [2,] 0.16535716
## [3,] 0.27220901
## [4,] 0.02090638
## [5,] 0.14349698
## [6,] 0.14415338
## [7,] 0.21251820
## traitBudSet2019:traitBudSet2019.units
## [1,] 0.6151109
## [2,] 0.6918221
## [3,] 0.6014122
## [4,] 0.6733291
## [5,] 0.6464305
## [6,] 0.7060680
## [7,] 0.5448681
## traitBudSet2020:traitBudSet2019.units
## [1,] 0.12974131
## [2,] 0.13162973
## [3,] 0.11622205
## [4,] 0.16906297
## [5,] 0.20330544
## [6,] 0.15740882
## [7,] 0.09603986
## traitBudSet2019:traitBudSet2020.units
## [1,] 0.12974131
## [2,] 0.13162973
## [3,] 0.11622205
## [4,] 0.16906297
## [5,] 0.20330544
## [6,] 0.15740882
## [7,] 0.09603986
## traitBudSet2020:traitBudSet2020.units
## [1,] 0.7475624
## [2,] 0.8964411
## [3,] 0.7071179
## [4,] 1.0604594
## [5,] 0.8480110
## [6,] 0.8080282
## [7,] 0.7930861# genetic correlation
bs_19_bs_20_pop_gencorr <- bs_19_bs_20_fam_blup_pop_corr$VCV[,"traitBudSet2019:traitBudSet2020.Population"]/sqrt(bs_19_bs_20_fam_blup_pop_corr$VCV[,"traitBudSet2019:traitBudSet2019.Population"] * bs_19_bs_20_fam_blup_pop_corr$VCV[,"traitBudSet2020:traitBudSet2020.Population"])
posterior.mode(bs_19_bs_20_pop_gencorr)## var1
## 0.322203HPDinterval(bs_19_bs_20_pop_gencorr)## lower upper
## var1 0.0106626 0.6950468
## attr(,"Probability")
## [1] 0.9499499Data
MCMC
## setting up dataframe for correlation based on the model blups
bs_19_bs_20_coefs <- data_frame(Trait = attr(colMeans(bs_19_bs_20_fam_blup_pop_corr$Sol), "names"),
Value = colMeans(bs_19_bs_20_fam_blup_pop_corr$Sol)) %>%
separate(Trait, c("Trait","Type","Population"), sep = "\\.", fill = "right") %>%
filter(Type == "Population") %>%
filter(Trait %in% c("traitBudSet2019", "traitBudSet2020")) %>%
select(-Type) %>%
spread(Trait, Value)
bs_19_bs_20_coefs <- merge(x=bs_19_bs_20_coefs,
y=exome_meta[,c("Pop","Region")],
all.x=T,
by.x="Population",
by.y="Pop")
# rename regions
bs_19_bs_20_coefs$Region <- revalue(bs_19_bs_20_coefs$Region, c("C"="Core", "M"="Margin","E"="Edge"))
# reorder regions
bs_19_bs_20_coefs$Region <- factor(bs_19_bs_20_coefs$Region,levels = c("Core", "Margin", "Edge"))
# merge climateNA data
bs_19_bs_20_coefs <- merge(x=bs_19_bs_20_coefs,
y=Population_climateNA[,c("Population","PC1","MAT")],
by="Population")
cor.test(bs_19_bs_20_coefs$traitBudSet2020,bs_19_bs_20_coefs$traitBudSet2019) #0.4646913 ##
## Pearson's product-moment correlation
##
## data: bs_19_bs_20_coefs$traitBudSet2020 and bs_19_bs_20_coefs$traitBudSet2019
## t = 9.6482, df = 338, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.3769606 0.5441582
## sample estimates:
## cor
## 0.4646913Family
budset_2019 <- budset_19
colnames(budset_2019)[colnames(budset_2019)=="BudSet"] <- "BudSet2019"
budset_2020 <- budset_20
colnames(budset_2020)[colnames(budset_2020)=="BudSet"] <- "BudSet2020"
bs_19_bs_20 <- merge(budset_2019,budset_2020[,c("Family","BudSet2020")])
bs_19_bs_20 <- bs_19_bs_20[!is.na(bs_19_bs_20$BudSet2020),]
bs_19_bs_20 <- bs_19_bs_20[!is.na(bs_19_bs_20$BudSet2019),]
# merge Region data
bs_19_bs_20 <- merge(x=bs_19_bs_20,
y=exome_meta[,c("Pop","Region")],
by.x="Population",
by.y="Pop")
# rename and reorder regions
bs_19_bs_20$Region <- plyr::revalue(bs_19_bs_20$Region,
c("C" = "Core",
"M" = "Margin",
"E" = "Edge"))
bs_19_bs_20$Region <- factor(bs_19_bs_20$Region,levels = c("Core", "Margin", "Edge"))
# merge climateNA data
bs_19_bs_20 <- merge(x=bs_19_bs_20,
y=Family_climateNA[,c("Population","PC1","MAT")],
by="Population")
bs_19_bs_20 <- distinct(bs_19_bs_20)Population
bs_19_bs_20$Population <- as.factor(bs_19_bs_20$Population)
bs_19_bs_20_pop <- bs_19_bs_20 %>%
select(Population, BudSet2019, BudSet2020, Region) %>%
group_by(Population) %>%
dplyr::summarise(BudSet2019 = mean(BudSet2019), BudSet2020 = mean(BudSet2020))
bs_19_bs_20_pop <- merge(x=bs_19_bs_20_pop,
y=Population_climateNA[,c("Population","PC1","MAT")],
by="Population")
bs_19_bs_20_pop <- merge(x=bs_19_bs_20_pop,
y=exome_meta[,c("Pop","Region")],
by.x="Population",
by.y="Pop")
# rename and reorder regions
bs_19_bs_20_pop$Region <- plyr::revalue(bs_19_bs_20_pop$Region,
c("C" = "Core",
"M" = "Margin",
"E" = "Edge"))
bs_19_bs_20_pop$Region <- factor(bs_19_bs_20_pop$Region,levels = c("Core", "Margin", "Edge"))
bs_19_bs_20_pop <- distinct(bs_19_bs_20_pop)MCMC plot
# plot
ggplot(bs_19_bs_20_coefs, aes(x = traitBudSet2020, y = traitBudSet2019, group = Population, color = PC1)) +
geom_point(aes(shape = Region), size=3, alpha = 0.7) +
scale_color_viridis(option = "A", direction = -1) +
labs(x = "Budset in 2020", y = "Budset in 2019") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) #+
# annotate("text", label = "Correlation between traits: 0.465", x = 1.1, y = -1.7)Family plot
# plot
ggplot(bs_19_bs_20, aes(x = BudSet2019, y = BudSet2020, group = Population, color = PC1)) +
geom_point(aes(shape = Region), size=3, alpha = 0.7) +
scale_color_viridis(option = "A", direction = -1) +
labs(x = "Budset in 2019", y = "Budset in 2020") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) #+
# annotate("text", label = "Correlation between traits: 0.465", x = 1.1, y = -1.7)Population plot
# plot
ggplot(bs_19_bs_20_pop, aes(x = BudSet2019, y = BudSet2020, group = Population, color = PC1)) +
geom_point(aes(shape = Region), size=3, alpha = 0.7) +
scale_color_viridis(option = "A", direction = -1) +
labs(x = "Budset in 2019", y = "Budset in 2020") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) 
Growth 2019-Growth 2020
Summary
gw_19_gw_20_fam_blup_pop_corr <- readRDS("./genetic_correlation_outputs/gw_19_gw_20_fam_blup_pop_corr")
summary(gw_19_gw_20_fam_blup_pop_corr)##
## Iterations = 10001:9999001
## Thinning interval = 1000
## Sample size = 9990
##
## DIC: 1522.016
##
## G-structure: ~us(trait):Population
##
## post.mean l-95% CI u-95% CI eff.samp
## traitGrowth2019:traitGrowth2019.Population 0.4206 0.2536 0.6211 9990
## traitGrowth2020:traitGrowth2019.Population 0.3128 0.1570 0.4808 9990
## traitGrowth2019:traitGrowth2020.Population 0.3128 0.1570 0.4808 9990
## traitGrowth2020:traitGrowth2020.Population 0.5557 0.3564 0.7908 9990
##
## R-structure: ~us(trait):units
##
## post.mean l-95% CI u-95% CI eff.samp
## traitGrowth2019:traitGrowth2019.units 0.5870 0.48919 0.6860 9990
## traitGrowth2020:traitGrowth2019.units 0.1105 0.04698 0.1713 10290
## traitGrowth2019:traitGrowth2020.units 0.1105 0.04698 0.1713 10290
## traitGrowth2020:traitGrowth2020.units 0.4347 0.36632 0.5113 9990
##
## Location effects: cbind(scale(Growth2019), scale(Growth2020)) ~ trait - 1
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## traitGrowth2019 -0.009875 -0.181653 0.171620 10328 0.905
## traitGrowth2020 -0.004009 -0.190909 0.200578 10313 0.957effectiveSize(gw_19_gw_20_fam_blup_pop_corr$VCV)## traitGrowth2019:traitGrowth2019.Population
## 9990.00
## traitGrowth2020:traitGrowth2019.Population
## 9990.00
## traitGrowth2019:traitGrowth2020.Population
## 9990.00
## traitGrowth2020:traitGrowth2020.Population
## 9990.00
## traitGrowth2019:traitGrowth2019.units
## 9990.00
## traitGrowth2020:traitGrowth2019.units
## 10290.24
## traitGrowth2019:traitGrowth2020.units
## 10290.24
## traitGrowth2020:traitGrowth2020.units
## 9990.00heidel.diag(gw_19_gw_20_fam_blup_pop_corr$VCV)##
## Stationarity start p-value
## test iteration
## traitGrowth2019:traitGrowth2019.Population passed 1 0.641
## traitGrowth2020:traitGrowth2019.Population passed 1 0.668
## traitGrowth2019:traitGrowth2020.Population passed 1 0.668
## traitGrowth2020:traitGrowth2020.Population passed 1 0.647
## traitGrowth2019:traitGrowth2019.units passed 1 0.505
## traitGrowth2020:traitGrowth2019.units passed 1 0.533
## traitGrowth2019:traitGrowth2020.units passed 1 0.533
## traitGrowth2020:traitGrowth2020.units passed 1 0.348
##
## Halfwidth Mean Halfwidth
## test
## traitGrowth2019:traitGrowth2019.Population passed 0.421 0.001907
## traitGrowth2020:traitGrowth2019.Population passed 0.313 0.001662
## traitGrowth2019:traitGrowth2020.Population passed 0.313 0.001662
## traitGrowth2020:traitGrowth2020.Population passed 0.556 0.002273
## traitGrowth2019:traitGrowth2019.units passed 0.587 0.000999
## traitGrowth2020:traitGrowth2019.units passed 0.111 0.000615
## traitGrowth2019:traitGrowth2020.units passed 0.111 0.000615
## traitGrowth2020:traitGrowth2020.units passed 0.435 0.000738head(gw_19_gw_20_fam_blup_pop_corr$VCV)## Markov Chain Monte Carlo (MCMC) output:
## Start = 10001
## End = 16001
## Thinning interval = 1000
## traitGrowth2019:traitGrowth2019.Population
## [1,] 0.3666177
## [2,] 0.3564076
## [3,] 0.3752535
## [4,] 0.4634167
## [5,] 0.3105347
## [6,] 0.5479918
## [7,] 0.6483959
## traitGrowth2020:traitGrowth2019.Population
## [1,] 0.4200191
## [2,] 0.3086909
## [3,] 0.3219415
## [4,] 0.3236853
## [5,] 0.2455000
## [6,] 0.3633362
## [7,] 0.2861508
## traitGrowth2019:traitGrowth2020.Population
## [1,] 0.4200191
## [2,] 0.3086909
## [3,] 0.3219415
## [4,] 0.3236853
## [5,] 0.2455000
## [6,] 0.3633362
## [7,] 0.2861508
## traitGrowth2020:traitGrowth2020.Population
## [1,] 0.9146002
## [2,] 0.8586466
## [3,] 0.6308588
## [4,] 0.5317278
## [5,] 0.4482986
## [6,] 0.5256747
## [7,] 0.5322276
## traitGrowth2019:traitGrowth2019.units
## [1,] 0.5804598
## [2,] 0.5559426
## [3,] 0.7001635
## [4,] 0.6080137
## [5,] 0.4862879
## [6,] 0.5352382
## [7,] 0.5864476
## traitGrowth2020:traitGrowth2019.units
## [1,] 0.1212451
## [2,] 0.1345896
## [3,] 0.1551620
## [4,] 0.1385535
## [5,] 0.1153481
## [6,] 0.1126667
## [7,] 0.1272492
## traitGrowth2019:traitGrowth2020.units
## [1,] 0.1212451
## [2,] 0.1345896
## [3,] 0.1551620
## [4,] 0.1385535
## [5,] 0.1153481
## [6,] 0.1126667
## [7,] 0.1272492
## traitGrowth2020:traitGrowth2020.units
## [1,] 0.5469751
## [2,] 0.3904507
## [3,] 0.5075336
## [4,] 0.4604836
## [5,] 0.4991465
## [6,] 0.4679778
## [7,] 0.4618191# genetic correlation
gw_19_gw_20_pop_gencorr <- gw_19_gw_20_fam_blup_pop_corr$VCV[,"traitGrowth2019:traitGrowth2020.Population"]/sqrt(gw_19_gw_20_fam_blup_pop_corr$VCV[,"traitGrowth2019:traitGrowth2019.Population"] * gw_19_gw_20_fam_blup_pop_corr$VCV[,"traitGrowth2020:traitGrowth2020.Population"])
posterior.mode(gw_19_gw_20_pop_gencorr)## var1
## 0.6761613HPDinterval(gw_19_gw_20_pop_gencorr)## lower upper
## var1 0.4524946 0.8213388
## attr(,"Probability")
## [1] 0.9499499Data
MCMC
## setting up dataframe for correlation based on the model blups
gw_19_gw_20_coefs <- data_frame(Trait = attr(colMeans(gw_19_gw_20_fam_blup_pop_corr$Sol), "names"),
Value = colMeans(gw_19_gw_20_fam_blup_pop_corr$Sol)) %>%
separate(Trait, c("Trait","Type","Population"), sep = "\\.", fill = "right") %>%
filter(Type == "Population") %>%
filter(Trait %in% c("traitGrowth2019", "traitGrowth2020")) %>%
select(-Type) %>%
spread(Trait, Value)
gw_19_gw_20_coefs <- merge(x=gw_19_gw_20_coefs,
y=exome_meta[,c("Pop","Region")],
all.x=T,
by.x="Population",
by.y="Pop")
# rename regions
gw_19_gw_20_coefs$Region <- revalue(gw_19_gw_20_coefs$Region, c("C"="Core", "M"="Margin","E"="Edge"))
# reorder regions
gw_19_gw_20_coefs$Region <- factor(gw_19_gw_20_coefs$Region,levels = c("Core", "Margin", "Edge"))
# merge climateNA data
gw_19_gw_20_coefs <- merge(x=gw_19_gw_20_coefs,
y=Population_climateNA[,c("Population","PC1","MAT")],
by="Population")
cor.test(gw_19_gw_20_coefs$traitGrowth2019,gw_19_gw_20_coefs$traitGrowth2020) #0.747941##
## Pearson's product-moment correlation
##
## data: gw_19_gw_20_coefs$traitGrowth2019 and gw_19_gw_20_coefs$traitGrowth2020
## t = 20.716, df = 338, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6970293 0.7913493
## sample estimates:
## cor
## 0.747941Family
heightGrowth_2019 <- heightGrowth_19
colnames(heightGrowth_2019)[colnames(heightGrowth_2019)=="Growth"] <- "Growth2019"
heightGrowth_2020 <- heightGrowth_20
colnames(heightGrowth_2020)[colnames(heightGrowth_2020)=="Growth"] <- "Growth2020"
gw_19_gw_20 <- merge(heightGrowth_2019,heightGrowth_2020[,c("Family","Growth2020")])
gw_19_gw_20 <- gw_19_gw_20[!is.na(gw_19_gw_20$Growth2019),]
gw_19_gw_20 <- gw_19_gw_20[!is.na(gw_19_gw_20$Growth2020),]
# merge Region data
gw_19_gw_20 <- merge(x=gw_19_gw_20,
y=exome_meta[,c("Pop","Region")],
by.x="Population",
by.y="Pop")
# rename and reorder regions
gw_19_gw_20$Region <- plyr::revalue(gw_19_gw_20$Region,
c("C" = "Core",
"M" = "Margin",
"E" = "Edge"))
gw_19_gw_20$Region <- factor(gw_19_gw_20$Region,levels = c("Core", "Margin", "Edge"))
# merge climateNA data
gw_19_gw_20 <- merge(x=gw_19_gw_20,
y=Family_climateNA[,c("Population","PC1","MAT")],
by="Population")
gw_19_gw_20 <- distinct(gw_19_gw_20)Population
gw_19_gw_20$Population <- as.factor(gw_19_gw_20$Population)
gw_19_gw_20_pop <- gw_19_gw_20 %>%
select(Population, Growth2019, Growth2020, Region) %>%
group_by(Population) %>%
dplyr::summarise(Growth2019 = mean(Growth2019), Growth2020 = mean(Growth2020))
gw_19_gw_20_pop <- merge(x=gw_19_gw_20_pop,
y=Population_climateNA[,c("Population","PC1","MAT")],
by="Population")
gw_19_gw_20_pop <- merge(x=gw_19_gw_20_pop,
y=exome_meta[,c("Pop","Region")],
by.x="Population",
by.y="Pop")
# rename and reorder regions
gw_19_gw_20_pop$Region <- plyr::revalue(gw_19_gw_20_pop$Region,
c("C" = "Core",
"M" = "Margin",
"E" = "Edge"))
gw_19_gw_20_pop$Region <- factor(gw_19_gw_20_pop$Region,levels = c("Core", "Margin", "Edge"))
gw_19_gw_20_pop <- distinct(gw_19_gw_20_pop)MCMC plot
# plot
ggplot(gw_19_gw_20_coefs, aes(x = traitGrowth2020, y = traitGrowth2019, group = Population, color = PC1)) +
geom_point(aes(shape = Region), size=3, alpha = 0.7) +
scale_color_viridis(option = "A", direction = -1) +
labs(x = "Height Growth in 2020", y = "Height Growth in 2019") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) #+
# annotate("text", label = "Correlation between traits: 0.748", x = 3.3, y = -1.2)Family plot
# plot
ggplot(gw_19_gw_20, aes(x = Growth2019, y = Growth2020, group = Population, color = PC1)) +
geom_point(aes(shape = Region), size=3, alpha = 0.7) +
scale_color_viridis(option = "A", direction = -1) +
labs(x = "Height Growth in 2019", y = "Height Growth in 2020") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) #+
# annotate("text", label = "Correlation between traits: 0.748", x = 3.3, y = -1.2)Population plot
# plot
ggplot(gw_19_gw_20_pop, aes(x = Growth2019, y = Growth2020, group = Population, color = PC1)) +
geom_point(aes(shape = Region), size=3, alpha = 0.7) +
scale_color_viridis(option = "A", direction = -1) +
labs(x = "Height Growth in 2019", y = "Height Growth in 2020") +
theme_bw() + theme(axis.text=element_text(size=14),
axis.title=element_text(size=14,face="bold"),
legend.title = element_text(color = "black", size = 14),
legend.text = element_text(color = "black", size = 13)) 