# Estimated QTL effects

## Overview

Teaching:10 min

Exercises:20 minQuestions

How do I find the estimated effects of a QTL on a phenotype?

Objectives

Obtain estimated QTL effects.

Plot estimated QTL effects.

The `scan1()`

function returns only LOD scores. To obtain estimated QTL effects, use the function `scan1coef()`

. This function takes a single phenotype and the genotype probabilities for a single chromosome and returns a matrix with the estimated coefficients at each putative QTL location along the chromosome.

For example, to get the estimated effects on chromosome 2 for the liver phenotype, we’d do the following:

```
c2eff <- scan1coef(pr[,"2"], iron$pheno[,"liver"])
```

The result is a matrix of 39 positions × 2 genotypes. An additional column contains the y-axis intercept values.

```
dim(c2eff)
```

```
[1] 39 3
```

```
head(c2eff)
```

```
SS SB BB
D2Mit379 86.47730 93.86458 104.9718
c2.loc39 86.34738 93.74505 105.3583
c2.loc40 86.26104 93.63473 105.6668
c2.loc41 86.21877 93.54212 105.8779
c2.loc42 86.21960 93.47498 105.9766
c2.loc43 86.26170 93.43885 105.9547
```

To plot the effects, use the function `plot_coef()`

. Use the argument `columns`

to indicate which coefficient columns to plot.

```
par(mar=c(4.1, 4.1, 1.1, 2.6), las=1)
col <- c("slateblue", "violetred", "green3")
plot(c2eff, map["2"], columns=1:3, col=col)
last_coef <- unclass(c2eff)[nrow(c2eff),] # pull out last coefficients
for(i in seq(along=last_coef))
axis(side=4, at=last_coef[i], names(last_coef)[i], tick=FALSE, col.axis=col[i])
```

The default is to provide phenotype averages for each genotype group. If instead you want additive and dominance effects, you can provide a square matrix of *contrasts*, as follows:

```
c2effB <- scan1coef(pr[,"2"], iron$pheno[,"liver"],
contrasts=cbind(mu=c(1,1,1), a=c(-1, 0, 1), d=c(-0.5, 1, -0.5)))
```

The result will then contain the estimates of `mu`

, `a`

, and `d`

.

```
dim(c2effB)
```

```
[1] 39 3
```

```
head(c2effB)
```

```
mu a d
D2Mit379 95.10455 9.247233 -1.239968
c2.loc39 95.15023 9.505449 -1.405183
c2.loc40 95.18752 9.702880 -1.552789
c2.loc41 95.21294 9.829582 -1.670822
c2.loc42 95.22372 9.878490 -1.748738
c2.loc43 95.21841 9.846482 -1.779555
```

Here’s a plot of the additive and dominance effects, which are in the second and third columns.

```
par(mar=c(4.1, 4.1, 1.1, 2.6), las=1)
plot(c2effB, map["2"], columns=2:3, col=col)
last_coef <- unclass(c2effB)[nrow(c2effB),2:3] # last two coefficients
for(i in seq(along=last_coef))
axis(side=4, at=last_coef[i], names(last_coef)[i], tick=FALSE, col.axis=col[i])
```

If you provide a kinship matrix to `scan1coef()`

, it fits a linear mixed model (LMM) to account for a residual polygenic effect. Here let’s use the kinship matrix from the LOCO method.

```
c2eff_pg <- scan1coef(pr[,"2"], iron$pheno[,"liver"], kinship_loco[["2"]])
dim(c2eff_pg)
```

```
[1] 39 3
```

```
head(c2eff_pg)
```

```
SS SB BB
D2Mit379 84.99427 93.95099 104.6940
c2.loc39 84.71083 93.87206 105.1537
c2.loc40 84.47392 93.80197 105.5385
c2.loc41 84.28851 93.74761 105.8280
c2.loc42 84.15801 93.71501 106.0058
c2.loc43 84.08450 93.70797 106.0625
```

Here’s a plot of the estimates.

```
par(mar=c(4.1, 4.1, 1.1, 2.6), las=1)
col <- c("slateblue", "violetred", "green3")
plot(c2eff_pg, map["2"], columns=1:3, col=col, ylab="Phenotype average")
last_coef <- unclass(c2eff_pg)[nrow(c2eff_pg),]
for(i in seq(along=last_coef))
axis(side=4, at=last_coef[i], names(last_coef)[i], tick=FALSE, col.axis=col[i])
```

You can also get estimated additive and dominance effects, using a matrix of contrasts.

```
c2effB_pg <- scan1coef(pr[,"2"], iron$pheno[,"liver"], kinship_loco[["2"]],
contrasts=cbind(mu=c(1,1,1), a=c(-1, 0, 1), d=c(-0.5, 1, -0.5)))
```

Here’s a plot of the results.

```
par(mar=c(4.1, 4.1, 1.1, 2.6), las=1)
plot(c2effB_pg, map["2"], columns=2:3, col=col)
last_coef <- unclass(c2effB_pg)[nrow(c2effB_pg),2:3]
for(i in seq(along=last_coef))
axis(side=4, at=last_coef[i], names(last_coef)[i], tick=FALSE, col.axis=col[i])
```

Another option for estimating the QTL effects is to treat them as random effects and calculate Best Linear Unbiased Predictors (BLUPs). This is particularly valuable for multi-parent populations such as the Collaborative Cross and Diversity Outbred mice, where the large number of possible genotypes at a QTL lead to considerable variability in the effect estimates. To calculate BLUPs, use `scan1blup()`

; it takes the same arguments as `scan1coef()`

, including
the option of a kinship matrix to account for a residual polygenic effect.

```
c2blup <- scan1blup(pr[,"2"], iron$pheno[,"liver"], kinship_loco[["2"]])
```

Here is a plot of the BLUPs (as dashed curves) alongside the standard estimates. Note that the BLUPs are centered at 0, while the coefficient estimates are centered at the phenotype average.

```
par(mar=c(4.1, 4.1, 1.1, 2.6), las=1)
col <- c("slateblue", "violetred", "green3")
ylim <- range(c(c2blup, c2eff))+c(-1,1)
plot(c2eff, map["2"], columns=1:3, col=col, ylab="Phenotype average", ylim=ylim,
xlab="Chr 2 position")
plot(c2blup, map["2"], columns=1:3, col=col, add=TRUE, lty=2)
last_coef <- unclass(c2eff)[nrow(c2eff),]
for(i in seq(along=last_coef))
axis(side=4, at=last_coef[i], names(last_coef)[i], tick=FALSE, col.axis=col[i])
```

The `scan1coef`

function can also provide estimated QTL effects for binary traits, with `model="binary"`

. (However, `scan1blup`

has not yet been implemented for binary traits.)

```
c2eff_bin <- scan1coef(pr[,"2"], bin_pheno[,"liver"], model="binary")
```

Here’s a plot of the effects. They’re a bit tricky to interpret, as they are basically log odds ratios.

```
par(mar=c(4.1, 4.1, 1.1, 2.6), las=1)
col <- c("slateblue", "violetred", "green3")
plot(c2eff_bin, map["2"], columns=1:3, col=col)
last_coef <- unclass(c2eff_bin)[nrow(c2eff_bin),] # pull out last coefficients
for(i in seq(along=last_coef))
axis(side=4, at=last_coef[i], names(last_coef)[i], tick=FALSE, col.axis=col[i])
```

Finally, to plot the raw phenotypes against the genotypes at a single putative QTL position, you can use the function `plot_pxg()`

. This takes a vector of genotypes as produced by the `maxmarg()`

function, which picks the most likely genotype from a set of genotype probabilities, provided it is greater than some specified value (the argument `minprob`

). Note that the “marg” in “maxmarg” stands for “marginal”, as this function is selecting the genotype at each position that has maximum marginal probability.

For example, we could get inferred genotypes at the chr 2 QTL for the liver phenotype (at 28.6 cM) as follows:

```
g <- maxmarg(pr, map, chr=2, pos=28.6, return_char=TRUE)
```

We use `return_char=TRUE`

to have `maxmarg()`

return a vector of character strings with the genotype labels.

We then plot the liver phenotype against these genotypes as follows:

```
par(mar=c(4.1, 4.1, 0.6, 0.6))
plot_pxg(g, iron$pheno[,"liver"], ylab="Liver phenotype")
```

```
par(mar=c(4.1, 4.1, 0.6, 0.6))
plot_pxg(g, iron$pheno[,"liver"], SEmult=2, swap_axes=TRUE, xlab="Liver phenotype")
```

## Challenge 1

In the code block above, we use

`swap_axes=TRUE`

to place the phenotype on the x-axis. We can use`SEmult=2`

to include the mean ± 2 SE intervals. 1) How would you decide which chromosome to plot? Discuss this with your neighbor, then write your responses in the collaborative document. 2) Calculate and plot the best linear unbiased predictors (blups) for spleen on the chromosome of your choice.## Solution to Challenge 3

## Key Points

.

.