mkdist copy.R.txt

#!/usr/bin/env Rscript
suppressPackageStartupMessages({ library(ape) })

args <- commandArgs(trailingOnly=TRUE)
if (length(args) < 3) {
  cat("usage: mkdist.R FASTA MODEL OUT_PREFIX\nmodels: RAW JC69 K80 TN93 HKY GTR\n"); quit(status=2)
}
fasta <- args[1]; model <- toupper(args[2]); outpref <- args[3]
dna <- read.dna(fasta, format="fasta")  # aligned COI

write_out <- function(dm, model, outpref) {
  dm <- as.matrix(dm)
  phylip <- paste0(outpref, "_", model, ".phylip")
  sink(phylip); cat(nrow(dm), "\n", sep=""); sink()
  safe <- function(x) substr(gsub("[^A-Za-z0-9_.-]", "_", x), 1, 10)
  rn <- safe(rownames(dm))
  write.table(cbind(rn, format(dm, scientific=FALSE)),
              file=phylip, append=TRUE, quote=FALSE, sep=" ",
              row.names=FALSE, col.names=FALSE)
  mega <- paste0(outpref, "_", model, ".mega.csv")
  write.table(dm, file=mega, sep=",", col.names=NA, quote=FALSE)
  cat("Wrote:\n", phylip, "\n", mega, "\n", sep="")
}

do_ape <- function(dna, ape_model, tag=model) {
  dm <- dist.dna(dna, model=ape_model, pairwise.deletion=FALSE, as.matrix=TRUE)
  write_out(dm, tag, outpref)
}

if (model %in% c("RAW","JC69","K80","TN93")) {
  do_ape(dna, model, model)
} else if (model %in% c("HKY","GTR")) {
  ok <- FALSE
  if (suppressWarnings(requireNamespace("phangorn", quietly=TRUE))) {
    suppressPackageStartupMessages(library(phangorn))
    phyd <- try(phyDat(dna, type="DNA"), silent=TRUE)
    if (!inherits(phyd, "try-error")) {
      dm <- try(dist.ml(phyd, model=model), silent=TRUE)
      if (!inherits(dm, "try-error")) { write_out(dm, model, outpref); ok <- TRUE }
    }
  }
  if (!ok) {
    if (model == "HKY") {
      message("HKY not supported by local phangorn::dist.ml(); falling back to ape F84.")
      do_ape(dna, "F84", "HKY_F84")
    } else {
      message("GTR not supported by local phangorn::dist.ml(); falling back to ape LOGDET.")
      do_ape(dna, "LOGDET", "GTR_LOGDET")
    }
  }
} else {
  stop("Unknown/unsupported model: ", model)
}