From 975571b956659b43cf1248fb06b04ca897b2e4e6 Mon Sep 17 00:00:00 2001 From: ayelet peres Date: Fri, 19 Jun 2026 20:59:12 -0400 Subject: [PATCH 1/5] update tigger reassign to use alakazam cpp function, included the confidence plot in the genotype, added depth per loci for the genotype inference, update the genotypeFasta function to include genes that were not seen in the genotype inference, included the genotyped_alleles function for the bayesian inference --- profiling/ProfileTigger.R | 15 --------------- 1 file changed, 15 deletions(-) delete mode 100644 profiling/ProfileTigger.R diff --git a/profiling/ProfileTigger.R b/profiling/ProfileTigger.R deleted file mode 100644 index 5f99e08..0000000 --- a/profiling/ProfileTigger.R +++ /dev/null @@ -1,15 +0,0 @@ -# Imports -library(tigger) -library(alakazam) -library(profvis) - -#### Load example data #### - -data("AIRRDb") -data("GermlineIGHV") - -#### Find novel alleles #### - -profvis({ - nv <- findNovelAlleles(AIRRDb, GermlineIGHV,nproc=6) -}) From 5df133e370c60d9823718bc425b9c5b55261c2c6 Mon Sep 17 00:00:00 2001 From: ayelet peres Date: Fri, 19 Jun 2026 21:10:25 -0400 Subject: [PATCH 2/5] Enhance genotype inference and allele reassignment features - Added `genotyped_alleles` parameter to `inferGenotypeBayesian` for including most likely alleles in the output. - Implemented locus-specific checks in `inferGenotype` and `reassignAlleles` to handle mixed loci warnings. - Updated `genotypeFasta` to support unseen alleles inclusion. - Enhanced tests for genotype inference and reassignment functions to validate new features. --- NAMESPACE | 3 + R/bayesian.R | 31 ++- R/functions.R | 308 ++++++++++++++++++++++++---- R/tigger.R | 6 +- tests/testthat/test_coreFunctions.R | 181 +++++++++++++++- 5 files changed, 489 insertions(+), 40 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 6d4c212..c999c5a 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -26,6 +26,7 @@ importFrom(alakazam,checkColumns) importFrom(alakazam,getAllele) importFrom(alakazam,getFamily) importFrom(alakazam,getGene) +importFrom(alakazam,getLocus) importFrom(alakazam,translateDNA) importFrom(doParallel,registerDoParallel) importFrom(dplyr,"%>%") @@ -57,6 +58,7 @@ importFrom(foreach,foreach) importFrom(foreach,registerDoSEQ) importFrom(graphics,plot) importFrom(gridExtra,arrangeGrob) +importFrom(gridExtra,grid.arrange) importFrom(gtools,ddirichlet) importFrom(iterators,icount) importFrom(lazyeval,interp) @@ -88,3 +90,4 @@ importFrom(tidyr,gather) importFrom(tidyr,spread) importFrom(tidyr,unnest) importFrom(utils,citation) +importFrom(utils,head) diff --git a/R/bayesian.R b/R/bayesian.R index 69e3504..6be84f2 100644 --- a/R/bayesian.R +++ b/R/bayesian.R @@ -52,6 +52,10 @@ #' assigned to \code{0}; e.g., the heterozygous case is #' \code{c(priors[1], priors[2], 0, 0)}. The prior for the #' homozygous distribution is fixed at \code{c(1, 0, 0, 0)}. +#' @param genotyped_alleles if \code{TRUE}, add a \code{genotyped_alleles} +#' column containing the most likely alleles based on +#' the highest Bayesian zygosity likelihood. Default +#' is \code{FALSE}. #' #' @return #' A \code{data.frame} of alleles denoting the genotype of the subject with the log10 @@ -70,6 +74,8 @@ #' \item \code{kt}: log10 likelihood that the \code{gene} is trizygous #' \item \code{kq}: log10 likelihood that the \code{gene} is quadrozygous. #' \item \code{k_diff}: log10 ratio of the highest to second-highest zygosity likelihoods. +#' \item \code{genotyped_alleles}: If requested, comma separated list of +#' alleles in the most likely genotype for the given \code{gene}. #' } #' #' @note @@ -97,7 +103,8 @@ inferGenotypeBayesian <- function(data, germline_db=NA, novel=NA, v_call="v_call", seq="sequence_alignment", find_unmutated=TRUE, - priors=c(0.6, 0.4, 0.4, 0.35, 0.25, 0.25, 0.25, 0.25, 0.25)){ + priors=c(0.6, 0.4, 0.4, 0.35, 0.25, 0.25, 0.25, 0.25, 0.25), + genotyped_alleles=FALSE){ # Visibility hack . <- NULL @@ -133,6 +140,17 @@ inferGenotypeBayesian <- function(data, germline_db=NA, novel=NA, } } + # Warn if the calls span more than one locus + if ("locus" %in% colnames(data) && length(allele_calls) == nrow(data)) { + loci <- as.character(data$locus) + } else { + loci <- getLocus(allele_calls, first=TRUE, strip_d=FALSE) + } + loci[is.na(loci)] <- "" + if (length(unique(loci[nzchar(loci)])) > 1) { + warning("Mixed loci detected. Repertoire depth and fractional thresholds are locus specific.") + } + # Find which rows' calls contain which genes gene_regex <- allele_calls %>% strsplit(",") %>% unlist() %>% getGene(strip_d=FALSE) %>% unique() %>% paste("\\*", sep="") @@ -266,6 +284,17 @@ inferGenotypeBayesian <- function(data, germline_db=NA, novel=NA, } } } + if (genotyped_alleles) { + # For each gene, keep the most likely alleles based on the highest + # Bayesian zygosity likelihood (kh, kd, kt, kq) + geno$genotyped_alleles <- apply(geno[, c("alleles", "kh", "kd", "kt", "kq")], 1, + function(row) { + m <- which.max(as.numeric(row[2:5])) + alleles <- unlist(strsplit(row[1], ",")) + m <- min(m, length(alleles)) + paste0(alleles[1:m], collapse=",") + }) + } rownames(geno) <- NULL return(geno) } diff --git a/R/functions.R b/R/functions.R index 0143961..5737520 100644 --- a/R/functions.R +++ b/R/functions.R @@ -1067,13 +1067,34 @@ inferGenotype <- function(data, germline_db=NA, novel=NA, v_call="v_call", } } + # Determine each call's locus so that fractional cutoffs use the + # locus-specific repertoire depth + if ("locus" %in% colnames(data) && length(allele_calls) == nrow(data)) { + call_loci <- as.character(data$locus) + } else { + call_loci <- getLocus(allele_calls, first=TRUE, strip_d=FALSE) + } + call_loci[is.na(call_loci)] <- "" + if (length(unique(call_loci[nzchar(call_loci)])) > 1) { + warning("Mixed loci detected. Repertoire depth and fractional thresholds are locus specific.") + } + # Find which rows' calls contain which genes - cutoff <- ifelse(gene_cutoff < 1, length(allele_calls)*gene_cutoff, gene_cutoff) gene_regex <- allele_calls %>% strsplit(",") %>% unlist() %>% getGene(strip_d=FALSE) %>% unique() %>% paste("\\*", sep="") gene_groups <- sapply(gene_regex, grep, allele_calls, simplify=FALSE) names(gene_groups) <- gsub("\\*", "", gene_regex, fixed=TRUE) - gene_groups <- gene_groups[sapply(gene_groups, length) >= cutoff] + gene_cutoffs <- sapply(names(gene_groups), function(g) { + if (gene_cutoff >= 1) { + gene_cutoff + } else { + g_locus <- getLocus(g, first=TRUE, strip_d=FALSE) + length(call_loci[call_loci == g_locus]) * gene_cutoff + } + }) + gene_groups <- gene_groups[sapply(seq_along(gene_groups), function(i) { + length(gene_groups[[i]]) >= gene_cutoffs[[i]] + })] gene_groups <- gene_groups[sortAlleles(names(gene_groups))] # Make a table to store the resulting genotype @@ -1166,11 +1187,23 @@ inferGenotype <- function(data, germline_db=NA, novel=NA, v_call="v_call", #' @param text_size point size of the plotted text. #' @param silent if \code{TRUE} do not draw the plot and just return the ggplot #' object; if \code{FALSE} draw the plot. +#' @param confidence_col name of a column in \code{genotype} holding a per-gene +#' confidence value to display as evidence, e.g. \code{"k_diff"} +#' from \link{inferGenotypeBayesian}. If \code{NULL} (default), +#' no confidence panel is drawn. The value is binned and shown on a +#' blue color scale, with the column name as the legend title. +#' @param confidence_breaks numeric vector of breaks used to bin the confidence +#' value into color groups. The default +#' \code{c(0, 1, 2, 3, 4, 5, 10, 20, 50, Inf)} matches the binning +#' used by RAbHIT for the haplotype lK panel. #' @param ... additional arguments to pass to ggplot2::theme. #' -#' @return A ggplot object defining the plot. +#' @return A ggplot object defining the plot. If \code{confidence_col} is supplied, a +#' \code{gridExtra} grob is returned instead, placing the confidence panel +#' next to the genotype. \code{confidence_col} is not combined with +#' \code{facet_by}; if both are given, \code{facet_by} is ignored. #' -#' @seealso \link{inferGenotype} +#' @seealso \link{inferGenotype}, \link{inferGenotypeBayesian} #' #' @examples #' # Plot genotype @@ -1183,11 +1216,26 @@ inferGenotype <- function(data, germline_db=NA, novel=NA, v_call="v_call", #' geno_sub <- rbind(genotype_a, genotype_b) #' plotGenotype(geno_sub, facet_by="SUBJECT", gene_sort="pos") #' +#' # Add a confidence evidence panel from a per-gene confidence column +#' geno_conf <- SampleGenotype +#' geno_conf$k_diff <- seq(0, 20, length.out=nrow(geno_conf)) +#' plotGenotype(geno_conf, confidence_col="k_diff") +#' #' @export plotGenotype <- function(genotype, facet_by=NULL, gene_sort=c("name", "position"), - text_size=12, silent=FALSE, ...) { + text_size=12, silent=FALSE, confidence_col=NULL, + confidence_breaks=c(0, 1, 2, 3, 4, 5, 10, 20, 50, Inf), ...) { # Check arguments gene_sort <- match.arg(gene_sort) + if (!is.null(confidence_col)) { + if (!confidence_col %in% colnames(genotype)) { + stop("confidence_col '", confidence_col, "' not found in genotype.") + } + if (!is.null(facet_by)) { + warning("confidence_col is not supported together with facet_by; ignoring facet_by.") + facet_by = NULL + } + } # Split genes' alleles into their own rows alleles = strsplit(genotype$alleles, ",") @@ -1202,8 +1250,8 @@ plotGenotype <- function(genotype, facet_by=NULL, gene_sort=c("name", "position" } # Set the gene order - geno2$gene = factor(geno2$gene, - levels=rev(sortAlleles(unique(geno2$gene), method=gene_sort))) + gene_levels = rev(sortAlleles(unique(geno2$gene), method=gene_sort)) + geno2$gene = factor(geno2$gene, levels=gene_levels) # Create the base plot p = ggplot(geno2, aes(x=!!rlang::sym("gene"), @@ -1228,10 +1276,54 @@ plotGenotype <- function(genotype, facet_by=NULL, gene_sort=c("name", "position" # Add additional theme elements p = p + do.call(theme, list(...)) + # Without a confidence column, return the genotype plot as is + if (is.null(confidence_col)) { + if (!silent) { plot(p) } + return(invisible(p)) + } + + # Bin the per-gene confidence value and draw it as a blue color panel beside the + # genotype. White marks NA/unscored genes; drop=FALSE keeps the full scale. + blues = c("#FFFFFF", "#F7FBFF", "#DEEBF7", "#C6DBEF", "#9ECAE1", "#6BAED6", + "#4292C6", "#2171B5", "#08519C", "#08306B") + bins = cut(suppressWarnings(as.numeric(genotype[[confidence_col]])), + confidence_breaks, include.lowest=TRUE, right=FALSE) + bin_levels = gsub(",", ", ", levels(bins)) + conf_levels = c("NA", bin_levels) + conf = data.frame(gene=factor(genotype$gene, levels=gene_levels), + confidence=factor(ifelse(is.na(bins), "NA", bin_levels[bins]), + levels=conf_levels)) + pc = ggplot(conf, aes(x=!!rlang::sym("gene"), + fill=!!rlang::sym("confidence"))) + + theme_bw() + + theme(axis.ticks=element_blank(), + axis.text=element_blank(), + panel.grid.major=element_blank(), + panel.grid.minor=element_blank(), + text=element_text(size=text_size)) + + geom_bar(position="fill") + + coord_flip() + xlab("") + ylab("") + + scale_fill_manual(name=confidence_col, + values=setNames(blues[seq_along(conf_levels)], conf_levels), + drop=FALSE) + + # Align the gene rows and place the confidence panel beside the genotype, with + # both legends moved to the right so the two panels sit next to each other + geno_grob = ggplotGrob(p + theme(legend.position="none")) + conf_grob = ggplotGrob(pc + theme(legend.position="none")) + conf_grob$heights = geno_grob$heights + p_legend = ggplotGrob(p) + pc_legend = ggplotGrob(pc) + leg_allele = p_legend$grobs[[which(p_legend$layout$name %in% c("guide-box", "guide-box-right"))[1]]] + leg_conf = pc_legend$grobs[[which(pc_legend$layout$name %in% c("guide-box", "guide-box-right"))[1]]] + panels = gridExtra::arrangeGrob(geno_grob, conf_grob, ncol=2, widths=c(0.85, 0.15)) + legends = gridExtra::arrangeGrob(leg_allele, leg_conf, ncol=1) + combined = gridExtra::arrangeGrob(panels, legends, ncol=2, widths=c(0.85, 0.15)) + # Plot - if (!silent) { plot(p) } + if (!silent) { gridExtra::grid.arrange(combined) } - invisible(p) + invisible(combined) } #' Return the nucleotide sequences of a genotype @@ -1246,6 +1338,10 @@ plotGenotype <- function(genotype, facet_by=NULL, gene_sort=c("name", "position" #' @param novel an optional \code{data.frame} containing putative #' novel alleles of the type returned by #' \link{findNovelAlleles}. +#' @param include_unseen if \code{TRUE}, include germline database alleles for +#' genes that are not present in \code{genotype}. For +#' genes present in \code{genotype}, include only the +#' genotyped alleles. #' #' @return A named vector of strings containing the germline nucleotide #' sequences of the alleles in the provided genotype. @@ -1257,7 +1353,7 @@ plotGenotype <- function(genotype, facet_by=NULL, gene_sort=c("name", "position" #' genotype_db <- genotypeFasta(SampleGenotype, SampleGermlineIGHV, SampleNovel) #' #' @export -genotypeFasta <- function(genotype, germline_db, novel=NA){ +genotypeFasta <- function(genotype, germline_db, novel=NA, include_unseen=FALSE){ if(!is.null(nrow(novel))){ # Extract novel alleles if any and add them to germline_db novel <- filter(novel, !is.na(!!rlang::sym("polymorphism_call"))) %>% @@ -1273,7 +1369,13 @@ genotypeFasta <- function(genotype, germline_db, novel=NA){ g_names <- names(germline_db) names(g_names) <- getAllele(names(germline_db), first = T, strip_d = T) - table_calls <- mapply(paste, genotype$gene, strsplit(genotype$alleles, ","), + allele_values <- genotype$alleles + if ("genotyped_alleles" %in% colnames(genotype)) { + use_genotyped <- !is.na(genotype$genotyped_alleles) & + nzchar(genotype$genotyped_alleles) + allele_values[use_genotyped] <- genotype$genotyped_alleles[use_genotyped] + } + table_calls <- mapply(paste, genotype$gene, strsplit(allele_values, ","), sep="*") table_calls_names <- unlist(table_calls) seq_names <- g_names[names(g_names) %in% table_calls_names] @@ -1285,6 +1387,11 @@ genotypeFasta <- function(genotype, germline_db, novel=NA){ paste(table_calls_names[not_found], collapse = ", ")) } + if (include_unseen) { + germline_genes <- getGene(names(germline_db), first=TRUE, strip_d=TRUE) + seqs <- c(germline_db[!germline_genes %in% genotype$gene], seqs) + } + return(seqs) } @@ -1321,10 +1428,34 @@ genotypeFasta <- function(genotype, germline_db, novel=NA){ #' (\code{"repertoire"}) assignments should be performed. #' Use of \code{"gene"} increases speed by minimizing required number of #' alignments, as gene level assignments will be maintained when possible. +#' @param trim_seq if \code{TRUE}, trim sample and germline sequences +#' to the segment boundaries before calculating Hamming +#' distance. Boundaries are determined from the segment +#' prefix of \code{v_call}, such as \code{v_*}, +#' \code{d_*}, or \code{j_*} columns. +#' @param overwrite if \code{TRUE}, replace \code{v_call} with reassigned +#' calls instead of writing a \code{*_call_genotyped} +#' column. +#' @param ignored_regex regular expression indicating characters to ignore +#' when comparing sequences. May also be \code{TRUE} to +#' ignore nothing (every position counts), as used for +#' D and J segments. +#' @param treat_multigene_as_uncalled if \code{TRUE}, sequences whose call +#' spans more than one gene are treated as uncalled and +#' realigned against the whole genotype rather than kept +#' at their first gene. Only applies when \code{keep_gene} +#' is \code{"gene"} or \code{"repertoire"}. +#' @param top_k maximum number of equally-best alleles to report per +#' sequence. \code{NULL} (default) keeps all ties. +#' @param top_by how to break ties when more than \code{top_k} alleles +#' are equally close. \code{"alphabetical"} keeps the +#' first \code{top_k} by name; \code{"mutation_count"} +#' keeps all ties. #' #' @return A modified input \code{data.frame} containing the best allele call from #' among the sequences listed in \code{genotype_db} in the -#' \code{v_call_genotyped} column. +#' \code{*_call_genotyped} column, or in \code{v_call} when +#' \code{overwrite=TRUE}. #' #' @examples #' # Extract the database sequences that correspond to the genotype @@ -1338,21 +1469,108 @@ genotypeFasta <- function(genotype, germline_db, novel=NA){ reassignAlleles <- function(data, genotype_db, v_call="v_call", seq="sequence_alignment", method="hamming", path=NA, - keep_gene=c("gene", "family", "repertoire")){ + keep_gene=c("gene", "family", "repertoire"), + trim_seq=FALSE, overwrite=FALSE, + ignored_regex="[\\.N-]", + treat_multigene_as_uncalled=FALSE, + top_k=NULL, top_by=c("alphabetical", "mutation_count")){ # Check arguments keep_gene <- match.arg(keep_gene) + top_by <- match.arg(top_by) + seg <- tolower(substr(v_call, 1, 1)) + if (!seg %in% c("v", "d", "j")) { + stop("Could not determine segment from call column: ", v_call, + ". Expected a column beginning with v, d, or j.") + } + output_col <- paste0(seg, "_call_genotyped") # Extract data subset and prepare output vector - v_sequences = as.character(data[[seq]]) - v_calls = getAllele(data[[v_call]], first=FALSE, strip_d=FALSE) - v_call_genotyped = rep("", length(v_calls)) + v_sequences <- as.character(data[[seq]]) + original_calls <- data[[v_call]] + # The call column typically has very few distinct values relative to the + # number of rows, so parse allele/gene/family on the unique calls and + # expand back via index instead of parsing the full-length vector. + uniq_calls <- unique(original_calls) + uniq_idx <- match(original_calls, uniq_calls) + uniq_alleles <- getAllele(uniq_calls, first=FALSE, strip_d=FALSE) + v_calls <- uniq_alleles[uniq_idx] + v_call_genotyped <- rep("", length(v_calls)) + has_call <- !is.na(v_calls) & nzchar(v_calls) + + if (trim_seq) { + germline_cols <- c(paste0(seg, "_germline_start"), + paste0(seg, "_germline_end")) + if (seq == "sequence_alignment") { + seq_cols <- germline_cols + } else { + seq_cols <- c(paste0(seg, "_sequence_start"), + paste0(seg, "_sequence_end")) + } + missing <- setdiff(unique(c(seq_cols, germline_cols)), colnames(data)) + if (length(missing) > 0) { + stop("Cannot trim sequences: missing columns ", + paste(missing, collapse=", ")) + } + v_sequences <- substr(v_sequences, data[[seq_cols[1]]], + data[[seq_cols[2]]]) + } + + mismatch_matrix <- function(samples, germlines, indices) { + # Map the ignored_regex to the equivalent alakazam Rcpp ignore-set + # (verified to reproduce getMutatedPositions bit-for-bit): + # "[\\.N-]" -> ignore gaps/N (V segments) + # TRUE -> ignore nothing (D/J, where ignored_regex is a logical + # and gregexpr("TRUE", ...) never matches a DNA sequence) + # Anything else falls back to the pure-R path below. + rcpp_ignore <- if (identical(ignored_regex, "[\\.N-]")) { + c(".", "N", "-") + } else if (isTRUE(ignored_regex)) { + character(0) + } else { + NULL + } + use_rcpp <- !is.null(rcpp_ignore) && + isTRUE(getOption("tigger.use_alakazam_rcpp_mismatch", TRUE)) && + exists("seqMismatchCountRcpp", envir=asNamespace("alakazam"), inherits=FALSE) && # this will be removed once alakazam updates + exists("seqMismatchMatrixRcpp", envir=asNamespace("alakazam"), inherits=FALSE) # this will be removed once alakazam updates + if (use_rcpp) { + ignore <- rcpp_ignore + if (!trim_seq) { + dist_mat <- get("seqMismatchMatrixRcpp", envir=asNamespace("alakazam"))( + samples, germlines, ignore=ignore) + } else { + dist_mat <- sapply(germlines, function(x) { + germline_seqs <- substr(rep(x, length(indices)), + data[[germline_cols[1]]][indices], + data[[germline_cols[2]]][indices]) + get("seqMismatchCountRcpp", envir=asNamespace("alakazam"))( + samples, germline_seqs, ignore=ignore) + }) + } + return(matrix(as.integer(dist_mat), nrow=length(samples), + ncol=length(germlines))) + } + + dists <- lapply(germlines, function(x) { + germline_seqs <- if (trim_seq) { + substr(rep(x, length(indices)), data[[germline_cols[1]]][indices], + data[[germline_cols[2]]][indices]) + } else { + x + } + sapply(getMutatedPositions(samples, germline_seqs, + ignored_regex=ignored_regex, + match_instead=FALSE), length) + }) + matrix(unlist(dists), ncol=length(germlines)) + } if (keep_gene == "gene") { - v = getGene(v_calls, first = TRUE, strip_d=FALSE) + v = getGene(uniq_alleles, first = TRUE, strip_d=FALSE)[uniq_idx] geno = getGene(names(genotype_db),strip_d=TRUE) names(geno) = names(genotype_db) } else if (keep_gene == "family") { - v <- getFamily(v_calls, first = TRUE, strip_d = FALSE) + v <- getFamily(uniq_alleles, first = TRUE, strip_d = FALSE)[uniq_idx] geno = getFamily(names(genotype_db),strip_d=TRUE) names(geno) = names(genotype_db) } else if (keep_gene == "repertoire") { @@ -1363,26 +1581,45 @@ reassignAlleles <- function(data, genotype_db, v_call="v_call", stop("Unknown keep_gene value: ", keep_gene) } + # Optionally treat multi-gene calls (e.g. "IGHV1-2*02,IGHV3-7*01") as + # uncalled so they are realigned against the whole genotype rather than + # kept at the first gene. Only meaningful for gene/repertoire grouping. + if (keep_gene %in% c("gene", "repertoire") && treat_multigene_as_uncalled) { + genes_full = getGene(uniq_alleles, first=FALSE, strip_d=FALSE)[uniq_idx] + genes_split = strsplit(genes_full, ",") + is_multigene = vapply(genes_split, + function(x) length(unique(trimws(x))) > 1, logical(1)) + } else { + is_multigene = rep(FALSE, length(v_calls)) + } + + # Optionally cap the number of equally-best alleles reported per sequence. + apply_top_k <- function(best_alleles) { + if (!is.null(top_k) && !is.na(top_k) && top_by == "alphabetical") { + best_alleles <- lapply(best_alleles, function(x) { + if (length(x) > top_k) head(sort(x), top_k) else x + }) + } + best_alleles + } + # keep_gene == FALSE # Find which genotype genes/families are homozygous and assign those alleles first hetero = unique(geno[which(duplicated(geno))]) homo = geno[!(geno %in% hetero)] homo_alleles = names(homo) names(homo_alleles) = homo - homo_calls_i = which(v %in% homo) + homo_calls_i = which(v %in% homo & !is_multigene & has_call) v_call_genotyped[homo_calls_i] = homo_alleles[v[homo_calls_i]] # Now realign the heterozygote sequences to each allele of that gene for (het in hetero){ - ind = which(v %in% het) + ind = which(v %in% het & !is_multigene & has_call) if (length(ind) > 0){ het_alleles = names(geno[which(geno == het)]) het_seqs = genotype_db[het_alleles] if(method == "hamming"){ - dists = lapply(het_seqs, function(x) - sapply(getMutatedPositions(v_sequences[ind], x, match_instead=FALSE), - length)) - dist_mat = matrix(unlist(dists), ncol = length(het_seqs)) + dist_mat <- mismatch_matrix(v_sequences[ind], het_seqs, ind) } else { stop("Only Hamming distance is currently supported as a method.") } @@ -1395,20 +1632,17 @@ reassignAlleles <- function(data, genotype_db, v_call="v_call", for (i in 1:nrow(dist_mat)) { best_match[[i]] = which(dist_mat[i, ]==min(dist_mat[i, ])) } - best_alleles = lapply(best_match, function(x) het_alleles[x]) + best_alleles = apply_top_k(lapply(best_match, function(x) het_alleles[x])) v_call_genotyped[ind] = unlist(lapply(best_alleles, paste, collapse=",")) } } # Now realign the gene-not-in-genotype calls to every genotype allele - hetero_calls_i = which(v %in% hetero) - not_called = setdiff(1:length(v), c(homo_calls_i, hetero_calls_i)) - if(length(not_called)>1){ + hetero_calls_i = which(v %in% hetero & !is_multigene & has_call) + not_called = setdiff(which(has_call), c(homo_calls_i, hetero_calls_i)) + if(length(not_called)>0){ if(method == "hamming"){ - dists = lapply(genotype_db, function(x) - sapply(getMutatedPositions(v_sequences[not_called], x, match_instead=FALSE), - length)) - dist_mat = matrix(unlist(dists), ncol = length(genotype_db)) + dist_mat <- mismatch_matrix(v_sequences[not_called], genotype_db, not_called) } else { stop("Only Hamming distance is currently supported as a method.") } @@ -1421,11 +1655,11 @@ reassignAlleles <- function(data, genotype_db, v_call="v_call", for (i in 1:nrow(dist_mat)) { best_match[[i]] = which(dist_mat[i, ]==min(dist_mat[i, ])) } - best_alleles = lapply(best_match, function(x) names(genotype_db[x])) + best_alleles = apply_top_k(lapply(best_match, function(x) names(genotype_db[x]))) v_call_genotyped[not_called] = unlist(lapply(best_alleles, paste, collapse=",")) } - if (all(v_call_genotyped == data[[v_call]])) { + if (all(v_call_genotyped == original_calls, na.rm=TRUE)) { msg <- ("No allele assignment corrections made.") if (all(v %in% homo) & length(hetero) > 0) { keep_opt <- eval(formals(reassignAlleles)$keep_gene) @@ -1436,7 +1670,11 @@ reassignAlleles <- function(data, genotype_db, v_call="v_call", warning(msg) } - data$v_call_genotyped <- v_call_genotyped + if (overwrite) { + data[[v_call]] <- v_call_genotyped + } else { + data[[output_col]] <- v_call_genotyped + } return(data) } diff --git a/R/tigger.R b/R/tigger.R index 078e93f..c481f5f 100644 --- a/R/tigger.R +++ b/R/tigger.R @@ -62,7 +62,7 @@ #' } #' #' @import ggplot2 -#' @importFrom alakazam getAllele getGene getFamily translateDNA DNA_COLORS checkColumns +#' @importFrom alakazam getAllele getGene getFamily getLocus translateDNA DNA_COLORS checkColumns #' @importFrom doParallel registerDoParallel #' @importFrom dplyr do n desc %>% #' glimpse distinct group_indices @@ -74,7 +74,7 @@ #' summarise transmute #' @importFrom foreach foreach %dopar% registerDoSEQ #' @importFrom graphics plot -#' @importFrom gridExtra arrangeGrob +#' @importFrom gridExtra arrangeGrob grid.arrange #' @importFrom gtools ddirichlet #' @importFrom iterators icount #' @importFrom lazyeval interp @@ -85,7 +85,7 @@ #' @importFrom stringi stri_length stri_detect_fixed stri_replace_all_regex #' stri_sub stri_sub<- stri_trans_toupper #' @importFrom tidyr gather spread unnest -#' @importFrom utils citation +#' @importFrom utils citation head NULL # Package loading actions diff --git a/tests/testthat/test_coreFunctions.R b/tests/testthat/test_coreFunctions.R index a69c4cc..0930131 100644 --- a/tests/testthat/test_coreFunctions.R +++ b/tests/testthat/test_coreFunctions.R @@ -65,6 +65,15 @@ test_that("Test findNovelAlleles",{ v_call="v_call", seq="sequence_alignment") expect_equivalent(geno_bay, geno_bay_airr) + expect_false("genotyped_alleles" %in% colnames(geno_bay_airr)) + + geno_bay_gt <- inferGenotypeBayesian(airr_db, + germline_db = germline_ighv, + novel = novel_df_airr, + v_call="v_call", + seq="sequence_alignment", + genotyped_alleles=TRUE) + expect_true("genotyped_alleles" %in% colnames(geno_bay_gt)) }) @@ -200,8 +209,178 @@ test_that("Test genotypeFasta",{ ) gtfa <- genotypeFasta(gt, germline_db) expect_equal(gtfa, germline_db[1:5]) + + gt$genotyped_alleles <- c("04", "01", "01", "09") + gtfa_unseen <- genotypeFasta(gt, germline_db, include_unseen=TRUE) + expect_equal(gtfa_unseen, germline_db[c(6, 1, 3, 4, 5)]) expect_error(genotypeFasta(gt, germline_db[-1]), regexp="IGHV1-2\\*04") -}) \ No newline at end of file +}) + +test_that("inferGenotype uses locus-specific fractional gene_cutoff", { + db <- data.frame( + v_call=c(rep("IGHV1-1*01", 10), rep("IGKV1-1*01", 10)), + sequence_alignment=rep("AAAA", 20), + locus=c(rep("IGH", 10), rep("IGK", 10)), + stringsAsFactors=FALSE + ) + + expect_warning( + geno <- inferGenotype(db, find_unmutated=FALSE, gene_cutoff=0.75), + regexp="Mixed loci detected" + ) + expect_equal(sort(geno$gene), c("IGHV1-1", "IGKV1-1")) + + db$locus <- NULL + expect_warning( + geno_from_call <- inferGenotype(db, find_unmutated=FALSE, gene_cutoff=0.75), + regexp="Mixed loci detected" + ) + expect_equal(sort(geno_from_call$gene), c("IGHV1-1", "IGKV1-1")) +}) + +test_that("reassignAlleles supports segment output, overwrite, and trimming", { + db_v <- data.frame( + v_call="IGHV1-1*01", + sequence_alignment="AAAT", + stringsAsFactors=FALSE + ) + genotype_v <- c("IGHV1-1*01"="AAAA", "IGHV1-1*02"="AAAT") + reassigned_v <- reassignAlleles(db_v, genotype_v, + v_call="v_call", + seq="sequence_alignment") + expect_equal(reassigned_v$v_call_genotyped, "IGHV1-1*02") + + overwritten_v <- reassignAlleles(db_v, genotype_v, + v_call="v_call", + seq="sequence_alignment", + overwrite=TRUE) + expect_equal(overwritten_v$v_call, "IGHV1-1*02") + expect_false("v_call_genotyped" %in% colnames(overwritten_v)) + + db_d <- data.frame( + d_call="IGHD1-1*01", + sequence_alignment="CCCA", + stringsAsFactors=FALSE + ) + genotype_d <- c("IGHD1-1*01"="CCCC", "IGHD1-1*02"="CCCA") + reassigned_d <- reassignAlleles(db_d, genotype_d, + v_call="d_call", + seq="sequence_alignment") + expect_equal(reassigned_d$d_call_genotyped, "IGHD1-1*02") + + db_j <- data.frame( + j_call="IGHJ1*01", + sequence_alignment="TTTA", + stringsAsFactors=FALSE + ) + genotype_j <- c("IGHJ1*01"="TTTT", "IGHJ1*02"="TTTA") + reassigned_j <- reassignAlleles(db_j, genotype_j, + v_call="j_call", + seq="sequence_alignment") + expect_equal(reassigned_j$j_call_genotyped, "IGHJ1*02") + + db_trim <- data.frame( + v_call="IGHV1-1*01", + sequence_alignment="GGAAAT", + v_germline_start=3, + v_germline_end=6, + stringsAsFactors=FALSE + ) + genotype_trim <- c("IGHV1-1*01"="GGAAAA", "IGHV1-1*02"="CCAAAT") + reassigned_trim <- reassignAlleles(db_trim, genotype_trim, + v_call="v_call", + seq="sequence_alignment", + trim_seq=TRUE) + expect_equal(reassigned_trim$v_call_genotyped, "IGHV1-1*02") + + expect_error( + reassignAlleles(db_v, genotype_v, trim_seq=TRUE), + regexp="missing columns" + ) +}) + +test_that("reassignAlleles Rcpp mismatch path matches fallback results", { + skip_if_not( + exists("seqMismatchCountRcpp", envir=asNamespace("alakazam"), inherits=FALSE) && + exists("seqMismatchMatrixRcpp", envir=asNamespace("alakazam"), inherits=FALSE), + "Alakazam Rcpp mismatch functions are not installed" + ) + + samples <- c("ACGT", "ACNT", "AC-T", "acgt", "AC.T") + germlines <- c(g1="ACGA", g2="ACGT", g3="TCGT") + fallback <- sapply(germlines, function(x) { + sapply(getMutatedPositions(samples, x, ignored_regex="[\\.N-]", + match_instead=FALSE), length) + }) + rcpp <- get("seqMismatchMatrixRcpp", envir=asNamespace("alakazam"))( + samples, germlines, ignore=c(".", "N", "-")) + expect_equal(unname(rcpp), unname(fallback)) + + data_trim <- data.frame( + sequence_alignment=c("GGAAAT", "TTACGT"), + v_germline_start=c(3, 3), + v_germline_end=c(6, 6), + stringsAsFactors=FALSE + ) + samples_trim <- substr(data_trim$sequence_alignment, + data_trim$v_germline_start, + data_trim$v_germline_end) + germlines_trim <- c(g1="GGAAAA", g2="TTACGT", g3="CCAAAT") + fallback_trim <- sapply(germlines_trim, function(x) { + ref <- substr(rep(x, nrow(data_trim)), data_trim$v_germline_start, + data_trim$v_germline_end) + sapply(getMutatedPositions(samples_trim, ref, + ignored_regex="[\\.N-]", + match_instead=FALSE), length) + }) + rcpp_trim <- sapply(germlines_trim, function(x) { + ref <- substr(rep(x, nrow(data_trim)), data_trim$v_germline_start, + data_trim$v_germline_end) + get("seqMismatchCountRcpp", envir=asNamespace("alakazam"))( + samples_trim, ref, ignore=c(".", "N", "-")) + }) + expect_equal(unname(rcpp_trim), unname(fallback_trim)) + + db <- data.frame( + v_call=c("IGHV1-1*01", "IGHV1-1*01", "IGHV1-1*01"), + sequence_alignment=c("AAAT", "AAGT", "AAAC"), + stringsAsFactors=FALSE + ) + genotype_db <- c("IGHV1-1*01"="AAAA", "IGHV1-1*02"="AAAT") + old_opt <- getOption("tigger.use_alakazam_rcpp_mismatch") + on.exit(options(tigger.use_alakazam_rcpp_mismatch=old_opt), add=TRUE) + + options(tigger.use_alakazam_rcpp_mismatch=FALSE) + fallback_db <- reassignAlleles(db, genotype_db, + v_call="v_call", + seq="sequence_alignment") + options(tigger.use_alakazam_rcpp_mismatch=TRUE) + rcpp_db <- reassignAlleles(db, genotype_db, + v_call="v_call", + seq="sequence_alignment") + expect_equal(rcpp_db$v_call_genotyped, fallback_db$v_call_genotyped) + + db_trim <- data.frame( + v_call=c("IGHV1-1*01", "IGHV1-1*01"), + sequence_alignment=c("GGAAAT", "TTAAAA"), + v_germline_start=c(3, 3), + v_germline_end=c(6, 6), + stringsAsFactors=FALSE + ) + genotype_trim <- c("IGHV1-1*01"="GGAAAA", "IGHV1-1*02"="CCAAAT") + options(tigger.use_alakazam_rcpp_mismatch=FALSE) + fallback_trim_db <- reassignAlleles(db_trim, genotype_trim, + v_call="v_call", + seq="sequence_alignment", + trim_seq=TRUE) + options(tigger.use_alakazam_rcpp_mismatch=TRUE) + rcpp_trim_db <- reassignAlleles(db_trim, genotype_trim, + v_call="v_call", + seq="sequence_alignment", + trim_seq=TRUE) + expect_equal(rcpp_trim_db$v_call_genotyped, + fallback_trim_db$v_call_genotyped) +}) From 3ab8d17358c6cc730e6c49758ef3cc3890f6c195 Mon Sep 17 00:00:00 2001 From: ayelet peres Date: Sat, 20 Jun 2026 13:01:51 -0400 Subject: [PATCH 3/5] Moved the confidence plot based on the bayesian inference to a new plot function. Added a variable to choose the allele column to plot for the genotyped_alleles column in the bayesian --- NAMESPACE | 1 + R/functions.R | 115 +++++++++++++++++++++++++++++++++----------------- 2 files changed, 77 insertions(+), 39 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index c999c5a..36c1f2a 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -12,6 +12,7 @@ export(inferGenotype) export(inferGenotypeBayesian) export(insertPolymorphisms) export(plotGenotype) +export(plotGenotypeConfidence) export(plotNovel) export(readIgFasta) export(reassignAlleles) diff --git a/R/functions.R b/R/functions.R index 5737520..79e01bd 100644 --- a/R/functions.R +++ b/R/functions.R @@ -1178,6 +1178,10 @@ inferGenotype <- function(data, germline_db=NA, novel=NA, v_call="v_call", #' #' @param genotype \code{data.frame} of alleles denoting a genotype, #' as returned by \link{inferGenotype}. +#' @param allele_col name of the column in \code{genotype} holding the +#' comma-separated alleles to plot. Defaults to +#' \code{"alleles"}; set to \code{"genotyped_alleles"} to plot +#' the most likely alleles from \link{inferGenotypeBayesian}. #' @param facet_by column name in \code{genotype} to facet the plot by. #' if \code{NULL}, then do not facet the plot. #' @param gene_sort string defining the method to use when sorting alleles. @@ -1187,23 +1191,11 @@ inferGenotype <- function(data, germline_db=NA, novel=NA, v_call="v_call", #' @param text_size point size of the plotted text. #' @param silent if \code{TRUE} do not draw the plot and just return the ggplot #' object; if \code{FALSE} draw the plot. -#' @param confidence_col name of a column in \code{genotype} holding a per-gene -#' confidence value to display as evidence, e.g. \code{"k_diff"} -#' from \link{inferGenotypeBayesian}. If \code{NULL} (default), -#' no confidence panel is drawn. The value is binned and shown on a -#' blue color scale, with the column name as the legend title. -#' @param confidence_breaks numeric vector of breaks used to bin the confidence -#' value into color groups. The default -#' \code{c(0, 1, 2, 3, 4, 5, 10, 20, 50, Inf)} matches the binning -#' used by RAbHIT for the haplotype lK panel. #' @param ... additional arguments to pass to ggplot2::theme. #' -#' @return A ggplot object defining the plot. If \code{confidence_col} is supplied, a -#' \code{gridExtra} grob is returned instead, placing the confidence panel -#' next to the genotype. \code{confidence_col} is not combined with -#' \code{facet_by}; if both are given, \code{facet_by} is ignored. +#' @return A ggplot object defining the plot. #' -#' @seealso \link{inferGenotype}, \link{inferGenotypeBayesian} +#' @seealso \link{inferGenotype}, \link{plotGenotypeConfidence} #' #' @examples #' # Plot genotype @@ -1216,46 +1208,34 @@ inferGenotype <- function(data, germline_db=NA, novel=NA, v_call="v_call", #' geno_sub <- rbind(genotype_a, genotype_b) #' plotGenotype(geno_sub, facet_by="SUBJECT", gene_sort="pos") #' -#' # Add a confidence evidence panel from a per-gene confidence column -#' geno_conf <- SampleGenotype -#' geno_conf$k_diff <- seq(0, 20, length.out=nrow(geno_conf)) -#' plotGenotype(geno_conf, confidence_col="k_diff") -#' #' @export plotGenotype <- function(genotype, facet_by=NULL, gene_sort=c("name", "position"), - text_size=12, silent=FALSE, confidence_col=NULL, - confidence_breaks=c(0, 1, 2, 3, 4, 5, 10, 20, 50, Inf), ...) { + text_size=12, silent=FALSE, allele_col="alleles", ...) { # Check arguments gene_sort <- match.arg(gene_sort) - if (!is.null(confidence_col)) { - if (!confidence_col %in% colnames(genotype)) { - stop("confidence_col '", confidence_col, "' not found in genotype.") - } - if (!is.null(facet_by)) { - warning("confidence_col is not supported together with facet_by; ignoring facet_by.") - facet_by = NULL - } + if (!allele_col %in% colnames(genotype)) { + stop("allele_col '", allele_col, "' not found in genotype.") } # Split genes' alleles into their own rows - alleles = strsplit(genotype$alleles, ",") + alleles = strsplit(genotype[[allele_col]], ",") geno2 = genotype r = 1 for (g in 1:nrow(genotype)){ for(a in 1:length(alleles[[g]])) { geno2[r, ] = genotype[g, ] - geno2[r, ]$alleles = alleles[[g]][a] + geno2[[allele_col]][r] = alleles[[g]][a] r = r + 1 } } # Set the gene order - gene_levels = rev(sortAlleles(unique(geno2$gene), method=gene_sort)) - geno2$gene = factor(geno2$gene, levels=gene_levels) + geno2$gene = factor(geno2$gene, + levels=rev(sortAlleles(unique(geno2$gene), method=gene_sort))) # Create the base plot p = ggplot(geno2, aes(x=!!rlang::sym("gene"), - fill=!!rlang::sym("alleles"))) + + fill=!!rlang::sym(allele_col))) + theme_bw() + theme(axis.ticks=element_blank(), axis.text.x=element_blank(), @@ -1276,14 +1256,71 @@ plotGenotype <- function(genotype, facet_by=NULL, gene_sort=c("name", "position" # Add additional theme elements p = p + do.call(theme, list(...)) - # Without a confidence column, return the genotype plot as is - if (is.null(confidence_col)) { - if (!silent) { plot(p) } - return(invisible(p)) + # Plot + if (!silent) { plot(p) } + + invisible(p) +} + + +#' Show a genotype with a confidence evidence panel +#' +#' \code{plotGenotypeConfidence} draws a genotype with \link{plotGenotype} and adds a +#' color panel beside it showing a per-gene confidence value, such as the \code{k_diff} +#' produced by \link{inferGenotypeBayesian}. +#' +#' @param genotype \code{data.frame} of alleles denoting a genotype, as +#' returned by \link{inferGenotypeBayesian}. +#' @param confidence_col name of the column in \code{genotype} holding the per-gene +#' confidence value, e.g. \code{"k_diff"}. The value is binned and +#' shown on a blue color scale, with the column name as the legend +#' title; white marks genes with no value. +#' @param allele_col name of the column in \code{genotype} holding the alleles to +#' plot, passed to \link{plotGenotype}. Set to +#' \code{"genotyped_alleles"} to plot the most likely alleles. +#' @param gene_sort string defining the method to use when sorting alleles, passed +#' to \link{plotGenotype}. +#' @param text_size point size of the plotted text. +#' @param confidence_breaks numeric vector of breaks used to bin the confidence +#' value into color groups. The default +#' \code{c(0, 1, 2, 3, 4, 5, 10, 20, 50, Inf)} matches the binning +#' used by RAbHIT for the haplotype lK panel. +#' @param silent if \code{TRUE} do not draw the plot and just return the grob; +#' if \code{FALSE} draw the plot. +#' @param ... additional arguments to pass to ggplot2::theme of the +#' genotype panel. +#' +#' @return A \code{gridExtra} grob combining the genotype plot with the confidence panel. +#' +#' @seealso \link{plotGenotype}, \link{inferGenotypeBayesian} +#' +#' @examples +#' # The Bayesian genotype carries a per-gene confidence (k_diff) and, optionally, +#' # a genotyped_alleles column +#' geno_bayesian <- inferGenotypeBayesian(AIRRDb, germline_db=SampleGermlineIGHV, +#' novel=SampleNovel, genotyped_alleles=TRUE) +#' plotGenotypeConfidence(geno_bayesian, confidence_col="k_diff", +#' allele_col="genotyped_alleles") +#' +#' @export +plotGenotypeConfidence <- function(genotype, confidence_col, allele_col="alleles", + gene_sort=c("name", "position"), text_size=12, + confidence_breaks=c(0, 1, 2, 3, 4, 5, 10, 20, 50, Inf), + silent=FALSE, ...) { + # Check arguments + gene_sort <- match.arg(gene_sort) + if (!confidence_col %in% colnames(genotype)) { + stop("confidence_col '", confidence_col, "' not found in genotype.") } + # Base genotype plot + p = plotGenotype(genotype, gene_sort=gene_sort, text_size=text_size, + allele_col=allele_col, silent=TRUE, ...) + # Bin the per-gene confidence value and draw it as a blue color panel beside the - # genotype. White marks NA/unscored genes; drop=FALSE keeps the full scale. + # genotype. White marks NA/unscored genes; drop=FALSE keeps the full scale. The + # gene order matches plotGenotype so the rows align. + gene_levels = rev(sortAlleles(unique(genotype$gene), method=gene_sort)) blues = c("#FFFFFF", "#F7FBFF", "#DEEBF7", "#C6DBEF", "#9ECAE1", "#6BAED6", "#4292C6", "#2171B5", "#08519C", "#08306B") bins = cut(suppressWarnings(as.numeric(genotype[[confidence_col]])), From 6695385eff07df73da1319576a5faaf71400f08b Mon Sep 17 00:00:00 2001 From: ayelet peres Date: Mon, 22 Jun 2026 09:35:37 -0400 Subject: [PATCH 4/5] Add strip_d parameter to genotypeFasta and reassignAlleles functions for gene annotation handling --- R/functions.R | 36 +++++++++++++++++++++++++++--------- 1 file changed, 27 insertions(+), 9 deletions(-) diff --git a/R/functions.R b/R/functions.R index 79e01bd..7890d9e 100644 --- a/R/functions.R +++ b/R/functions.R @@ -1379,6 +1379,10 @@ plotGenotypeConfidence <- function(genotype, confidence_col, allele_col="alleles #' genes that are not present in \code{genotype}. For #' genes present in \code{genotype}, include only the #' genotyped alleles. +#' @param strip_d if \code{TRUE} (default) remove the "D" from the end of +#' gene annotations (denoting a duplicate gene in the locus) +#' when matching genotype alleles to \code{germline_db}. If +#' \code{FALSE}, alleles are matched exactly. #' #' @return A named vector of strings containing the germline nucleotide #' sequences of the alleles in the provided genotype. @@ -1390,7 +1394,8 @@ plotGenotypeConfidence <- function(genotype, confidence_col, allele_col="alleles #' genotype_db <- genotypeFasta(SampleGenotype, SampleGermlineIGHV, SampleNovel) #' #' @export -genotypeFasta <- function(genotype, germline_db, novel=NA, include_unseen=FALSE){ +genotypeFasta <- function(genotype, germline_db, novel=NA, include_unseen=FALSE, + strip_d=TRUE){ if(!is.null(nrow(novel))){ # Extract novel alleles if any and add them to germline_db novel <- filter(novel, !is.na(!!rlang::sym("polymorphism_call"))) %>% @@ -1402,9 +1407,9 @@ genotypeFasta <- function(genotype, germline_db, novel=NA, include_unseen=FALSE) } } - genotype$gene <- getGene(genotype$gene, first = T, strip_d = T) + genotype$gene <- getGene(genotype$gene, first = T, strip_d = strip_d) g_names <- names(germline_db) - names(g_names) <- getAllele(names(germline_db), first = T, strip_d = T) + names(g_names) <- getAllele(names(germline_db), first = T, strip_d = strip_d) allele_values <- genotype$alleles if ("genotyped_alleles" %in% colnames(genotype)) { @@ -1425,7 +1430,7 @@ genotypeFasta <- function(genotype, germline_db, novel=NA, include_unseen=FALSE) } if (include_unseen) { - germline_genes <- getGene(names(germline_db), first=TRUE, strip_d=TRUE) + germline_genes <- getGene(names(germline_db), first=TRUE, strip_d=strip_d) seqs <- c(germline_db[!germline_genes %in% genotype$gene], seqs) } @@ -1488,6 +1493,15 @@ genotypeFasta <- function(genotype, germline_db, novel=NA, include_unseen=FALSE) #' are equally close. \code{"alphabetical"} keeps the #' first \code{top_k} by name; \code{"mutation_count"} #' keeps all ties. +#' @param strip_d if \code{TRUE} (default) remove the "D" from the end of +#' gene annotations (denoting a duplicate gene in the locus) +#' when grouping the \code{genotype_db} alleles by gene. If +#' \code{FALSE}, the "D" is kept, so the genotype grouping +#' and the sequence calls are matched consistently (use +#' together with \code{genotypeFasta(strip_d=FALSE)}). +#' @param reassign_uncalled if \code{TRUE} (default), sequences whose call is +#' empty or \code{NA} are also realigned against the whole +#' genotype. If \code{FALSE}, they are left unassigned. #' #' @return A modified input \code{data.frame} containing the best allele call from #' among the sequences listed in \code{genotype_db} in the @@ -1510,7 +1524,8 @@ reassignAlleles <- function(data, genotype_db, v_call="v_call", trim_seq=FALSE, overwrite=FALSE, ignored_regex="[\\.N-]", treat_multigene_as_uncalled=FALSE, - top_k=NULL, top_by=c("alphabetical", "mutation_count")){ + top_k=NULL, top_by=c("alphabetical", "mutation_count"), + strip_d=TRUE, reassign_uncalled=TRUE){ # Check arguments keep_gene <- match.arg(keep_gene) top_by <- match.arg(top_by) @@ -1604,11 +1619,11 @@ reassignAlleles <- function(data, genotype_db, v_call="v_call", if (keep_gene == "gene") { v = getGene(uniq_alleles, first = TRUE, strip_d=FALSE)[uniq_idx] - geno = getGene(names(genotype_db),strip_d=TRUE) + geno = getGene(names(genotype_db),strip_d=strip_d) names(geno) = names(genotype_db) } else if (keep_gene == "family") { v <- getFamily(uniq_alleles, first = TRUE, strip_d = FALSE)[uniq_idx] - geno = getFamily(names(genotype_db),strip_d=TRUE) + geno = getFamily(names(genotype_db),strip_d=strip_d) names(geno) = names(genotype_db) } else if (keep_gene == "repertoire") { v <- rep(v_call, length(v_calls)) @@ -1674,9 +1689,12 @@ reassignAlleles <- function(data, genotype_db, v_call="v_call", } } - # Now realign the gene-not-in-genotype calls to every genotype allele + # Now realign the gene-not-in-genotype calls to every genotype allele. + # By default only sequences that carry a call are realigned; set + # reassign_uncalled=TRUE to also realign empty/NA-call sequences (legacy). hetero_calls_i = which(v %in% hetero & !is_multigene & has_call) - not_called = setdiff(which(has_call), c(homo_calls_i, hetero_calls_i)) + candidate_rows = if (reassign_uncalled) seq_along(v) else which(has_call) + not_called = setdiff(candidate_rows, c(homo_calls_i, hetero_calls_i)) if(length(not_called)>0){ if(method == "hamming"){ dist_mat <- mismatch_matrix(v_sequences[not_called], genotype_db, not_called) From 48df6553cd681fcdaecc340c1ac85bb093d63c41 Mon Sep 17 00:00:00 2001 From: ayelet peres Date: Mon, 22 Jun 2026 17:48:27 -0400 Subject: [PATCH 5/5] Refactor plotGenotypeConfidence function to improve confidence visualization with color bins and tile geometry --- R/functions.R | 26 ++++++++++++++++---------- 1 file changed, 16 insertions(+), 10 deletions(-) diff --git a/R/functions.R b/R/functions.R index 7890d9e..67b6966 100644 --- a/R/functions.R +++ b/R/functions.R @@ -1317,20 +1317,21 @@ plotGenotypeConfidence <- function(genotype, confidence_col, allele_col="alleles p = plotGenotype(genotype, gene_sort=gene_sort, text_size=text_size, allele_col=allele_col, silent=TRUE, ...) - # Bin the per-gene confidence value and draw it as a blue color panel beside the - # genotype. White marks NA/unscored genes; drop=FALSE keeps the full scale. The - # gene order matches plotGenotype so the rows align. + # Confidence plot gene_levels = rev(sortAlleles(unique(genotype$gene), method=gene_sort)) blues = c("#FFFFFF", "#F7FBFF", "#DEEBF7", "#C6DBEF", "#9ECAE1", "#6BAED6", "#4292C6", "#2171B5", "#08519C", "#08306B") bins = cut(suppressWarnings(as.numeric(genotype[[confidence_col]])), confidence_breaks, include.lowest=TRUE, right=FALSE) bin_levels = gsub(",", ", ", levels(bins)) + bins_chr = gsub(",", ", ", as.character(bins), fixed=TRUE) conf_levels = c("NA", bin_levels) + conf_colors = setNames(blues[seq_along(conf_levels)], conf_levels) conf = data.frame(gene=factor(genotype$gene, levels=gene_levels), - confidence=factor(ifelse(is.na(bins), "NA", bin_levels[bins]), - levels=conf_levels)) - pc = ggplot(conf, aes(x=!!rlang::sym("gene"), + confidence=factor(ifelse(is.na(bins_chr), "NA", bins_chr), + levels=conf_levels)) + + pc = ggplot(conf, aes(x = 1, y=!!rlang::sym("gene"), fill=!!rlang::sym("confidence"))) + theme_bw() + theme(axis.ticks=element_blank(), @@ -1338,10 +1339,15 @@ plotGenotypeConfidence <- function(genotype, confidence_col, allele_col="alleles panel.grid.major=element_blank(), panel.grid.minor=element_blank(), text=element_text(size=text_size)) + - geom_bar(position="fill") + - coord_flip() + xlab("") + ylab("") + - scale_fill_manual(name=confidence_col, - values=setNames(blues[seq_along(conf_levels)], conf_levels), + geom_tile(width=1, height=0.9, show.legend=TRUE) + + scale_x_continuous(expand=c(0, 0)) + + scale_y_discrete(drop=FALSE) + + xlab("") + ylab("") + + scale_fill_manual(name=confidence_col, values=conf_colors, + limits=conf_levels, breaks=conf_levels, + guide=guide_legend(override.aes=list( + fill=unname(conf_colors), + colour=NA,alpha=1)), drop=FALSE) # Align the gene rows and place the confidence panel beside the genotype, with