---
title: "Linkage map construction with ASMap"
author: "Julian Taylor and David Butler"
bibliography: ref.bib
link-citations: true
nocite: |
  @tb17
output:
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 2
    mathjax: null
vignette: >
  %\VignetteIndexEntry{Linkage map construction with ASMap}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  comment = NA,
  collapse = FALSE,
  fig.width = 7,
  fig.height = 4.5,
  fig.align = "center",
  out.width = "100%",
  dpi = 96
)
library(ASMap)
```

This vignette introduces the construction of a linkage map with ASMap, the
arguments that most influence the outcome, and the facilities provided for
assessing the result. It is intended to be read once, in order, by someone
meeting the package for the first time.

> **Scope.** This is a short introduction. Complete documentation — nine
> further articles covering every function, the technical basis of the MSTmap
> algorithm, and a worked example carrying an unconstructed 3,023 marker barley
> backcross through to a diagnosed seven linkage group map — is published at
> <https://drj001.github.io/ASMap/> and listed at the end of this document.
> Individual functions are documented in the help pages.

Construction is performed by the MSTmap algorithm [@mst08], which clusters
markers into linkage groups and determines an optimal order within each group in
a single pass. This distinguishes it from the two-stage procedures common
elsewhere, in which a non-exhaustive initial ordering is followed by an
exhaustive search within a sliding marker window. The returned object is an
R/qtl `cross` object [@br14], so the facilities of that package remain available
throughout.

## The package in outline

The functions divide into three stages of a map construction workflow.

| Stage | Functions | Purpose |
|---|---|---|
| Pre-construction | `statGen()`, `profileGen()`, `genClones()`, `fixClones()`, `statMark()`, `profileMark()`, `pValue()`, `pullCross()` | Diagnose and curate genotypes and markers before a map is built, and set problematic markers aside |
| Construction | `mstmap()` | Cluster markers into linkage groups, order them, and estimate genetic distances |
| Post-construction | `pushCross()`, `heatMap()`, `alignCross()`, `breakCross()`, `mergeCross()`, `combineMap()`, `subsetCross()`, `quickEst()` | Reintroduce markers, display and interpret the map, and manipulate its structure |

The two pre- and post-construction stages are what distinguish a usable workflow
from a single call to an ordering algorithm. Poor genotypes and uninformative
markers degrade an order, and are more cheaply removed before construction than
diagnosed afterwards.

## Marker data

The package supplies five datasets. `mapDHf` is an unconstructed doubled haploid
wheat marker set, in the layout the construction function requires: **markers in
rows, genotypes in columns**, with marker names in the `rownames` and genotype
names in the `names`.

> The datasets are not lazy-loaded, and therefore require an explicit `data()`
> call. Its omission results in `object 'mapDHf' not found`.

```{r look}
data(mapDHf)

dim(mapDHf)
mapDHf[1:5, 1:6]
```

## Construction

A single call clusters the markers into linkage groups, orders the markers
within each group, and estimates the genetic distances between them.

```{r build}
map <- mstmap(mapDHf, pop.type = "DH", dist.fun = "kosambi", p.value = 1e-06)

nchr(map)
round(chrlen(map), 1)
```

The linkage group summary printed above is emitted by the MSTmap algorithm
itself rather than by R, and appears independently of the `trace` argument.

```{r plotmap, fig.height = 5}
plotMap(map)
```

Where the markers already reside in a `cross` object, the same generic accepts
it directly. Two combinations of arguments recur throughout the documentation:

```r
# cluster and order from scratch, disregarding existing linkage groups
mstmap(cross, bychr = FALSE, p.value = 1e-12)

# re-order markers within existing linkage groups, without splitting them
mstmap(cross, bychr = TRUE, p.value = 2)
```

A `p.value` greater than unity suppresses clustering entirely, which is the
mechanism underlying the second of these.

## The choice of `p.value`

The `p.value` argument sets the threshold at which markers are separated into
distinct linkage groups, and is the single argument that most affects the
result. Its appropriate value depends strongly on the size of the population,
larger populations requiring a smaller value. The relationship may be inspected
directly, so that the choice need not be made by trial alone.

```{r pvalue, fig.height = 5}
pValue(dist = c(25, 30, 35, 40), pop.size = 100:500, map.function = "kosambi")
```

For a population of approximately 300 individuals, a `p.value` of `1e-12`
corresponds to a threshold of some 30 cM before two markers are placed in a
common linkage group.

## Diagnosing genotypes

The diagnostic functions operate on a `cross` object, and are therefore applied
either to a map obtained elsewhere or to an existing map that is about to be
reconstructed. `statGen()` summarises each genotype by its number of
crossovers, its number of double crossovers, and its proportion of missing
values.

```{r statgen}
data(mapDH)
sg <- statGen(mapDH, bychr = FALSE, stat.type = c("xo", "dxo", "miss"))

range(sg$xo)
sort(sg$dxo, decreasing = TRUE)[1:5]
```

A double crossover across a short interval is improbable under normal
recombination, so genotypes accumulating them are candidates for genotyping
error. One genotype, `DH208`, carries seven.

`profileGen()` presents the same statistics graphically, which exposes an
isolated offender more readily than a sorted table does. Supplying an expected
recombination rate to `xo.lambda` marks those genotypes whose crossover count
departs significantly from it.

```{r profilegen, fig.height = 4}
profileGen(mapDH, bychr = FALSE, stat.type = c("xo", "dxo", "miss"),
           id = "Genotype", xo.lambda = 25, layout = c(1, 3), lty = 2)
```

A single genotype is flagged here, `DH218`, whose 47 crossovers are the most in
the population and well above the expectation of 25. Derivation of an
appropriate `xo.lambda` for a given population is set out in the
[diagnostics article](https://drj001.github.io/ASMap/articles/diagnostics.html).

> `profileGen()` returns its statistics invisibly, including a logical vector
> named `xo.lambda` that marks the genotypes it has flagged. This pairs directly
> with `subsetCross()`, which removes the offending individuals.

Genotypes that are unexpectedly similar, whether through inadvertent
duplication or genuine clonal relationship, are identified separately by
`genClones()` and may be combined by `fixClones()`.

## Diagnosing markers and intervals

`statMark()` summarises the markers and the intervals between them, returning a
component for each.

```{r statmark}
sm <- statMark(mapDH, stat.type = c("marker", "interval"))

colnames(sm$marker)
colnames(sm$interval)
```

`profileMark()` profiles any combination of these statistics along the linkage
groups, drawing each group in a distinct colour. Setting `crit.val = "bonf"`
annotates those markers whose segregation distortion remains significant after
a Bonferroni adjustment.

```{r profilemark, fig.height = 4.5, warning = FALSE}
profileMark(mapDH, stat.type = c("seg.dist", "prop"), crit.val = "bonf",
            layout = c(1, 3), type = "l", cex = 0.5)
```

Markers exhibiting distorted segregation or extreme allele proportions are
candidates for withdrawal before construction, which is considered below.

## Heat maps

`heatMap()` displays the pairwise estimated recombination fractions and the
pairwise LOD scores of linkage in a single figure, each with its own legend.
Consistent heat within a linkage group indicates strong linkage between
neighbouring markers; blocks of shared heat between two groups suggest groups
that may warrant merging.

```{r heat, fig.height = 7, warning = FALSE}
heatMap(mapDH, lmax = 50)
```

## Setting markers aside and reintroducing them

Uninformative markers are better withheld during construction than deleted.
`pullCross()` removes markers of a nominated type from the map but retains them
within the `cross` object, and `pushCross()` returns them afterwards. The
supported types are `"co.located"`, `"seg.distortion"` and `"missing"`.

```{r pull}
mapDH2 <- pullCross(mapDH, type = "missing", pars = list(miss.thresh = 0.03))

sum(nmar(mapDH2))
names(mapDH2)
nrow(mapDH2$missing$table)
```

Nine markers have been withdrawn to a new `missing` element of the object, from
which they can be restored once a stable map has been obtained.

> **The thresholds are deliberately reversed.** Markers are *pulled* when the
> statistic falls below the threshold and *pushed back* when it exceeds one.
> Supplying the same threshold to both calls therefore selects nothing, and
> `pushCross()` fails with `There are no markers to push back with these
> characteristics`. The threshold given to `pushCross()` must be the more
> permissive of the two.

```{r push}
mapDH3 <- pushCross(mapDH2, type = "missing", pars = list(miss.thresh = 0.1))

sum(nmar(mapDH3))
names(mapDH3)
```

All nine markers have been reinstated and the `missing` element has been
removed, leaving an object of the original size.

## Further documentation

| Article | Subject |
|---|---|
| [How the MSTmap algorithm works](https://drj001.github.io/ASMap/articles/mstmap-algorithm.html) | The minimum spanning tree formulation, clustering and marker ordering |
| [Constructing a linkage map](https://drj001.github.io/ASMap/articles/constructing-a-map.html) | Both construction functions and the full set of MSTmap parameters |
| [Pulling and pushing markers](https://drj001.github.io/ASMap/articles/pulling-and-pushing.html) | Setting problematic markers aside and reintroducing them |
| [Diagnosing genotypes and markers](https://drj001.github.io/ASMap/articles/diagnostics.html) | Profiling of individuals, markers and intervals; genetic clones |
| [Heat maps](https://drj001.github.io/ASMap/articles/heat-maps.html) | Interpretation and tuning of the recombination fraction and LOD display |
| [Manipulating linkage maps](https://drj001.github.io/ASMap/articles/manipulating-maps.html) | Breaking, merging, subsetting and combining maps |
| [Worked example I: construction](https://drj001.github.io/ASMap/articles/worked-example-construction.html) | A barley backcross, from raw markers to a constructed map |
| [Worked example II: refinement](https://drj001.github.io/ASMap/articles/worked-example-refinement.html) | Reintroduction of markers, merging of groups, fine mapping |
| [Notes on the MSTmap algorithm](https://drj001.github.io/ASMap/articles/algorithm-notes.html) | Distance calculations, `mvest.bc` and `detectBadData` |

## Citing ASMap

The package is described in the *Journal of Statistical Software* [@tb17], which
should be cited in published work. The current entry is given by
`citation("ASMap")`.

## References