Back to index.html.

GWAS using the ssGBLUP framework

The single-step GBLUP methodology has been extended to support genome-wide association studies (GWAS) as a post-mixed-model analysis. This is an iterative procedure with two steps: (1) prediction of GEBVs with ssGBLUP and (2) prediction of SNP effects based on the GEBVs. Here, we call this method ssGWAS. Wang et al. (2012) described the detailed algorithm, and POSTGSF90 implements this method along with additional procedures.

Algorithm

To describe the ssGWAS algorithm, we use the notation from Wang et al. (2012). The notation may differ from that used in previous chapters. Do not confuse these symbols with the symbols used earlier.

We assume that \(\mathbf{a}\) is a vector of breeding values for genotyped animals and that \(\mathbf{u}\) is a vector of SNP marker effects. They are related by the following equation. \[ \mathbf{a}=\mathbf{Zu} \] Their variances can be written as \[ \mathrm{var}(\mathbf{a}) = \mathbf{G}\sigma_a^2 \quad\text{and}\quad \mathrm{var}(\mathbf{Zu}) = \mathbf{ZDZ}'\sigma_u^2 \] where \(\mathbf{D}\) is a diagonal matrix of weights accounting for the variances of SNP markers. This matrix is usually assumed to be \(\mathbf{I}\) in regular ssGBLUP. The two variances above are identical, so we can derive \[ \mathrm{var}(\mathbf{a}) = \mathrm{var}(\mathbf{Zu}) = \mathbf{G}\sigma_a^2 = \mathbf{ZDZ}'\sigma_u^2 \] and this formula implies \[ \mathbf{G}=\mathbf{ZDZ}'\frac{\sigma_u^2}{\sigma_a^2}=\mathbf{ZDZ}'\lambda \] where \(\lambda=\sigma_u^2/\sigma_a^2\). According to the definition of \(\mathbf{G}\), the variance ratio can be \[ \lambda=\frac{\sigma_u^2}{\sigma_a^2}=\frac{1}{2\sum_j p_j(1-p_j)}. \] The prediction of breeding values, \(\hat{\mathbf{a}}\), is calculated with ssGBLUP. The SNP effects, \(\hat{\mathbf{u}}\), are also predicted using best prediction. \[ \begin{aligned} \hat{\mathbf{u}} &= \mathrm{cov}(\mathbf{u},\mathbf{a}')\left[\mathrm{var}(\mathbf{a})\right]^{-1}\hat{\mathbf{a}}\\ &= \lambda\mathbf{DZ}'\mathbf{G}^{-1}\hat{\mathbf{a}}\\ &= \mathbf{DZ}'(\mathbf{ZDZ}')^{-1}\hat{\mathbf{a}} \end{aligned} \] Using the predicted SNP effects, we can assign weights to markers. The current default weight is defined as in Wang et al. (2012):

\[ w_i = 2 p_i q_i \hat{u}_i^2. \]

This may be used as a weight in \(\mathbf{D}=\mathrm{diag}(w_i)\) after appropriate scaling (see Wang et al. 2012). We restart the entire process with this new \(\mathbf{D}\) and can repeat the procedure until reasonable results are obtained. However, this weighting scheme is not the most accurate and can lead to unstable iterations. A well-tested weighting method is VanRaden’s “nonlinearA” (2008), which is likely to become the default. \[ w_i = 1.125^{ \frac{|\hat{u}_i|}{\sqrt{var\hat{u}}} -2} \]

The SNP effects can be used to obtain an indirect prediction of GEBVs based on \(\mathbf{a} = \mathbf{Zu}\). At the end of this section, we demonstrate this approach.

Programs

We will use several programs to perform ssGWAS.

1. RENUMF90: general preparation of the data set.
2. PREGSF90: quality control and creation of cleaned genotypes.
3. BLUPF90+: prediction of \(\mathbf{a}\) with a weighted \(\mathbf{G}\).
4. POSTGSF90: prediction of \(\mathbf{u}\) and a new weight \(\mathbf{D}\).
5. PREDF90 (optional): computation of GEBVs based on the predicted SNP effects and genotypes.

The first two programs are needed only once, during the initial preparation. If your data are well prepared, you only need to run BLUPF90+ and POSTGSF90 repeatedly. In this case, it may be easier to use BLUPF90+ directly to create \(\mathbf{G}\) and calculate the GEBV solutions in the same program.

Numerical example

Files

We will use numerical examples similar to those used in a summer course at UGA in 2012. I made some corrections to the files. The modified package is available in the author’s repository.

• lab_gwas_mod2019.zip at https://github.com/masuday/data/blob/master/tutorial/examples/genomic_analysis/

Download the zip package and extract it to your computer. It has several files as follows.

• phenotypes.txt = observation data file
• pedigree = pedigree file
• marker.geno.clean = marker file
• chrmap.txt = chromosome map file
• w = the default weight file
• renum.par = instruction file for RENUMF90

The package also has two Bash scripts to automate renumbering and ssGWAS. I do not explain how they work—you can try the scripts after reading the following sections.

• run_renum.sh = for renumbering, and for generating the parameter files
• run_ssgwas.sh = for iterative ssGWAS with BLUPF90 and POSTGSF90

The files were created on Linux (or macOS), and some Windows editors may not display their content correctly because of the line-ending format. You may see additional files in the directory, but they will not be used in this section. We will start with RENUMF90 to prepare the data set.

Preparation

First, we run RENUMF90 to generate the required files. After renumbering, we can see renf90.par, a template parameter file. We make two copies of this template: one for BLUPF90+ and one for POSTGSF90. Here, we use param_ssgwas1a.txt as the parameter file for BLUPF90+ and param_ssgwas1b.txt for POSTGSF90. The parameter files are available in the author’s GitHub repository: https://github.com/masuday/data/tree/master/tutorial.

BLUPF90+ will generate \(\mathbf{G}\) using weights supplied in a weight file (see the previous section). The ZIP package contains a weight file (w), and the option in the parameter file refers to this file. The complete parameter file for BLUPF90+ is as follows.

#
# This is for BLUPF90+
#
DATAFILE
renf90.dat
NUMBER_OF_TRAITS
          1
NUMBER_OF_EFFECTS
          2
OBSERVATION(S)
 1
WEIGHT(S)

EFFECTS:
 2      1 cross
 3  15800 cross
RANDOM_RESIDUAL VALUES
0.50000
RANDOM_GROUP
    2
RANDOM_TYPE
add_animal
FILE
renadd02.ped
(CO)VARIANCES
0.50000
OPTION SNP_file marker.geno.clean
OPTION saveGInverse
OPTION weightedG w

POSTGSF90 needs to read \(\mathbf{G}^{-1}\) generated by BLUPF90+. It also needs a chromosome map file because the program predicts each SNP effect. The following is the complete parameter file for POSTGSF90.

#
# This is for POSTGSF90
#
DATAFILE
renf90.dat
NUMBER_OF_TRAITS
           1
NUMBER_OF_EFFECTS
           2
OBSERVATION(S)
 1
WEIGHT(S)

EFFECTS:
2      1 cross
3  15800 cross
RANDOM_RESIDUAL VALUES
0.50000
RANDOM_GROUP
    2
RANDOM_TYPE
add_animal
FILE
renadd02.ped
(CO)VARIANCES
0.50000
OPTION SNP_file marker.geno.clean
OPTION readGInverse
OPTION weightedG w
OPTION chrinfo chrmap.txt

Iterative run and Manhattan plot

First, we run BLUPF90 with the parameter file param_ssgwas1a.txt, followed by POSTGSF90 with the different parameter file param_ssgwas1b.txt. POSTGSF90 creates several files.

• snp_sol = Solutions for SNP markers (that is, \(\mathbf{u}\)). It has seven columns.

  1. Trait
  2. Effect
  3. SNP number
  4. Chromosome
  5. Position
  6. SNP solution (\(\hat{u}_i\)). From this column, you can compute weights yourself.
  7. Weight (\(d_i\)) (see the official manual for what is the default method to get this value)

• chrsnp = Used in the Manhattan plot. It has six columns.

  1. Trait
  2. Effect
  3. % of variance explained by adjacent SNPs
  4. SNP number
  5. Chromosome
  6. Position

• snp_pred = Used in the program to predict GEBV based on the SNP effects. A user does not have to know the exact format.

• Sft1e2.gnuplot and Sft1e2.R = Scripts to draw Manhattan plots with Gnuplot and R. Each script reads chrsnp.

For the next round, we need the seventh column in snp_sol. This column should be extracted and saved in the file w, which is updated with the new values. Once w is updated, rerun BLUPF90+ and POSTGSF90. Repeat the process as needed. In this example, we repeat it three times.

This iterative process can be automated using a script. On Linux and macOS, Bash shell scripts are commonly used. On Windows, a traditional batch file is useful, but additional software may be needed to extract the seventh column of the file. We do not explain these techniques here. You can find Bash scripts in the same directory as the data files in the ZIP package. The scripts show you how to run the programs repeatedly.

After 3 rounds, we can draw the Manhattan plot with the script. Using R, you will see the following figure. According to this graph, there is a large SNP effect on chromosome 1.

Here, we arbitrarily repeat the process three times. Although two or three iterations may be sufficient, the optimal number of iterations for prediction depends on the data. In some cases, only the first round (no iterations) may produce the best result. Find a suitable procedure for your data.

Indirect prediction of GEBV based on SNP effects

Once you find the SNP effects, another program PREDF90 can calculate indirect GEBV from the SNP effects based on the equation \(\mathbf{\hat{a}} = \mathbf{Z\hat{u}}\). This program reads the marker file and snp_pred, generated with POSTGSF90. Run the program, and see the following message.

 name of genotype file?

Unlike other programs, it requires the name of the marker file. Type the name of the genotype file here. The program immediately calculates the indirect GEBVs for genotyped animals.

PREDF90 creates a file of SNP predictions. It has three columns.

1. Animal ID
2. Call rate
3. Indirect prediction of GEBV (based on \(\mathbf{\hat{a}} = \mathbf{Z\hat{u}}\)).

In other words, the prediction for an animal is the sum of the estimated SNP effects, weighted according to its genotypes, from the overall ssGBLUP analyses. This approach is especially useful when there are many young genotyped animals but many of them should not be included in the ssGBLUP analysis because they are part of an interim evaluation between two full genomic evaluations or are excluded for other reasons.

Back to index.html.