options(scipen = 999) if(!require("optparse")) {install.packages("optparse", repos = "http://cran.us.r-project.org")} if(!require("BAMMtools")) {install.packages("BAMMtools", repos = "http://cran.us.r-project.org")} library(optparse) library(BAMMtools) opt_list = list( make_option(c("-e", "--eventdata"), type = "character", default = NULL, help = "name of the event data file", metavar = "character"), make_option(c("-t", "--tree"), type = "character", default = NULL, help = "name of the tree file", metavar = "character"), make_option(c("-c", "--controlfile"), type = "character", default = NULL, help = "name of the control file", metavar = "character"), make_option(c("-b", "--burnin"), type = "double", default = 0, help = "proportion of chain length to discard as burnin", metavar = "double") ) opt_parser = OptionParser(option_list = opt_list) opt = parse_args(opt_parser) if(is.null(opt$eventdata)) { print_help(opt_parser) stop("Name of the event data file must be supplied.n", call. = FALSE) } if(is.null(opt$tree)) { print_help(opt_parser) stop("Name of the tree file must be supplied.n", call. = FALSE) } if(is.null(opt$controlfile)) { print_help(opt_parser) stop("Name of the control file must be supplied.n", call. = FALSE) } if(is.null(opt$burnin)) { opt$burnin <- 0 } ctl_file <- readLines(opt$controlfile) line <- grep("expectedNumberOfShifts", ctl_file) expNumShifts <- as.numeric(strsplit(ctl_file[line], " = ")[[1]][2]) getCredShiftSet <- function(tree, eventdata, expectedNumberOfShifts, burnin) { file_prefix <- ifelse(grepl("renumbered", eventdata), strsplit(eventdata, "event_data_renumbered.txt")[[1]][1], strsplit(eventdata, "event_data.txt")[[1]][1]) ed <- getEventData(tree, eventdata, burnin, nsamples = 500) css <- credibleShiftSet(ed, expectedNumberOfShifts, threshold = 5, set.limit = 0.95) output <- data.frame(Rank = 1:(css$number.distinct), Probability = css$frequency, Cumulative = css$cumulative, `Core shifts` = sapply(css$shiftnodes, function(x) length(x)), stringsAsFactors = F) write.table(output, paste0(file_prefix, "credshiftsets.txt"), quote = F, row.names = F, sep = "\t") } getCredShiftSet(opt$tree, opt$eventdata, expNumShifts, opt$burnin)