# Workshop on Genomic Evaluation

Marie Pégard

► To cite this version:

Marie Pégard. Workshop on Genomic Evaluation. École thématique. France. 2020, pp.107. hal-04731328

HAL Id: hal-04731328

https://hal.science/hal-04731328v1

Submitted on 10 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

Horizon 2020 of European Union

This project has received funding from the European Union's Horizon 2020
Programme for Research & Innovation under grant agreement n727312

<figure><img src="image_3.png" /></figure>

EUC
LEG

# WorkShop Genomic Evaluation

05/11/2020

Marie PEGARD

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>
Introduction

<figure><img src="image_6.png" /></figure>
A brief presentation of myself

* 2015 - 2018 → PhD : New models for implementation of genome-wide evaluation in black poplar breeding program. - INRAE UMR BioForA

* 2019 → Post-Doctoral Fellowship : Studying the genetic diversity management in the genomic selection context for two species : the maritime pine and the black poplar. - INRAE UMR Biogeco

* Since January → PhD on the EUCLEG projet in the WP5 INRE URP3F Lusignan

* My email : marie.pegard@inrae.fr

<figure><img src="image_7.png" /></figure>
# Workshop Goals

This workshop :

✔ Takes place in the WP5 of the EUCLEG project.

✔ Provide knowledge and tools for breeders

✔ An user-friendly software (Progeno) was developed

✔ R programming language linked to the Progeno program in command line.

This workshop covers :

✔ data preparation

✔ QTL detection

✔ genomic evaluation

✔ genomic selection.

⇒ Methods easy to implement and robust

<figure><img src="image_8.png" /></figure>
# Workshop Organization

<figure><img src="image_9.png" /></figure>

## Theoretical part

<figure><img src="image_10.png" /></figure>

* Data Observation and understanding

<figure><img src="image_11.png" /></figure>

* GWAS: Genome Wide Association Studies

<figure><img src="image_12.png" /></figure>

* Genomic selection

<figure><img src="image_13.png" /></figure>

* R

<figure><img src="image_14.png" /></figure>

* R markdown

<figure><img src="image_15.png" /></figure>

## Overview of the practical part

<figure><img src="image_16.png" /></figure>

* Presentation of the dataset

<figure><img src="image_17.png" /></figure>

* What are we going to do in this workshop?

<figure><img src="image_18.png" /></figure>

* How the practical part will take place ?

<figure><img src="image_19.png" /></figure>

## Practical part

<figure><img src="image_20.png" /></figure>
Data observation and understanding

<figure><img src="image_21.png" /></figure>
# Observation and understanding of data

A good understanding of an experimental design and the collected data observation :

* Is part of the experiment.

* Essential pre-requisite for the field data preparation or analysis

I will provide **guidelines** on :

* How to observe and analyse field data

* How to prepare the data

* How to fit a linear model to adjust the phenotypic data before GWAS, GWE or GS analysis.

<figure><img src="image_22.png" /></figure>
# Why prepare and observe data ?

In genome-wide association studies :

<figure><img src="image_23.png" /></figure>

* Determine the **relationships between phenotypes and genotypes**

<figure><img src="image_24.png" /></figure>

* Rare phenotypes and their associated genotypes are of great interest

**Rare phenotypes** → the result of technical errors or bad models

**False Rare phenotypes** ⇒ leads to the detection of false QTLs or inaccurate prediction

To avoid this ⇒ have the same rigour to the analysis of phenotypic and genotypic data.

<figure><img src="image_25.png" /></figure>
An overview of the process

<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /></figure>
# Phenotypic Data observation

Raw description of our example:

Let us take the following example: Eucleg's Soy data for 1 trial over 2 years.

There is 15 columns, including :

* Field information : Plot, Row, Colum, Trial, Year

* Genetic information : EntryNo, EUCLEGID, Accession name

* Trait of interest : Seed yield, Protein content

* Explanatory variables : R1 (beginning of flowering), R2 (full flowering), R5 (start of seed filling), R8 (maturity)

For the example we are going to look at the protein content trait only

<figure><img src="image_28.png" /></figure>
# Response variable protein content

Descriptive statistics summarize data to facilitate the exploration :

<figure><img src="image_29.png" /></figure>

* Include measures of centrality

<figure><img src="image_30.png" /></figure>

* the mean : 42.34

<figure><img src="image_31.png" /></figure>

* the median : 42.13

<figure><img src="image_32.png" /></figure>

* Include measures of dispersion

<figure><img src="image_33.png" /></figure>

* the deviation : 1.88

<figure><img src="image_34.png" /></figure>

* the range 36.9, 49.48

<figure><img src="image_35.png" /></figure>
## Response variable protein content

I look at the distribution and the phenotypic mean of the traits depending on the year.
The distribution looks normal for the two years.

<figure><img src="image_36.png" /></figure>

<figure><img src="image_37.png" /></figure>
Explanatory variable R8

Descriptive statistics summarize data to facilitate the exploration :

* Include measures of centrality

the mean : 122.19

the median : 123

* Include measures of dispersion

the deviation : 16.79

the range 78, 166

<figure><img src="image_38.png" /></figure>
Explanatory variable R8
---

The distribution is not normal but bimodal.

<figure><img src="image_39.png" /></figure>

<figure><img src="image_40.png" /></figure>
<figure><img src="image_41.png" /></figure>

Correlation between Response and Explanatory variable

<figure><img src="image_42.png" /></figure>

<figure><img src="image_43.png" /></figure>
# Genotyping data - missing values

An important aspect in genotyping data is missing data:

<figure><img src="image_44.png" /></figure>

* too much missing data ⇌ to remove markers or individuals

<figure><img src="image_45.png" /></figure>

* reasonable number of missing data ⇌ several solutions are possible

<figure><img src="image_46.png" /></figure>

* imputation software (BEAGLE, FImpute, Alphalmpute) can be used, or

<figure><img src="image_47.png" /></figure>

* imputation with the mean or allelic frequency

<figure><img src="image_48.png" /></figure>
# Genotyping data - missing values

There is no fixed threshold for determining whether the number of missing data is reasonable or not, it depends on the population of interest:

<figure><img src="image_49.png" /></figure>

* Access to pedigree and full genotyping in the parents ~ imputation of 80% or even 90%

<figure><img src="image_50.png" /></figure>

* Complete independent individuals ~ beyond 1% or 2% of the population, the bias brought about by imputation can be critical.

<figure><img src="image_51.png" /></figure>
# Genotyping data

In the following example there is :

* 360 phenotyped individuals

* 357 genotyped individuals

* 226798 SNPs

<figure><img src="image_52.png" /></figure>
# Visualisation of missing data

In the present situation, the number of missing value is really low (0.211%)

⇒ I replaced the missing value by the minor allelic frequency

<figure><img src="image_53.png" /></figure>

<figure><img src="image_54.png" /></figure>
# Visualisation of the Minor allelic frequency

<figure><img src="image_55.png" /><figcaption>Histogram of freqp</figcaption></figure>

<figure><img src="image_56.png" /></figure>
SNPs trimming

We remove SNPs with a MAF of less than 0.01 (62,863 SNPs are concerned)

At this frequency a very large number of individuals are needed to be able to estimate
their effect correctly.

To let them bias the analysis would be a mistake.

<figure><img src="image_57.png" /></figure>
<figure><img src="image_58.png" /><figcaption>Histogram of freqp</figcaption></figure>

<figure><img src="image_59.png" /></figure>
# Visualization Marker density along the genome

Observation of the distribution of markers along the chromosomes.

There is 99.2041376 % of the markers on the chromosomes, the rest on the scaffolds.

When the physical position of the markers is available, it can be interesting to look at their distribution.

Are there areas of the genome that are absent or less well covered by genotyping?

<figure><img src="image_60.png" /></figure>
Visualization Marker density along the genome

<figure><img src="image_61.png" /></figure>

<figure><img src="image_62.png" /></figure>
Linkage disequilibrium

The Linkage disequilibrium is the Non random association between alleles from different (linked) loci.

A quick way to do this is to calculate a partial correlation squared.

When we have the physical (or genetic map) we can represent it according to the
distance as here with chromosome 1

<figure><img src="image_63.png" /></figure>
Linkage disequilibrium - visualisation

<figure><img src="image_64.png" /></figure>

<figure><img src="image_65.png" /></figure>
# Genetic structure

The genetic structure can be assessed with several softwares and methods:

* STRUCTURE, ADMIXTURE,

* The R package Adegenet (Discriminant Analysis of Principal Components (DAPC))

If a structure in the population is observed. It is always interesting to represent it graphically.

The structure can be used later for:

* QTL detection

* data adjustment

* to optimize the creation of training and validation populations.

<figure><img src="image_66.png" /></figure>
# Visualisation - Genetic structure

<figure><img src="image_67.png" /></figure>

<figure><img src="image_68.png" /></figure>
# Visualisation - Pedigree

In some case, the pedigree information is available but not the genotypes

The pedigree is a kind of genetic structure

## Name P1 P2
## 1: EUC_GM_009 EUC_GM_239 <NA>
## 2: EUC_GM_050 EUC_GM_179 <NA>
## 3: EUC_GM_055 EUC_GM_204 <NA>
## 4: EUC_GM_056 EUC_GM_328 <NA>
## 5: EUC_GM_072 EUC_GM_328 <NA>
## 6: EUC_GM_097 EUC_GM_604 EUC_GM_615

<figure><img src="image_69.png" /></figure>
# Understanding the data

We carried out the first steps of observing the available data

<table><thead><tr><th colspan="2">Understanding the Response variables</th></tr></thead><tbody><tr><td>Graphs</td><td>Descriptive Statistics</td></tr></tbody></table>

<table><thead><tr><th colspan="2">Understanding the explanatory variables</th></tr></thead><tbody><tr><td>Graphs</td><td>Cross-table</td></tr><tr><td></td><td>Freq</td></tr></tbody></table>

<table><thead><tr><th colspan="2">Understanding the Genetic structure</th></tr></thead><tbody><tr><td>Graphs</td><td>Descriptive Statistics</td></tr><tr><td></td><td>Freq</td></tr></tbody></table>

<figure><img src="image_70.png" /></figure>
# Understanding the data

## What did we learn from our observations

1. The response variables have a distribution close to normal

2. The explanatory variable R8 is not correlated with the response variable Proteincontent but has it with the variable Seedyield

3. The genotyping data are of good quality

4. We have replaced the NA and removed the low frequency markers

5. A genetic structure was detected in the data

6. Not all phenotyped individuals are genotyped

7. ...

<figure><img src="image_71.png" /></figure>
Understanding the data

The next step concerns the outliers detection and the data transformation

There is no need to transform our response variable in our example.

In all case, It is necessary to be careful about data transformation

<figure><img src="image_72.png" /></figure>

<figure><img src="image_73.png" /></figure>
## Outliers

Outliers can be detected with a boxplot.

The point after the brackets are outliers.

In the true worlds when there is a continuum like here, it is not outliers.

<figure><img src="image_74.png" /></figure>

<figure><img src="image_75.png" /></figure>
## Linear models

Linear models are a good way of extracting the components of the phenotype and remove all undesired effects

Fit a linear model

ANOVA

Model diagnostics

<figure><img src="image_76.png" /></figure>
Let us specify the writing of the linear model for any individual $i$ ($i = 1, \dots, n$) of a sample of size $n$:

$$Y_i = \sum_{j=1}^{p} \beta_j X_i^j + U_i$$

* $Y_i$ is the true response random variable

* $U_i$ is the true error random variable, assumed to be $N(0, \sigma^2)$ and independent

* $\beta_j$ are coefficients, unknown parameters, to be estimated: generally random effects

* $X_i^j$ are the values of the explanatory variables: generally fixed effects

<figure><img src="image_77.png" /></figure>
Linear models and Phenotypes

Matrix notation

$$Y = X\beta + U$$

* Y and U are random vectors

* X is a matrix n × p

* β is the vector p parameters.

<figure><img src="image_78.png" /></figure>
In our example

The model :

$$Y_{jk} = \mu + s_k + g_j + e_{jk}$$

* Y<sub>jk</sub> : phenotype of the <i>j</i><sup>th</sup> genotype in the <i>k</i><sup>th</sup> spatial location (row and column)

* μ : overall mean

* s<sub>k</sub> : the effect of the <i>k</i><sup>th</sup> spatial location

* g<sub>j</sub> : the genetic effect of the <i>j</i><sup>th</sup> genotype

* e<sub>jk</sub> : the residual

<figure><img src="image_79.png" /></figure>
In our example

The model :

$$
Y_{jk} = \mu + s_k + g_j + e_{jk}
$$

We are going to integrate in our model:

► the intercept

a genetic effect: via the EUCLEDID code, the genetic structure or molecular
markers

a spatial effect: plot (xy coordinates), calculated with a Bi-spline model

the residual

<figure><img src="image_80.png" /></figure>
Bi-splines

✓ A B-spline is a linear combination of positive splines

✓ B-splines are the generalisation of the Bézier curves

✓ The shape of the basic functions is determined by the position of the nodes.

✓ The curve is within the convex envelope of the control points.

<figure><img src="image_81.png" /></figure>

<figure><img src="image_82.png" /></figure>
Linear model and Phenotypic adjustment

For GWAS or GS cases, the only parts of the phenotype we are interested in are :

* the genetic part

* the residual part.

We want the phenotypes to be cleaned of ground and environmental effects

<figure><img src="image_83.png" /></figure>
Linear model and Phenotypic adjustment

It is possible to do this year by year or to combine years.

Let's take the example of the year 2018 for the protein content.

I will use an “EuclegID” effect to represent the genetic effect

I use the R package breedR to do it.

<figure><img src="image_84.png" /></figure>
# Linear model and Phenotypic adjustment

## BreedR function

```r
resA1 <- remlf90(
  fixed = Proteincontent ~ 1,
  random = ~ EUCLEGID, #,
  spatial = list(model = 'splines',
                coord = phenoSoy[phenoSoy$Year == 2018,
                                  c('Row', 'Column')]),
  data = phenoSoy[phenoSoy$Year == 2018,],
  method = 'em'
)
```

<figure><img src="image_85.png" /></figure>
Linear model and Phenotypic adjustment

## Formula: Proteincontent ~ 0 + Intercept + EUCLEGID + spatial

## Data: phenoSoy[phenoSoy$Year == 2018, ]

## AIC BIC logLik

## 1484 unknown -739

##

## Parameters of special components:

## spatial: n.knots: 12 9

##

## Variance components:

## Estimated variances

## EUCLEGID 2.47300

## spatial 0.03121

## Residual 0.23300

##

<figure><img src="image_86.png" /></figure>
Linear model and Phenotypic adjustment

Spatial effet :

<figure><img src="image_87.png" /></figure>

<figure><img src="image_88.png" /></figure>
# Linear model and Phenotypic adjustment

## Residuals

<figure><img src="image_89.png" /></figure>

<figure><img src="image_90.png" /></figure>
# Linear model and Phenotypic adjustment

## Fitted values vs Residuals

<figure><img src="image_91.png" /></figure>

<figure><img src="image_92.png" /></figure>
Linear model and Phenotypic adjustment

Heritability

## [1] 0.9034747

<figure><img src="image_93.png" /></figure>
Linear model and Phenotypic adjustment - 2019

We are doing the same for the second year.

Heritability

## [1] 0.8198598

<figure><img src="image_94.png" /></figure>
Linear model and Phenotypic adjustment - 2019

It is clear that the field effect is different between the two years.

<figure><img src="image_95.png" /></figure>

<figure><img src="image_96.png" /></figure>
Combined

Now We make a model that combines the two years. This will allow a better estimate of
the genetic value of individuals because several years of phenotyping will be cumulated.

The model :

$$
Y_{ijk} = \mu + s_{kl} + g_j + a_i + e_{ijkl}
$$

Y_{ijk} : phenotype of the j^{th} genotype in the i^{th} year of the k^{th} spatial location (row and column)

µ : overall mean

s<sub>k</sub> : the effect of the k<sup>th</sup> spatial location

g<sub>j</sub> : the genetic effect of the j<sup>th</sup> genotype

a<sub>i</sub> : the year effect

e<sub>ijk</sub> : the residual

<figure><img src="image_97.png" /></figure>
Combined

Glob_res <- remlf90(
fixed = Proteincontent ~ 1 + Year,
random = ~ EUCLEGID,
generic = list(sp18 = list(inc.sp18,
                      breedR:::get_structure(sp18)),
            sp19 = list(inc.sp19,
                      breedR:::get_structure(sp19))
),
data = phenoSoy,
method = 'em'
)

<figure><img src="image_98.png" /></figure>
Combined

## Formula: Proteincontent ~ 0 + Intercept + Year + EUCLEGID

## Data: phenoSoy

## AIC BIC logLik

## 3023 unknown -1508

##

## Parameters of special components:

##

##

## Variance components:

## Estimated variances

## EUCLEGID 2.415000

## sp18 0.002286

## sp19 0.338500

## Residual 0.890300

<figure><img src="image_99.png" /></figure>
Combined

<figure><img src="image_100.png" /></figure>

<figure><img src="image_101.png" /></figure>
Combined

<figure><img src="image_102.png" /></figure>
# What if we integrate genomic information

This is equivalent to making a genomic evaluation.

A basic formula is :

$$y = \mu + Xg + \epsilon$$

<figure><img src="image_103.png" /></figure>

g is a vector (p × 1) of GEBV,

<figure><img src="image_104.png" /></figure>

X a design matrix linking observations to individuals.

GEBVs follow a normal distribution, with a covariance between effects modelled through an additive relationship matrix that is derived from markers for GEBV or from pedigree in classical pedigree-based BLUPs.

I will give the details when we discuss GWAS and genomic selection.

<figure><img src="image_105.png" /></figure>
# What if we integrate genomic information

## In breedR

```r
Glob_res_lopez <- remlf90(
  fixed = Proteincontent ~ 1 + Year,
  # random = ~ 1, # If you need random effects
  generic = list(genetic = list(inc.H,HX),
    sp18 = list(inc.sp18,
      breedR:::get_structure(sp18)),
    sp19 = list(inc.sp19,
      breedR:::get_structure(sp19))
  ),
  data = phenoSoy,
  method = 'em'
)
```

<figure><img src="image_106.png" /></figure>
## What if we integrate genomic information

### Formula: Proteincontent ~ 0 + Intercept + Year

### Data: phenoSoy

### AIC BIC logLik

### 2876 unknown -1434

###

### Parameters of special components:

###

###

### Variance components:

### Estimated variances

### generic_genetic 2.46900

### sp18 0.01584

### sp19 0.06144

### Residual 0.95260

<figure><img src="image_107.png" /></figure>
# What if we integrate genomic information

## Spatial effect

<figure><img src="image_108.png" /></figure>

<figure><img src="image_109.png" /></figure>
# What if we integrate genomic information

## [1] 0.7056544

<figure><img src="image_110.png" /></figure>

<figure><img src="image_111.png" /></figure>
# Summary

We have just seen that the model combining the data from the two rating years gives better results.

For the following steps, two solutions are possible:

<figure><img src="image_112.png" /></figure>

* Either we use the combined model without the modular data.

<figure><img src="image_113.png" /></figure>

* Either we use the combined model with the molecular data.

In both cases, the part due to the year and the spatial effect will be removed from the phenotype, leaving only the genetic and residual part.

We are ready to perform GWAS and/or GS

<figure><img src="image_114.png" /></figure>
GWAS

<figure><img src="image_115.png" /></figure>
# Introduction

The GWAS can be used in several situations.

Firstly, when seeking to determine the genetic architecture of a trait of ecological or agronomic interest.

It is also useful in the context of Marker Assisted Breeding and Genomic Selection.

<figure><img src="image_116.png" /></figure>
Introduction

GWAS aims at determine the association between the phenotype and genotype of
individuals in a population by the way of linear regression.

Some parameters influence the GWAS results :

- LD : Non random association between alleles from different (linked) loci

- Population structure : Allele frequencies vary across sub-populations and cause long-range disequilibrium, i.e. between unlinked loci

<figure><img src="image_117.png" /></figure>
Introduction

Several methods have been developed over the years to perform GWAS :

✔️ GRAMMAR (Aulchenko et al., 2007)

✔️ P3D (Zhang et al., 2010) & EMMAX (Kang et al., 2010)

✔️ FaST-LMM (Lippert et al., 2011)

✔️ GEMMA (Zhou & Stephens, 2012)

✔️ MixABEL (Aulchenko, et al.)

✔️ CMLM (Zhang et al., 2010)

✔️ ECMLM (Li et al., 2014)

✔️ MLMM (Segura et al., 2012)

✔️ FarmCPU (Liu et al., 2016)

✔️ ...

<figure><img src="image_118.png" /></figure>
Introduction

A very good illustration in the paper of Liu et al, 2016 shows the differences between the methods.

<figure><img src="image_119.png" /></figure>

<figure><img src="image_120.png" /></figure>
We are going to test the MLMM with the genomic relationship matrix.

As the genomic relationship matrix includes information on both population structure and relatedness, it is in general not useful to consider admixture information as fixed effects covariates (Astle and Balding 2009).

This method is a stepwise model regression with a forward inclusion and a backward elimination.

This limits the space of the model to be explored and makes the method computationally efficient.

<figure><img src="image_121.png" /></figure>
# MLMM - Proteincontent

Graph showing the variance explained by the SNPs at each stage :

* In red the residual variance

* In green the genetic variance

* In blue the variance explained by the QTLs detected.

<figure><img src="image_122.png" /></figure>

<figure><img src="image_123.png" /></figure>
# MLMM - Proteincontent

Representation of the QTLs detected with different step of the MLMM.

<figure><img src="image_124.png" /></figure>

<figure><img src="image_125.png" /></figure>
Genomic Selection

<figure><img src="image_126.png" /></figure>
# Introduction

Genomic selection (GS) is the direct descendant of marker-assisted selection (SAM)

<figure><img src="image_127.png" /><figcaption>The genomic selection</figcaption></figure>

<figure><img src="image_128.png" /></figure>
Accuracy of genomic selection

The accuracy of genomic prediction : the correlation between the true breeding value (TBV) and its estimate (GEBV : Genomic Estimated Breeding Value).

If the observed values are phenotypes, the accuracy of prediction : $\frac{predictive\-ability}{\sqrt{h^2}}$
($r_{GEBV, TBV} = r_{GEBV, Phenotypes}/\sqrt{h^2}$) (Falconer, 1981).

<figure><img src="image_129.png" /></figure>
# Accuracy of genomic selection

<figure><img src="image_130.png" /></figure>

<figure><img src="image_131.png" /></figure>
## Other quality criteria

The predicting ability (or accuracy) is the most common criteria

Other complementary criteria for prediction quality are :

* the slope : should be close to 1

* the intercept : close to zero (but it is not a critical biais)

These two parameter allow to observe model's deficiencies :

* wrong variance partition

* incomplete model,

* nonrandom choice of individuals for training and validation population

<figure><img src="image_132.png" /></figure>
# Cross-validation

## Cross validation & Test Set

<figure><img src="image_133.png" /></figure>

<figure><img src="image_134.png" /></figure>
# GBLUP and other models

Two main groups can be distinguished :

* approaches that estimate an additive effect associated with each marker

* approaches that directly give the additive value of individuals.

In recent years, several methods have been developed :

* derivatives of BLUP methods (Henderson, 1975),

* Bayesian methods and non-parametric methods (Gianola and van Kaam, 2008; Neves et al., 2012 ...).

Roughly, the different strategies are thought to accommodate implicitly different
underlying genetic architectures of quantitative traits:

* from simpler infinitesimal-like architectures

* to those more complex with a variety of heterogeneous effects across genes.

<figure><img src="image_135.png" /></figure>
The formula to estimate GEBV : $y = \mu + Xg + \epsilon$

<figure><img src="image_136.png" /></figure>

$g$ is a vector ($p \times 1$) of GEBV,

<figure><img src="image_137.png" /></figure>

$X$ a design matrix linking observations to individuals.

GEBVs follow a normal distribution, with a **covariance** between effects modelled through an **additive relationship matrix** that is derived from :

* markers for GEBV (GRM)

* from pedigree in classical pedigree-based BLUPs (NRM)

<figure><img src="image_138.png" /></figure>
# GRM matrice

There are several methods for calculating this matrix G (or GRM) :

* the first methods uses a similarity index calculated for each locus : This method is based on the assumption that all identical alleles (IBS) are all identical by descent (IBD).

* This formula can be corrected by taking into account the probability of an allele being IBS for a locus by using the founder's genotypes.

* Another method is to correct the matrix G by twice the allele frequency of the minor allele as in equation (VanRaden 2007, Habier 2008).

<figure><img src="image_139.png" /></figure>
GS example

We will see an example of genomic predictions.

Here I varied the size of the training population.

I divided the data as follows:

<figure><img src="image_140.png" /></figure>

20% in test

<figure><img src="image_141.png" /></figure>

80% in cross-validation

<figure><img src="image_142.png" /></figure>

Training: 20%, 50% or 70%.

<figure><img src="image_143.png" /></figure>

Validation: 80%, 50%, 70%.

**Just as a reminder** : The phenotypes of the individuals in the validation and test population are masked and the model must predict them.

<figure><img src="image_144.png" /></figure>
GS example

We will look at the quality of the prediction with the following parameters:

<figure><img src="image_145.png" /></figure>

Predictive ability

<figure><img src="image_146.png" /></figure>

Accuracy

<figure><img src="image_147.png" /></figure>

the Ranking

<figure><img src="image_148.png" /></figure>

the slope

10 repetitions have been performed

<figure><img src="image_149.png" /></figure>
# GS example

## The predicting ability in validation

<figure><img src="image_150.png" /></figure>

<figure><img src="image_151.png" /></figure>
# GS example

## In the test Set

<figure><img src="image_152.png" /></figure>

<figure><img src="image_153.png" /></figure>
# GS example

## The ranking in validation

<figure><img src="image_154.png" /></figure>

<figure><img src="image_155.png" /></figure>
GS example

Then in Test Set

<figure><img src="image_156.png" /></figure>

<figure><img src="image_157.png" /></figure>
# GS example

## The accuracy in validation

<figure><img src="image_158.png" /></figure>

<figure><img src="image_159.png" /></figure>
# GS example

## Then in Test Set

<figure><img src="image_160.png" /></figure>

<figure><img src="image_161.png" /></figure>
# GS example

## The slope in validation

<figure><img src="image_162.png" /></figure>

<figure><img src="image_163.png" /></figure>
# GS example

## Then in Test Set

<figure><img src="image_164.png" /></figure>

<figure><img src="image_165.png" /></figure>
To go further

This example is very simple, to go further in your analysis you can :

Test several different marker sets

selected according to the LD

Around the QTLs

According to the recombination rate

Optimising the choice of individuals of the training population

Using non infinitesimal models

Using GWAS data in your predictions

Multiple-trait

...

<figure><img src="image_166.png" /></figure>
R, Rstudio and R Markdown

<figure><img src="image_167.png" /></figure>
Quick introduction

In recent years R has become the statistical programming language :

+ statisticians

+ the most widely used generic environment for analysis of high-throughput genomic
data.

R's core strength are :

+ the literally thousands of packages freely available

+ There is a good chance that for any given task there already is a package that can
do the job at hand

<figure><img src="image_168.png" /></figure>
Quick introduction

More specifically for genome-wide association studies (GWAS) or genomic Selection (GS), there are hundreds of packages available for the various analytical steps.

There are packages for :

* importing a wide range of data formats,

* preprocessing data,

* to perform quality control tests,

* to run the analysis per se,

* A large number of new algorithms and methods are published and at the same time released as an R package,

<figure><img src="image_169.png" /></figure>
Quick introduction

R is **free** (released under the GNU public license).

It is free to use, free to modify, and an open source.

R is also **platform independent**.

Scripts will *generally* run on any operating system without changes, only a small number of packages are not available on all platforms.

Since R is a scripted language it is very easy to essentially assemble various packages, add some personalized routines, and chain-link it all into a full analysis pipeline all the way from raw data to final report.

<figure><img src="image_170.png" /></figure>
<figure><img src="image_171.png" /></figure>

<figure><img src="image_172.png" /></figure>

<figure><img src="image_173.png" /></figure>
<figure><img src="image_174.png" /></figure>

<figure><img src="image_175.png" /></figure>

<figure><img src="image_176.png" /></figure>
R markdown

```r
## Why prepare and Observe data ?

In genome-wide association studies :
* Determine the **relationships between phenotypes and genotypes**
* Rare phenotypes and their associated genotypes are of great interest

**Rare phenotypes** $ \$\text{Rightarrow}$ the result of technical errors or bad models

**False Rare phenotypes** $ \$\text{Rightarrow}$ leads to the detection of false QTLs or inaccurate prediction

To avoid this $ \$\text{Rightarrow}$ bring the same *rigour* to the analysis of phenotypic and genotypic data.

## An overview of the process

The first step is data observation.

Let us take the following example: Eucleg's Soy data for 1 trial over 2 years.

We will apply the following process steps

[[First step] (Processi.png) width=70%]

## Phenotypic Data observation

```r
""[r data-server-loading-vcf]

genome <- read("DataSet/euclegnewblasted sorted.vcf", stringsAsFactors = F) read the file and load it
ped <- read("DataSet/Pedigree Soja.txt", header=T, na strings = "")

pheno18 <- read("DataSet/IFVCONS 2018.txt", header = T, na strings = "NA", dec = ",", stringsAsFactors = F)
pheno18 Trial <- "IFVCONS" # Create a column and fill it with the character string "IFVCONS"
pheno18 Year <- 2018 # Create a column and fill it with 2018

pheno19 <- read("DataSet/IFVCONS 2019.txt", header = T, na strings = "NA", dec = ",", stringsAsFactors = F)
pheno19 Trial <- "IFVCONS"
pheno19 Year <- 2019

genoet struc : read("DataSet/Genetic structure.txt", header = T, na strings = "NA", stringsAsFactors = F)

phenoSoy <- rbind(pheno18, pheno19) #Concatenates data by lines

rm(pheno18, pheno19)
```

<figure><img src="image_177.png" /></figure>
Overview of the practical part

<figure><img src="image_178.png" /></figure>
# Presentation of the dataset

Soybean Dataset field experiment of IFVCNS in 2018 and 2019

There is 15 columns, including :

* Field information : Plot, Raow, Colum, Trial, Year

* Genetic information : EntryNo, EUCLEGID, Accessionname

* Trait of interest : Seed yield, Protein content

* Explanatory variables : R1 (beginning of flowering), R2 (full flowering), R5 (start of seed filling), R8 (maturity)

The genomic data :

* 360 phenotyped individuals

* 357 genotyped individuals

* 226798 SNPs

<figure><img src="image_179.png" /></figure>
# What are we going to do in this workshop?

* Observation and understanding of data

* Dataset Visualization and Preparation

* Relationship matrix construction

* NRM matrice

* GRM matrice

* Hybrid matrice

* Heritability estimation

* Estimation of genetic correlation with a multiple-trait model

* Phenotypic adjustment

* GWAS

* Genomic Selection

* GBLUP

* Q-GBLUP

<figure><img src="image_180.png" /></figure>
# How the practical part will take place ?

The different steps :

1. download the following files from the platform:

<figure><img src="image_181.png" /></figure>

complete dataset file

<figure><img src="image_182.png" /></figure>

the slides

<figure><img src="image_183.png" /></figure>

the workshop file

2. Go to the Rstudio de progeno (procedure in the following slides)

3. Create a working folder

4. Download the data from your computer to progeno or use the database to extract and import them into Rstudio.

5. Start the tutorial

6. Try on your own data or use the sample data

<figure><img src="image_184.png" /></figure>
Pratical part

<figure><img src="image_185.png" /></figure>
step 2

<figure><img src="image_186.png" /></figure>

<figure><img src="image_187.png" /></figure>

<figure><img src="image_188.png" /></figure>
Step 3

<figure><img src="image_189.png" /></figure>

Click here

6

<figure><img src="image_190.png" /></figure>

<figure><img src="image_191.png" /></figure>

Click here

Then here

Click
here

<figure><img src="image_192.png" /></figure>

<figure><img src="image_193.png" /></figure>
Step 4

<figure><img src="image_194.png" /></figure>

Click here

<figure><img src="image_195.png" /></figure>

Click here and
chose the
download files

<figure><img src="image_196.png" /></figure>
Acknowledgements

<figure><img src="image_197.png" /></figure>
Acknowledgements

* Bernadette Julier and Philippe Barre

* Hilde Muylle, Jonas Asper and Aamir Saleem

* Steven Maenhout

* Cloe Paul-Victor

<figure><img src="image_198.png" /></figure>
Let's GO !

<figure><img src="image_199.png" /></figure>
just few last point

If you need help please email me : marie.pegard@inrae.fr

I may have forgot tu write some explaination to understand how works a function use
?function_name (example ?summary )

Take your time to test and look the results. Even if you don't have the time to do everything !

At the end I will give you the scripts that I used to generate the pdf and the slides.

I will also provide a script with “all-inclusive” function to performe genomic prediction.

<figure><img src="image_200.png" /></figure>