2 Examples
We will demonstrate the basic functionality of the package using some example data sets. The first example will be for a real life data set for a relatively simple doubled haploid population. To highlight some specifics of more complex populations and crosses, we will use simulated data.
2.1 Steptoe Morex
The first example is the Steptoe Morex data, described by (Hayes et al. 1993). This is a population of 150 barley doubled haploid lines with parentage Steptoe / Morex. The population was developed for the North American Barley Genome Mapping Project by the Oregon State University Barley Breeding Program. For a full description see the of the program (currently offline).
The map and genotypic file for this population are included in the package.
Both the map and the genotypic file should be in a tab-delimited format. The map file should be a file without a header. It should contain three columns, corresponding to marker name, chromosome number and position on the chromosome in centiMorgan.
## Read the map and display the first rows.
map <- read.table(system.file("extdata/SxM", "SxM_map.txt", package = "statgenIBD"))
head(map)
#> V1 V2 V3
#> 1 plc 1 0.0
#> 2 glx 1 18.7
#> 3 wg789a 1 24.1
#> 4 abg380 1 33.5
#> 5 abc158 1 40.2
#> 6 ksua1a 1 44.2The genotypic file should a file with a header containing the marker names. The first column should contain the genotypes, starting with the parents, in this case Morex and Steptoe. The name of this column is irrelevant and can even be an empty string as in this example.
The parents are assumed to be inbred lines, without missing scores. For this Steptoe x Morex example, the markers for Morex are all 1, Steptoe all 2. The original SNP data can also be used. Missing marker scores in the offspring should have value “-”. More general, for populations with SNP scores a and b for the parents, the offspring can be scored as a (or a/a), b (or b/b), a/b, a/-, -/b, or -.
## Read the genotypic file and display the first rows and columns.
geno <- read.table(system.file("extdata/SxM", "SxM_geno.txt", package = "statgenIBD"),
header = TRUE)
head(geno[, 1:5])
#> plc glx wg789a abg380 abc158
#> Morex 1 1 1 1 1
#> Steptoe 2 2 2 2 2
#> dh001 2 1 1 1 1
#> dh002 2 2 2 2 2
#> dh003 2 2 2 2 2
#> dh004 2 2 2 2 1Using the map and genotypic files we can compute the IBD
probabilities using calcIBD. This function computes, per
genotype and marker, the probability that it descended from either
parent. If this information is known for the marker from the input, the
probabilities will simply be 0 for one of the parents and 1 for the
other. When a marker score is missing from the input, the probabilities
are calculated using the information of the other markers on the
chromosome.
## Compute IBD probabilities for Steptoe Morex.
SxMIBD <- calcIBD(popType = "DH",
markerFile = system.file("extdata/SxM", "SxM_geno.txt",
package = "statgenIBD"),
mapFile = system.file("extdata/SxM", "SxM_map.txt",
package = "statgenIBD"))
## Print summary.
summary(SxMIBD)
#> population type: DH
#> Number of evaluation points: 116
#> Number of individuals: 150
#> Parents: Morex SteptoeThe output of the calcIBD function is an object of class
IBDprob. This object is a list that consists
of five elements:
| element | class | description |
|---|---|---|
| map | data.frame |
The map for the population. |
| markers | array |
A three dimensional array with the IBD probabilities. The dimensions of this array are #genotypes x #markers x #parents. The array contains per combination of genotypes and marker the probabilities that it descended from either of the parents. |
| popType | character |
The population type. |
| parents | character |
The parents. |
The summary function can be used to get a short
description of the content of a IBDprob object.
Looking at genotype dh001 and marker
plc we can see from the genotypic file printed above
that it has a value of 2, equal to that of Steptoe. If we now look at
the output of calcIBD it shows a value of 1 for Steptoe and
a value of 0 for Morex.
For the combination of marker abg313b and genotype dh005 the genotypic file contained a missing value. Looking at the output we can see that the probability that this marker came from Morex is about 0.67 and the probability that it came from Steptoe is about 0.33.
2.2 Visualizing results
The calculated IBD probabilities can be visualized using the
plot function. Several types of plot can be made using this
function. Besides creating a plot, all plots also return a
ggplot object that can be further customized.
- a plot of the structure of the population for which the IBD probabilities were computed
- a plot of the genetic map
- a plot of the computed IBD probabilities per parent for a single genotype
- a plot of the computed IBD probabilities across the genome for all genotypes
- a plot of the coverage of each parent across the genome
- a plot of the total coverage of each parent
All plots are described in more detail in the sections below.
2.2.1 Pedigree plot
A plot of the structure of the population for which the IBD
probabilities were computed can be made specifying
plotType = "pedigree". This plot is fully determined by the
population type of the population for which the probabilies were
computed.
2.2.2 Genetic map plot
It is also possible to plot the genetic map by specifying
plotType = "map". This will show the length of the
chromosomes and the position of the markers on the chromosomes.
2.2.3 Single genotype plot
plotType = "singleGeno" will generate a plot of the
probabilities per parent for a selected genotype. The plot for
dh005 is made below. Note that it contains a lot of
white space since the markers in the file are spread widely along the
genome. In the next section we will show how to compute probabilities on
a denser grid leading to a plot with a lot less white space.
2.2.4 Plot of IBD probabilities for all genotypes.
When specifying plotType = "allGeno" a plot is made of
all genotypes colored per evaluation position according to the parent
with the highest probability for that combination of genotype and
position. Darker colors indicate a higher probability. Dashed lines
indicate the start of a new chromosome.
2.2.5 Plot the coverage of each parent across the genome.
When specifying plotType = "meanProbs" a plot is made of
the mean probabilities of each parent for each evaluation position. The
plot can be restricted to one or more chromosomes by specifying
chr.
2.2.6 Plot the total coverage of each parent
A plot of the total coverage of each parent can be made by specifying
plotType = "totalCovarage". This will create a bar plot of
the coverage per parent over the genome. The plot can be restricted to
one or more chromosomes by specifying chr.
2.3 Evaluation positions
When calling the calcIBD function without extra
parameters, the only positions for which the calculations are made are
the positions present in the map. There are two ways to change this. A
first option is to specify the parameter evalDist. Setting
this to a value of e.g. 5, assures that extra evaluation positions are
added in such a way that the maximum distance between evaluation points
will be 5 cM.
Adding the extra evaluation points can be done on a grid by setting
grid = TRUE. Doing this defines a grid of evaluation points
from the beginning of each chromosome to the end with a distance between
the evaluation points of evalDist. The original marker
positions will be ignored in this case. To use evaluation points on a
regular grid is important for haplotype construction and is efficient
for detection of QTLs (Van Eeuwijk et al. 2009).
## Compute IBD probabilities for Steptoe Morex.
## Add extra evaluation positions on dense grid.
SxMIBD_Ext_grid <- calcIBD(popType = "DH",
markerFile = system.file("extdata/SxM", "SxM_geno.txt",
package = "statgenIBD"),
mapFile = system.file("extdata/SxM", "SxM_map.txt",
package = "statgenIBD"),
evalDist = 1,
grid = TRUE)
## Print summary.
summary(SxMIBD_Ext_grid)
#> population type: DH
#> Number of evaluation points: 1138
#> Number of individuals: 150
#> Parents: Morex SteptoeAs the summary shows, the number of evaluation points has increased to 1138. This will also show when we visualize the results. The positions of the crossovers on the different chromosomes are clearly visible now.
When setting grid = FALSE extra evaluation positions
will be included between existing marker positions. These extra
evaluation points will be spread evenly. If we look at the map for the
Steptoe Morex example, we see that the distance between the first two
markers on chromosome 1, plc and glx,
is 18.7. To assure a maximum distance between evaluation points of 5, 3
extra evaluation points will be added. Since they will be spread evenly,
the distance between them will be 18.7/3 = 4.68. The extra evaluation
points will be named EXT_chr_pos, e.g the extra evaluation point at
position 4.68 of chromosome 1, will be named EXT_1_4.68.
## Compute IBD probabilities for Steptoe Morex.
## Add extra evaluation positions between existing markers.
SxMIBD_Ext <- calcIBD(popType = "DH",
markerFile = system.file("extdata/SxM", "SxM_geno.txt",
package = "statgenIBD"),
mapFile = system.file("extdata/SxM", "SxM_map.txt",
package = "statgenIBD"),
evalDist = 5,
grid = FALSE)
## Print summary.
summary(SxMIBD_Ext)
#> population type: DH
#> Number of evaluation points: 290
#> Number of individuals: 150
#> Parents: Morex SteptoeAs the summary shows, the number of evaluation point has increased to 290, from 116 in the original calculation with only the map positions. Looking at the first rows of the map in the output, we see that indeed 3 extra evaluation points have been added between plc and glx.
## Show first rows of map in output.
head(SxMIBD_Ext$map)
#> chr pos
#> plc 1 0.00
#> EXT_1_4.68 1 4.68
#> EXT_1_9.35 1 9.35
#> EXT_1_14.02 1 14.02
#> glx 1 18.70
#> EXT_1_21.4 1 21.40A second way of changing the position for which the computations are
made, is by specifying them in a data.frame. This
data.frame should have two columns, “chr” and “pos” and can
be specified in the calcIBD function using the
evalPos parameter.
In the package a very simple example of such a
data.frame for the Steptoe Morex data is included as a .txt
file . In it three evaluation positions are specified for each of the
chromosomes.
## Read the evalPos file and display the first rows.
evalPos <- read.table(system.file("extdata/SxM", "SxM_eval.txt", package = "statgenIBD"),
header = TRUE)
head(evalPos)
#> chr pos
#> 1 1 0
#> 2 1 20
#> 3 1 40
#> 4 2 0
#> 5 2 20
#> 6 2 40When using a data.frame with evaluation position, IBD
calculations are made only for the positions in the
data.frame. The information in the map file is still used
in those computations, but for the marker positions themselves the
evaluations are not done. The evaluation points will be named
EVAL_chr_pos, e.g the first evaluation point in the file, at position 0
of chromosome 1, will be named EVAL_1_0.
SxMIBD_evalPos <- calcIBD(popType = "DH",
markerFile = system.file("extdata/SxM", "SxM_geno.txt",
package = "statgenIBD"),
mapFile = system.file("extdata/SxM", "SxM_map.txt",
package = "statgenIBD"),
evalPos = evalPos)
## Print summary.
summary(SxMIBD_evalPos)
#> population type: DH
#> Number of evaluation points: 21
#> Number of individuals: 150
#> Parents: Morex SteptoeAs the summary shows, the number of evaluation points now decreased to 21, 3 for each of the 7 chromosomes.
In case both evalPos and evalDist are
specified the former takes prevalence and evalDist is
ignored.
2.4 Extracting value for markers of interest
Often only the subset of markers is of interest for further analysis,
e.g. for QTL Mapping of for fitting a multi-QTL model. To extract the
computed probabilities for such a marker subset from a
calcIBD object we can use the getProbs
function.
## Extract marker probabilities for markers plc and ABG053.
SxM_probs <- getProbs(SxMIBD, markers = c("plc", "ABG053"))
head(SxM_probs)
#> geno plc_Morex plc_Steptoe ABG053_Morex ABG053_Steptoe
#> 1 dh001 0 1 0 1
#> 2 dh002 0 1 1 0
#> 3 dh003 0 1 1 0
#> 4 dh004 0 1 1 0
#> 5 dh005 0 1 1 0
#> 6 dh006 0 1 1 0The output of the getProbs function is a
data.frame with a column geno containing the
genotype and for both markers the probability that it descended from
Morex and that it descended from Steptoe.
2.5 Write and read tab-delimited *.txt files
The results of IBD calculations can be written to a tab-delimited
.txt file using writeIBDs. Details about this file
format can be found in the accompanying IBD
file format vignette.
Likewise, IBD files can be read using the readIBDs
method.
2.6 Write results to Flapjack format
The results of IBD calculations can be written to Flapjack (Milne et al. 2010)
format using writeFlapjack. Two files will be written that
can be imported directly into Flapjack, a map file and a genotypic file.
For the genotypic file the probabilities from the IBD calculations are
converted to combinations of parents. If for a genotype x marker
combination the probability that it comes from a particular parent is
high enough (higher than \(0.85 / number \; of
\; parents\)), that parent is included in the output for that
genotype x marker combination.
## Write results to Flapjack format.
writeFlapjack(SxMIBD_Ext,
outFileMap = "map.txt",
outFileGeno = "geno.txt")When opened in Flapjack the results look like this:
