Parallel_quartetTable_byTree.R

Parallel_quartetTable_byTree<-function (trees, taxonnames = NULL, 
                                         epsilon = 0, random = 0,
                                         parallel=FALSE,RAM_Gigs=1.0) {
  if (random < 0) 
    stop("Parameter 'random' must be non-negative.")
  random = round(random)
  if (is.null(taxonnames)) {
    taxonnames = trees[[1]]$tip.label
    message("Using taxa that appear on first gene tree.")
  }
  taxonnames = sort(taxonnames)
  nt = length(trees)
  N = length(taxonnames)
  if (random > 0) 
    M = random
  else M = choose(N, 4)
  Q = matrix(0, M, N + 4)
  qnames = c("12|34", "13|24", "14|23", "1234")
  colnames(Q) = c(taxonnames, qnames)
  warnMissing = 0
  if (random == 0) {
    m = 0
    for (i in 1:(N - 3)) {
      for (j in (i + 1):(N - 2)) {
        for (k in (j + 1):(N - 1)) {
          for (l in (k + 1):N) {
            m = m + 1
            Q[m, c(i, j, k, l)] = 1
          }
        }
      }
    }
  }
  else {
    i = 1
    while (i <= random) {
      q = sample(N, size = 4, replace = FALSE)
      row = integer(N)
      row[q] = 1
      j = 1
      while (j < i) {
        if (identical(Q[j, 1:N], row)) 
          j = i + 1
        else j = j + 1
      }
      if (j == i) {
        Q[i, 1:N] = row
        i = i + 1
      }
    }
  }
  
  
  
  message("Counting occurrences of displayed quartets for ", 
          M, " four-taxon subsets of ", N, " taxa across ", nt, 
          " gene trees.")
  start_time <- Sys.time()
 if(parallel==TRUE){
   options(future.globals.maxSize = RAM_Gigs * 1024^3) #1024^3 = ca. 1G RAM
   fopts <- options()[c("future.globals.maxSize")]
   array.out<-future_sapply(1:nt,FUN=function(numt,Q) {
     ## Set 'fopts' options (here fopts is a global variable)
       options(fopts)
      ## Show that those options are set
       cat(sprintf("test = %s\n", getOption("test")))
    t = trees[[numt]]
    if (is.null(t$edge.length)) {
      t = compute.brlen(t, 1)
    }
    else {
      if (sum(t$edge.length < 0) > 0) 
        stop("Error: Negative branch length in tree")
    }
    zeros = which(t$edge.length <= epsilon)
    t$edge.length[] = 1
    t$edge.length[zeros] = 0
    d = cophenetic.phylo(t)
    for (m in 1:M) {
      tax = as.character(taxonnames[which(Q[m, 1:N] == 
                                            1)])
      if (all(tax %in% colnames(d))) {
        a = d[tax[1], tax[2]] + d[tax[3], tax[4]]
        b = d[tax[1], tax[3]] + d[tax[2], tax[4]]
        c = d[tax[1], tax[4]] + d[tax[2], tax[3]]
        z = sort(c(a, b, c))
        if (z[1] == z[2]) {
          Q[m, "1234"] = Q[m, "1234"] + 1
        }
        else {
          if (z[1] == a) {
            Q[m, "12|34"] = Q[m, "12|34"] + 1
          }
          else {
            if (z[1] == b) {
              Q[m, "13|24"] = Q[m, "13|24"] + 1
            }
            else {
              Q[m, "14|23"] = Q[m, "14|23"] + 1
            }
          }
        }
      }
      else warnMissing = 1
    }
    return(Q[,c("12|34", "13|24", "14|23", "1234")])
    # return(Q)
  },Q=Q,simplify = "array")
 }else{
   array.out<-sapply(1:nt,FUN=function(numt,Q) {
     t = trees[[numt]]
     if (is.null(t$edge.length)) {
       t = compute.brlen(t, 1)
     }
     else {
       if (sum(t$edge.length < 0) > 0) 
         stop("Error: Negative branch length in tree")
     }
     zeros = which(t$edge.length <= epsilon)
     t$edge.length[] = 1
     t$edge.length[zeros] = 0
     d = cophenetic.phylo(t)
     for (m in 1:M) {
       tax = as.character(taxonnames[which(Q[m, 1:N] == 
                                             1)])
       if (all(tax %in% colnames(d))) {
         a = d[tax[1], tax[2]] + d[tax[3], tax[4]]
         b = d[tax[1], tax[3]] + d[tax[2], tax[4]]
         c = d[tax[1], tax[4]] + d[tax[2], tax[3]]
         z = sort(c(a, b, c))
         if (z[1] == z[2]) {
           Q[m, "1234"] = Q[m, "1234"] + 1
         }
         else {
           if (z[1] == a) {
             Q[m, "12|34"] = Q[m, "12|34"] + 1
           }
           else {
             if (z[1] == b) {
               Q[m, "13|24"] = Q[m, "13|24"] + 1
             }
             else {
               Q[m, "14|23"] = Q[m, "14|23"] + 1
             }
           }
         }
       }
       else warnMissing = 1
     }
     return(Q[,c("12|34", "13|24", "14|23", "1234")])
     # return(Q)
   },Q=Q,simplify = "array")
 }
  
  array.summed<-rowSums(array.out, dims = 2)
  Q<-cbind(Q[,-which(names(as.data.frame(Q))%in%qnames)],array.summed)
  # Q<-array.out
  
  if (warnMissing == 1) 
    warning("Some taxa missing from some trees.")
  if ((N > 4) && (sum(rowSums(Q[, qnames, drop = FALSE]) == 
                      0) > 0)) 
    warning("Some 4-taxon set not present on any tree.")
  if (sum(Q[, "1234"]) > 0) 
    warning("Some quartets unresolved.")
  current_time = Sys.time()
  elapsedTime = as.numeric(difftime(current_time, start_time, 
                                    units = "mins"))
  # if (elapsedTime > 15) {
    message("Time to process quartets on gene trees was ", 
            elapsedTime, " minutes.")
  # }
  return(Q)
}