R/haplotype.R

In pegas: Population and Evolutionary Genetics Analysis System

## haplotype.R (2023-11-21)

##   Haplotype Extraction, Frequencies, and Networks

## Copyright 2009-2023 Emmanuel Paradis, 2013 Klaus Schliep

## This file is part of the R-package `pegas'.
## See the file ../DESCRIPTION for licensing issues.

getAllCombs <- function(n) {

    foo <- function(set, i) {
        if (j > K) return()
        for (z in set) {
            v[i] <<- z
            if (i == n) {
                res[, j] <<- v
                j <<- j + 1L

            } else {
                foo(set[set != z], i + 1)
            }
        }
    }

    res <- matrix(NA_integer_, n, K <- factorial(n))
    j <- 1L
    v <- integer(n)
    foo(1:n, 1L)
    res
}

rmst <- function(d, B = NULL, stop.criterion = NULL, iter.lim = 1000,
                 quiet = FALSE)
{
    updateMat <- function(mat, rep) {
        i <- match(mat, MAT)
        Nnew <- length(inew <- which(is.na(i)))
        if (Nnew) {
            TAB <<- c(TAB, rep(1L, Nnew))
            MAT <<- c(MAT, mat[inew])
            NOnew <<- 0L
            if (Nnew < n - 1) {
                i <- i[-inew]
                TAB[i] <<- TAB[i] + 1L
            }
        } else {
            TAB[i] <<- TAB[i] + 1L
            NOnew <<- NOnew + 1L
        }
    }

    ## do MST and return vector of mode character with the labelled-links
    foo <- function(d) {
        d <<- as.dist(d)
        rnt <<- rnt <- mst(d)
        m <- rnt[, 1:2]
        dm <- dim(m)
        mat <- attr(rnt, "labels")[m]
        dim(mat) <- dm
        mat <- apply(mat, 1, sort)
        paste(mat[1, ], mat[2, ], sep = "\r")
    }

    if (!is.matrix(d) && !inherits(d, "dist"))
        stop("'d' must be a matrix or a 'dist' object")
    D <- as.matrix(d)
    n <- nrow(D)
    randomize <- TRUE
    if (n < 6 && is.null(B)) {
        B <- factorial(n)
        COMBS <- getAllCombs(n)
        randomize <- FALSE
        if (!quiet)
            cat("Less than 6 observations: doing complete enumeration.\n")
    }
    MAT <- character()
    TAB <- integer()
    NOnew <- 0L # number of iterations without new links
    if (is.null(B)) {
        rep <- 1L
        if (is.null(stop.criterion)) stop.criterion <- ceiling(sqrt(n))
        repeat {
            if (rep >= iter.lim) {
                B <- rep
                break
            }
            ri <- sample(n)
            d <- D[ri, ri]
            mat <- foo(d)
            updateMat(mat, rep)
            if (!quiet)
                cat("\rIteration:", rep, "  Number of links:", length(MAT))
            if (NOnew >= stop.criterion)  {
                B <- rep
                break
            }
            rep <- rep + 1L
        }
    } else {
        for (rep in 1:B) {
            if (!quiet) cat("\rIteration:", rep)
            ri <- if (randomize) sample(n) else COMBS[, rep]
            d <- D[ri, ri]
            mat <- foo(d)
            updateMat(mat, rep)
        }
    }
    if (!quiet) cat("\n")
    if (length(TAB) == n - 1) {
        attr(rnt, "weight") <- rep(B, n - 1)
        return(rnt)
    }
    links <- names(TAB) <- MAT

    ## define the backbone-network:
    labs <- rownames(D)
    forest <- 1:n
    k <- 0L
    done <- FALSE
    alt.logi <- !logical(length(TAB))
    ## the backbone-network is built by checking the links in the order
    ## of pairs of observations, and including the links which were
    ## found during the iterations of RMST.
    for (i in 1:(n - 1)) {
        labi <- labs[i]
        for (j in (i + 1):n) {
            tmp <- paste(sort(c(labi, labs[j])), collapse = "\r")
            m <- match(tmp, links)
            if (is.na(m)) next
            if (forest[i] == forest[j]) next
            alt.logi[m] <- FALSE # so this is not an alternative link
            k <- k + 1L
            rnt[k, 1L] <- i
            rnt[k, 2L] <- j
            if (k == n - 1) {
                done <- TRUE
                break
            }
            forest[forest == forest[j]] <- forest[i]
        }
        if (done) break
    }
    rnt[, 3L] <- D[rnt[, 1:2]]

    ## An alternative is to use the last MST and recode the haplotypes
    ## in the same order than in the input distance matrix:
    ##labs <- rownames(D)
    ##recode <- match(attr(rnt, "labels"), labs)
    ##rnt[, 1L] <- recode[rnt[, 1L]]
    ##rnt[, 2L] <- recode[rnt[, 2L]]

    attr(rnt, "labels") <- labs
    alt <- unlist(strsplit(links[alt.logi], "\r"))
    alt <- match(alt, labs)
    alt <- matrix(alt, length(alt)/2, 2, byrow = TRUE)
    alt <- t(apply(alt, 1, sort))
    alt <- cbind(alt, D[alt])
    colnames(alt) <- c("", "", "step")
    attr(rnt, "alter.links") <- alt
    names(TAB) <- gsub("\r", "--", names(TAB))
    attr(rnt, "weight") <- TAB
    rnt
}

msn <- function(d)
{
    getIandJ <- function(ij, n) {
        ## assumes a lower triangle, so i > j
        ## n must be > 1 (not checked)
        ## ij must be <= (n - 1)*n/2 (not checked too)
        j <- 1L
        N <- n - 1L
        while (ij > N) {
            j <- j + 1L
            N <- N + n - j
        }
        i <- n - (N - ij)
        c(j, i) # return the smaller index first
    }
    if (is.matrix(d)) d <- as.dist(d)
    MST <- mst(d)
    n <- attr(d, "Size")
    if (n < 3) {
        warning("less than 3 observations: function mst() was used")
        return(MST)
    }
    m <- matrix(NA_real_, 0, 3)
    forest <- 1:n
    ud <- sort(unique(d)) # unique distances in increasing order

    i <- 1L
    while (length(unique(forest)) > 1) {
        delta <- ud[i]
        for (iud in which(d == delta)) {
            p <- getIandJ(iud, n)
            f1 <- forest[p[1L]]
            f2 <- forest[p[2L]]
            m <- rbind(m, c(p, delta))
            if (f1 != f2) forest[forest == f2] <- f1
        }
        i <- i + 1L
    }
    ## define the alternative links wrt the MST:
    MST.str <- paste(MST[, 1], MST[, 2], sep = "\r")
    m.str <- paste(m[, 1], m[, 2], sep = "\r")
    k <- match(m.str, MST.str)
    nak <- is.na(k)
    if (any(nak)) {
        alt <- m[nak, , drop = FALSE]
        colnames(alt) <- c("", "", "step")
        m <- m[!nak, , drop = FALSE]
        attr(m, "alter.links") <- alt
    }
    colnames(m) <- c("", "", "step")
    attr(m, "labels") <- attr(d, "Labels")
    class(m) <- "haploNet"
    m
}

mst <- function(d)
{
    if (is.matrix(d)) d <- as.dist(d)
    n <- attr(d, "Size")
    if (n < 2) stop("less than 2 observations in distance matrix")
    m <- .Call(mst_C, d, n)
    colnames(m) <- c("", "", "step")
    attr(m, "labels") <- attr(d, "Labels")
    class(m) <- "haploNet"
    m
}

haplotype <- function(x, ...) UseMethod("haplotype")

haplotype.DNAbin <- function(x, labels = NULL, strict = FALSE, trailingGapsAsN = TRUE, ...)
{
    nms.x <- deparse(substitute(x))
    if (is.list(x)) x <- as.matrix(x)
    n <- nrow(x)
    if (!n) {
        class(x) <- c("haplotype", "DNAbin")
        attr(x, "index") <- list()
        attr(x, "from") <- nms.x
        warning("empty DNAbin object")
        return(x)
    }
    if (trailingGapsAsN) x <- latag2n(x)
    segs <- seg.sites(x, strict = strict, trailingGapsAsN = FALSE)
    if (length(segs)) {
        xb <- x[, segs]
        s <- ncol(xb)
        res <- .C(haplotype_DNAbin, xb, n, s, integer(n), c(0L, 0L),
                  as.integer(strict), NAOK = TRUE)
        if (res[[5]][1]) warning("some sequences of different lengths were assigned to the same haplotype")
        if (res[[5]][2]) warning("some sequences were not assigned to the same haplotype because of ambiguities")
        h <- res[[4]]
        u <- h == 0
        h[u] <- i <- which(u)
    } else {
        i <- 1L
        h <- rep.int(1L, n)
        warning("no segregating site detected with these options")
    }
    obj <- x[i, ]
    if (is.null(labels))
        labels <- as.character(as.roman(seq_along(i)))
    rownames(obj) <- labels
    class(obj) <- c("haplotype", "DNAbin")
    attr(obj, "index") <- lapply(i, function(x) which(h == x))
    attr(obj, "from") <- nms.x
    obj
}

haplotype.character <- function(x, labels = NULL, ...)
{
    nms.x <- deparse(substitute(x))
    if (!is.matrix(x))
        stop("x must be a matrix")
    h <- apply(x, 1, paste, collapse = "\r")
    hnotdup <- !duplicated(h)
    obj <- x[which(hnotdup), , drop = FALSE]
    hnotdup <- h[hnotdup]
    N <- nrow(obj)
    if (is.null(labels))
        labels <- as.character(as.roman(1:N))
    rownames(obj) <- labels
    class(obj) <- c("haplotype", "character")
    attr(obj, "index") <- lapply(1:N, function(i) which(h == hnotdup[i]))
    attr(obj, "from") <- nms.x
    obj
}

haplotype.numeric <- function(x, labels = NULL, ...)
    haplotype.character(x = x, labels = labels, ...)

haploFreq <- function(x, fac, split = "_", what = 2, haplo = NULL)
{
    if (missing(fac)) {
        fac <- strsplit(rownames(x), split)
        fac <- factor(sapply(fac, function(xx) xx[what]))
    } else {
        if (length(fac) != nrow(x))
            stop("number of elements in 'fac' not the same than number of sequences")
        if (!is.factor(fac)) fac <- factor(fac) # added 2021-02-10
    }

    if (is.null(haplo)) haplo <- haplotype(x)
    h.index <- attr(haplo, "index")
    l <- nlevels(fac)
    if (l == 1) {
        res <- sapply(h.index, length)
        dim(res) <- c(length(res), 1L)
    } else {
        res <- sapply(h.index, function(xx) tabulate(fac[xx], l))
        res <- t(res)
    }
    rownames(res) <- rownames(haplo)
    colnames(res) <- levels(fac)
    res
}

diffHaplo <- function(h, a = 1, b = 2, strict = FALSE, trailingGapsAsN = TRUE)
{
    x <- h[c(a, b), ]
    s <- seg.sites(x, strict = strict, trailingGapsAsN = trailingGapsAsN)
    data.frame(pos = s, toupper(t(as.character(x[, s]))))
}

.TempletonProb <- function(j, S, b = 2, r = 1)
{
    br <- b * r
    P <- numeric(max(j))
    L_jm <- function(q, j, m) {
        jm1 <- j - 1
        qonbr <- q/br
        (2*q)^jm1 * (1 - q)^(2*m + 1) * (1 - qonbr) *
            (2 - q*(br + 1)/br)^jm1 *
                (1 - 2*q*(1 - qonbr))
    }
    for (i in seq_along(P)) {
        M <- S - i
        denom <- integrate(L_jm, 0, 1, j = i, m = M)$value
        ## eq.7 from Templeton et al. 1992:
        out <- integrate(function(q) q*L_jm(q, j = i, m = M), 0, 1)$value/denom
        P[i] <- 1 - out
    }
    cumprod(P)[j]
}

haploNet <- function(h, d = NULL, getProb = TRUE)
{
    if (!inherits(h, "haplotype"))
        stop("'h' must be of class 'haplotype'")
    freq <- sapply(attr(h, "index"), length)
    n <- length(freq) # number of haplotypes
    link <- matrix(0, 0, 3)
    if (is.null(d)) {
        ## d <- dist.dna(h, "N", pairwise.deletion = TRUE)
        d <- .C(distDNA_pegas, h, as.integer(n), ncol(h), numeric(n * (n - 1)/2), NAOK = TRUE)[[4]]
        attr(d, "Size") <- n
        attr(d, "Labels") <- rownames(h)
        attr(d, "Diag") <- attr(d, "Upper") <- FALSE
        class(d) <- "dist"
    }
    d <- as.matrix(d)
    d[col(d) >= row(d)] <- NA # put NA's in the diag and above-diag elts
    one2n <- seq_len(n)
    dimnames(d) <- list(one2n, one2n)
    step <- 1
    gr <- one2n
    altlink <- matrix(0, 0, 3) # alternative links
    while (length(unique(gr)) > 1) {
        newLinks <- which(d == step, TRUE)
        if (length(newLinks)) {
            del <- NULL
            for (i in seq_len(nrow(newLinks))) {
                h1 <- newLinks[i, 1]
                h2 <- newLinks[i, 2]
                a <- gr[h1]
                b <- gr[h2]
                ## if both nodes are already in the same
                ## subnet, then count the link as alternative
                if (a == b) {
                    del <- c(del, i)
                    if (d[h1, h2] <= step)
                        altlink <- rbind(altlink, c(newLinks[i, ], step))
                } else gr[which(gr == b)] <- a
            }
            if (!is.null(del)) newLinks <- newLinks[-del, , drop = FALSE]
            newLinks <- cbind(newLinks, rep(step, nrow(newLinks)))
            link <- rbind(link, newLinks)
        }
        step <- step + 1
    }
    Probs <- if (getProb) .TempletonProb(link[, 3], ncol(h)) else NA
    link <- cbind(link, Probs)
    dimnames(link) <- list(NULL, c("", "", "step", "Prob"))
    attr(link, "freq") <- freq
    attr(link, "labels") <- rownames(h)
    if (nrow(altlink)) {
        Probs <- if (getProb) .TempletonProb(altlink[, 3], ncol(h)) else NA
        altlink <- cbind(altlink, Probs)
        dimnames(altlink) <- dimnames(link)
        attr(link, "alter.links") <- altlink
    }
    class(link) <- "haploNet"
    link
}

print.haploNet <- function(x, ...)
{
    cat("Haplotype network with:\n")
    cat(" ", length(attr(x, "labels")), "haplotypes\n")
    N <- n <- nrow(x)
    altlinks <- attr(x, "alter.links")
    if (!is.null(altlinks)) N <- N + nrow(altlinks)
    cat(" ", N, if (N > 1) "links\n" else "link\n")
    cat("  link lengths between", x[1, 3], "and", x[n, 3], "step(s)\n\n")
    cat("Use print.default() to display all elements.\n")
}

circle <- function(x, y, size, col = "black", pie = NULL, bg = NULL)
{
    if (is.null(pie)) {
        symbols(x, y, circles = size / 2, inches = FALSE,
                add = TRUE, fg = col, bg = bg)
    } else {
        OPTS <- get("plotHaploNetOptions", envir = .PlotHaploNetEnv)
        op <- par(fg = OPTS$pie.inner.segments.color)
        floating.pie.asp(x, y, pie, radius = size / 2, col = bg)
        par(op)
    }
}

.squarePegas <- function(x, y, size, pie = NULL)
{
    l <- size * sqrt(pi) / 4
    left <- x - l
    bottom <- y - l
    right <- x + l
    top <- y + l
    res <- list(c(left, bottom, right, top))
    if (is.null(pie)) return(res)
    n <- length(pie)
    pie <- pie / sum(pie)
    m <- matrix(NA_real_, 4, n)
    area <- l * l * 4
    TOP <- TRUE
    y2 <- top
    x1 <- left
    x2 <- right
    for (i in 1:n) {
        if (TOP) {
            y1 <- y2 - pie[i] * area / (x2 - x1)
        } else {
            x2 <- pie[i] * area / (y2 - y1) + x1
        }
        m[, i] <- c(x1, y1, x2, y2)
        if (TOP) {
            y2 <- y1
            y1 <- bottom
        } else {
            x1 <- x2
            x2 <- right
        }
        TOP <- !TOP
    }
    c(res, list(m))
}

square <- function(x, y, size, col, pie = NULL, bg = NULL)
{
    XY <- .squarePegas(x, y, size, pie = pie)
    OPTS <- get("plotHaploNetOptions", envir = .PlotHaploNetEnv)
    border <- OPTS$pie.inner.segments.color
    if (!is.null(pie)) {
        for (i in seq_along(pie)) {
            s <- XY[[2]][, i]
            rect(s[1], s[2], s[3], s[4], col = bg[i], border = border)
        }
        bg <- NULL # for the call to rect() below
    }
    s <- XY[[1]]
    rect(s[1], s[2], s[3], s[4], col = bg, border = border)
}

.rotateSquares <- function(XY)
{
    ROT <- pi / 4
    foo <- function(rec) {
        x <- rec[c(1, 3, 3, 1)] - x0
        y <- rec[c(2, 2, 4, 4)] - y0
        tmp <- rect2polar(x, y)
        o <- polar2rect(tmp$r, tmp$angle + ROT)
        o$x <- o$x + x0
        o$y <- o$y + y0
        o
    }
    rec <- XY[[1]]
    x0 <- (rec[1] + rec[3]) / 2
    y0 <- (rec[2] + rec[4]) / 2
    res <- list(foo(rec))
    if (length(XY) > 1)
        res[[2]] <- apply(XY[[2]], 2, foo)
    res
}

diamond <- function(x, y, size, col, pie = NULL, bg = NULL)
{
    XY <- .squarePegas(x, y, size, pie = pie)
    XY <- .rotateSquares(XY)
    OPTS <- get("plotHaploNetOptions", envir = .PlotHaploNetEnv)
    border <- OPTS$pie.inner.segments.color
    if (!is.null(pie)) {
        for (i in seq_along(pie)) {
            xy <- XY[[2]][[i]]
            polygon(xy$x, xy$y, col = bg[i], border = border)
        }
        bg <- NULL # for the call to polygon() below
    }
    xy <- XY[[1]]
    polygon(xy$x, xy$y, border = border, col = bg)
}

.drawSymbolsHaploNet <- function(xx, yy, size, col, bg, shape, pie)
{
    nopie <- is.null(pie)
    if (identical(shape, "circles") && nopie) {
        ## the only case with vectorised arguments
        symbols(xx, yy, circles = size / 2, inches = FALSE,
                add = TRUE, fg = col, bg = bg)
    } else {
        n <- length(xx) # nb of nodes
        size <- rep(size, length.out = n)
        col <- rep(col, length.out = n)
        bg <- if (!nopie && is.function(bg)) bg(ncol(pie)) else rep(bg, length.out = ifelse(nopie, n, ncol(pie)))
        ## 'bg' should always be a vector of colours
        shape <- rep(shape, length.out = n)
        for (i in 1:n)
            switch(shape[i],
                   "circles" = circle,
                   "squares" = square,
                   "diamonds" = diamond)(xx[i], yy[i], size[i], col[i], pie[i, ], if (nopie) bg[i] else bg)
    }
}

.drawAlternativeLinks <- function(xx, yy, altlink, threshold, show.mutation, scale.ratio, size)
{
    s <- altlink[, 3] >= threshold[1] & altlink[, 3] <= threshold[2]
    if (!any(s)) return(NULL)
    xa0 <- xx[altlink[s, 1]]
    ya0 <- yy[altlink[s, 1]]
    xa1 <- xx[altlink[s, 2]]
    ya1 <- yy[altlink[s, 2]]
    segments(xa0, ya0, xa1, ya1, col = "grey", lty = 2)
    if (show.mutation) {
        n <- length(xx)
        .labelSegmentsHaploNet(xx, yy, altlink[s, 1:2, drop = FALSE],
                               altlink[s, 3, drop = FALSE], size,
                               1, "black", show.mutation, scale.ratio)
    }
}

.mutationRug <- function(x0, y0, x1, y1, n, deltaRadii, space = 0.05, length = 0.2)
{
    ## convert inches into user-coordinates:
    l <- xinch(length)
    sp <- xinch(space)
    for (i in seq_along(x0)) {
        xstart <- seq(-(sp*(n[i] - 1))/2, by = sp, length.out = n[i])
        ystart <- rep(-l/2, n[i])
        xstop <- xstart
        ystop <- -ystart
        ## rotation if the line is not horizontal
        if (y0[i] != y1[i]) {
            theta <- atan2(y1[i] - y0[i], x1[i] - x0[i])
            tmpstart <- rect2polar(xstart, ystart)
            tmpstop <- rect2polar(xstop, ystop)
            Rstart <- tmpstart$r
            Rstop <- tmpstop$r
            THETAstart <- tmpstart$angle + theta
            THETAstop <- tmpstop$angle + theta
            xy.start <- polar2rect(Rstart, THETAstart)
            xy.stop <- polar2rect(Rstop, THETAstop)
            xstart <- xy.start$x
            ystart <- xy.start$y
            xstop <- xy.stop$x
            ystop <- xy.stop$y
        }
        ## translation:
        if (deltaRadii[i] == 0) {
            xm <- (x0[i] + x1[i]) / 2
            ym <- (y0[i] + y1[i]) / 2
           } else {
               ## see explanations below
               li <- sqrt((x0[i] - x1[i])^2 + (y0[i] - y1[i])^2)
               w0 <- 1 / (li - deltaRadii[i] / 2)
               w1 <- 1 / (li + deltaRadii[i] / 2)
               sumw <- w0 + w1
               xm <- (x0[i] * w0 + x1[i] * w1) / sumw
               ym <- (y0[i] * w0 + y1[i] * w1) / sumw
        }
        xstart <- xstart + xm
        ystart <- ystart + ym
        xstop <- xstop + xm
        ystop <- ystop + ym
        segments(xstart, ystart, xstop, ystop)
    }
}

.mutationDots <- function(x0, y0, x1, y1, n, deltaRadii, lwd, col.link,
                          space = 0.1, diameter = 0.05)
{
    ## convert inches into user-coordinates:
    l <- xinch(diameter)
    sp <- xinch(space)
    for (i in seq_along(x0)) {
        x <- seq(-(sp*(n[i] - 1))/2, by = sp, length.out = n[i])
        y <- numeric(n[i])
        ## rotation if the line is not horizontal
        if (y0[i] != y1[i]) {
            theta <- atan2(y1[i] - y0[i], x1[i] - x0[i])
            tmp <- rect2polar(x, y)
            xy <- polar2rect(tmp$r, tmp$angle + theta)
            x <- xy$x
            y <- xy$y
        }
        ## translation:
        if (deltaRadii[i] == 0) {
            xm <- (x0[i] + x1[i]) / 2
            ym <- (y0[i] + y1[i]) / 2
           } else {
               ## see explanations below
               li <- sqrt((x0[i] - x1[i])^2 + (y0[i] - y1[i])^2)
               w0 <- 1 / (li - deltaRadii[i] / 2)
               w1 <- 1 / (li + deltaRadii[i] / 2)
               sumw <- w0 + w1
               xm <- (x0[i] * w0 + x1[i] * w1) / sumw
               ym <- (y0[i] * w0 + y1[i] * w1) / sumw
        }
        x <- x + xm
        y <- y + ym
        symbols(x, y, circles = rep(lwd/15, length(x)), inches = FALSE,
                add = TRUE, fg = col.link, bg = col.link)
    }
}

.labelSegmentsHaploNet <- function(xx, yy, link, step, size, lwd, col.link, method, scale.ratio)
{
### method: the way the segments are labelled
### 1: small segments
### 2: Klaus's method
### 3: show the number of mutations on each link

    l1 <- link[, 1]
    l2 <- link[, 2]
    switch(method, {
        .mutationRug(xx[l1], yy[l1], xx[l2], yy[l2], step, size[l2] - size[l1])
    }, {
        .mutationDots(xx[l1], yy[l1], xx[l2], yy[l2], step, size[l2] - size[l1],
                      lwd, col.link, space = 0.1, diameter = 0.05)
###        ld1 <- step
###        ld2 <- step * scale.ratio
###        for (i in seq_along(ld1)) {
###            pc <- ((1:ld1[i]) / (ld1[i] + 1) * ld2[i] + size[l1[i]]/2) / (ld2[i] + (size[l1[i]] + size[l2[i]])/2)
###            xr <- pc * (xx[l2[i]] - xx[l1[i]]) +  xx[l1[i]]
###            yr <- pc * (yy[l2[i]] - yy[l1[i]]) +  yy[l1[i]]
###            symbols(xr, yr, circles = rep(lwd/15, length(xr)), inches = FALSE,
###                    add = TRUE, fg = col.link, bg = col.link)
###        }
    }, {
        x0 <- xx[l1]
        x1 <- xx[l2]
        y0 <- yy[l1]
        y1 <- yy[l2]
        x <- (x0 + x1) / 2
        y <- (y0 + y1) / 2
        deltaRadii <- size[l2] - size[l1]
        if (any(s <- deltaRadii != 0)) {
### if there are links between symbols of different size: need to place the
### annotation in the middle of the "visible" part of the link. We do that
### with a weighting mean of the coordinates (simpler to calculate than
### trogonometric formulas). We do that only for the different sized symbols.
            s <- which(s)
            ## lengths between each pair of nodes:
            l <- sqrt((x0[s] - x1[s])^2 + (y0[s] - y1[s])^2)
            ## difference in radii (keep in mind that 'size' are diameters,
            ## so we need to divide by 2):
            deltaRadii <- deltaRadii[s] / 2
            ## weights for the mean:
            w0 <- 1 / (l - deltaRadii)
            w1 <- 1 / (l + deltaRadii)
            sumw <- w0 + w1
            x[s] <- (x0[s] * w0 + x1[s] * w1) / sumw
            y[s] <- (y0[s] * w0 + y1[s] * w1) / sumw
        }
        BOTHlabels(step, NULL, x, y, c(0.5, 0.5), "rect", NULL, NULL,
                   NULL, NULL, "black", "lightgrey", FALSE, NULL, NULL)
    })
}

replot <- function(xy = NULL, col.identifier = "purple", ...)
{
    on.exit(return(list(x = xx, y = yy)))

    Last.phn <- get("last_plot.haploNet", envir = .PlotHaploNetEnv)
    xx <- Last.phn$xx
    yy <- Last.phn$yy
    l1 <- Last.phn$link[, 1]
    l2 <- Last.phn$link[, 2]
    step <- Last.phn$step
    lwd <- Last.phn$lwd
    size <- Last.phn$size
    col <- Last.phn$col
    bg <- Last.phn$bg
    lty <- Last.phn$lty
    shape <- Last.phn$shape
    col.link <- Last.phn$col.link
    labels <- Last.phn$labels
    asp <- Last.phn$asp
    pie <- Last.phn$pie
    labels <- Last.phn$labels
    show.mutation <- Last.phn$show.mutation
    scale.ratio <- Last.phn$scale.ratio
    altlink <- Last.phn$alter.links
    threshold <- Last.phn$threshold

    if (is.character(labels)) {
        font <- Last.phn$font
        cex <- Last.phn$cex
    }

    getMove <- function() {
        p1 <- locator(1)
        if (is.null(p1)) return("done")
        ## find the closest node to p1...
        i <- which.min(sqrt((p1$x - xx)^2 + (p1$y - yy)^2))
        points(xx[i][rep(1L, 3L)], yy[i][rep(1L, 3L)],
               pch = 10, cex = 3:5, col = col.identifier)
        p2 <- locator(1)
        ## ... and change its coordinates:
        xx[i] <<- p2$x
        yy[i] <<- p2$y
    }

    if (is.null(xy)) {
        cat("\n    *** You are about to edit your haplotype network ***\n\n1. Click once on the node to be moved (it will be visually identified).\n2. Click a second time where to place it.\n3. Repeat 1. and 2. as many times as you want.\n4. Right-click to exit.\n\nNOTES:\n* The whole network is redrawn after each node move.\n* This function returns the final coordinates of the nodes so it is\n  recommended to save them (e.g., xy <- replot()).\n* The saved coordinates can later be used non-interactively\n  (e.g., replot(xy) or plot(net, xy = xy)).\n")
        res <- getMove()
    } else {
        xx <- xy$x
        yy <- xy$y
    }

    repeat {
        .setWindow4haploNet(xx, yy, size, asp, ...)
        segments(xx[l1], yy[l1], xx[l2], yy[l2], lwd = lwd,
                 lty = lty, col = col.link)
        if (show.mutation)
            .labelSegmentsHaploNet(xx, yy, cbind(l1, l2), step, size, lwd,
                                   col.link, as.numeric(show.mutation), scale.ratio)
        if (!is.null(altlink) && !identical(as.numeric(threshold), 0))
            .drawAlternativeLinks(xx, yy, altlink, threshold, show.mutation, scale.ratio, size)
        .drawSymbolsHaploNet(xx, yy, size, col, bg, shape, pie)
        if (is.character(labels))
            text(xx, yy, labels, font = font, cex = cex)
        assign("tmp", list(xx, yy), envir = .PlotHaploNetEnv)
        eval(quote(last_plot.haploNet[1:2] <- tmp), envir = .PlotHaploNetEnv)
        if (!is.null(xy)) break
        res <- getMove()
        if (identical(res, "done")) break
    }
}


## set the (empty) graphic window before plotting the network
.setWindow4haploNet <- function(xx, yy, size, asp, ...)
{
    dots <- list(...)
    noxlim <- noylim <- TRUE
    if (length(dots)) {
        ## find the bounding box unless xlim and/or ylim given
        if ("xlim" %in% names(dots)) noxlim <- FALSE
        if ("ylim" %in% names(dots)) noylim <- FALSE
    }
    half.size <- size / 2
    if (noxlim) dots$xlim <- c(min(xx - half.size), max(xx + half.size))
    if (noylim) dots$ylim <- c(min(yy - half.size), max(yy + half.size))
    dots$x <- NA
    dots$type <- dots$bty <- "n"
    dots$ann <- dots$axes <- FALSE
    dots$asp <- asp
    do.call(plot, dots)
}

plot.haploNet <- function(x, size = 1, col, bg, col.link, lwd, lty,
             shape = "circles", pie = NULL, labels, font, cex, col.lab,
             scale.ratio, asp = 1, legend = FALSE, fast = FALSE,
             show.mutation, threshold = c(1, 2), xy = NULL, ...)
{
    ## map options to arguments
    OPTS <- get("plotHaploNetOptions", envir = .PlotHaploNetEnv)
    if (missing(col)) col <- OPTS$haplotype.outer.color
    if (missing(bg)) {
        bg <- if (is.null(pie)) OPTS$haplotype.inner.color else OPTS$pie.colors.function
    }
    if (missing(col.link)) col.link <- OPTS$link.color
    if (missing(lwd)) lwd <- OPTS$link.width
    if (missing(lty)) lty <- OPTS$link.type
    if (missing(labels)) labels <- OPTS$labels
    if (missing(font)) font <- OPTS$labels.font
    if (missing(cex)) cex <- OPTS$labels.cex
    if (missing(col.lab)) col.lab <- OPTS$labels.color
    if (missing(scale.ratio)) scale.ratio <- OPTS$scale.ratio
    if (missing(show.mutation)) show.mutation <- OPTS$show.mutation

    SHAPES <- c(c = "circles", s = "squares", d = "diamonds")
    shape <- SHAPES[substr(tolower(shape), 1L, 1L)]
    ## par(xpd = TRUE) # normalement plus utile...
    link <- x[, 1:2, drop = FALSE]
    l1 <- x[, 1]
    l2 <- x[, 2]
    ld <- x[, 3] * scale.ratio

    tab <- tabulate(link)
    n <- length(tab)
    show.mutation <- as.integer(show.mutation)

    ## adjust 'ld' wrt the size of the symbols:
    size <- rep(size, length.out = n)
    ld <- ld + (size[l1] + size[l2])/2

    if (!is.null(xy)) {
        xx <- xy$x
        yy <- xy$y
        fast <- TRUE
    } else {
        if (n < 3) {
            xx <- c(-ld, ld)/2
            yy <- rep(0, 2)
            fast <- TRUE
        } else {
            xx <- yy <- angle <- theta <- r <- numeric(n)
            avlb <- !logical(length(ld))

            H <- vector("list", n) # the list of hierarchy of nodes...

            ## the recursive function to allocate quadrants
            foo <- function(i) {
                j <- integer() # indices of the haplotypes linked to 'i'
                for (ii in 1:2) { # look at both columns
                    ll <- which(link[, ii] == i & avlb)
                    if (length(ll)) {
                        newj <- link[ll, -ii]
                        r[newj] <<- ld[ll]
                        j <- c(j, newj)
                        avlb[ll] <<- FALSE
                    }
                }
                if (nlink <- length(j)) {
                    H[[i]] <<- j
                    start <- theta[i] - angle[i]/2
                    by <- angle[i]/nlink
                    theta[j] <<- seq(start, by = by, length.out = nlink)
                    angle[j] <<- by
                    xx[j] <<- r[j] * cos(theta[j]) + xx[i]
                    yy[j] <<- r[j] * sin(theta[j]) + yy[i]
                    for (ii in j) foo(ii)
                }
            }

            ## start with the haplotype with the most links:
            central <- which.max(tab)
            angle[central] <- 2*pi
            foo(central)
        }
    }

    if (!fast) {
        fCollect <- function(i) {
            ## find all nodes to move simultaneously
            ii <- H[[i]]
            if (!is.null(ii)) {
                j <<- c(j, ii)
                for (jj in ii) fCollect(jj)
            }
        }

        ## NOTE: moving this function outside of the body of plot.haploNet() is not more efficient
        energy <- function(xx, yy) {
            ## First, check line crossings
            nlink <- length(l1)
            ## rounding makes it work better (why???)
            x0 <- round(xx[l1])
            y0 <- round(yy[l1])
            x1 <- round(xx[l2])
            y1 <- round(yy[l2])
            ## compute all the slopes and intercepts:
            beta <- (y1 - y0)/(x1 - x0)
            if (any(is.na(beta))) return(Inf)
            intp <- y0 - beta*x0
            for (i in 1:(nlink - 1)) {
                for (j in (i + 1):nlink) {
                    ## in case they are parallel:
                    if (beta[i] == beta[j]) next
                    ## if both lines are vertical:
                    if (abs(beta[i]) == Inf && abs(beta[j]) == Inf) next
                    ## now find the point where both lines cross
                    if (abs(beta[i]) == Inf) { # in case the 1st line is vertical...
                        xi <- x0[i]
                        yi <- beta[j]*xi + intp[j]
                    } else if (abs(beta[j]) == Inf) { # ... or the 2nd one
                        xi <- x0[j]
                        yi <- beta[i]*xi + intp[i]
                    } else {
                        xi <- (intp[j] - intp[i])/(beta[i] - beta[j])
                        yi <- beta[i]*xi + intp[i]
                    }
                    xi <- round(xi) # rounding as above
                    yi <- round(yi)

                    if (x0[i] < x1[i]) {
                        if (xi <= x0[i] || xi >= x1[i]) next
                    } else {
                        if (xi >= x0[i] || xi <= x1[i]) next
                    }

                    if (y0[i] < y1[i]) {
                        if (yi <= y0[i] || yi >= y1[i]) next
                    } else {
                        if (yi >= y0[i] || yi <= y1[i]) next
                    }

                    ## if we reach here, the intersection point is on
                    ## the 1st segment, check if it is on the 2nd one
                    if (x0[j] < x1[j]) {
                        if (xi <= x0[j] || xi >= x1[j]) next
                    } else {
                        if (xi >= x0[j] || xi <= x1[j]) next
                    }

                    if (y0[i] < y1[j]) {
                        if (yi <= y0[j] || yi >= y1[j]) next
                    } else {
                        if (yi >= y0[j] || yi <= y1[j]) next
                    }
                    return(Inf)
                }
            }
            D <- dist(cbind(xx, yy))
            sum(1/c(D)^2, na.rm = TRUE)
        }

        Rotation <- function(rot, i, beta) {
            ## rot: indices of the nodes to rotate
            ## i: index of the node connected to 'rot' (= fixed rotation point)
            xx.rot <- xx[rot] - xx[i]
            yy.rot <- yy[rot] - yy[i]
            theta <- atan2(yy.rot, xx.rot) + beta
            h <- sqrt(xx.rot^2 + yy.rot^2)
            new.xx[rot] <<- h*cos(theta) + xx[i]
            new.yy[rot] <<- h*sin(theta) + yy[i]
        }

        OptimizeRotation <- function(node, rot) {
            ## test the direction 1st
            inc <- pi/90
            Rotation(rot, node, inc)
            new.E <- energy(new.xx, new.yy)
            if (new.E >= E) {
                inc <- -inc
                Rotation(rot, node, inc)
                new.E <- energy(new.xx, new.yy)
            }

            while (new.E < E) {
                xx <<- new.xx
                yy <<- new.yy
                E <<- new.E
                Rotation(rot, node, inc)
                new.E <- energy(new.xx, new.yy)
            }
        }

        E <- energy(xx, yy)
        new.xx <- xx
        new.yy <- yy

        nextOnes <- NULL
        for (i in H[[central]][-1]) {
            ## collect the nodes descending from 'i':
            j <- NULL # j must be initialized before calling fCollect
            fCollect(i)
            rot <- c(i, j) # index of the nodes that will be moved
            OptimizeRotation(central, rot)
            nextOnes <- c(nextOnes, i)
        }

        while (!is.null(nextOnes)) {
            newNodes <- nextOnes
            nextOnes <- NULL
            for (i in newNodes) {
                if (is.null(H[[i]])) next
                for (j in H[[i]]) {
                    fCollect(j)
                    rot <- j
                    OptimizeRotation(i, rot)
                    nextOnes <- c(nextOnes, rot)
                }
            }
        }
    }

    .setWindow4haploNet(xx, yy, size, asp, ...)
    segments(xx[l1], yy[l1], xx[l2], yy[l2], lwd = lwd,
             lty = lty, col = col.link)

    ## draw alternative links
    altlink <- attr(x, "alter.links")
    if (!is.null(altlink) && !identical(as.numeric(threshold), 0))
        .drawAlternativeLinks(xx, yy, altlink, threshold, show.mutation, scale.ratio, size)

    if (show.mutation) {
        if (show.mutation == 1 && all(x[, 3] < 1)) {
            warning("link lengths < 1: cannot use the default for 'show.mutation' (changed to 3)")
            show.mutation <- 3
        }
        .labelSegmentsHaploNet(xx, yy, link, x[, 3], size, lwd, col.link,
                               show.mutation, scale.ratio)
    }

    .drawSymbolsHaploNet(xx, yy, size, col, bg, shape, pie)

    if (labels) {
        labels <- attr(x, "labels") # for export below
        text(xx, yy, labels, font = font, cex = cex, col = col.lab)
    }
    if (legend[1]) {
        if (is.logical(legend)) {
            cat("Click where you want to draw the legend")
            xy <- unlist(locator(1))
            cat("\nThe coordinates x = ", xy[1], ", y = ", xy[2], " are used\n", sep = "")
        } else {
            if (!is.numeric(legend) || length(legend) < 2)
                stop("wrong coordinates of legend")
            xy <- legend
        }
        if (length(SZ <- unique(size)) > 1) {
            SZ <- unique(c(min(SZ), floor(median(SZ)), max(SZ)))
            SHIFT <- max(SZ) / 2
            vspace <- strheight(" ")
            if (any(shape == "circles")) {
                for (sz in SZ) {
                    seqx <- seq(-sz / 2, sz / 2, length.out = 100)
                    seqy <- sqrt((sz /2)^2 - seqx^2)
                    seqx <- seqx + xy[1] + SHIFT
                    seqy <- xy[2] + seqy - SHIFT
                    lines(seqx, seqy)
                    text(seqx[100], seqy[100], sz, adj = c(0.5, 1.1))
                }
                xy[2] <- xy[2] - SHIFT - 2 * vspace # update 'y'
            }
            if (any(shape == "squares") || any(shape == "diamonds")) {
                sqrtPIon4 <- sqrt(pi) / 4
                ## center of the squares:
                orig.x <- xy[1] + max(SZ) * sqrtPIon4
                orig.y <- xy[2] - max(SZ) * sqrtPIon4
                for (sz in SZ) {
                    TMP <- sz * sqrtPIon4
                    lines(orig.x + c(-TMP, -TMP, TMP, TMP),
                          orig.y + c(0, TMP, TMP, 0))
                    text(orig.x + TMP, orig.y, sz, adj = c(0.5, 1.1))
                }
                xy[2] <- xy[2] - SHIFT - 2 * vspace # update
            }
        }
        if (!is.null(pie)) {
            nc <- ncol(pie)
            ##co <- if (length(bg) == 1 && bg == "white") rainbow(nc) else rep(bg, length.out = nc)
            co <- if (is.function(bg)) bg(nc) else rep(bg, length.out = nc)
            w <- diff(par("usr")[3:4]) / 40
            TOP <- seq(xy[2], by = -w, length.out = nc)
            BOTTOM <- TOP + diff(TOP[1:2]) * 0.9
            LEFT <- rep(xy[1], nc)
            RIGHT <- LEFT + w
            rect(LEFT, BOTTOM, RIGHT, TOP, col = co)
            text(RIGHT + strwidth("n"), (TOP + BOTTOM) /2, colnames(pie), adj = 0)
        }
    }
    assign("last_plot.haploNet",
           list(xx = xx, yy = yy, link = link, step = x[, 3], size = size,
                col = col, bg = bg, lwd = lwd, lty = lty, shape = shape,
                col.link = col.link,labels = labels, font = font, cex = cex,
                asp = asp, pie = pie, show.mutation = show.mutation,
                scale.ratio = scale.ratio, alter.links = altlink,
                threshold = threshold),
           envir = .PlotHaploNetEnv)
}

plot.haplotype <- function(x, xlab = "Haplotype", ylab = "Number", ...)
{
    barplot(sapply(attr(x, "index"), length), xlab = xlab, ylab = ylab,
            names.arg = rownames(x), ...)
}

sort.haplotype <-
    function(x, decreasing = ifelse(what == "frequencies", TRUE, FALSE), what = "frequencies", ...)
{
    oc <- oldClass(x)
    from <- attr(x, "from")
    what <- match.arg(what, c("frequencies", "labels"))
    idx <- attr(x, "index")
    o <- switch(what,
                frequencies = order(sapply(idx, length), decreasing = decreasing),
                labels = order(rownames(x), decreasing = decreasing))
    x <- x[o, ]
    attr(x, "index") <- idx[o]
    class(x) <- oc
    attr(x, "from") <- from
    x
}

subset.haplotype <- function(x, minfreq = 1, maxfreq = Inf, maxna = Inf,
                             na = c("N", "?"), ...)
{
    oc <- oldClass(x)
    from <- attr(x, "from")
    idx <- attr(x, "index")
    f <- sapply(idx, length)
    s <- f <= maxfreq & f >= minfreq
    if (is.finite(maxna)) {
        na <- tolower(na)
        all <- c("r", "m", "w", "s", "k", "y", "v", "h", "d", "b", "n", "-", "?")
        if (identical(na, "all")) na <- all
        if (identical(na, "ambiguous")) na <- all[1:11]
        freq <- if (maxna < 1) FALSE else TRUE
        foo <- function(x) sum(base.freq(x, freq, TRUE)[na])
        count.na <- numeric(n <- nrow(x))
        for (i in seq_len(n)) count.na[i] <- foo(x[i, ]) # cannot use apply
        s <- s & count.na <= maxna
    }
    x <- x[s, ]
    attr(x, "index") <- idx[s]
    class(x) <- oc
    attr(x, "from") <- from
    x
}

summary.haplotype <- function(object, ...)
    setNames(lengths(attr(object, "index")), rownames(object))

print.haplotype <- function(x, ...)
{
    n <- (d <- dim(x))[1]
    DF <- summary.haplotype(x)
    cat("\nHaplotypes extracted from:", attr(x, "from"), "\n\n")
    cat("    Number of haplotypes:", n, "\n")
    cat("         Sequence length:", d[2], "\n\n")
    cat("Haplotype labels and frequencies:\n\n")
    if (n <= 40) print(DF)
    else {
        print(DF[1:40])
        cat("...\n(use summary() to print all)\n")
    }
}

"[.haplotype" <- function(x, ...)
{
    cls <- class(x)
    if (length(cls) > 1) cls <- cls[-1]
    y <- NextMethod("[")
    class(y) <- cls
    y
}

as.phylo.haploNet <- function(x, quiet = FALSE, ...)
{
    if (!is.null(attr(x, "alter.links")) && !quiet)
        warning("some links (edges) were dropped because of reticulations")

    foo <- function(xx) {
        NODES <<- c(NODES, xx)
        W <- which(mat == xx, TRUE)
        for (i in 1:nrow(W)) {
            k <- W[i, 1]
            if (DONE[k]) next
            DONE[k] <<- TRUE
            link <- mat[k, ]
            if (link[1] != xx) link <- rev(link)
            edge[ie, ] <<- link
            el[ie] <<- x[k, 3]
            ie  <<- ie + 1L
            neigh <- link[2]
            if (deg[neigh] == 1L) TIPS <<- c(TIPS, neigh) else foo(neigh)
        }
    }

    LABS <- attr(x, "labels")
    n <- length(LABS) # number of haplotype nodes
    mat <- x[, 1:2]
    deg <- tabulate(mat, n) # get the degree of each node
    deg1 <- deg == 1
    ntip <- sum(deg1)

    TIPS <- integer()
    NODES <- integer()
    edge <- mat
    edge[] <- 0L
    el <- numeric(nrow(mat))
    ie <- 1L
    DONE <- logical(nrow(mat))
    ROOT <- which(!deg1)[1]
    foo(ROOT)
    TIPNODE <- c(TIPS, NODES)
    edge <- match(edge, TIPNODE)
    dim(edge) <- dim(mat)
    phy <- list(edge = edge, edge.length = el, tip.label = LABS[TIPS],
                Nnode = length(NODES), node.label = LABS[NODES])
    class(phy) <- "phylo"
    phy
}

as.evonet.haploNet <- function(x, ...)
{
    res <- as.phylo.haploNet(x, quiet = TRUE)
    alt <- attr(x, "alter.links")
    if (!is.null(alt)) {
        LABS <- attr(x, "labels")
        o <- match(LABS, c(res$tip.label, res$node.label))
        alt <- o[alt[, 1:2]]
        dim(alt) <- c(length(alt)/2, 2)
        res$reticulation <- alt
        class(res) <- c("evonet", "phylo")
    }
    res
}

if (getRversion() >= "2.15.1") utils::globalVariables(c("network", "network.vertex.names<-"))
as.network.haploNet <- function(x, directed = FALSE, altlinks = TRUE, ...)
{
    res <- x[, 1:2]
    if (altlinks) res <- rbind(res, attr(x, "alter.links")[, 1:2])
    res <- network(res, directed = directed, ...)
    network.vertex.names(res) <- attr(x, "labels")
    res
}

if (getRversion() >= "2.15.1") utils::globalVariables("graph.edgelist")
as.igraph.haploNet <- function(x, directed = FALSE, use.labels = TRUE,
                               altlinks = TRUE, ...)
{
    y <- x[, 1:2]
    if (altlinks) y <- rbind(y, attr(x, "alter.links")[, 1:2])
    y <-
        if (use.labels) matrix(attr(x, "labels")[y], ncol = 2)
        else y - 1L
    graph.edgelist(y, directed = directed, ...)
}

## approximate function giving the cophenetic distance on a network
## -> can be used to find nodes with short distances to most nodes, eg:
## rowSums(cophenetic.haploNet(net))
cophenetic.haploNet <- function(x)
{
    LABS <- attr(x, "labels")
    phy <- as.phylo.haploNet(x, quiet = TRUE)
    res <- dist.nodes(phy)
    o <- match(LABS, c(phy$tip.label, phy$node.label))
    res <- res[o, o]
    alt <- attr(x, "alter.links")
    if (!is.null(alt)) {
        for (i in 1:nrow(alt)) {
            a <- alt[i, 1]
            b <- alt[i, 2]
            res[a, b] <- res[b, a] <- alt[i, 3]
        }
    }
    dimnames(res) <- list(LABS, LABS)
    res
}

labels.haploNet <- function(object, ...) attr(object, "labels")

haplotype.loci <- function(x, locus = 1:2, quiet = FALSE, compress = TRUE,
                           check.phase = TRUE, ...)
{
    x <- x[, attr(x, "locicol")[locus]]
    nloc <- ncol(x)
    n <- nrow(x)

    if (check.phase) {
        ## drop the rows with at least one unphased genotype:
        s <- apply(is.phased(x), 1, all)
        if (any(!s)) {
            x <- x[s, ]
            warning(paste("dropping", sum(!s), "individual(s) out of", n, "due to unphased genotype(s)"))
            n <- nrow(x)
        }
    }

### NOTE: trying to find identical rows first does not speed calculations

    ## initialise (works for all levels of ploidy)
    nh <- .checkPloidy(x[, 1, drop = FALSE]) # the number of haplotypes
    names(nh) <- NULL
    res <- matrix("", nloc, nh * n)
    class(x) <- "data.frame" # drop "loci"

### unlist(lapply()) is a bit faster than filling the matrix with 'for'
### thanks to use.names = FALSE in unlist(). The old code is:
###    y <- matrix(NA_integer_, n, nloc)
###    for (j in seq_len(nloc)) y[, j] <- as.integer(x[[j]])
    y <- unlist(lapply(x, as.character), FALSE, FALSE)
    dim(y) <- c(n, nloc)

    mysplit <- function(x)
        unlist(strsplit(x, "|", fixed = TRUE, useBytes = TRUE), FALSE, FALSE)

    k <- 1:nh
    for (i in seq_len(n)) { # loop along each individual
        if (!quiet && !(i %% 100)) cat("\rAnalysing individual no.", i, "/", n)
        tmp <- mysplit(y[i, ])
        dim(tmp) <- c(nh, nloc) # arrange the alleles in a matrix...
        res[, k] <- t(tmp)      # faster than using matrix(, byrow = TRUE)
        k <- k + nh
    }

### experimental code to run the above loop with parallel
### foo <- function(i) {
###     tmp <- mysplit(y[i, ])
###     dim(tmp) <- c(nh, nloc)
###     t(tmp)
### }
### res <- mclapply(seq_len(n), foo)
### res <- unlist(res, FALSE, FALSE)
### dim(res) <- c(nloc, n * nh)

    if (!compress) {
        rownames(res) <- colnames(x)
        return(res)
    }

    nc <- ncol(res)
    h <- .Call(unique_haplotype_loci, res, nloc, nc)
    u <- h == 0
    if (all(u)) freq <- rep(1L, nc) else {
        i <- which(u)
        res <- res[, i, drop = FALSE]
        h[u] <- i
        freq <- tabulate(h)
        freq <- freq[freq > 0]
    }

    if (!quiet) cat("\rAnalysing individual no.", n, "/", n, "\n")

    rownames(res) <- colnames(x)
    class(res) <- "haplotype.loci"
    attr(res, "freq") <- freq
    res
}

plot.haplotype.loci <- function(x, ...)
{
    y <- attr(x, "freq")
    names(y) <- apply(x, 2, paste, collapse = "-")
    barplot(y, ...)
}

dist.hamming <- function(x)
{
    n <- nrow(x)
    if (n < 2) stop("less than two haplotypes")
    d <- numeric(n * (n - 1)/2)
    k <- 1L
    for (i in 1:(n - 1)) {
        for (j in (i + 1):n) {
            d[k] <- sum(x[i, ] != x[j, ])
            k <- k + 1L
        }
    }
    attr(d, "Size") <- n
    attr(d, "Labels") <- rownames(x)
    attr(d, "Diag") <- attr(d, "Upper") <- FALSE
    attr(d, "call") <- match.call()
    attr(d, "method") <- "N"
    class(d) <- "dist"
    d
}

dist.haplotype.loci <- function(x)
{
    n <- ncol(x)
    if (n < 2) stop("less than two haplotypes")
    d <- numeric(n * (n - 1) / 2)
    k <- 1L
    for (i in 1:(n - 1)) {
        for (j in (i + 1):n) {
            d[k] <- sum(x[, i] != x[, j])
            k <- k + 1L
        }
    }
    attr(d, "Size") <- n
    attr(d, "Labels") <- colnames(x)
    attr(d, "Diag") <- attr(d, "Upper") <- FALSE
    attr(d, "call") <- match.call()
    attr(d, "method") <- "N"
    class(d) <- "dist"
    d
}

LD <- function(x, locus = c(1, 2), details = TRUE)
{
    if (length(locus) != 2)
        stop("you must specify two loci to compute linkage disequilibrium")
    hap <- haplotype.loci(x, locus = locus, quiet = TRUE)
    alleles1 <- unique(hap[1, ])
    alleles2 <- unique(hap[2, ])
    k <- length(alleles1)
    m <- length(alleles2)
    nij <- matrix(0L, k, m, dimnames = list(alleles1, alleles2))
    nij[t(hap)] <- attr(hap, "freq")
    N <- sum(nij)
    pij <- nij / N
    pi <- rowSums(pij)
    qj <- colSums(pij)
    eij <- pi %o% qj * N
    pi <- rep(pi, ncol(pij))
    qj <- rep(qj, each = nrow(pij))
    D <- pij - pi * qj
    rij <- D / sqrt(pi * (1 - pi) * qj * (1 - qj))
    df <- (k - 1) * (m - 1)
    T2 <- df * N * sum(rij^2) / (k * m)
    res <- c("T2" = T2, "df" = df, "P-val" = pchisq(T2, df, lower.tail = FALSE))
    if (details) {
        res <- list(nij, eij, rij, 2 * sum(nij * log(nij / eij)),
                    2 * sum((nij - eij)^2 / eij), res)
        names(res) <- c("Observed frequencies", "Expected frequencies", "Correlations among alleles",
                       "LRT (G-squared)", "Pearson's test (chi-squared)", "T2")
    }
    res
}

LD2 <- function(x, locus = c(1, 2), details = TRUE)
{
    if (length(locus) != 2)
        stop("you must specify two loci to compute linkage disequilibrium")
    x <- x[, attr(x, "locicol")[locus]]
    if (any(.checkPloidy(x) != 2))
        stop("linkage disequilibrium with unphased genotypes works only for diploid data")
    n <- nrow(x)
    s <- summary(x)
    alleles1 <- names(s[[1]]$allele)
    alleles2 <- names(s[[2]]$allele)
    genotypes1 <- names(s[[1]]$genotype)
    genotypes2 <- names(s[[2]]$genotype)

    ## get allele proportions:
    P <- lapply(s, "[[", "allele")
    P <- rapply(P, function(x) x/sum(x), how = "replace")

    ## compute HW disequilibrium coefficients:
    G <- lapply(s, "[[", "genotype")
    G <- rapply(G, function(x) x/sum(x), how = "replace")
    DHW <- vector("list", 2)
    for (i in 1:2) {
        alls <- if (i == 1) alleles1 else alleles2
        tmp <- G[[i]][paste(alls, alls, sep = "/")] # frequencies of all possible homozygotes
        tmp[is.na(tmp)] <- 0 # for the unobserved
        names(tmp) <- alls
        DHW[[i]] <- tmp - P[[i]]^2
    }

    Delta <- r <- numeric()

    x1 <- as.integer(x[, 1L, drop = TRUE])
    x2 <- as.integer(x[, 2L, drop = TRUE])

    for (a in alleles1) {
        Pa <- P[[1]][a]
        Da <- DHW[[1]][a]
        for (b in alleles2) {
            Pb <- P[[2]][b]
            Db <- DHW[[2]][b]

            ## find the homozygote genotypes:
            i <- match(paste(a, a, sep = "/"), genotypes1)
            j <- match(paste(b, b, sep = "/"), genotypes2)

            n1 <- n2 <- n3 <- n4 <- 1e100 # initialise
            ## if the homozygotes were not observed set the n's accordingly,
            ## else find them in the respective columns:
            if (is.na(i)) n1 <- n2 <- 0 else Ho1 <- x1 == i
            if (is.na(j)) n1 <- n3 <- 0 else Ho2 <- x2 == j

            ## count the homozygotes at both loci:
            if (as.logical(n1)) n1 <- sum(Ho1 & Ho2)

            ## the next command will find the genotypes with the 1st allele (homozygote OR heterozygote):
            tmp1 <- grep(paste0("^", a, "/|/", a, "$"), genotypes1)
            if (!is.na(i)) tmp1 <- tmp1[tmp1 != i] # remove the homozygote genotype if needed
            ## and get the heterozygotes:
            if (length(tmp1)) He1 <- x1 %in% tmp1 else n3 <- n4 <- 0

            ## ... and the same for the 2nd allele:
            tmp2 <- grep(paste0("^", b, "/|/", b, "$"), genotypes2)
            if (!is.na(j)) tmp2 <- tmp2[tmp2 != j]
            if (length(tmp2)) He2 <- x2 %in% tmp2 else n2 <- n4 <- 0

            ## if homozygotes were observed for each locus:
            if (as.logical(n2)) n2 <- sum(Ho1 & He2)
            if (as.logical(n3)) n3 <- sum(Ho2 & He1)
            if (as.logical(n4)) n4 <- sum(He1 & He2)

            delta <- (2 * n1 + n2 + n3 + 0.5 * n4)/n - 2 * Pa * Pb
            Delta <- c(Delta, delta)
            r <- c(r, delta^2 / ((Pa * (1 - Pa) + Da) * (Pb * (1 - Pb) + Db)))
        }
    }
    k <- length(alleles1)
    m <- length(alleles2)
    df <- (k - 1) * (m - 1)
    T2 <- df * n * sum(r) / (k * m)
    res <- c("T2" = T2, "df" = df, "P-val" = pchisq(T2, df, lower.tail = FALSE))
    if (details) {
        dim(Delta) <- c(k, m)
        dimnames(Delta) <- list(alleles1, alleles2)
        res <- list(Delta = Delta, T2 = res)
    }
    res
}

LDscan <- function(x, ...) UseMethod("LDscan")

LDscan.DNAbin <- function(x, quiet = FALSE, what = c("r", "Dprime"), ...)
{
    what <- match.arg(what)
    if (!quiet) cat("Scanning haplotypes... ")
    ss <- seg.sites(x)
    nloci <- length(ss)
    hap <- t(as.character(x[, ss]))
    hap[hap == "n"] <- NA_character_
    if (!quiet) cat("done.\n")
    .LD <- function (hap, loc1, loc2) {
        nij <- table(hap[loc1, ], hap[loc2, ])
        if (any(dim(nij) != 2)) return(NA_real_)
        N <- sum(nij)
        pij <- nij/N
        pi <- rep(rowSums(pij), 2)
        qj <- rep(colSums(pij), each = 2)
        D <- pij - pi * qj
        if (what == "Dprime") {
            D <- D[1]
            denom <- if (D > 0) min(-pi[1] * qj[3], -pi[2] * qj[1]) else max(-pi[1] * qj[1], -pi[2] * qj[2])
            return(D/denom)
        }
        rij <- D/sqrt(pi * (1 - pi) * qj * (1 - qj))
        abs(rij[1])
    }
    M <- nloci * (nloci - 1) / 2
    ldx <- numeric(M)
    k <- 0L
    for (i in 1:(nloci - 1)) {
        for (j in (i + 1):nloci) {
            k <- k + 1L
            ldx[k] <- .LD(hap, i, j)
            if (!quiet) cat("\r", round(100 * k / M), "%")
        }
    }
    if (!quiet) cat("\n")
    class(ldx) <- "dist"
    attr(ldx, "Size") <- nloci
    attr(ldx, "Labels") <- names(x)
    attr(ldx, "Diag") <- attr(ldx, "Upper") <- FALSE
    attr(ldx, "call") <- match.call()
    ldx
}

LDscan.loci <- function(x, depth = NULL, quiet = FALSE, what = c("r", "Dprime"), ...)
{
    what <- match.arg(what)
    nloci <- length(attr(x, "locicol"))
    if (!quiet) cat("Scanning haplotypes... ")
    hap <- haplotype.loci(x, seq_len(nloci), TRUE, FALSE, FALSE)
    if (!quiet) cat("done.\n")
    .LD <- function (hap, loc1, loc2) {
        nij <- table(hap[loc1, ], hap[loc2, ])
        N <- sum(nij)
        pij <- nij/N
        pi <- rep(rowSums(pij), 2)
        qj <- rep(colSums(pij), each = 2)
        D <- pij - pi * qj
        if (what == "Dprime") {
            D <- D[1]
            denom <-
                if (D > 0) min(-pi[1] * qj[3], -pi[2] * qj[1])
                else max(-pi[1] * qj[1], -pi[2] * qj[2])
            return(D/denom)
        }
        rij <- D/sqrt(pi * (1 - pi) * qj * (1 - qj))
        abs(rij[1])
    }
    if (is.null(depth)) {
        M <- nloci * (nloci - 1) / 2
        ldx <- numeric(M)
        k <- 0L
        for (i in 1:(nloci - 1)) {
            for (j in (i + 1):nloci) {
                k <- k + 1L
                ldx[k] <- .LD(hap, i, j)
                if (!quiet) cat("\r", round(100 * k / M), "%")
            }
        }
        if (!quiet) cat("\n")
        class(ldx) <- "dist"
        attr(ldx, "Size") <- nloci
        attr(ldx, "Labels") <- names(x)
        attr(ldx, "Diag") <- attr(ldx, "Upper") <- FALSE
        attr(ldx, "call") <- match.call()
    } else {
        depth <- as.integer(depth)
        ldx <- vector("list", length(depth))
        k <- 0L
        for (o in depth) {
            vec <- numeric(nloci - o)
            j <- 0L
            for (i in 1:(nloci - o)) {
                j <- j + 1L
                vec[j] <- .LD(hap, i, i + o)
                if (!quiet) cat("\rScanning at depth", o, ":\t", round(100 * j / (nloci - o)), "%")
            }
            if (!quiet) cat("\n")
            k <- k + 1L
            ldx[[k]] <- vec
        }
        names(ldx) <- depth
    }
    ldx
}

LDmap <- function(d, POS = NULL, breaks = NULL, col = NULL, border = NA,
                  angle = 0, asp = 1, cex = 1, scale.legend = 0.8, ...)
{
    if (!is.null(POS) & angle != 0) {
        warning("'angle' set to 0 since 'POS' is provided")
        angle <- 0
    }
    if (is.matrix(d)) d <- as.dist(d)
    nloci <- attr(d, "Size")
    n <- nloci - 1
    m <- nloci * n /2
    nl <- if (is.null(col)) 10 else length(col) # Nb of colour levels
    if (is.null(breaks)) {
        rgd <- range(d, na.rm = TRUE)
        breaks <- seq(rgd[1], rgd[2], length.out = nl + 1)
    } else {
        nl <- length(breaks) - 1
        if (!is.null(col))
            col <- rep(col, length.out = nl)
    }
    if (is.null(col))
        col <- colorRampPalette(c("lightyellow", "red"))(nl)

    co <- col[cut(d, breaks, include.lowest = TRUE)]

    x.lab.loci <- 1:nloci - 0.75
    y.lab.loci <- 1:nloci - 0.25

    use.rect <- angle == 45

    ## assume angle = 45
    yb <- unlist(mapply(":", 1:n, n, USE.NAMES = FALSE))
    yt <- yb + 1
    xl <- unlist(mapply(rep, 0:(n - 1), n:1, USE.NAMES = FALSE))
    xr <- xl + 1

    if (!use.rect) {
        xx <- rbind(xl, xl, xr, xr, rep(NA, m))
        yy <- rbind(yb, yt, yt, yb, rep(NA, m))
        dim(xx) <- dim(yy) <- NULL
        tmp <- rect2polar(xx, yy)
        new.angle <- tmp$angle + 2 * pi * (angle - 45)/360
        xy <- polar2rect(tmp$r, new.angle)
        X <- range(xy$x, na.rm = TRUE)
        Y <- range(xy$y, na.rm = TRUE)
        ## adjust the coordinates of the labels of the loci:
        tmp <- rect2polar(x.lab.loci, y.lab.loci)
        new.angle <- tmp$angle + 2 * pi * (angle - 45)/360
        tmp <- polar2rect(tmp$r, new.angle)
        x.lab.loci <- tmp$x
        y.lab.loci <- tmp$y
    } else {
        X <- Y <- c(0, nloci)
    }

    if (!is.null(POS)) Y[1] <- -0.1 * Y[2]

    plot.default(X, Y, "n", axes = FALSE, xlab = "", ylab = "",
                 asp = asp, ...)
    if (use.rect) rect(xl, yb, xr, yt, col = co, border = border)
    else polygon(xy, col = co, border = border)

    psr <- par("usr")

    if (!is.null(POS)) {
        f <- function(x, mn, mx) {
            x <- x - mn # translate to the origin
            x <- x / mx # rescale to 1
            (psr[2] - psr[1]) * x + psr[1]
        }
        POS <- as.double(POS)
        mn <- POS[1]
        mx <- POS[nloci] - mn
        AT <- pretty(POS)
        at <- f(AT, mn, mx)
        segments(psr[1], Y[1], psr[2], Y[1], lwd = 2)
        segments(at, Y[1] - 1, at, Y[1])
        text(at, Y[1], AT, adj = c(0.5, 2), cex = cex)
        segments(f(POS, mn, mx), Y[1], x.lab.loci, y.lab.loci, lty = 3)
        mtext("Position", 1, 1.5, cex = cex)
    }

    x1 <- psr[1] + scale.legend
    y1 <- seq(psr[4] - scale.legend * nl, psr[4] - scale.legend,
              by = scale.legend)
    y2 <- y1 + scale.legend
    rect(psr[1], y1, x1, y2, col = col, border = border)
    text(x1, c(y1[1], y2), round(breaks, 3), adj = -0.1, xpd = TRUE,
         cex = cex)
    text(x.lab.loci, y.lab.loci, attr(d, "Labels"),
         srt = angle - 90, adj = 0, cex = cex)
}

all.equal.haploNet <- function(target, current, use.steps = TRUE, ...)
{
    nt1 <- sQuote(deparse(substitute(target)))
    nt2 <- sQuote(deparse(substitute(current)))
    if (identical(target, current)) return(TRUE)
    ## function to build a list to make comparisons easier
    foo <- function(x) {
        mat <- x[, 1:2]
        step <- x[, 3]
        if (!is.null(attr(x, "alter.links"))) {
            mat <- rbind(mat, attr(x, "alter.links")[, 1:2])
            step <- c(step, attr(x, "alter.links")[, 3])
        }
        dm <- dim(mat)
        mat <- attr(x, "labels")[mat]
        dim(mat) <- dm
        mat <- t(apply(mat, 1, sort))
        list(mat = mat, step = step)
    }
    ## function to arrange print of links:
    bar <- function(x) gsub("\r", "--", x)
    ## another one to print comparison of link lengths:
    bar2 <- function(x, y, z) paste0(bar(x), " (", y, ", ", z, ")")

    msg <- NULL

    labs1 <- attr(target, "labels")
    labs2 <- attr(current, "labels")
    comp12 <- is.na(match(labs1, labs2))
    comp21 <- is.na(match(labs2, labs1))
    if (all(comp12) && all(comp21))
        return(paste0("No common label between ", nt1, " and ", nt2))
    else {
        if (any(comp12))
            msg <- c(msg, paste0("Labels in ", nt2, " not in ", nt1, ":"),
                     paste(labs1[comp12], sep = ", "))
        if (any(comp21))
            msg <- c(msg, paste0("Labels in ", nt1, " not in ", nt2, ":"),
                     paste(labs2[comp21], sep = ", "))
    }

    X1 <- foo(target)
    X2 <- foo(current)
    links1 <- paste(X1$mat[, 1], X1$mat[, 2], sep = "\r")
    links2 <- paste(X2$mat[, 1], X2$mat[, 2], sep = "\r")
    comp12 <- match(links1, links2)
    if (length(links1) == length(links2)) {
        if (anyNA(comp12)) {
            comp21 <- match(links2, links1)
            msg <- c(msg, "Number of links equal",
                     paste0("Links in ", nt1, " not in ", nt2, ":"),
                     bar(links1[is.na(comp12)]),
                     paste0("Links in ", nt2, " not in ", nt1, ":"),
                     bar(links2[is.na(comp21)]))
        }
        if (use.steps) {
            tmp <- X2$step[comp12]
            test <- X1$step != tmp
            if (anyNA(test)) test[is.na(test)] <- FALSE
            if (any(test))
                msg <- c(msg, "Link lengths different (in target, in current):",
                         bar2(links1[test], X1$step[test], tmp[test]))
        }
    } else {
        msg <- c(msg, "Number of links different.")
        comp21 <- match(links2, links1)
        if (anyNA(comp12))
            msg <- c(msg, paste0("Links in ", nt1, " not in ", nt2, ":"),
                     bar(links1[is.na(comp12)]))
        if (anyNA(comp21))
            msg <- c(msg, paste0("Links in ", nt2, " not in ", nt1, ":"),
                     bar(links2[is.na(comp21)]))
        if (use.steps) {
            if (length(links1) > length(links2)) {
                tmp1 <- X1$step[!is.na(comp12)]
                tmp2 <- X2$step[comp21]
            } else {
                tmp1 <- X1$step[comp12]
                tmp2 <- X2$step[!is.na(comp21)]
            }
            test <- tmp1 != tmp2
            if (any(test))
                msg <- c(msg,
                         paste0("Links identical but of lengths different (in ", nt1, ", in ", nt2, "):"),
                         bar2(links1[tmp1][test], tmp1[test], tmp2[test]))
        }
    }
    if (is.null(msg)) TRUE else msg
}

.diffHaplo.loci <- function(h, a = 1, b = 2)
{
    s <- which(h[, a] != h[, b])
    data.frame(pos = s, h[s, c(a, b)])
}

utils::globalVariables(c("mutations.arrow.color", "mutations.arrow.type", "mutations.cex", "mutations.font", "mutations.frame.background", "mutations.frame.border", "mutations.sequence.color", "mutations.sequence.end", "mutations.sequence.length", "mutations.sequence.width", "mutations.text.color"))

mutations <- function(haploNet, link, x, y, data = NULL, style = "table",
                      POS, SEQLEN, ...)
{
    m <- haploNet[, 1:2, drop = FALSE]
    labs <- labels.haploNet(haploNet)
    if (!is.null(attr(haploNet, "alter.links")))
        m <- rbind(m, attr(haploNet, "alter.links")[, 1:2, drop = FALSE])
    if (missing(link)) {
        cat("Link is missing: select one below\n")
        cat(paste0(1:nrow(m), ": ", labs[m[, 1]], "-", labs[m[, 2]]), sep = "\n")
        cat("\nEnter a link number: ")
        link <- as.integer(readLines(n = 1))
    }
    if (link < 1 || link > nrow(m) || is.na(link)) stop("wrong value")
    Last.phn <- get("last_plot.haploNet", envir = .PlotHaploNetEnv)
    nodes <- m[link, 1:2]
    xx.link <- sum(Last.phn$xx[nodes]) / 2
    yy.link <- sum(Last.phn$yy[nodes]) / 2
    if (missing(x) || missing(y)) {
        cat("Coordinates are missing: click where you want to place the annotations:")
        xy <- locator(1)
        x <- xy$x
        y <- xy$y
        cat("\nThe coordinates x = ", x, ", y = ", y, " are used\n", sep = "")
    }
    style <- match.arg(style, c("table", "sequence"))
    if (is.null(data)) {
        data <- attr(haploNet, "data")
        if (is.null(data)) stop("no data attached to network")
    }
    ## for the moment only two main classes are accepted:
    dnabin <- inherits(data, "DNAbin")
    if (!dnabin && !inherits(data, "haplotype.loci"))
        stop("data must have the class \"DNAbin\" or \"haplotype.loci\"")
    if (!dnabin) { # => class(data) == "haplotype.loci"
        if (missing(POS)) stop("'POS' should be given")
        if (style == "sequence" && missing(SEQLEN))
            stop("'SEQLEN' should be given")
    }
    ## get/check the options
    OPTS <- getHaploNetOptions()
    dots <- list(...)
    options.names <- switch(style, "table" = {
        c("mutations.cex", "mutations.font", "mutations.frame.background",
          "mutations.frame.border", "mutations.text.color",
          "mutations.arrow.color", "mutations.arrow.type")
    }, "sequence" = {
        c("mutations.arrow.color", "mutations.arrow.type",
          "mutations.sequence.color", "mutations.sequence.end",
          "mutations.sequence.length", "mutations.sequence.width")
    })
    mapply(assign, options.names, OPTS[options.names],
           MoreArgs = list(envir = environment()))
    if (length(dots)) {
        if (any(i <- names(dots) %in% options.names)) {
            dots <- dots[i]
            mapply(assign, names(dots), dots,
                   MoreArgs = list(envir = environment()))
        }
    }
    ## get the mutations
    if (dnabin) {
        mu <- diffHaplo(data, labs[m[link, 1]], labs[m[link, 2]])
    } else {
        mu <- .diffHaplo.loci(data, labs[m[link, 1]], labs[m[link, 2]])
        mu$pos <- POS[mu$pos]
    }
    nmu <- nrow(mu)
    switch(style, "table" = {
        strg <- paste0(mu[[1]], ": ", mu[[2]], "|", mu[[3]])
        maxw <- max(strwidth(strg))
        strh <- max(strheight(strg))
        spc <- strh * 0.5
        ystrings <- seq(y - strh/2, by = -1.5 * strh, length.out = nmu)
        xleft <- x - spc
        ybottom <- ystrings[nmu] - strh/2 - spc
        xright <- x + maxw + spc
        ytop <- ystrings[1] + strh/2 + spc
        rect(xleft, ybottom, xright, ytop, col = mutations.frame.background,
             border = mutations.frame.border)
        text(x + maxw, ystrings, strg, adj = 1, cex = mutations.cex,
             font = mutations.font, col = mutations.text.color)
        if (xx.link <= xleft) {
            x0 <- xleft
        } else {
            if (xx.link >= xright) {
                x0 <- xright
            } else {
                x0 <- (xleft + xright) / 2
            }
        }
        if (yy.link <= ybottom) {
            y0 <- ybottom
        } else {
            if (yy.link >= ytop) {
                y0 <- ytop
            } else {
                y0 <- (ybottom + ytop) / 2
            }
        }
    }, "sequence" = {
        WIDTH <- diff(par("usr")[1:2])
        if (dnabin) SEQLEN <- ncol(data)
        W <- WIDTH * mutations.sequence.length # real width of the segment
        x2 <- x + W
        segments(x, y, x2, y, lwd = mutations.sequence.width,
                 col = mutations.sequence.color, lend = mutations.sequence.end)
        ## text(x2, y, " 1 kb", adj = 0)
        xx <- x + W * mu[[1]] / SEQLEN
        segments(xx, y - 0.5, xx, y + 0.5)
        x0 <- if (abs(x - xx.link) <= abs(x2 - xx.link)) x - WIDTH / 250 else x2 + WIDTH / 250
        y0 <- y
    })
    fancyarrows(x0, y0, xx.link, yy.link, col = mutations.arrow.color,
                code = mutations.arrow.type)
}