R/stage4-tools.R

In ropensci/phylotaR: Automated Phylogenetic Sequence Cluster Identification from 'GenBank'

#' @name seeds_blast
#' @title BLAST seed sequences
#' @description Runs all-v-all blast for seed sequences.
#' @param sqs All seed sequences to be BLASTed
#' @template ps
#' @family run-private
#' @return blast res data.frame
seeds_blast <- function(sqs, ps) {
  info(lvl = 2, ps = ps, "BLASTing [", length(sqs@ids), " sqs]")
  dbfl <- "seeds-db.fa"
  outfl <- "seeds-db-blastout.txt"
  file.path(ps[["wd"]], "blast", dbfl)
  blastdb_gen(sqs, dbfl = dbfl, ps = ps)
  blast_res <- blastn_run(dbfl = dbfl, outfl = outfl, ps = ps)
  blast_res
}

#' @name clstrs_join
#' @title Join clusters for merging
#' @description Uses seed sequence BLAST results and IDs to join clusters
#' identified as sisters into single clusters. Resulting object is of joined
#' clusters, merging is required to reformat the clusters for subsequent
#' analysis.
#' @param blast_res Seed sequence BLAST results
#' @param seed_ids Seed sequence IDs
#' @param all_clstrs List of all clusters
#' @template ps
#' @family run-private
#' @return list of joined clusters
clstrs_join <- function(blast_res, seed_ids, all_clstrs, ps) {
  join <- function(x) {
    pull <- seed_ids %in% x[["sids"]]
    jnd_clstr <- all_clstrs[pull]
    nms <- slotNames(jnd_clstr[[1]])
    clstr <- lapply(nms, function(nm) {
      unlist(lapply(jnd_clstr, function(cl) slot(cl, nm)))
    })
    names(clstr) <- nms
    clstr[["typ"]] <- "merged"
    clstr[["seed"]] <- x[["seed"]]
    # ensure no dups seqs in joined cluster
    pull <- !duplicated(clstr[["sids"]])
    clstr[["sids"]] <- clstr[["sids"]][pull]
    clstr[["txids"]] <- clstr[["txids"]][pull]
    clstr
  }
  pull <- blast_res[["query.id"]] != blast_res[["subject.id"]] &
    blast_res[["qcovs"]] > ps[["mncvrg"]]
  blast_res <- blast_res[pull, ]
  clstr_list <- blast_clstr(blast_res = blast_res)
  info(lvl = 2, ps = ps, "Identified [", length(clstr_list), "] clusters")
  lapply(clstr_list, join)
}

#' @name clstrs_merge
#' @title Merge joined clusters
#' @description Takes a list of joined clusters and computes each
#' data slot to create a single merged cluster. txdct is required for
#' parent look-up.
#' @param jnd_clstrs List of joined clusters
#' @param txdct Taxonomic dictionary
#' @return list of ClstrRecs
#' @family run-private
clstrs_merge <- function(jnd_clstrs, txdct) {
  mrg_clstrs <- vector("list", length = length(jnd_clstrs))
  for (i in seq_along(jnd_clstrs)) {
    cl <- jnd_clstrs[[i]]
    prnt <- parent_get(id = cl[["txids"]], txdct = txdct)
    nsqs <- length(cl[["sids"]])
    ntx <- length(unique(cl[["txids"]]))
    clstrrec <- new("ClstrRec",
      sids = cl[["sids"]], txids = cl[["txids"]],
      nsqs = nsqs, ntx = ntx, typ = "merged", prnt = prnt,
      seed = cl[["seed"]]
    )
    mrg_clstrs[[i]] <- clstrrec
  }
  mrg_clstrs
}

#' @name clstrs_renumber
#' @title Renumber cluster IDs
#' @description Returns a ClstrArc with ID determined by the number
#' of sequences in each cluster.
#' @param clstrrecs List of clusters
#' @return ClstrArc
#' @family run-private
clstrs_renumber <- function(clstrrecs) {
  nsqs <- vapply(clstrrecs, function(x) length(x@sids), numeric(1))
  ord <- order(nsqs, decreasing = TRUE)
  clstrrecs <- clstrrecs[ord]
  for (i in seq_along(clstrrecs)) {
    clstrrecs[[i]]@id <- as.integer(i - 1)
  }
  clstrarc_gen(clstrrecs)
}