1 Input

1.1 Parameters

Code
printParams(params)
Table 1.1: Input parameters.

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)$table
Code
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
input_summary <- db %>%
   group_by(!!!rlang::syms(unique(c("sample_id",params$cloneby, "locus")))) %>%
   summarize(n=n(), .groups="drop") %>%
       pivot_wider(names_from=locus, values_from=n)

input_summary_table <- eetable(input_summary) 
input_summary_table$table
Table 1.2: Input data summary.

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

Code
caption <- dtnplot$enchantr$html_caption

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
# This chunk runs when all dist_nearest values are NA
threshold <- NA
threshold_summary <- data.frame(mean_threshold = NA)
cat("All `dist_nearest` values are NA. Skipping threshold analysis.")
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)
Table 2.1: Summary of threshold values, p-value from Hartigans<U+2019> dip statistic (HDS) test with <0.05 indicating significant bimodality. Mean threshold is: 0.1 File can be found here: all_reps_threshold-summary.tsv
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

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. <a href='ggplots/M5_dtnthrplot.RData'>ggplot file: M5_dtnthrplot.RData</a>

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

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. <a href='ggplots/M4_dtnthrplot.RData'>ggplot file: M4_dtnthrplot.RData</a>

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

Code
sessionInfo()
## 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