Construction of norms of reaction - Case study : Maritime pine (Pinus pinsater)
Victor Papin, Laurent Bouffier, Leopoldo Sanchez-Rodriguez
B4EST International conference : Data integration for prediction and breeding - June, 23th 2022
The approach presented here is under development in the framework of my thesis and will soon be published in a scientific article. It is presented here as an example but should be discussed and validated with us before being reused and reproduced.
Introduction : Norms of reaction : definition and modeling
What is a norm of reaction ?
For a genotype, the norm of reaction describes all phenotypes that it can produce within a range of environments.
We refere here to the phenotypic plasticity
What is the point of norms of reaction ?
Conventional and genomic selection focuses mainly on final growth traits (height, circumference, straightness). These phenotypic parameters are very integrative: they do not take into account the evolutions during the long trees growth, in connection with the various environments encountered.
The modeling of individual norms of reaction allows to describe phenotypic evolution of trees across a range of environments. Genetic parameters can be then declined for each environment considered.
Therefore, breeders can apply a selection strategy that explicitly takes into account the environmental effect (ex: Specialized varieties for a single type of environment, or generalist varieties with high genetic values in all environments)
How to construct norms of reaction in a forest tree context ?
To construct norms of reaction, repeated phenotypic measurments across environmentsfor each individual are needed. But in some contexts, such as for maritime pine, it is very difficult to evaluate the same individual in different environments. Clonal experimental design are rather rare.
Fortunately, we can use longitudinal data based on wood. Wood is a record of tree anatomical reaction to environmental changes. During wood formation, the cambium reacts to environmental changes by adjusting the anatomy of the newly formed xylem cells. This continuous modification leaves a permanent anatomical trace in the wood that follows seasonal changes and weather events (= rings in a temperate climate)
Densitometric profiles continuously describe the trees growth over the years. Environmental conditions being variable within and between years, we have repeated phenotypic measurements of the same individuals in different environments.
Several approches are possible to construct norms of reaction from densitometric profils. Here we describe probably the simpliest one : norms of reaction based on annual growth. In this approach we will consider each ring as a unit of measurement. For a given ring, the growth of the tree can be summarized by a general phenotypic variable (ring width, ring surface, ring mean density). At the same time, this same ring corresponding to one year can also be associated with the global climate fo the year. From then on, we can model the evolution of the ring phenotypic variable as a function of the annual climate index.
The statistical tool : the random regression model
Here, we want to model individual norms of reaction for maritime pines, using the phenotypic trajectories previously described. The statistical tool we’ll use is the random regression model.
Random regression models integrate three types of data :
Phenotypic data : phenotypic trajectories we want to model
Environmental/Climatic data : envrionmental parameters characterizing each environments (each year in this case)
Genetic data : to describre relationships between indviduals (using pedigree or genomic information)
1. Phenotypic data : exprimental design and measurments
a. Experimental design
In this design, we have two sites (Cestas & Escource) located the Landes forest in the south west of France. Experimental devices were installed in 1996.
Each site includes 150 half-sib families with 35 individual per family.
We have a complete block plan, with one indiviudal of each family in each block.
b. Phenotypic measurments
In each site 25 families were selected as representative of the phenotypic diversity. And in each families, 13 individuals were cored.
That’s mean we have densitometric profiles for 325 individuals per site.
phenotypes <- read.table("Phenotypes_cleaned.txt", header=TRUE)
phenotypes[1:10,c("Stem_id","Site","Position","Ligne","Year","Ring_width","Ring_surface","Ring_mean_density")]## Stem_id Site Position Ligne Year Ring_width Ring_surface Ring_mean_density
## 1 3966538 CESTAS 12 1 2000 12.5222 492.6190 0.3274944
## 2 3966538 CESTAS 12 1 2001 13.6906 1666.0035 0.3410641
## 3 3966538 CESTAS 12 1 2002 9.1186 1763.0522 0.3214382
## 4 3966538 CESTAS 12 1 2003 7.7724 1915.2083 0.3261796
## 5 3966538 CESTAS 12 1 2004 5.5118 1588.1976 0.3778939
## 6 3966562 CESTAS 36 1 1999 14.6558 674.7905 0.2911754
## 7 3966562 CESTAS 36 1 2000 14.8336 2057.2162 0.2848890
## 8 3966562 CESTAS 36 1 2001 11.8618 2639.8710 0.2925545
## 9 3966562 CESTAS 36 1 2002 8.0264 2287.7883 0.2627189
## 10 3966562 CESTAS 36 1 2003 6.4008 2114.5511 0.2897420c. First visualizations
Some plots from direct phenotypic data availiable
2. Climatic data : characterization of each year
a. Available weather data
weather <- read.table("Landes_weather_data.csv", header=TRUE, sep=";")
weather[1:10,c("YEAR","MONTH","DAY","Rainfall","Temperature")]## YEAR MONTH DAY Rainfall Temperature
## 1 2000 1 1 0.0 9.0
## 2 2000 1 2 0.0 8.1
## 3 2000 1 3 0.0 4.3
## 4 2000 1 4 0.0 6.7
## 5 2000 1 5 0.0 8.2
## 6 2000 1 6 0.0 5.1
## 7 2000 1 7 0.0 8.5
## 8 2000 1 8 7.0 10.0
## 9 2000 1 9 1.0 6.5
## 10 2000 1 10 0.5 1.1
b. Global aridity for each year
For each year, we can calculate a global level of aridity using the Martonne index :
Years.aridity <- NULL
for(year in 2000:2019){
weather.year <- weather[weather$YEAR==year,]
ARIDITY <- NULL
for(month in 3:10){
weather.year.month <- weather.year[weather.year$MONT==month,]
aridity <- (12 *sum(weather.year.month$Rainfall)) / (mean(weather.year.month$Temperature) + 10)
ARIDITY <- c(ARIDITY, aridity)
}
Years.aridity <- c(Years.aridity, mean(ARIDITY))
}## YEAR MEAN_ARIDITY
## 1 2000 39.22
## 2 2001 34.16
## 3 2002 27.52
## 4 2003 27.54
## 5 2004 29.64
## 6 2005 22.44
## 7 2006 30.85
## 8 2007 30.67
## 9 2008 39.96
## 10 2009 28.59
## 11 2010 23.82
## 12 2011 16.66
## 13 2012 39.02
## 14 2013 45.78
## 15 2014 28.53
## 16 2015 22.66
## 17 2016 30.68
## 18 2017 32.72
## 19 2018 32.24
## 20 2019 32.32
c. Link to our phenotypic file
phenotypes$Aridity <- Results.Martonne$MEAN_ARIDITY[match(phenotypes$Year, Results.Martonne$YEAR)]
phenotypes <- phenotypes[phenotypes$Year %in% c(2004:2019),]
phenotypes[c(1:10),c("Stem_id","Site","Year","Aridity","Ring_width","Ring_surface","Ring_mean_density")]## Stem_id Site Year Aridity Ring_width Ring_surface Ring_mean_density
## 6470 3953442 ESCOURCE 2004 29.64 8.8392 2561.7833 0.3175980
## 8088 3953442 ESCOURCE 2005 22.44 5.5626 1863.8343 0.2996103
## 8089 3953442 ESCOURCE 2006 30.85 6.0706 2255.9083 0.3472942
## 8090 3953442 ESCOURCE 2007 30.67 6.9342 2860.1347 0.3391853
## 8091 3953442 ESCOURCE 2008 39.96 5.1308 2310.7645 0.3304337
## 9388 3953442 ESCOURCE 2009 28.59 3.9116 1872.7908 0.3364061
## 9389 3953442 ESCOURCE 2010 23.82 2.5654 1280.4599 0.3101849
## 9390 3953442 ESCOURCE 2011 16.66 0.5080 258.4614 0.3501913
## 9391 3953442 ESCOURCE 2012 39.02 0.9652 495.5437 0.3306175
## 10688 3953442 ESCOURCE 2013 45.78 2.5654 1345.5577 0.3289281phenotypes.mean.year <- aggregate(phenotypes$Ring_surface, by=list(phenotypes$Year), FUN=mean)
colnames(phenotypes.mean.year) <- c("Year","Mean.Ring.Surface")
phenotypes.mean.year$Year <- as.numeric(as.character(phenotypes.mean.year$Year))
p1 <- ggplot(data=phenotypes.mean.year, aes(x=Year, y=Mean.Ring.Surface))+geom_smooth()+geom_point() ; p1+ ggtitle("Mean evolution of ring surface across years")
3. Genetic data : relationship between individuals
a. Pedigree information and A matrix
pedigree <- phenotypes[,c("Stem_id","Mum_name","Dad_name")]
head(unique(pedigree))## Stem_id Mum_name Dad_name
## 6470 3953442 F1.2522 F1.0729
## 6475 3953479 F1.2438 0
## 6479 3953481 F1.1678 F1.0777
## 6484 3953510 0251-7 0
## 6489 3953534 3601-1 0
## 6495 3953536 0044-4 0#library(AGHmatrix)
Amat <- Amatrix(pedigree, ploidy=2, w =0)## Verifying conflicting data...
## Organizing data...
## To organize the data in a fast way wasn't possible...
## Trying to organize in a slow (naive) way...
## Your data was chronologically organized with success.
## Processing a large pedigree data... It may take a couple of minutes...
## Constructing matrix A using ploidy = 2
## Completed! Time = 0.4981667 minutesheatmap(as.matrix(Amat), col = coul)
b. Genomic information and G matrix
Genomic_data <- read.table("Genomic_data1.txt", header=T)
Genomic_data[1:5,1:4]## AX.117412740 AX.117420102 AX.117432316 AX.117435752
## 3955008 1 1 0 2
## 3954794 1 1 0 2
## 3954591 2 0 0 2
## 3955305 1 1 1 2
## 3955321 2 2 1 2#Gmat <- Gmatrix(Genomic_data, ploidy=2, w=0)
heatmap(as.matrix(Gmat[1:641,1:641]), col = coul)
c. Comparison between the two matrices :
4. The random regression model
a. Definition
b. Launch of the model with Wombat
The program we use to run our models is wombat (K. Meyer). Free download here : http://didgeridoo.une.edu.au/km/wombat.php
The distribution format is an executable :
After several steps of formatting the data, we need to do a description of the model in the PARAMETER FILE :
Fixed_effects <- c("Polymix","Mean_trajectory")
Random_effetcs <- c("Genetic","Permanent_envrionement")
Polynome_order <- 2
trait <- "Ring_surface"
relationship <- "genomic"
parameter.file <- write.parameter.file(fix.eff = Fixed_effects,
ran.eff = Random_effects,
poly.ord = 2,
trait = "Ring_surface",
relationship = "genomic")
write.table(parameter.file, "wombat.par")No matter how the file is generated (directly written in Notepad++, using R functions..), it must have this final format :
To finally launch the model from R :
#shell("wombat.exe")
c. Outputs of the model
Model fitting : convergence, BIC & R²
Convergence <- read.table("SumEstimates.out", skip=54, nrows=5, header=F, sep="\t")
print(Convergence)## V1
## 1 Convergence criteria for last 3 iterates
## 2 Change in log likelihood = 0.045578 0.001650 0.000166
## 3 Change in parameter vector = 0.002353 0.000483 0.000152
## 4 Norm of gradient vector = 0.2786 0.2779 0.2778
## 5 Newton decrement = -43.7533 -36.1847 -36.2005BIC <- read.table("SumEstimates.out", skip=14, nrows=1, header=F, sep="\t")
print(BIC)## V1
## 1 -1/2 BIC = -68457.044 Penalty factor = 4.585R² (determination coefficient) is the percentage of variance explained by the model. R² = (G_var + Env_var) / (G_var + Env_var + Residual_var). As it is not directly given by Wombat, we need to calculate it :
BestPoint <- read_table2("BestPoint", col_names=F,cols(),skip=c(1))
kable(BestPoint[,c(1:15)])| X1 | X2 | X3 | X4 | X5 | X6 | X7 | X8 | X9 | X10 | X11 | X12 | X13 | X14 | X15 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 132817.0 | 179719.54 | 351074.30 | 316028.30 | 371407.9703 | 532066.035 | 330827.1 | 627008.2 | 730472.8 | 461956.6 | 341916.4 | 412502.5 | 516851.7 | 313739 | 366831.1 |
| 360564.9 | 52158.67 | -17532.08 | 20231.67 | -356.8617 | 6723.318 | NA | NA | NA | NA | NA | NA | NA | NA | NA |
| 1357091.4 | 146583.08 | -26288.49 | 16034.41 | -4777.5818 | 28909.985 | NA | NA | NA | NA | NA | NA | NA | NA | NA |
## [1] 132817.0 179719.5 351074.3 316028.3 371408.0 532066.0 330827.1 627008.2
## [9] 730472.8 461956.6 341916.4 412502.5 516851.7 313739.0 366831.1## [1] 99275.98 123415.86 130849.44 169768.58 174056.55 182359.40 210077.73
## [8] 213030.36 214173.39 218535.33 220087.68 222090.73 229001.62 241564.62
## [15] 277194.67## [1] 480703.6 504879.3 521766.6 626301.0 638644.4 662733.2 743238.2 751651.5
## [9] 754891.9 767160.9 771487.5 777037.6 795920.3 830307.9 951475.3
Genetic parameters : heritability and genetic correlations
RanRegRatios <- read.table("RanRegVarRatios.dat", skip=1)
colnames(RanRegRatios) <- c("level","rrc","animal","animal.se","subject","subject.se")
RanRegR_animal <- RanRegRatios[,c("rrc","animal","animal.se")]
RanRegR_animal$composante <- "animal"
colnames(RanRegR_animal) <- c("Index_variable","value","se","composante")
head(RanRegR_animal)## Index_variable value se composante
## 1 179 0.121 0.052 animal
## 2 220 0.121 0.044 animal
## 3 229 0.135 0.046 animal
## 4 273 0.183 0.049 animal
## 5 278 0.160 0.043 animal
## 6 288 0.166 0.044 animal
CorAll <- read.table("RanRegCorrAll.dat")
CorAll.tmp <- CorAll[CorAll$V1=="animal", c("V5","V6","V8")]
CorAll.tmp.wide <- spread(CorAll.tmp, V5, V8)
row.names(CorAll.tmp.wide) <- CorAll.tmp.wide[,1]
CorAll.tmp.wide <- CorAll.tmp.wide[,-1]
CorAll.tmp.wide <- as.matrix(CorAll.tmp.wide)
corrplot(CorAll.tmp.wide , method="circle", col=COL2('RdBu', 200) , addCoef.col = "white",number.cex = 0.8, tl.col="black", tl.srt=45, diag=TRUE )
Genetic values got from the random regression model
sln <- read.delim(paste("RnSoln_animal.dat"), sep="", header=F)
sln <- sln[-1,-c(6,7,8)]
colnames(sln) <- c("Run_No.", "Original_ID","Tr", "Solution", "Std.Error")
sln$Tr <- as.numeric(as.character(sln$Tr))
sln$Solution <- as.numeric(as.character(sln$Solution))
sln$Original_ID <- as.numeric(as.character(sln$Original_ID))
sln$Run_No. <- as.numeric(as.character(sln$Run_No.))
PolynomOrder <- max(sln$Tr)-1
sln[1:3,2:4]## Original_ID Tr Solution
## 2 3953442 1 -180.17125
## 3 3953442 2 11.25321
## 4 3953442 3 84.40917Exploitation of genetic values and reconstruction of trajectories
GEBVs at each time point were be obtained following to Mrode (2014). For an individual j in the environment t, the GEBVs can be obtained by :
where ϕt is the row vector of the matrix of Legendre polynomials of order 3.
“Legendre” and “stdtime” functions were adapted from Mrode (2005) by Gota Morota :
legendre(2, gengler = F)## [,1] [,2] [,3]
## [1,] 0.7071068 0.000000 0.000000
## [2,] 0.0000000 1.224745 0.000000
## [3,] -0.7905694 0.000000 2.371708##Given time points covariate and order of fit for Legendre polynomials, return matrix 'M' containing the polynomials of standardized time.
stdtime(rrc, 2)## [,1] [,2] [,3]
## [1,] 1 -1.00000000 1.00000000
## [2,] 1 -0.70397112 0.49557534
## [3,] 1 -0.63898917 0.40830716
## [4,] 1 -0.32129964 0.10323346
## [5,] 1 -0.28519856 0.08133822
## [6,] 1 -0.21299639 0.04536746
## [7,] 1 0.06859206 0.00470487
## [8,] 1 0.10469314 0.01096065
## [9,] 1 0.11913357 0.01419281
## [10,] 1 0.17689531 0.03129195
## [11,] 1 0.19855596 0.03942447
## [12,] 1 0.22743682 0.05172751
## [13,] 1 0.33574007 0.11272140
## [14,] 1 0.55956679 0.31311499
## [15,] 1 1.00000000 1.00000000Phi1 <- stdtime(rrc, 2) %*% t(legendre(2, gengler = F))
Phi1## [,1] [,2] [,3]
## [1,] 0.7071068 -1.22474487 1.58113883
## [2,] 0.7071068 -0.86218502 0.38479070
## [3,] 0.7071068 -0.78259871 0.17781604
## [4,] 0.7071068 -0.39351009 -0.54572977
## [5,] 0.7071068 -0.34929547 -0.59765890
## [6,] 0.7071068 -0.26086624 -0.68297103
## [7,] 0.7071068 0.08400777 -0.77941084
## [8,] 0.7071068 0.12822239 -0.76457394
## [9,] 0.7071068 0.14590823 -0.75690821
## [10,] 0.7071068 0.21665162 -0.71635404
## [11,] 0.7071068 0.24318039 -0.69706608
## [12,] 0.7071068 0.27855208 -0.66788686
## [13,] 0.7071068 0.41119593 -0.52322715
## [14,] 0.7071068 0.68532655 -0.04795201
## [15,] 0.7071068 1.22474487 1.58113883ghat.t.y2 <- t(apply(sln, 1, function (x) Phi1 %*% x))
colnames(ghat.t.y2) <- rrc
ghat.t.y2[1:8,]## 179 220 229 273 278 288
## 3953442 -24.840029 -101.47648 -115.95519 -171.06180 -175.44785 -182.92874
## 3953479 161.760571 96.32177 91.42761 104.30537 108.83986 119.29361
## 3953481 -6.338911 28.90119 35.74957 63.31198 65.73666 70.09960
## 3953510 144.524728 116.60699 109.42754 73.25658 69.31735 61.72032
## 3953534 -173.112801 -174.88331 -176.63174 -189.55665 -191.31218 -194.89107
## 3953536 -123.745343 -167.24215 -177.49780 -227.00447 -232.28515 -242.46169
## 3953553 16.130109 -13.16064 -23.48202 -86.44974 -94.42266 -110.55975
## 3953564 11.551102 -73.33041 -91.35434 -171.27069 -179.16337 -194.01721
## 327 332 334 342 345 349
## 3953442 -193.99201 -193.13386 -192.63536 -189.73490 -188.26641 -185.97854
## 3953479 171.92746 179.37609 182.36407 194.29758 198.74103 204.61527
## 3953481 80.29625 80.74694 80.86883 81.01550 80.92731 80.68569
## 3953510 37.73726 35.54958 34.74381 31.94073 31.07139 30.07490
## 3953534 -208.45593 -210.00511 -210.60486 -212.87536 -213.66857 -214.67167
## 3953536 -275.22381 -278.39372 -279.58348 -283.87216 -285.27797 -286.97238
## 3953553 -172.32491 -179.69375 -182.58379 -193.77353 -197.80189 -203.01597
## 3953564 -237.70359 -241.39114 -242.72957 -247.27523 -248.63664 -250.15097
## 364 395 456
## 3953442 -173.94025 -130.60653 36.96463
## 3953479 225.78022 261.20864 264.58856
## 3953481 78.47936 66.98433 13.45373
## 3953510 28.11877 34.40904 99.43898
## 3953534 -217.80202 -220.24452 -201.89004
## 3953536 -291.37123 -289.28117 -229.38985
## 3953553 -220.75113 -245.81757 -228.56106
## 3953564 -252.62283 -240.04411 -131.85987
To facilitate the visualization and interpretation of trajectories, we proceed to a clustering of curves :
#library(kml)
traj <- as.matrix(tab[,-ncol(tab)])
obj <- clusterLongData(traj, idAll=row.names(traj), time=rrc, varNames = "Largeur")
kml(obj,nbClusters=2:10,nbRedrawing=1,toPlot='none')## ~ Fast KmL ~
## *********S#choice(obj, typeGraph = "png")
affectation <- getClusters(obj, nbCluster=5, clusterRank = 1, asInteger = FALSE)
affectation <- data.frame(affectation)
affectation$ind <- row.names(tab)
tablong$group <- affectation$affectation[match(tablong$ind, affectation$ind)]
p <- ggplot(data=tablong, aes(x=control.variable ,y=gBLUP, group=ind))+geom_smooth(se=F,aes(color=group), size=0.7)+xlab("Level of aridity (Martonne index)")+ggtitle("Evolution of the indiviudal genetic value across environments") ; p+theme(plot.title = element_text(hjust = 0.5))