Immcantation - enchantR
Find Threshold
Updated: Fri Jan 16 19:05:38 2026
1 Input
1.2 Data Sources
Code
# Read repertoires,
# keep only needed columns
selected_columns <- unique(c("sample_id", "sequence_id",
params$cloneby,
params$crossby,
"v_call", "j_call", "junction",
"cell_id", params$singlecell, "locus"))
db <- readInput(params[['input']], col_select=selected_columns)
input_sizes <- db %>%
count(input_file)
input_sizes <- input_sizes %>% rename("sequences" = "n")
eetable(input_sizes)$tableCode
input_size <- nrow(db)
is_heavy <- isHeavyChain(db[['locus']])
singlecell <- params$singlecell
# Check for single cell label
if (!is.null(singlecell)) {
if (singlecell %in% colnames(db)) {
db[[ singlecell ]] <- as.logical(db[[ singlecell ]] )
} else {
stop("`", singlecell, "` is not a valid field in `db`.")
}
} else {
singlecell <- "single_cell"
db[[singlecell]] <- F
if ("cell_id" %in% colnames(db)) {
message("Setting `singlecell` using `cell_id`.")
db[[singlecell]][!db[['cell_id']] %in% c(NA, '')] <- T
}
}
na_single <- is.na(db[[singlecell]])
if (sum(na_single)>0) {
warning(sum(na_single), " sequences are missing single cell information. Using single_cell=F")
db[[singlecell]][na_single] <- FALSE
}
bulk_heavy <- is_heavy & !db[[singlecell]]
sc_heavy <- is_heavy & db[[singlecell]]
bulk_dtn <- data.frame()
sc_dtn <- data.frame()
subsample <- params[['subsample']]
if (subsample %in% c(NA, 0, NULL)) {
subsample <- NULL
}Number of sequences loaded: 12307.
Code
2 Calculate Threshold for Clonal Relationships
The shazam package provides methods to find an appropriate distance threshold to determine clonal relationships for each dataset (distToNearest and findThreshold) by generating a distribution of distances between each sequence and its closest non-identical neighbor. Typically, the distance-to-nearest distribution for a repertoire is bimodal. The first mode (on the left) represents sequences that have at least one clonal relative in the dataset, while the second mode (on the right) is representative of the sequences that do not have any clonal relatives in the data. A reasonable threshold will separate these two modes of the distribution. The distance may be selected manually, or findThreshold can be used to automatically select a threshold.
Specifying the cross argument to distToNearest forces distance calculations to be performed across groups, such that the nearest neighbor of each sequence will always be a sequence in a different group.
Code
if (any(db[[singlecell]] == F) && sum(bulk_heavy)>0) {
message("Using ",sum(bulk_heavy), " bulk heavy chain sequences.")
# Add distances to nearest neighbor within a "cloneby" group
bulk_dtn <- distToNearest(db[bulk_heavy,] ,
fields=params$cloneby,
sequenceColumn="junction",
vCallColumn="v_call", jCallColumn="j_call",
model="ham", first=FALSE, VJthenLen=TRUE, normalize="len",
nproc=params$nproc,
subsample = subsample)
bulk_dtn <- distToNearest(bulk_dtn,
fields=NULL,
cross=params$crossby,
sequenceColumn="junction",
vCallColumn="v_call", jCallColumn="j_call",
model="ham", first=FALSE, VJthenLen=TRUE, normalize="len",
nproc=params$nproc,
subsample = subsample)
} Code
if (any(db[[singlecell]] == T) && sum(sc_heavy)>0) {
message("Using ",sum(sc_heavy), " sc heavy chain sequences.")
# Create unique cell_id
# db[['cell_id']] <- paste0(db[['sample_id']],db[['cell_id']])
# Add distances to nearest neighbor within a "cloneby" group
sc_dtn <- distToNearest(db[sc_heavy,] ,
fields=params$cloneby,
sequenceColumn="junction",
vCallColumn="v_call", jCallColumn="j_call",
model="ham", first=FALSE, VJthenLen=TRUE, normalize="len",
nproc=params$nproc,
cellIdColumn = "cell_id",
locusColumn = "locus",
onlyHeavy = TRUE,
subsample = subsample)
# Add distances to nearest neighbor across a "crossby" group
sc_dtn <- distToNearest(sc_dtn,
fields=NULL,
cross=params$crossby,
sequenceColumn="junction",
vCallColumn="v_call", jCallColumn="j_call",
model="ham", first=FALSE, VJthenLen=TRUE, normalize="len",
nproc=params$nproc,
cellIdColumn = "cell_id",
locusColumn = "locus",
onlyHeavy = TRUE,
subsample = subsample)
} 2.1 distToNearest plot
Code
caption <- paste0("Distribution of the distance to the nearest sequence within ",
paste(params$cloneby, collapse=","),", displayed by sample_id. When ",
"cross ",paste(params$crossby, collapse=","),
" distances are available, they are shown in an inverse y-axis.")
if (!is.null(subsample)) {
caption <- paste0(caption," Subsampling requested: ", subsample, ".")
}
wrap_formula <- as.formula(paste0(
c(
singlecell,
"~",
paste(unique(c("sample_id", params$clone_by)), collapse = "+")
),collapse="")
)
db <- bind_rows(bulk_dtn,sc_dtn)
dtnplot <- NULL
if (any(!is.na(db[['dist_nearest']]))) {
dtnplot <- ggplot(db) +
geom_histogram(aes(dist_nearest, fill=sample_id),binwidth = 0.01) +
scale_x_continuous(
breaks = seq(0, 1, 0.1)
)
} else {
cat("All `dist_nearest` values are NA.")
warning("All `dist_nearest` values are NA.")
}
if ("cross_dist_nearest" %in% c(colnames(bulk_dtn),colnames(sc_dtn))) {
if (any(!is.na(db[['cross_dist_nearest']]))) {
dtnplot <- dtnplot +
geom_histogram(
aes(x=cross_dist_nearest, y=-(..count..), fill=sample_id),
binwidth=0.01, position="identity") +
#scale_y_continuous(labels = abs) +
ylab("count")
} else {
cat("All `cross_dist_nearest` values are NA.")
warning("All `cross_dist_nearest` values are NA.")
}
}
if (is.null(dtnplot)) {
cat("All distance values are NA.")
}Code
dtnplot <- dtnplot +
facet_wrap(wrap_formula, scales = "free_y", ncol=2, labeller = "label_both") +
expand_limits(x = 0, y =0) +
theme_enchantr() +
theme(legend.position = "bottom")
dtnplot <- eeplot(dtnplot,
outdir=params$outdir,
file=knitr::opts_current$get('dtnplot'),
caption=caption
)
ggplotly(dtnplot + theme(panel.spacing=unit(2, 'lines'), legend.position="none"))Figure 2.1: Distribution of the distance to the nearest sequence within subject_id, displayed by sample_id. When cross subject_id distances are available, they are shown in an inverse y-axis. ggplot file: dtnplot.RData
2.2 Find threshold
findThreshold uses the distribution of distances calculated in the previous step to determine an appropriate threshold for the dataset.
2.2.1 Threshold(s) summary table
Code
Code
cross_dist_nearest <- NULL
if ("cross_dist_nearest" %in% colnames(db)) {
if (any(!is.na(db[['cross_dist_nearest']]))) {
cross_dist_nearest <- "cross_dist_nearest"
} else {
warning("All `cross_dist_nearest` values are NA.")
}
}
threshold <- findThresholdDb(db,
distanceColumn="dist_nearest",
crossDistanceColumn=cross_dist_nearest,
method=params$findthreshold_method,
model=params$findthreshold_model,
edge=params$findthreshold_edge,
cutoff = params$findthreshold_cutoff,
spc = params$findthreshold_spc,
nproc=params$nproc,
fields=params$cloneby,
subsample=subsample)
threshold_summary <- gmmSummary(threshold) %>%
mutate(mean_threshold=mean(threshold, na.rm=TRUE))
mean_threshold <- round(threshold_summary$mean_threshold[1],2)
tab_caption <- paste0("Summary of threshold values, p-value from Hartigans’ dip statistic (HDS) test with <0.05 indicating significant bimodality. Mean threshold is: ", mean_threshold)
if (!is.null(subsample)) {
tab_caption <- paste0(tab_caption," Subsampling requested: ", subsample, ".")
}
tab <- eetable(threshold_summary,
outdir=params$outdir,
file=paste(params$outname, "threshold-summary", sep="_"),
caption=tab_caption)
tab$table %>%
DT::formatRound(columns=setdiff(colnames(threshold_summary)[-1],"model"),
digits=3)Code
caption_list <- NULL
tmp <- lapply(threshold, function(thr) {
this_thr <- NA
if (!is.null(attributes(thr$GmmThreshold))) {
this_thr <- round(attributes(thr$GmmThreshold)$threshold,2)
}
cat ("\n\n### ",paste(c(params$crossby, " ", thr[['fields']], ", estimated threshold:", this_thr),collapse=" "),"\n\n")
fn <- paste(c(thr[['fields']],"dtnthrplot"),collapse="_")
dtnthrcaption <- paste0("Distribution of the distance to the nearest sequence within ",
paste(params$cloneby, collapse=",")," (",
paste(thr[['fields']],collapse="_")
,"). When ",
"cross ",paste(params$crossby, collapse=","),
" distances are available, they are shown in an inverse y-axis.")
p <- thr[['plot']] +
scale_x_continuous(
breaks = seq(0, 1, 0.1)
)
p <- eeplot(p,
outdir=params$outdir,
file=fn,
caption=dtnthrcaption
)
print(p)
caption_list <<- c(caption_list, p$enchantr$html_caption)
})2.2.2 subject_id M5 , estimated threshold: 0.06
Figure 2.2: Distribution of the distance to the nearest sequence within subject_id (M5). When cross subject_id distances are available, they are shown in an inverse y-axis. ggplot file: M5_dtnthrplot.RData
2.2.3 subject_id M4 , estimated threshold: 0.14
Figure 2.3: Distribution of the distance to the nearest sequence within subject_id (M4). When cross subject_id distances are available, they are shown in an inverse y-axis. ggplot file: M4_dtnthrplot.RData
3 Save
Code
# pass <- db[['collapse_pass']]
if (!is.null(params$outname)) {
output_fn <- paste(params$outname,"threshold-pass.tsv", sep="_")
} else {
output_fn <- sub(".tsv$", "_threshold-pass.tsv", basename(params$input))
}
if (!is.null(params$log)) {
log_fn <- paste0(params$log,".txt")
} else {
log_fn <- sub("threshold-pass.tsv$", "command_log.txt", basename(output_fn))
}
# write_rearrangement(db , file=output_fn)
tables_dir <- file.path(params$outdir,"tables")
if (!dir.exists(tables_dir)) {
dir.create(tables_dir, recursive = T)
}
cat(threshold_summary[['mean_threshold']][1],
file=file.path(tables_dir,
sub("threshold-pass.tsv$", "threshold-mean.tsv" ,output_fn)),
append=F
)Code
cat("START> FindThreshold", file=log_fn, append=F)
cat(paste0("\nFILE> ",basename(params$input)), file=log_fn, append=T)
cat(paste0("\nOUTPUT> ",basename(output_fn)), file=log_fn, append=T)
cat(paste0("\nPASS> ",nrow(db)), file=log_fn, append=T)
cat(paste0("\nFAIL> ",input_size-nrow(db)), file=log_fn, append=T)4 Software versions
## R version 4.4.3 (2025-02-28)
## Platform: x86_64-redhat-linux-gnu
## Running under: Fedora Linux 40 (Container Image)
##
## Matrix products: default
## BLAS/LAPACK: FlexiBLAS OPENBLAS-OPENMP; LAPACK version 3.12.0
##
## locale:
## [1] C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] R.utils_2.13.0 R.oo_1.27.0 R.methodsS3_1.8.2 plotly_4.10.4
## [5] tidyr_1.3.1 shazam_1.2.0 alakazam_1.3.0 ggplot2_3.5.2
## [9] airr_1.5.0 dplyr_1.1.4 DT_0.33 enchantr_0.1.21
##
## loaded via a namespace (and not attached):
## [1] bitops_1.0-9 gridExtra_2.3
## [3] rlang_1.1.6 magrittr_2.0.3
## [5] clue_0.3-66 GetoptLong_1.0.5
## [7] ade4_1.7-23 matrixStats_1.5.0
## [9] compiler_4.4.3 png_0.1-8
## [11] vctrs_0.6.5 stringr_1.5.1
## [13] pkgconfig_2.0.3 shape_1.4.6.1
## [15] crayon_1.5.3 fastmap_1.2.0
## [17] XVector_0.46.0 labeling_0.4.3
## [19] ggraph_2.2.1 Rsamtools_2.22.0
## [21] rmarkdown_2.29 tzdb_0.5.0
## [23] UCSC.utils_1.2.0 bit_4.6.0
## [25] purrr_1.0.4 xfun_0.52
## [27] zlibbioc_1.52.0 cachem_1.1.0
## [29] seqinr_4.2-36 GenomeInfoDb_1.42.3
## [31] jsonlite_2.0.0 progress_1.2.3
## [33] DelayedArray_0.32.0 BiocParallel_1.40.1
## [35] tweenr_2.0.3 parallel_4.4.3
## [37] prettyunits_1.2.0 cluster_2.1.8
## [39] R6_2.6.1 bslib_0.9.0
## [41] stringi_1.8.7 RColorBrewer_1.1-3
## [43] diptest_0.77-1 jquerylib_0.1.4
## [45] GenomicRanges_1.58.0 bookdown_0.42
## [47] knitr_1.50 Rcpp_1.0.14
## [49] SummarizedExperiment_1.36.0 iterators_1.0.14
## [51] readr_2.1.5 IRanges_2.40.1
## [53] Matrix_1.7-3 igraph_2.1.4
## [55] tidyselect_1.2.1 rstudioapi_0.17.1
## [57] abind_1.4-8 yaml_2.3.10
## [59] viridis_0.6.5 doParallel_1.0.17
## [61] codetools_0.2-20 lattice_0.22-6
## [63] tibble_3.2.1 Biobase_2.66.0
## [65] withr_3.0.2 evaluate_1.0.3
## [67] polyclip_1.10-7 circlize_0.4.16
## [69] Biostrings_2.74.1 pillar_1.10.2
## [71] MatrixGenerics_1.18.1 KernSmooth_2.23-26
## [73] foreach_1.5.2 stats4_4.4.3
## [75] generics_0.1.3 vroom_1.6.5
## [77] S4Vectors_0.44.0 hms_1.1.3
## [79] munsell_0.5.1 scales_1.3.0
## [81] glue_1.8.0 lazyeval_0.2.2
## [83] tools_4.4.3 data.table_1.17.0
## [85] GenomicAlignments_1.42.0 graphlayouts_1.2.2
## [87] tidygraph_1.3.1 grid_4.4.3
## [89] ape_5.8-1 crosstalk_1.2.1
## [91] colorspace_2.1-1 nlme_3.1-167
## [93] GenomeInfoDbData_1.2.13 ggforce_0.4.2
## [95] cli_3.6.4 S4Arrays_1.6.0
## [97] viridisLite_0.4.2 ComplexHeatmap_2.18.0
## [99] gtable_0.3.6 sass_0.4.10
## [101] digest_0.6.37 BiocGenerics_0.52.0
## [103] SparseArray_1.6.2 ggrepel_0.9.6
## [105] rjson_0.2.23 htmlwidgets_1.6.4
## [107] farver_2.1.2 memoise_2.0.1
## [109] htmltools_0.5.8.1 lifecycle_1.0.4
## [111] httr_1.4.7 GlobalOptions_0.1.2
## [113] bit64_4.6.0-1 MASS_7.3-64