Vectorise get_consensus

The two apply() calls over the consensus matrix cost ~4 ms per call
regardless of how many sequences are being consensed, which dominates when
the caller consenses thousands of small groups. max.col(ties.method =
"first") reproduces which.max exactly, so tied columns still resolve
deterministically. Verified byte-identical over 400 random alignments across
every combination of excluded_characters, min_freq and fill.

drop = FALSE keeps a single remaining character class from collapsing the
matrix to a vector, which apply() could not handle.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

Author: SamT123

R/apply_substitutions.R

@@ -7,36 +7,36 @@
 #' @return Modified sequence with substitutions applied
 #'
 #' @export
-apply_substitutions = function(sequence, substitutions, type = "aa") {
-  alphabet = switch(
+apply_substitutions <- function(sequence, substitutions, type = "aa") {
+  alphabet <- switch(
     type,
     "aa" = c(Biostrings::AA_STANDARD, "X", "-"),
     "nt" = c(Biostrings::DNA_BASES, "N", "-"),
     stop("type must be 'aa' or 'nt'")
   )
 
-  new_sequence = sequence
+  new_sequence <- sequence
   for (s in substitutions) {
-    fat = convert_substitution(s, alphabet)
-    old_char = stringr::str_sub(new_sequence, fat$at, fat$at)
+    fat <- convert_substitution(s, alphabet)
+    old_char <- stringr::str_sub(new_sequence, fat$at, fat$at)
     if (!old_char == fat$from) {
       # fmt: skip
       stop("sub = ", s, "; from_char mismatch, expected: ", fat$from, "; got ", old_char)
     }
 
-    stringr::str_sub(new_sequence, fat$at, fat$at) = fat$to
+    stringr::str_sub(new_sequence, fat$at, fat$at) <- fat$to
   }
 
   new_sequence
 }
 
-convert_substitution = function(s, alphabet) {
-  invalid = function(s) {
+convert_substitution <- function(s, alphabet) {
+  invalid <- function(s) {
     stop("Invalid substitution? ", s)
   }
-  from = stringr::str_sub(s, 1, 1)
-  to = stringr::str_sub(s, -1, -1)
-  at = tryCatch(
+  from <- stringr::str_sub(s, 1, 1)
+  to <- stringr::str_sub(s, -1, -1)
+  at <- tryCatch(
     as.integer(stringr::str_sub(s, 2, -2)),
     warning = function(msg) invalid(s)
   )

R/clean_sequences.R

@@ -29,7 +29,7 @@
 #'
 #' @importFrom stringr str_replace_all
 #' @importFrom Biostrings AA_STANDARD
-clean_sequences = function(
+clean_sequences <- function(
   sequences,
   type,
   alphabet = NULL,
@@ -40,26 +40,26 @@ clean_sequences = function(
     stop("type must be either 'aa' or 'nt', got '", type, "'")
   }
 
-  sequences = toupper(sequences)
+  sequences <- toupper(sequences)
 
   if (is.null(replacement_character)) {
-    replacement_character = c(aa = "X", nt = "N")[[type]]
+    replacement_character <- c(aa = "X", nt = "N")[[type]]
   }
 
   if (is.null(alphabet)) {
-    alphabet = list(
+    alphabet <- list(
       aa = c(Biostrings::AA_STANDARD, "-"),
       nt = c("A", "T", "C", "G", "-")
     )[[type]]
   }
 
-  regex = paste0(
+  regex <- paste0(
     "[^",
     paste(alphabet, collapse = ""),
     "]"
   )
 
-  sequences = setNames(
+  sequences <- setNames(
     str_replace_all(
       sequences,
       regex,

R/fasta.R

@@ -20,9 +20,9 @@
 #'
 #' @importFrom readr read_lines
 #' @export
-fast_fasta = function(path, noisy = F) {
-  lines = readr::read_lines(path, skip_empty_rows = FALSE)
-  name_lines = which(substr(lines, 1, 1) == '>')
+fast_fasta <- function(path, noisy = F) {
+  lines <- readr::read_lines(path, skip_empty_rows = FALSE)
+  name_lines <- which(substr(lines, 1, 1) == '>')
   if (noisy) {
     message("Loaded ", length(name_lines), " sequences")
   }
@@ -32,25 +32,25 @@ fast_fasta = function(path, noisy = F) {
       all(name_lines == seq(1, length(lines), 2))
   ) {
     # case: single-line sequences
-    sequences = lines[-name_lines]
-    names(sequences) = lines[name_lines]
+    sequences <- lines[-name_lines]
+    names(sequences) <- lines[name_lines]
   } else {
     # case: multi-line sequences
-    sequences = rep(NA, length(name_lines))
-    names(sequences) = rep(NA, length(name_lines))
+    sequences <- rep(NA, length(name_lines))
+    names(sequences) <- rep(NA, length(name_lines))
 
-    name_lines = c(name_lines, length(lines) + 1)
+    name_lines <- c(name_lines, length(lines) + 1)
 
     for (i in seq_along(name_lines[-1])) {
-      names(sequences)[[i]] = lines[[name_lines[[i]]]]
-      sequences[[i]] = paste(
+      names(sequences)[[i]] <- lines[[name_lines[[i]]]]
+      sequences[[i]] <- paste(
         lines[(name_lines[[i]] + 1):(name_lines[[i + 1]] - 1)],
         collapse = ""
       )
     }
   }
 
-  names(sequences) = substring(names(sequences), 2)
+  names(sequences) <- substring(names(sequences), 2)
   sequences
 }
 
@@ -77,26 +77,26 @@ fast_fasta = function(path, noisy = F) {
 #'
 #' @importFrom readr write_lines
 #' @export
-write_fast_fasta = function(seqs, names = NULL, path, line_width = NULL) {
+write_fast_fasta <- function(seqs, names = NULL, path, line_width = NULL) {
   if (is.null(names)) {
-    names = names(seqs)
+    names <- names(seqs)
   }
 
   stopifnot(length(seqs) == length(names))
 
   if (!is.null(line_width)) {
     # Split sequences into chunks of line_width characters
-    seqs = vapply(
+    seqs <- vapply(
       seqs,
       function(seq) {
         if (nchar(seq) <= line_width) {
           return(seq)
         }
-        seq_chunks = character(ceiling(nchar(seq) / line_width))
+        seq_chunks <- character(ceiling(nchar(seq) / line_width))
         for (i in seq_along(seq_chunks)) {
-          start_pos = (i - 1) * line_width + 1
-          end_pos = min(i * line_width, nchar(seq))
-          seq_chunks[i] = substr(seq, start_pos, end_pos)
+          start_pos <- (i - 1) * line_width + 1
+          end_pos <- min(i * line_width, nchar(seq))
+          seq_chunks[i] <- substr(seq, start_pos, end_pos)
         }
         paste(seq_chunks, collapse = "\n")
       },

R/get_consensus.R

@@ -41,22 +41,15 @@ get_consensus <- function(
 ) {
   sequences <- toupper(sequences)
   cm <- Biostrings::consensusMatrix(sequences)
-  cm <- cm[!rownames(cm) %in% excluded_characters, ]
+  cm <- cm[!rownames(cm) %in% excluded_characters, , drop = FALSE]
 
-  get_aa <- function(x) {
-    names(x)[which.max(x)]
-  }
+  # ties.method = "first" matches which.max, so a tied column resolves the same
+  # way on every run.
+  winners <- max.col(t(cm), ties.method = "first")
+  counts <- cm[cbind(winners, seq_len(ncol(cm)))]
+  highest_frequencies <- counts / colSums(cm)
 
-  get_freq <- function(x) {
-    max(x) / sum(x)
-  }
-
-  highest_frequencies <- setNames(
-    apply(cm, 2, get_freq),
-    apply(cm, 2, get_aa)
-  )
-
-  aas <- names(highest_frequencies)
+  aas <- rownames(cm)[winners]
   aas[!is.finite(highest_frequencies) | highest_frequencies < min_freq] <- fill
   paste0(aas, collapse = "")
 }

R/get_substitutions.R

@@ -36,7 +36,7 @@
 #' @importFrom purrr map2 pmap
 #' @importFrom stringr str_split
 #' @export
-get_substitutions = function(
+get_substitutions <- function(
   sequence_1,
   sequence_2,
   position_map = 1:nchar(sequence_1),
@@ -58,12 +58,12 @@ get_substitutions = function(
     }
   }
 
-  lengths = nchar(c(sequence_1, sequence_2))
+  lengths <- nchar(c(sequence_1, sequence_2))
 
   if (!all(lengths == min(lengths))) {
     message('Trimming to shortest length: ', min(lengths))
-    sequence_1 = substr(sequence_1, 1, min(lengths))
-    sequence_2 = substr(sequence_2, 1, min(lengths))
+    sequence_1 <- substr(sequence_1, 1, min(lengths))
+    sequence_2 <- substr(sequence_2, 1, min(lengths))
   }
 
   if (length(sequence_1) == 1 & length(sequence_2) == 1 & simplify) {
@@ -75,21 +75,21 @@ get_substitutions = function(
     ))
   }
 
-  sequence_2_uniques = unique(sequence_2)
-  idxs = match(sequence_2, sequence_2_uniques)
+  sequence_2_uniques <- unique(sequence_2)
+  idxs <- match(sequence_2, sequence_2_uniques)
 
   # Store original names from sequence_2
   sequence_2_names <- names(sequence_2)
   has_names <- !is.null(sequence_2_names) && !all(is.na(sequence_2_names))
 
-  unique_substitutions = get_substitutions_CASE_MULTIPLE(
+  unique_substitutions <- get_substitutions_CASE_MULTIPLE(
     sequence_1,
     sequence_2_uniques,
     position_map,
     exclude
   )
 
-  all_substitutions = unique_substitutions[idxs]
+  all_substitutions <- unique_substitutions[idxs]
 
   # Preserve names from sequence_2 in output (only for multi-sequence case)
   if (has_names) {
@@ -100,13 +100,13 @@ get_substitutions = function(
 }
 
 
-get_substitutions_CASE_SINGLE = function(
+get_substitutions_CASE_SINGLE <- function(
   sequence_1,
   sequence_2,
   position_map = 1:max(nchar(sequence_1), nchar(sequence_2)),
   exclude = c()
 ) {
-  min_len = min(
+  min_len <- min(
     stringr::str_length(sequence_1),
     stringr::str_length(sequence_2)
   )
@@ -122,7 +122,7 @@ get_substitutions_CASE_SINGLE = function(
   sequence_1_split <- stringr::str_split(sequence_1, '')[[1]][1:min_len]
   sequence_2_split <- stringr::str_split(sequence_2, '')[[1]][1:min_len]
 
-  diffs = unlist(
+  diffs <- unlist(
     purrr::pmap(
       list(
         s1 = sequence_1_split,
@@ -139,14 +139,14 @@ get_substitutions_CASE_SINGLE = function(
 
   # Convert NULL to character(0) for consistent return type
   if (is.null(diffs)) {
-    diffs = character(0)
+    diffs <- character(0)
   }
 
   diffs
 }
 
 
-get_substitutions_CASE_MULTIPLE = function(
+get_substitutions_CASE_MULTIPLE <- function(
   sequence_1,
   sequence_2,
   position_map = 1:nchar(sequence_1),
@@ -157,17 +157,17 @@ get_substitutions_CASE_MULTIPLE = function(
     return(rep(list(character(0)), length(sequence_2)))
   }
 
-  diffs = rep(list(NULL), length(sequence_2))
+  diffs <- rep(list(NULL), length(sequence_2))
 
-  sequence_1_split = stringr::str_split(sequence_1, pattern = "")[[1]]
-  sequence_2_split = Biostrings::AAStringSet(sequence_2)
+  sequence_1_split <- stringr::str_split(sequence_1, pattern = "")[[1]]
+  sequence_2_split <- Biostrings::AAStringSet(sequence_2)
 
   for (i in seq_along(sequence_1_split)) {
-    sequence_2_pos_i = as.character(Biostrings::subseq(sequence_2_split, i, i))
-    incl = sequence_1_split[[i]] != sequence_2_pos_i &
+    sequence_2_pos_i <- as.character(Biostrings::subseq(sequence_2_split, i, i))
+    incl <- sequence_1_split[[i]] != sequence_2_pos_i &
       (!(sequence_2_pos_i %in% exclude | sequence_1_split[[i]] %in% exclude))
 
-    diffs[incl] = purrr::map2(
+    diffs[incl] <- purrr::map2(
       diffs[incl],
       paste0(
         sequence_1_split[[i]],
@@ -179,7 +179,7 @@ get_substitutions_CASE_MULTIPLE = function(
   }
 
   # Convert NULL to character(0) for consistent return type
-  diffs = lapply(diffs, function(x) {
+  diffs <- lapply(diffs, function(x) {
     if (is.null(x)) character(0) else x
   })
 

R/mafft_align.R

@@ -22,7 +22,7 @@
 #'
 #' @note Requires MAFFT to be installed and available in PATH
 #' @export
-mafft_align = function(unaligned_sequences, reference_sequence, noisy = F) {
+mafft_align <- function(unaligned_sequences, reference_sequence, noisy = F) {
   # Input validation
   if (length(reference_sequence) != 1) {
     stop(
@@ -36,28 +36,28 @@ mafft_align = function(unaligned_sequences, reference_sequence, noisy = F) {
     stop("MAFFT is not installed or not available in PATH")
   }
 
-  NA_sequence_locations = is.na(unaligned_sequences)
+  NA_sequence_locations <- is.na(unaligned_sequences)
 
-  temp_mafft_folder = tempdir()
+  temp_mafft_folder <- tempdir()
 
-  reference_seq_file = tempfile(
+  reference_seq_file <- tempfile(
     pattern = "reference",
     fileext = ".fasta",
     tmpdir = temp_mafft_folder
   )
-  unaligned_seqs_file = tempfile(
+  unaligned_seqs_file <- tempfile(
     pattern = "sequences",
     fileext = ".fasta",
     tmpdir = temp_mafft_folder
   )
-  aligned_seqs_file = tempfile(
+  aligned_seqs_file <- tempfile(
     pattern = "aligned",
     fileext = ".fasta",
     tmpdir = temp_mafft_folder
   )
 
-  unique_unaligned_sequences = unique(unaligned_sequences)
-  all_to_unique_sequence_map = match(
+  unique_unaligned_sequences <- unique(unaligned_sequences)
+  all_to_unique_sequence_map <- match(
     unaligned_sequences,
     unique_unaligned_sequences
   )
@@ -88,12 +88,12 @@ mafft_align = function(unaligned_sequences, reference_sequence, noisy = F) {
     )
   )
 
-  unique_aligned_sequences = toupper(fast_fasta(aligned_seqs_file)[-1])
-  all_aligned_sequences = unique_aligned_sequences[all_to_unique_sequence_map]
+  unique_aligned_sequences <- toupper(fast_fasta(aligned_seqs_file)[-1])
+  all_aligned_sequences <- unique_aligned_sequences[all_to_unique_sequence_map]
 
-  all_aligned_sequences[NA_sequence_locations] = NA
+  all_aligned_sequences[NA_sequence_locations] <- NA
 
-  names(all_aligned_sequences) = names(unaligned_sequences)
+  names(all_aligned_sequences) <- names(unaligned_sequences)
 
   all_aligned_sequences
 }