1 Input parameters

Code
printParams(params)
Table 1.1: Input parameters.

2 Read repertoires

Code
# Read repertoire
db <- readInput(params[['input']], col_select = NULL)

if (params$species =="auto") {
    if ("species" %in% colnames(db)) {
        species <- unique(db[['species']])
        if (length(species)>1) {
            stop("Multiple species detected. Only one allowed.")
        }
    } else  {
        stop("Can't detect species. The column `species` does not exist in `db`.")
    }
} else {
    species <- params$species
}

# 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
   }

}

if (!"locus" %in% colnames(db)) {
    db[['locus']] <- getLocus(db[['v_call']])
}
heavy_chains <- isHeavyChain(db[['locus']])

if ("clone_id" %in% colnames(db)) {
    if (params$force) {
        # Reset if force
        warning("Overwritting clone_id.")
        db$clone_id <- NULL
    }   
    # Reset always if clone_id exists
    if ("clone_size_count" %in% colnames(db)) {    
        warning("Overwritting clone_size_count.")    
        db$clone_size_count <- NULL
    }
    if ("clone_size_freq" %in% colnames(db)) {    
        warning("Overwritting clone_size_freq.")    
        db$clone_size_freq <- NULL
    }
}
Code
input_size <- nrow(db)
input_sizes <- db %>%
    group_by(!!!rlang::syms(unique(c("input_file")))) %>%
    summarize(input_size=n())

# Track input size of the expected output groups, for the
# end of report summary
input_sizes_byoutput <- db %>%
    group_by(!!!rlang::syms(unique(params$outputby))) %>%
    mutate(
        num_input_files = length(unique(input_file)),
        input_files = paste(unique(input_file), collapse=",")
    ) %>%
    ungroup() %>%
    group_by(!!!rlang::syms(unique(c(params$outputby, "input_file", "num_input_files", "input_files")))) %>%
    summarize(
        input_size = n()
    )

Number of sequences loaded: 12297. Number of heavy chain sequences loaded: 12297.

Code
input_samples_summary <- db %>%
   group_by(sample_id, subject_id, tissue) %>%
   summarize(size = n()) %>%
   arrange(subject_id)

eetable(input_samples_summary)$table
Table 2.1: Input samples summary.

2.1 Sequences per locus

Code
input_locus_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) %>%
   rowwise() %>%
   mutate(Total = sum(!!!rlang::syms(unique(db[['locus']]))))

total <- data.frame(list("sample_id"="Total", 
                    t(input_locus_summary %>%
   select(!!!rlang::syms(c(unique(db[['locus']]), "Total"))) %>%
   colSums(na.rm = T))))

input_locus_summary <- bind_rows(input_locus_summary, total)

tab_caption <- paste0("Input data. Number of sequences in each ", 
                   paste("sample", params$cloneby, sep=", "),
                   " and locus."
      )
eetable(input_locus_summary)$table
Table 2.2: Input data. Number of sequences in each sample, subject_id and locus.
Code
cat("## Sequences per c_call\n")

2.2 Sequences per c_call

Code
input_c_call_summary <- db %>%
   group_by(!!!rlang::syms(unique(c("sample_id",params$cloneby, "c_call")))) %>%
   summarize(n=n(), .groups="drop") %>%
       pivot_wider(names_from=c_call, values_from=n)

tab_caption <- paste0("Input data. Number of sequences in each ", 
                   paste("sample", params$cloneby, sep=", "),
                   " and c_call"
      )
eetable(input_c_call_summary)$table
Table 2.3: Input data. Number of sequences in each sample, subject_id and c_call
Code
cat("## Sequences per constant region\n")

2.3 Sequences per constant region

Code
# Create c_gene column only for plotting (will be removed before storing the df again)
db <- db %>% filter(locus %in% c("IGH", "IGK", "IGL")) %>%
    mutate(c_gene = alakazam::getGene(c_call, first=TRUE))

input_cgene_summary <- db %>%
   group_by(!!!rlang::syms(unique(c("sample_id",params$cloneby, "c_gene")))) %>%
   summarize(n=n(), .groups="drop") %>%
       pivot_wider(names_from=c_gene, values_from=n)

tab_caption <- paste0("Input data. Number of sequences in each ", 
                   paste("sample", params$cloneby, sep=", "),
                   " and c_gene"
      )
eetable(input_cgene_summary)$table
Table 2.4: Input data. Number of sequences in each sample, subject_id and c_gene

3 Clonal assignment

The input data contains clonal assignments, and are being reused.

4 Clone size distribution

Find your clone sizes table here. Most real datasets, will have most clones of size 1 (one sequence). Straight sequence count as a mesure of the size of the clones is not the best measure to compare clone size between samples due to possible disproportionate sampling. See the Clonal abundance section.

Description of terms:

  • clone_size_count: Clone size as sequence counts. In a sample (sample_id), the number of heavy chain sequences with the same clone_id.

  • clone_size_freq: Clone size as percent of the repertoire. clone_size_count divided by the number of heavy chain sequences in the sample (sample_id).

4.1 Number of clones (heavy chain, incl. singletons)

Code
# Add clone_size 
clone_sizes <- countClones(
            db %>% filter(isHeavyChain(locus)), # Keep heavy chains only
            groups=unique(c("sample_id", params$cloneby)))
Code
tmp <- eetable(clone_sizes, 
               caption=NULL, 
               outdir=params$outdir, file="clone_sizes_table")

db <- db %>%
   left_join(clone_sizes) %>%
   rename(
      clone_size_count = seq_count,
      clone_size_freq = seq_freq
   ) %>% 
   mutate_at(vars(starts_with("clone_size")), round,2)
Code
num_clones_table <- db %>%
   filter(isHeavyChain(locus)) %>% # Keep heavy chains only
   group_by(sample_id) %>%
   mutate(sequences=n()) %>%
   group_by(sample_id) %>%
   mutate(
      number_of_clones=length(unique(clone_id)),
   ) %>%
   group_by(!!!rlang::syms(unique(c("sample_id","sequences", params$cloneby, "number_of_clones")))) %>%
   summarize_at(vars(starts_with("clone_size")), list("min"=min, "median"=median, "max"=max)) %>%
   mutate_at(vars(starts_with("clone_size")), round, 2)

tab_caption <- "Summary of the number of clones, and clone size, per sample. Includes singletons (clone_size == 1)."
tab <- eetable(num_clones_table, caption=tab_caption, outdir=params$outdir, file="num_clones_table")
tab$table
Table 4.1: Summary of the number of clones, and clone size, per sample. Includes singletons (clone_size == 1). File can be found here: num_clones_table.tsv

4.2 Clone size distribution

Code
caption <- "Clone size distribution. Size is measured as number of heavy chain sequences belonging to the same clone."
clone_size_plot <- ggplot(clone_sizes, aes(x=seq_count, color=sample_id, fill=sample_id))+
    geom_bar() + theme_enchantr() +
    facet_wrap(~sample_id, scales = "free_y", ncol=3) +
    xlab("Clone size (Number of sequences per clone)")
clone_size_plot <- eeplot(clone_size_plot, 
       outdir=params$outdir, 
       file=knitr::opts_current$get('clone-size'),
       caption=caption
       )
ggplotly(clone_size_plot + theme(panel.spacing=unit(2, 'lines'), legend.position="right"))

Figure 4.1: Clone size distribution. Size is measured as number of heavy chain sequences belonging to the same clone. ggplot file: clone_size_plot.RData

Code
caption <- clone_size_plot$enchantr$html_caption

4.3 Clone size distribution without singletons

Code
caption <- "Clone size distribution, excluding singletons (subset to clone size > 1). Size is measured as number of heavy chain sequences belonging to the same clone."
clone_size_atleast2 <- ggplot(clone_sizes %>% filter(seq_count>1), aes(x=seq_count, color=sample_id, fill=sample_id))+
    geom_bar() + theme_enchantr() +
    facet_wrap(~sample_id, scales = "free_y", ncol=3) +
    xlab("Clone size (Sequences per clone)")
clone_size_atleast2 <- eeplot(clone_size_atleast2, 
       outdir=params$outdir, 
       file=knitr::opts_current$get('clone-size_atleast2'),
       caption=caption
       )
ggplotly(clone_size_atleast2 + theme(panel.spacing=unit(2, 'lines'), legend.position="right"))

Figure 4.2: Clone size distribution, excluding singletons (subset to clone size > 1). Size is measured as number of heavy chain sequences belonging to the same clone. ggplot file: clone_size_atleast2.RData

Code
caption <- clone_size_atleast2$enchantr$html_caption
Code
cat("Only singletons detected: there aren't clones of size>1.\n\n")

5 Clonal abundance

Code
# set empty default abundance
a <- new("AbundanceCurve")

Clonal abundance is the size of each clone (as a fraction of the entire repertoire). To correct for the different number of sequences in each of the samples, estimateAbundance estimates the clonal abundance distribution along with confidence intervals on these clone sizes using bootstrapping. 200 random bootstrap samples were taken, with size the number of sequences in the sample with less sequences (N). The y-axis shows the clone abundance (i.e., the size as a percent of the repertoire) and the x-axis is a rank of each clone, where the rank is sorted by size from larger (rank 1, left) to smaller (right). The shaded areas are confidence intervals.

Code
# calculate the rank-abundance curve
a <- estimateAbundance(db %>% filter(isHeavyChain(locus)),
                       group = "sample_id", min_n = params$min_n)
if (nrow(a@abundance)==0) {
    cat("\nAll groups failed to pass the threshold min_n=",params$min_n,". Skipping clonal abundance report.\n\n")
}
Code
# annotate
a@abundance <- a@abundance %>%
   left_join( db %>% 
                 select(any_of(c(unique(c("sample_id", 
                                                "subject_id",
                                                params$cloneby)))
                                       )
                  ) %>%
                 distinct(),
              
              )
p <- plotAbundanceCurve(a, annotate="depth", silent = T)
Code
tab <- eetable(a@abundance, file = "clonal_abundance", 
               outdir=params$outdir, show_max=10,
               caption="Example 10 lines of the clonal abundance file.")
tab$table %>%
        DT::formatRound(columns=c("p","p_sd","lower","upper"),digits=3)
Table 5.1: Example 10 lines of the clonal abundance file. File can be found here: clonal_abundance.tsv

5.1 Abundance plot by sample

Code
abundanceSample <- p + facet_wrap(~ subject_id + sample_id, ncol=3)
abundanceSample <- eeplot(abundanceSample, 
                     outdir=params$outdir, 
                     file=knitr::opts_current$get('abundanceSample'),
                     caption=paste0("Clonal abundance plot per `sample_id`. Clonal abundance was calculated with a sample of N=",
                      a@n,
                      " sequences with ",
                      a@nboot,
                      " boostrapping repetitions."),
                     a)
abundanceSample +
    theme(legend.position = "none")
Clonal abundance plot per `sample_id`. Clonal abundance was calculated with a sample of N=65 sequences with 200 boostrapping repetitions. <a href='ggplots/abundanceSample.RData'>ggplot file: abundanceSample.RData</a>

Figure 5.1: Clonal abundance plot per sample_id. Clonal abundance was calculated with a sample of N=65 sequences with 200 boostrapping repetitions. ggplot file: abundanceSample.RData

Code
caption <- abundanceSample$enchantr$html_caption  

6 Diversity

The clonal abundance distribution can be characterized using diversity statistics. Diversity scores (D) are calculated using the generalized diversity index (Hill numbers), which covers many different measures of diversity in a single function with a single varying parameter, the diversity order q.

The function alphaDiversity resamples the sequences (200 random bootstrapping events, with the number of sequences in the sample with less sequences (N)) and calculates diversity scores (D) over a interval of diversity orders (q). The diversity (D) is shown on the y-axis and the x-axis is the parameter q. - q = 0 corresponds to Species Richness - q = 1 corresponds to Shannon Entropy - q = 2 corresponds to Simpson Index

Inspection of this figure is useful to determine whether any difference in diversity between two repertoires depends on the statistic used or if it is a universal property.

The clonal diversity \(D\) of the repertoire was calculated according to the general formula of Hill Diversity numbers:

\[ \begin{aligned} ^{q}D = \left( \sum_{i=1}^Rp_i^q \right)^{1/(1-q)} \end{aligned} \]

where:

  • \(p_i\) is the proportion of unique sequences belonging to clone \(i\).
  • \(q\) are the values of the different diversity numbers.
  • \(R\) is the Richness, the number of different clones in the sample.

At \(q=1\) the function is undefined and the limit to zero equals the exponential of the Shannon Entropy:

\[ \begin{aligned} ^{1}D = exp \left( \sum_{i=1}^Rp_i ln(p_i) \right) \end{aligned} \]

The intuition about the different Hill Diversity values is the following:

  • At \(q=0\) the diversity index equals the number of clones in the sample.
  • At \(q=1\) the diversity index is the geometric mean of the clones in the sample, weighted by their proportion in the sample.
  • At \(q>1\) more weight is given to the clones with higher proportions in the sample.

6.1 Diversity curves

The following table shows the summary of the diversity calculations per sample.

Code
# generate the Hill diversity curve
d <- alphaDiversity(db %>% filter(isHeavyChain(locus)), 
                    group = "sample_id", 
                    min_n=params$min_n)

d@diversity <- d@diversity %>%
   left_join( db %>% 
                 select(any_of(c(unique(c("sample_id", 
                                                "subject_id", 
                                                params$cloneby)))
                                       )
                  ) %>%
                 distinct(),
              
              )

diversitySample <- plotDiversityCurve(d, silent = T, annotate="depth")

tab <- eetable(d@diversity, file = "clonal_diversity", 
               outdir=params$outdir, show_max=10,
               caption="Example 10 lines of the clonal diversity file.")
tab$table %>%
        DT::formatRound(columns=c("q","d","d_sd","d_lower", "d_upper", "e", "e_lower", "e_upper"),digits=3)
Table 6.1: Example 10 lines of the clonal diversity file. File can be found here: clonal_diversity.tsv

6.2 Diversity plot by sample

Code
# plot duplicated cells
diversitySample <- diversitySample + 
    geom_vline(xintercept = c(0,1,2), color = "grey50", linetype = "dashed") +
    facet_wrap(~sample_id + subject_id, ncol=3, scales = "free_x")
diversitySample <- eeplot(diversitySample, 
                     outdir=params$outdir, 
                     file=knitr::opts_current$get('label'),
                     caption=paste0("Clonal diversity per `sample_id`.Clonal diversity was calculated with a sample of N=",
                      d@n,
                      " sequences with ",
                      a@nboot,
                      " boostrapping repetitions."),
                     d)
diversitySample +
    theme(legend.position = "none")
Clonal diversity per `sample_id`.Clonal diversity was calculated with a sample of N=65 sequences with 200 boostrapping repetitions. <a href='ggplots/diversitySample.RData'>ggplot file: diversitySample.RData</a>

Figure 6.1: Clonal diversity per sample_id.Clonal diversity was calculated with a sample of N=65 sequences with 200 boostrapping repetitions. ggplot file: diversitySample.RData

Code
caption <- diversitySample$enchantr$html_caption

6.3 Diversity at different q

  • q=0 corresponds to Richness (number of different clones in the sample).
  • q=1 corresponds to Shannon diversity.
  • q=2 corresponds to Simpson diversity.
Code
# plot duplicated cells
div_data <- d@diversity %>% dplyr::filter(q==0 | q==1 | q==2)

div_at_q <- ggplot(div_data, aes(x=q, y=d, color=sample_id)) +
    geom_point() + theme_enchantr() +
    scale_x_continuous(breaks= c(0,1,2))
    
div_at_q <- eeplot(div_at_q, 
                     outdir=params$outdir, 
                     file=knitr::opts_current$get('label'),
                     caption="Diversity at different q values.",
                     div_at_q)
ggplotly(div_at_q)

(#fig:div_at_q)Diversity at different q values. ggplot file: div_at_q.RData

Code
caption <- div_at_q$enchantr$html_caption
Code
cat("# Mutation frequency\n")

7 Mutation frequency

Code
db <- shazam::observedMutations(db,
                                sequenceColumn = "sequence_alignment",
                                germlineColumn = "germline_alignment_d_mask",
                                frequency = T,
                                combine = T,
                                nproc = params$nproc)
Code
cat("Showing mutation frequency per c_gene only if the `c_call` column is present.\n")

Showing mutation frequency per c_gene only if the c_call column is present.

Code
cat("The mutation frequency per sequence is stored in the final dataframes in the `mu_freq` column.\n")

The mutation frequency per sequence is stored in the final dataframes in the mu_freq column.

Code
db_cgene <- db %>% filter(!is.na(c_gene), !(is.na(mu_freq)))
db_cgene$mu_freq <- as.numeric(db_cgene$mu_freq)

if (nrow(db_cgene) > 0) {
    mufreqSample <- ggplot(db_cgene, aes(x=c_gene, y=mu_freq, fill=c_gene)) + 
                    geom_boxplot() +
                    facet_wrap(~ subject_id + sample_id, ncol=3, scales = "free") +
                    theme_enchantr() +
                    theme(axis.text.x = element_text(angle = 45, hjust = 1))

    mufreqSample <- eeplot(mufreqSample, 
                        outdir=params$outdir, 
                        file='mutation_frequency_sample',
                        caption=paste0("Mutation frequency per C gene per `sample_id`."),
                        mufreqSample)
    print(mufreqSample +
        theme(legend.position = "none"))
    caption <- mufreqSample$enchantr$html_caption  
}
Mutation frequency per C gene per `sample_id`. <a href='ggplots/mutation_frequency_sample.RData'>ggplot file: mutation_frequency_sample.RData</a>

Figure 7.1: Mutation frequency per C gene per sample_id. ggplot file: mutation_frequency_sample.RData

Code
db_locus <- db %>% filter(!is.na(locus), !is.na(mu_freq))
db_locus$mu_freq <- as.numeric(db_locus$mu_freq)

if (nrow(db_locus) > 0) {
    mufreqSampleLocus <- ggplot(db_locus, aes(x=locus, y=mu_freq, fill=locus)) + 
                    geom_boxplot() +
                    facet_wrap(~ subject_id + sample_id, ncol=3, scales = "free") +
                    theme_enchantr() +
                    theme(axis.text.x = element_text(angle = 45, hjust = 1))

    mufreqSampleLocus <- eeplot(mufreqSampleLocus, 
                        outdir=params$outdir, 
                        file='mutation_frequency_sample_locus',
                        caption=paste0("Mutation frequency per locus per `sample_id`."),
                        mufreqSampleLocus)
    print(mufreqSampleLocus +
        theme(legend.position = "none"))
    caption <- mufreqSampleLocus$enchantr$html_caption 
}

8 Final repertoires and tables

Summary tables:

Find your final repertoires here:

9 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] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] plotly_4.10.4         ComplexHeatmap_2.18.0 enchantr_0.1.21      
##  [4] dowser_2.3            scoper_1.3.0          shazam_1.2.0         
##  [7] alakazam_1.3.0        ggplot2_3.5.2         airr_1.5.0           
## [10] tidyr_1.3.1           dplyr_1.1.4           DT_0.33              
## 
## loaded via a namespace (and not attached):
##   [1] RColorBrewer_1.1-3          rstudioapi_0.17.1          
##   [3] jsonlite_2.0.0              shape_1.4.6.1              
##   [5] magrittr_2.0.3              farver_2.1.2               
##   [7] rmarkdown_2.29              GlobalOptions_0.1.2        
##   [9] fs_1.6.6                    zlibbioc_1.52.0            
##  [11] vctrs_0.6.5                 memoise_2.0.1              
##  [13] Rsamtools_2.22.0            ggtree_3.14.0              
##  [15] htmltools_0.5.8.1           S4Arrays_1.6.0             
##  [17] progress_1.2.3              SparseArray_1.6.2          
##  [19] gridGraphics_0.5-1          sass_0.4.10                
##  [21] KernSmooth_2.23-26          bslib_0.9.0                
##  [23] htmlwidgets_1.6.4           cachem_1.1.0               
##  [25] GenomicAlignments_1.42.0    igraph_2.1.4               
##  [27] lifecycle_1.0.4             iterators_1.0.14           
##  [29] pkgconfig_2.0.3             Matrix_1.7-3               
##  [31] R6_2.6.1                    fastmap_1.2.0              
##  [33] GenomeInfoDbData_1.2.13     MatrixGenerics_1.18.1      
##  [35] clue_0.3-66                 digest_0.6.37              
##  [37] aplot_0.2.5                 colorspace_2.1-1           
##  [39] patchwork_1.3.0             S4Vectors_0.44.0           
##  [41] crosstalk_1.2.1             GenomicRanges_1.58.0       
##  [43] labeling_0.4.3              phylotate_1.3              
##  [45] httr_1.4.7                  polyclip_1.10-7            
##  [47] abind_1.4-8                 compiler_4.4.3             
##  [49] bit64_4.6.0-1               withr_3.0.2                
##  [51] doParallel_1.0.17           BiocParallel_1.40.1        
##  [53] viridis_0.6.5               ggforce_0.4.2              
##  [55] MASS_7.3-64                 DelayedArray_0.32.0        
##  [57] rjson_0.2.23                tools_4.4.3                
##  [59] ape_5.8-1                   quadprog_1.5-8             
##  [61] glue_1.8.0                  nlme_3.1-167               
##  [63] cluster_2.1.8               ade4_1.7-23                
##  [65] generics_0.1.3              seqinr_4.2-36              
##  [67] gtable_0.3.6                tzdb_0.5.0                 
##  [69] data.table_1.17.0           hms_1.1.3                  
##  [71] tidygraph_1.3.1             XVector_0.46.0             
##  [73] BiocGenerics_0.52.0         ggrepel_0.9.6              
##  [75] foreach_1.5.2               pillar_1.10.2              
##  [77] markdown_2.0                stringr_1.5.1              
##  [79] vroom_1.6.5                 yulab.utils_0.2.0          
##  [81] circlize_0.4.16             tweenr_2.0.3               
##  [83] treeio_1.30.0               lattice_0.22-6             
##  [85] bit_4.6.0                   tidyselect_1.2.1           
##  [87] Biostrings_2.74.1           knitr_1.50                 
##  [89] gridExtra_2.3               bookdown_0.42              
##  [91] IRanges_2.40.1              SummarizedExperiment_1.36.0
##  [93] stats4_4.4.3                xfun_0.52                  
##  [95] graphlayouts_1.2.2          Biobase_2.66.0             
##  [97] diptest_0.77-1              matrixStats_1.5.0          
##  [99] stringi_1.8.7               UCSC.utils_1.2.0           
## [101] lazyeval_0.2.2              ggfun_0.1.8                
## [103] yaml_2.3.10                 evaluate_1.0.3             
## [105] codetools_0.2-20            ggraph_2.2.1               
## [107] tibble_3.2.1                ggplotify_0.1.2            
## [109] cli_3.6.4                   munsell_0.5.1              
## [111] jquerylib_0.1.4             Rcpp_1.0.14                
## [113] GenomeInfoDb_1.42.3         png_0.1-8                  
## [115] parallel_4.4.3              readr_2.1.5                
## [117] prettyunits_1.2.0           bitops_1.0-9               
## [119] phangorn_2.12.1             viridisLite_0.4.2          
## [121] tidytree_0.4.6              scales_1.3.0               
## [123] purrr_1.0.4                 crayon_1.5.3               
## [125] GetoptLong_1.0.5            rlang_1.1.6                
## [127] fastmatch_1.1-6