In rnabioco/djvdj: A collection of single-cell V(D)J tools

# Chunk opts
knitr::opts_chunk$set(
  collapse   = TRUE,
  comment    = "#>",
  warning    = FALSE,
  message    = FALSE,
  fig.width  = 8,
  fig.height = 3
)

This vignette provides detailed examples for quantifying repertoire diversity. For the examples shown below, we use data for splenocytes from BL6 and MD4 mice collected using the 10X Genomics scRNA-seq platform. MD4 B cells are monoclonal and specifically bind hen egg lysozyme.

library(djvdj)
library(Seurat)
library(ggplot2)

# Load GEX data
data_dir <- system.file("extdata/splen", package = "djvdj")

gex_dirs <- c(
  BL6 = file.path(data_dir, "BL6_GEX/filtered_feature_bc_matrix"),
  MD4 = file.path(data_dir, "MD4_GEX/filtered_feature_bc_matrix")
)

so <- gex_dirs |>
  Read10X() |>
  CreateSeuratObject() |>
  AddMetaData(splen_meta)

# Add V(D)J data to object
vdj_dirs <- c(
  BL6 = system.file("extdata/splen/BL6_BCR", package = "djvdj"),
  MD4 = system.file("extdata/splen/MD4_BCR", package = "djvdj")
)

so <- so |>
  import_vdj(vdj_dirs, define_clonotypes = "cdr3_gene")

Calculating diversity

To calculate repertoire diversity and store the results in the object meta.data, the calc_diversity() function can be used. This function is designed to specifically work with the R package abdiv. The diversity metric can be selected by passing the name of the function to the method argument. Any alpha diversity function from the abdiv package that takes species counts as input can be used. Be sure to read the documentation for the function you are using to ensure it is appropriate for your analysis.

In this example we are calculating the Shannon entropy for BL6 and MD4 samples.

so_vdj <- so |>
  calc_diversity(
    data_col    = "clonotype_id",
    cluster_col = "sample",
    method      = abdiv::shannon
  )

Estimations of species diversity are influenced by sample size. One approach to address this is to equalize the number of cells present in each cluster. The downsample argument will randomly sample cells so each sample being tested has the same number of cells as the smallest cluster. The bootstrapped standard error can also be calculated by setting the number of bootstrap samples with the n_boots argument.

so_vdj <- so |>
  calc_diversity(
    data_col    = "clonotype_id",
    cluster_col = "sample",
    downsample  = TRUE,
    n_boots     = 50
  )

Diversity metrics can also be calculated for a specific chain. To do this, the column passed to the data_col argument must contain per-chain data, such as CDR3 amino acid or nucleotide sequences. In this example diversity is calculated based only on heavy chain CDR3 sequences.

so_vdj <- so |>
  calc_diversity(
    data_col    = "cdr3_nt",
    cluster_col = "sample",
    chain       = "IGH"
  )

Plotting diversity

The plot_diversity() function will create plots summarizing repertoire diversity for each sample. A named list of functions can also be passed to plot multiple metrics. Two metrics for measuring diversity are the Simpson index and Shannon entropy. Both of these metrics are influenced by species richness (number of unique sequences) and evenness (relative abundance of sequences). Pielou's index will specifically measure species evenness. For these metrics, maximally diverse samples will return a value of 1.

As expected, BL6 B cells have a very diverse repertoire, while MD4 cells have a restricted repertoire.

div_fns <- list(
  "simpson" = abdiv::simpson,
  "shannon" = abdiv::shannon,
  "pielou evenness" = abdiv::pielou_e
)

so |>
  plot_diversity(
    data_col    = "clonotype_id",
    cluster_col = "sample",
    method      = div_fns
  )

Diversity plots can also be separated based on an additional grouping variable such as treatment group (e.g. pacebo vs drug) or disease status (e.g. healthy vs disease). This will generate boxplots with each point representing a label present in the cluster_col column. In this example we have 3 BL6 and 3 MD4 samples, so there should be 5 points shown for each boxplot.

so |>
  plot_diversity(
    data_col    = "clonotype_id",
    cluster_col = "sample",
    group_col   = "orig.ident",
    method      = div_fns
  )

Additional arguments are provided to adjust plot aesthetics. The plot_colors parameter can be used to modify colors, the panel_nrow and panel_scales arguments will adjust the plot scales and number of rows used to arrange plots.

so |>
  plot_diversity(
    data_col    = "clonotype_id",
    cluster_col = "sample",
    group_col   = "orig.ident",
    method      = div_fns,
    plot_colors = c(BL6 = "#3182bd", MD4 = "#fec44f"),
    panel_nrow  = 2
  )

plot_diversity() returns a ggplot object that can be modified with ggplot2 functions such as ggplot2::theme(). Plots can be further adjusted by passing aesthetic parameters directly to ggplot2, e.g. alpha, linetype, color, etc.

so |>
  plot_diversity(
    data_col    = "clonotype_id",
    cluster_col = "sample",
    method      = div_fns,

    alpha    = 0.5,         # parameters to pass to ggplot2
    linetype = 2,
    color    = "black"    
  ) +
  theme(strip.text = element_text(face = "bold"))

Rarefaction curves

Another approach to ensure differences in sample size are not having an undue influence on diversity results is to plot rarefaction curves. This method involves calculating species diversity for different sized samples generated by randomly downsampling each cluster. By default the bootstrapped 95% confidence interval will also be plotted.

Calculations used to generate rarefaction curves are performed using the iNEXT package. There are three diversity calculations that can be specified with the method argument:

'richness', species richness, this is equivalent to the calculation performed by abdiv::richness()
'shannon', the exponential of Shannon entropy
'invsimpson', the inverse Simpson index, this is equivalent to the calculation performed by abdiv::invsimpson()

so |>
  plot_rarefaction(
    data_col    = "clonotype_id",
    cluster_col = "orig.ident",
    method      = c("richness", "shannon", "invsimpson"),
    plot_colors = c("#3182bd", "#fec44f")
  )

If the 95% confidence interval is not desired, set n_boots to 0. In this example we also plot a separate line for each BL6 and MD4 sample.

so |>
  plot_rarefaction(
    data_col    = "clonotype_id",
    cluster_col = "sample",
    method      = c("richness", "shannon", "invsimpson"),
    n_boots     = 0
  )