## This code is part of the polenta package
## © C. Heibl 2017, F.-S. Krah (last update 2017-11-08)
#' @title Ultra-Large Multiple Sequence Alignment with PASTA
#' @description Provides a complete reimplementation of the PASTA algorithm
#' (Mirarab, Nguyen, and Warnow 2014) in R.
#' @param seqs An object of class \code{\link{DNAbin}} or \code{\link{AAbin}}
#' containing unaligned sequences of DNA or amino acids.
#' @param gt \emph{Currently unused.}
#' @param k An integer giving the size of cluster in which the dataset is split.
#' @param bootstrap An integer giving the number of bootstrap replicates.
#' @param method A character string choosing a method of the alignment program;
#' see \code{\link{mafft}} for possible options.
#' @param exec A character string giving the path to the alignment program
#' executable.
#' @param ncore An integer giving the number of cores to use in parallel mode.
#' @return An object of class \code{"\link[=polentaDNA-class]{polentaDNA}"}.
#' @seealso \code{\link{extractMSA}} for extractiong the multiple sequence
#' alignment of a \code{"\link[=polentaDNA-class]{polentaDNA}"} object.
#' @importFrom ape del.gaps dist.dna dist.aa
#' @importFrom igraph as_edgelist
#' @importFrom ips mafft mafft.merge
#' @export
polenta <- function(seqs, gt, k = 200, bootstrap = 100,
method = "auto", exec, ncore){
## remove gaps from aligned sequences
## ----------------------------------
if (is.matrix(seqs)) {
seqs <- del.gaps(seqs)
}
## less than k species will be aligned with MAFFT-LINSI
## ----------------------------------------------------
if (length(seqs) <= k){
cat(length(seqs), "species will be aligned with MAFFT L-INS-i\n")
seqs <- guidance(seqs, ncore = ncore,
bootstrap = bootstrap,
method = method, msa.exec = exec)
## more than k species will be aligned with PASTA
## ----------------------------------------------
} else {
## This is a quick hack to get an inital guide tree
## Should be replaced by the method used by Mirarab and Warnow
## or perhaps a hybrid with taxonomy.
if (missing(gt)){
gt <- mafft(seqs, method = "auto")
if (inherits(gt, "DNAbin")){
gt <- dist.dna(gt, model = "F81")
} else {
gt <- dist.aa(gt)
}
gt <- nj(gt)
}
## split dataset in subsets of size <= k
## -------------------------------------
subtrees <- centroidDecomposition(gt, k = k)
subtrees <- lapply(subtrees, function(z) z$tip.label)
names(subtrees) <- paste0("S", seq_along(subtrees))
## alignment of subtrees
## ---------------------
foo <- function(seqs, taxa){
guidance(seqs[taxa], ncore = ncore,
bootstrap = bootstrap,
method = method, msa.exec = exec)
}
seqs <- lapply(subtrees, foo, seqs = seqs)
names(seqs) <- names(subtrees)
## special case: there are only two subMSAa and they will be simply merged
## with out transitivity merging
if (length(seqs) == 2){
seqlist <- extractMSA(seqs)
seqlist <- list(mafft.merge(seqlist, method = method, exec = exec))
names(seqlist) <- paste(names(seqs), collapse = "-")
seqs <- reappendScores(names(seqs), merged = seqlist, scored = seqs)
} else {
## compute spanning tree of subsets
## --------------------------------
st <- spanningTree(gt, subtrees)
# save.image("devworkspace.rda")
## do profile-alignment
## --------------------
e <- igraph::as_edgelist(st)
merger <- function(seqlist, index, exec){
mafft.merge(seqlist[index], exec = exec)
}
seqlist <- extractMSA(seqs)
seqlist <- apply(e, 1, merger, seqlist = seqlist, exec = exec)
names(seqlist) <- paste(e[, 1], e[, 2], sep = "-")
## reappend scores to merged alignments
## ------------------------------------
seqs <- apply(e, 1, reappendScores, merged = seqlist, scored = seqs)
names(seqs) <- paste(e[, 1], e[, 2], sep = "-")
## do transitivity merging
## -----------------------
# load("devworkspace.rda")
## calculate pairings for transitivity merging
## this is probably very inefficient
vertex.set <- strsplit(names(seqs), "-")
pairings <- function(z){
obj <- list(); meta <- list()
for (i in 1:(length(z) - 1)){
zz <- sapply(z, intersect, x = z[[i]])
zz <- sapply(zz, length)
p <- which(zz > 0)
p <- p[p > i][1]
if (is.na(p)) next
p <- c(i, p)
meta <- c(meta, list(sort(unique(unlist(z[p])))))
obj <- c(obj, list(p))
}
attr(obj, "vertices") <- meta
obj
}
while (length(seqs) > 1){
p <- pairings(vertex.set)
seqs <- lapply(p, transitivityMerge, x = seqs)
vertex.set <- attr(p, "vertices")
}
seqs <- seqs[[1]]
## next steps
## - put pairing() in its on file
## - testing
## - parallelisation
## - make pairing more efficient
# save.image("devworkspace.rda")
}
}
seqs
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.