Skip to content

seqMismatch functions rename #139

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
ssnn-airr merged 2 commits into from
Sep 28, 2026
Merged

seqMismatch functions rename #139

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
98e6978
Refactor seqMismatch functions: rename Rcpp exports and update docume…
ayeletperes Sep 28, 2026
c688b75
Merge branch 'master' into mismatch-wrappers
ayeletperes Sep 28, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 6 additions & 6 deletions R/RcppExports.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -78,16 +78,16 @@ nonsquareDistRcpp <- function(seq, indx, dist_mat) {
.Call(`_alakazam_nonsquareDistRcpp`, seq, indx, dist_mat)
}

seqMismatchCountRcpp <- function(samples, germlines, ignore = as.character( c("N", ".", "-"))) {
.Call(`_alakazam_seqMismatchCountRcpp`, samples, germlines, ignore)
seqMismatchCountRcpp <- function(samples, germlines, ignore = as.character( c("N", "-", ".", "?")), count_trailing = FALSE) {
.Call(`_alakazam_seqMismatchCountRcpp`, samples, germlines, ignore, count_trailing)
}

seqMismatchMatrixRcpp <- function(samples, germlines, ignore = as.character( c("N", ".", "-"))) {
.Call(`_alakazam_seqMismatchMatrixRcpp`, samples, germlines, ignore)
seqMismatchMatrixRcpp <- function(samples, germlines, ignore = as.character( c("N", "-", ".", "?")), count_trailing = FALSE) {
.Call(`_alakazam_seqMismatchMatrixRcpp`, samples, germlines, ignore, count_trailing)
}

seqMismatchPositionsRcpp <- function(samples, germlines, ignore = as.character( c("N", ".", "-"))) {
.Call(`_alakazam_seqMismatchPositionsRcpp`, samples, germlines, ignore)
seqMismatchPositionsRcpp <- function(samples, germlines, ignore = as.character( c("N", "-", ".", "?")), count_trailing = FALSE) {
.Call(`_alakazam_seqMismatchPositionsRcpp`, samples, germlines, ignore, count_trailing)
}

countSeqsWithInvalidBases_rcpp <- function(seqs) {
Expand Down
59 changes: 58 additions & 1 deletion R/Sequence.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1376,4 +1376,61 @@ seqMismatchPositions <- function(samples, germlines, ignore=c("N", ".", "-")) {
#' @export
nonsquareDist <- function(seq, indx, dist_mat=getDNAMatrix()) {
nonsquareDistRcpp(seq, indx, dist_mat)
}
}

#' Count or locate mismatches between sample and germline sequences
#'
#' \code{seqMismatchCount} counts Hamming-style mismatches between paired sample
#' and germline sequences, \code{seqMismatchMatrix} counts them between every
#' sample and every germline, and \code{seqMismatchPositions} returns the
#' mismatch positions of paired sequences.
#'
#' @param samples character vector of sample sequences.
#' @param germlines character vector of germline sequences. For
#' \code{seqMismatchCount} and \code{seqMismatchPositions},
#' a single germline is recycled across all samples.
#' @param ignore vector of characters to ignore, in either sequence.
#' Default is to ignore \code{c("N", "-", ".", "?")}.
#' @param count_trailing if \code{TRUE}, sample positions past the end of a
#' shorter germline count as mismatches, so a germline
#' gains nothing from ending early. If \code{FALSE},
#' sequences are compared only through the length of the
#' shorter one.
#'
#' @return \code{seqMismatchCount}: an integer vector of mismatch counts.
#' \code{seqMismatchMatrix}: an integer matrix of mismatch counts, samples
#' in rows and germlines in columns.
#' \code{seqMismatchPositions}: a list of integer vectors of 1-based
#' mismatch positions.
#'
#' @details Comparisons are case-insensitive. A missing (\code{NA}) sample or
#' germline gives \code{NA}.
#'
#' @examples
#' seqMismatchCount(c("ATGGC", "ATGGN"), "ATGGC")
#' seqMismatchMatrix(c("ATGGC", "ATGGN"), c("ATGGC", "ATGGG"))
#' seqMismatchPositions("ATGGCA", "ATGG")
#'
#' # A germline that ends early is not rewarded for it
#' seqMismatchMatrix("ATGGCA", c(full="ATGGCC", short="ATGG"))
#' seqMismatchMatrix("ATGGCA", c(full="ATGGCC", short="ATGG"), count_trailing=TRUE)
#'
#' @export
seqMismatchCount <- function(samples, germlines, ignore=c("N", "-", ".", "?"),
count_trailing=FALSE) {
seqMismatchCountRcpp(samples, germlines, ignore=ignore, count_trailing=count_trailing)
}

#' @rdname seqMismatchCount
#' @export
seqMismatchMatrix <- function(samples, germlines, ignore=c("N", "-", ".", "?"),
count_trailing=FALSE) {
seqMismatchMatrixRcpp(samples, germlines, ignore=ignore, count_trailing=count_trailing)
}

#' @rdname seqMismatchCount
#' @export
seqMismatchPositions <- function(samples, germlines, ignore=c("N", "-", ".", "?"),
count_trailing=FALSE) {
seqMismatchPositionsRcpp(samples, germlines, ignore=ignore, count_trailing=count_trailing)
}
27 changes: 15 additions & 12 deletions src/RcppExports.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -74,41 +74,44 @@ BEGIN_RCPP
END_RCPP
}
// seqMismatchCountRcpp
IntegerVector seqMismatchCountRcpp(CharacterVector samples, CharacterVector germlines, CharacterVector ignore);
RcppExport SEXP _alakazam_seqMismatchCountRcpp(SEXP samplesSEXP, SEXP germlinesSEXP, SEXP ignoreSEXP) {
IntegerVector seqMismatchCountRcpp(CharacterVector samples, CharacterVector germlines, CharacterVector ignore, bool count_trailing);
RcppExport SEXP _alakazam_seqMismatchCountRcpp(SEXP samplesSEXP, SEXP germlinesSEXP, SEXP ignoreSEXP, SEXP count_trailingSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
Rcpp::traits::input_parameter< CharacterVector >::type samples(samplesSEXP);
Rcpp::traits::input_parameter< CharacterVector >::type germlines(germlinesSEXP);
Rcpp::traits::input_parameter< CharacterVector >::type ignore(ignoreSEXP);
rcpp_result_gen = Rcpp::wrap(seqMismatchCountRcpp(samples, germlines, ignore));
Rcpp::traits::input_parameter< bool >::type count_trailing(count_trailingSEXP);
rcpp_result_gen = Rcpp::wrap(seqMismatchCountRcpp(samples, germlines, ignore, count_trailing));
return rcpp_result_gen;
END_RCPP
}
// seqMismatchMatrixRcpp
IntegerMatrix seqMismatchMatrixRcpp(CharacterVector samples, CharacterVector germlines, CharacterVector ignore);
RcppExport SEXP _alakazam_seqMismatchMatrixRcpp(SEXP samplesSEXP, SEXP germlinesSEXP, SEXP ignoreSEXP) {
IntegerMatrix seqMismatchMatrixRcpp(CharacterVector samples, CharacterVector germlines, CharacterVector ignore, bool count_trailing);
RcppExport SEXP _alakazam_seqMismatchMatrixRcpp(SEXP samplesSEXP, SEXP germlinesSEXP, SEXP ignoreSEXP, SEXP count_trailingSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
Rcpp::traits::input_parameter< CharacterVector >::type samples(samplesSEXP);
Rcpp::traits::input_parameter< CharacterVector >::type germlines(germlinesSEXP);
Rcpp::traits::input_parameter< CharacterVector >::type ignore(ignoreSEXP);
rcpp_result_gen = Rcpp::wrap(seqMismatchMatrixRcpp(samples, germlines, ignore));
Rcpp::traits::input_parameter< bool >::type count_trailing(count_trailingSEXP);
rcpp_result_gen = Rcpp::wrap(seqMismatchMatrixRcpp(samples, germlines, ignore, count_trailing));
return rcpp_result_gen;
END_RCPP
}
// seqMismatchPositionsRcpp
List seqMismatchPositionsRcpp(CharacterVector samples, CharacterVector germlines, CharacterVector ignore);
RcppExport SEXP _alakazam_seqMismatchPositionsRcpp(SEXP samplesSEXP, SEXP germlinesSEXP, SEXP ignoreSEXP) {
List seqMismatchPositionsRcpp(CharacterVector samples, CharacterVector germlines, CharacterVector ignore, bool count_trailing);
RcppExport SEXP _alakazam_seqMismatchPositionsRcpp(SEXP samplesSEXP, SEXP germlinesSEXP, SEXP ignoreSEXP, SEXP count_trailingSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
Rcpp::traits::input_parameter< CharacterVector >::type samples(samplesSEXP);
Rcpp::traits::input_parameter< CharacterVector >::type germlines(germlinesSEXP);
Rcpp::traits::input_parameter< CharacterVector >::type ignore(ignoreSEXP);
rcpp_result_gen = Rcpp::wrap(seqMismatchPositionsRcpp(samples, germlines, ignore));
Rcpp::traits::input_parameter< bool >::type count_trailing(count_trailingSEXP);
rcpp_result_gen = Rcpp::wrap(seqMismatchPositionsRcpp(samples, germlines, ignore, count_trailing));
return rcpp_result_gen;
END_RCPP
}
Expand Down Expand Up @@ -163,9 +166,9 @@ static const R_CallMethodDef CallEntries[] = {
{"_alakazam_seqDistRcpp", (DL_FUNC) &_alakazam_seqDistRcpp, 3},
{"_alakazam_pairwiseDistRcpp", (DL_FUNC) &_alakazam_pairwiseDistRcpp, 2},
{"_alakazam_nonsquareDistRcpp", (DL_FUNC) &_alakazam_nonsquareDistRcpp, 3},
{"_alakazam_seqMismatchCountRcpp", (DL_FUNC) &_alakazam_seqMismatchCountRcpp, 3},
{"_alakazam_seqMismatchMatrixRcpp", (DL_FUNC) &_alakazam_seqMismatchMatrixRcpp, 3},
{"_alakazam_seqMismatchPositionsRcpp", (DL_FUNC) &_alakazam_seqMismatchPositionsRcpp, 3},
{"_alakazam_seqMismatchCountRcpp", (DL_FUNC) &_alakazam_seqMismatchCountRcpp, 4},
{"_alakazam_seqMismatchMatrixRcpp", (DL_FUNC) &_alakazam_seqMismatchMatrixRcpp, 4},
{"_alakazam_seqMismatchPositionsRcpp", (DL_FUNC) &_alakazam_seqMismatchPositionsRcpp, 4},
{"_alakazam_countSeqsWithInvalidBases_rcpp", (DL_FUNC) &_alakazam_countSeqsWithInvalidBases_rcpp, 1},
{"_alakazam_fastDist_rcpp", (DL_FUNC) &_alakazam_fastDist_rcpp, 1},
{"_alakazam_countAASeqsWithInvalidChars_rcpp", (DL_FUNC) &_alakazam_countAASeqsWithInvalidChars_rcpp, 1},
Expand Down
51 changes: 37 additions & 14 deletions src/RcppSeqMismatch.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,8 @@ inline bool isIgnored(char c, const std::vector<char>& ignore) {

inline int countMismatches(std::string sample,
std::string germline,
const std::vector<char>& ignore) {
const std::vector<char>& ignore,
bool count_trailing) {
int sample_len = sample.length();
int germline_len = germline.length();
int len = std::min(sample_len, germline_len);
Expand All @@ -61,13 +62,23 @@ inline int countMismatches(std::string sample,
}
}

// Sample positions past the end of the germline cannot match it
if (count_trailing) {
for (int i = len; i < sample_len; i++) {
if (!isIgnored(sample[i], ignore)) {
count++;
}
}
}

return count;
}


inline std::vector<int> findMismatchPositions(std::string sample,
std::string germline,
const std::vector<char>& ignore) {
const std::vector<char>& ignore,
bool count_trailing) {
int sample_len = sample.length();
int germline_len = germline.length();
int len = std::min(sample_len, germline_len);
Expand All @@ -87,15 +98,25 @@ inline std::vector<int> findMismatchPositions(std::string sample,
}
}

// Sample positions past the end of the germline cannot match it
if (count_trailing) {
for (int i = len; i < sample_len; i++) {
if (!isIgnored(sample[i], ignore)) {
pos.push_back(i + 1);
}
}
}

return pos;
}


// seqMismatchCount
// Internal; the documented, exported interface is seqMismatchCount() in R/Sequence.R.
// [[Rcpp::export]]
IntegerVector seqMismatchCountRcpp(CharacterVector samples,
CharacterVector germlines,
CharacterVector ignore=CharacterVector::create("N", ".", "-")) {
CharacterVector germlines,
CharacterVector ignore=CharacterVector::create("N", "-", ".", "?"),
bool count_trailing=false) {
int n = samples.length();
int m = germlines.length();

Expand All @@ -122,18 +143,19 @@ IntegerVector seqMismatchCountRcpp(CharacterVector samples,
std::string sample = as<std::string>(samples[i]);
std::string germline = as<std::string>(germlines[germline_i]);

counts[i] = countMismatches(sample, germline, ignore_chars);
counts[i] = countMismatches(sample, germline, ignore_chars, count_trailing);
}

return counts;
}


// seqMismatchMatrix
// Internal; the documented, exported interface is seqMismatchMatrix() in R/Sequence.R.
// [[Rcpp::export]]
IntegerMatrix seqMismatchMatrixRcpp(CharacterVector samples,
CharacterVector germlines,
CharacterVector ignore=CharacterVector::create("N", ".", "-")) {
CharacterVector germlines,
CharacterVector ignore=CharacterVector::create("N", "-", ".", "?"),
bool count_trailing=false) {
int n = samples.length();
int m = germlines.length();

Expand All @@ -157,7 +179,7 @@ IntegerMatrix seqMismatchMatrixRcpp(CharacterVector samples,
}

std::string germline = as<std::string>(germlines[j]);
counts(i, j) = countMismatches(sample, germline, ignore_chars);
counts(i, j) = countMismatches(sample, germline, ignore_chars, count_trailing);
}
}

Expand All @@ -169,11 +191,12 @@ IntegerMatrix seqMismatchMatrixRcpp(CharacterVector samples,
}


// seqMismatchPositions
// Internal; the documented, exported interface is seqMismatchPositions() in R/Sequence.R.
// [[Rcpp::export]]
List seqMismatchPositionsRcpp(CharacterVector samples,
CharacterVector germlines,
CharacterVector ignore=CharacterVector::create("N", ".", "-")) {
CharacterVector germlines,
CharacterVector ignore=CharacterVector::create("N", "-", ".", "?"),
bool count_trailing=false) {
int n = samples.length();
int m = germlines.length();

Expand All @@ -200,7 +223,7 @@ List seqMismatchPositionsRcpp(CharacterVector samples,
std::string sample = as<std::string>(samples[i]);
std::string germline = as<std::string>(germlines[germline_i]);

positions[i] = wrap(findMismatchPositions(sample, germline, ignore_chars));
positions[i] = wrap(findMismatchPositions(sample, germline, ignore_chars, count_trailing));
}

return positions;
Expand Down
16 changes: 16 additions & 0 deletions tests/testthat/testSequence.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -196,6 +196,22 @@ test_that("seqMismatchPositions", {
expect_equal(obs[[3]], integer(0))
})

test_that("seqMismatch ignore default and count_trailing", {
# "?" is ignored by default, as in seqEqual and collapseDuplicates
expect_equal(seqMismatchCount("AC?T", "ACGT"), 0L)
expect_equal(seqMismatchCount("AC?T", "ACGT", ignore="N"), 1L)

# A germline that ends early is compared only through its length...
samples <- c("ACGTAA", "ACGTAA")
expect_equal(seqMismatchCount(samples, c("ACGTCC", "ACGT")), c(2L, 0L))
# ...unless sample positions past its end count as mismatches
expect_equal(seqMismatchCount(samples, c("ACGTCC", "ACGT"), count_trailing=TRUE), c(2L, 2L))
expect_equal(seqMismatchCount("ACGTN.", "ACGT", count_trailing=TRUE), 0L)
expect_equal(unname(seqMismatchMatrix("ACGTAA", c("ACGTCC", "ACGT"), count_trailing=TRUE)),
matrix(c(2L, 2L), nrow=1))
expect_equal(seqMismatchPositions("ACGTAA", "ACGT", count_trailing=TRUE)[[1]], c(5L, 6L))
})

#### pairwiseDist ####

test_that("pairwiseDist Nucleotide", {
Expand Down
Loading