## This code is part of the megaptera package
## © C. Heibl 2014 (last update 2017-03-28)
#' @export
speciesConsensus <- function(megProj, spec){
md5 <- spec[3]
n <- spec[2]
spec <- spec[1]
## PARAMETERS
## -----------
gene <- megProj@locus@sql
acc.tab <- paste("acc", gsub("^_", "", gene), sep = "_")
align.exe <- megProj@align.exe
max.bp <- megProj@params@max.bp
min.identity <- megProj@locus@min.identity
min.coverage <- megProj@locus@min.coverage
logfile <- paste0("log/", gene, "stepF.log")
## Read species alignments
## -----------------------
obj <- dbReadDNA(spec, x = megProj, tab.name = acc.tab,
max.bp = max.bp,
min.identity = min.identity,
min.coverage = min.coverage)
if (!is.matrix(obj)) obj <- mafft(obj, exec = megProj@align.exe)
## Species consensus sequences
## ---------------------------
obj <- specCons(obj, log = logfile)
obj <- list(obj)
names(obj) <- spec
class(obj) <- "DNAbin"
## Write MSA to database
## ---------------------
dbWriteMSA(megProj, obj, n = n, md5 = md5)
cat(paste0("\n>", paste(gsub(" ", "_", names(obj)))), file = "aaa.fas", append = TRUE)
cat(paste0("\n", seqinr::c2s(unlist(as.character(obj)))), file = "aaa.fas", append = TRUE)
# z <- read.FASTA("aaa.fas")
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.