Example simulated data
As a first example for performing IBD calculations and QTL mapping for a multi-parent population we use a relatively simple simulated data set. The example contains simulated data for two F4DH populations. The population type F4DH is cross between two parents, followed by 3 generations of selfings, followed by a DH generation, see statgenIBD for details.
For the first population the parents where A and B, for the second the parents where A and C. The first population consists of 100 individuals, the second of 80 individuals. This is a simple example of a NAM population, having parent A as central parent. The data is simulated with three QTLs, on chromosome 1, 2, and 3. All necessary data for this population is available in the package.
Before doing QTL detection we compute IBD probabilities on a grid of
positions along the genome. This can be done using the
calcIBDMPP function in the statgenMPP package.
To perform IBD calculations a marker file is required for each of the
populations. These files should be a tab-delimited file with first
column ID identifying the genotype. The following columns should contain
marker information. The first rows should contain the parents used in
the cross. As an example, the file for the first cross, AxB, starts like
this:
| ID | M1_1 | M1_2 | M1_3 | M1_4 |
|---|---|---|---|---|
| A | 1 | 2 | 2 | 2 |
| B | 2 | 2 | 2 | 1 |
| AxB0001 | 1 | 2 | 2 | 1 |
| AxB0002 | 2 | 2 | 2 | 2 |
In this example markers M1_1 and M1_4 are segregating in the AxB
cross, M1_2 and M1_3 not.
A map file is also required. This should also be a tab-delimited file
with three columns, “marker name”, “chromosome” and “position”. The map
file cannot contain a header and has to be identical for all
crosses.
Phenotypic data can be added as a data.frame when
computing IBD probabilities. Such a data.frame should have
a first column “genotype” and all other columns have to be numerical.
For the simulated NAM population the data.frame with
phenotypic data starts like this:
| genotype | yield |
|---|---|
| AxB0001 | 9.89 |
| AxB0002 | 6.55 |
| AxB0003 | 7.90 |
| AxB0004 | 4.46 |
| AxB0005 | 5.21 |
| AxB0006 | 5.27 |
The phenotypic data only contains a single trait, yield. When
performing IBD calculations and specifying phenotypic data, the
phenotypic data will be combined with computed IBD probabilities. For
this, the genotype specified in the ID column in the marker file(s) will
be matched with the genotype in the genotype column in the
data.frame with phenotypic information. Phenotypic data for
all crosses can either be added from a single data.frame
containing phenotypic data for all crosses, or a list of
data.frames each containing phenotypic data for a single
cross.
## Specify files containing markers.
# One file for each of the two crosses.
markerFiles <- c(system.file("extdata/multipop", "AxB.txt",
package = "statgenMPP"),
system.file("extdata/multipop", "AxC.txt",
package = "statgenMPP"))
## Specify file containing map.
# Both crosses use the same map file.
mapFile <- system.file("extdata/multipop", "mapfile.txt",
package = "statgenMPP")
## Read phenotypic data
phenoDat <- read.delim(system.file("extdata/multipop", "AxBxCpheno.txt",
package = "statgenMPP"))
# Check contents.
head(phenoDat)
#> genotype yield
#> 1 AxB0001 9.89
#> 2 AxB0002 6.55
#> 3 AxB0003 7.90
#> 4 AxB0004 4.46
#> 5 AxB0005 5.21
#> 6 AxB0006 5.27
## Perform IBD calculations.
ABCMPP <- calcIBDMPP(crossNames = c("AxB", "AxC"),
markerFiles = markerFiles,
pheno = phenoDat,
popType = "F4DH",
mapFile = mapFile,
evalDist = 5)With calcIBDMPP IBD probabilities are computed on a grid
for each of the crosses separately and then combined into a single
output object. The value of evalDist can be used to specify
the maximum distance (in cM) between two evaluation points on the grid.
The exact distance depends on the length of the chromosomes. The output
is stored in an object of class gDataMPP (genomic Data for
Multi Parent Populations) in which information about map, markers, and
phenotypic data is combined. With the summary function we
can get some insight in this information.
## Print summary
summary(ABCMPP)
#> map
#> Number of markers: 95
#> Number of chromosomes: 5
#>
#> markers
#> Number of markers: 95
#> Number of genotypes: 180
#> Parents: A, B, C
#> pheno
#> Number of traits: 1
#> Traitnames: yield
#> Number of genotypes: 180
#>
#> crosses
#> AxB:100
#> AxC: 80To get a further idea about the population and the computed IBD
probabilities we can visualize the results. First we have a look at
structure of the pedigree of the population using
plotType = "pedigree" to get a general idea of what the
design looks like.
Next we look at the genetic map using
plotType = "genMap". This will display the genetic map of
the population showing the length of each of the chromosomes and
indicating the positions where the IBD probabilities were calculated.
Optionally it is possible to highlight one or more markers using
highlight argument.
## Plot genetic map.
# Highlight marker on chromosome 3 at position 40.
plot(ABCMPP, plotType = "genMap", highlight = "EXT_3_40")
Finally we visualize the computed IBD probabilities across the genome
for a selected genotype using plotType = "singleGeno". This
plot will show the IBD probabilities for all parents for all positions
on the genome for the selected genotype.
## Plot IBD probabilities for genotype AxB0001.
plot(ABCMPP, plotType = "singleGeno", genotype = "AxB0001")
Using the computed IBD probabilities we can now do the actual QTL
Mapping using the selQTLMPP function. First we perform a
SQM by setting maxCofactors = 0. Our trait of interest,
yield, is specified in the trait argument of the
function.
## Perform Single-QTL Mapping.
ABCSQM <- selQTLMPP(MPPobj = ABCMPP,
trait = "yield",
maxCofactors = 0)The results of the SQM can be plotted using the plot
function with plotType = "QTLProfile".
Already from the SQM the position of the three simulated QTLs is
quite clear. We can get an even better result using MQM. For that, again
we use the selQTLMPP function. We don’t specify the
maxCofactors to let the algorithm determine the number of
cofactors based on the threshold and QTLwindow. As long as new markers
are found with a \(-log10(p)\) value
above threshold and the maximum number of cofactors is not
reached, a new round of scanning is done. In the new round of scanning
the marker with the highest \(-log10(p)\) value from the previous round
is added to the cofactors. Based on the profile plot for the SQM the
threshold is set to 3. The QTLwindow is not
specified and therefore left at its default value of 10cM.
It is possible to include a polygenic term in the model to control
for background genetic information. To do this a list of chromosome
specific kinship matrices can either be computed by the
selQTLMPP function or specified by the user. In the former
case it is enough to specify computeKin = TRUE. When doing
so the kinship matrix is computed by averaging \(Z Z^t\) over all markers, where \(Z\) is the genotype \(\times\) parents matrix with IBD
probabilities for the marker. A list of precomputed chromosome specific
kinship matrices can be specified in K. Note that adding a
kinship matrix to the model increases the computation time a lot,
especially for large populations. It is advisable to use parallel
computing when doing so (see Parallel
computing).
## Perform Multi-QTL Mapping.
# Compute kinship matrices.
ABCMQM_kin <- selQTLMPP(MPPobj = ABCMPP,
trait = "yield",
threshold = 3,
computeKin = TRUE)Next we can visualize the results. First we plot the positions of the
QTLs found on the genetic map. This will produce a plot similar the the
genetic map plot we have seen before, but now the QTLs will be
highlighted. This plot can be made by specifying
plotType = "QTLRegion"
As for the SQM we can also plot the QTL profile. The QTLs found will be highlighted in the profile in red.
It is also possible to plot the size of the parental effects for each
of the QTLs found. Positive effects of a parent on the trait will be
indicated by shades of red, negative effects by shades of blue. The
stronger the color, the larger the effect for the specific parent is.
This plot can be made using plotType = "parEffs".
Also a combined plot of the QTL profile and the parental can be made.
In this plot the two previous plots are plotted above each other with
the chromosomes and positions aligned to allow for easily getting an
overview of which effect belongs to which QTL in the QTL profile. This
plot can be made using plotType = "QTLProfileExt".
Finally a plot of the confidence intervals around the estimated
parental effects can be made for each of the QTLs found. This plot can
be made using plotType = "parCIs". The plot will show per
QTL the effects for each of the parents, as well as the confidence
intervals around those effects (calculated as \(2 \times se\)). Confidence intervals that
don’t contain the origin, are assumed to be significant and are shown in
red for positive effects and in blue for negative effects. All other
intervals are shown in black.
A summary of the QTL-analysis gives a short overview containing the total number of markers and the number of QTLs found. Also for all QTL their position on the chromosome is shown as well as the nearest marker on the original map, the explained variance and the effects and the standard errors of all parents.
## Print summary
summary(ABCMQM)
#> Trait analysed: yield
#>
#> Data are available for 95 markers.
#> Threshold: 3
#>
#> Number of QTLs: 3
#>
#> evalPos chr pos mrkNear minlog10p varExpl eff_A eff_B eff_C se_eff_A se_eff_B
#> EXT_1_25 1 25 M1_3 16.45 0.242 -1.417 0.859 0.558 0.145 0.239
#> EXT_2_75 2 75 M2_8 8.59 0.217 0.833 -1.327 0.495 0.166 0.238
#> EXT_3_65 3 65 M3_7 11.69 0.287 0.984 0.560 -1.544 0.139 0.219
#> se_eff_C
#> 0.212
#> 0.274
#> 0.214From the output of selQTLMPP the p-Values and effects
for all markers can be extracted. They are stored in a
data.table within the output object. The example below
shows how to extract them. The output will contain columns evalPos, chr
and pos with name, chromosome number and evaluation position and columns
pValue and eff_par1, eff_par2, eff_par.. with the effects of all parents
for that marker.
## Extract results of QTL mapping.
ABCMQMres <- ABCMQM$GWAResult$yield
head(ABCMQMres[, 1:8])
#> NULLIt is also possible to only extract the markers that are either QTLs
or within the window of one of the selected QTLs, as specified when
calling the selQTLMPP function. As for the full results,
this information in stored in the output as a data.table
that can be extracted as shown below. The columns in this
data.table are identical to those in the full results
except for an additional column at the end, snpStatus, that shows
whether a marker is a QTL or within the window of a QTL.
