PlatypusV3 AgedCNS vignette

6.1 Checking vgm output

Now we can analyze the VDJ repertoire data before integrating VDJ and GEX information. This may be useful if only VDJ libraries have been sequenced without the accompanying gene expression data. We first start by examining the structure of vgm[[1]]

print(colnames(vgm[[1]]))

If GEX data has been integrated, this dataframe also contains dimensional reduction coordinates and cluster id. If trim.and.align was set to TRUE, the columns VDJ_sequence_nt_trimmed to VJ_trimmed_ref contain aligned and reference sequences. The clonotype_frequency column takes input from the clonotyping output of cellranger, which is saved int the clonotype_id_10x column.

This object contains both cells with less than 2 V(D)J chains as well as with more than 2 To show how this is done the VDJ_cgene column is taken as an example

print(unique(vgm[[1]]$VDJ_cgene))

Apart from classical isotypes, we can see an empty (““) element. This is for cells, for which only a VJ chain is available. Secondly we can see 3 examples of cells containing more than one VDJ chains. These are delimited by”;“. Importantly, this format is maintained throughout all columns.

Both cells with less and more than 2 chains can be cumbersome for analysis. Platypus V3 functions generally support such cells. To check their quantity and, if needed, exclude these, the following line can be used

cat(" Cell count by number of VDJ chains")
print(table(vgm[[1]]$Nr_of_VDJ_chains))
cat("\n Cell count by number of VJ chains")
print(table(vgm[[1]]$Nr_of_VJ_chains))

#Subset the VGM matrix to only include cells with 1 VDJ and 1VJ chain
#vgm[[1]] <- subset(vgm[[1]], Nr_of_VJ_chains == 1 & Nr_of_VJ_chains == 1)

6.2 Changing the clonotype strategy

Often paired nucleotide CDRH3 and CDRL3 clonotyping may not be the best strategy given somatic hypermutation may occur in the CDR3 region. Therefore, there could be highly similar clones that likely bind the same antigen that are officially part of different clonal families. To address this, we have added a function that allows for various heuristic clonotyping strategies. This involves clonotyping by identical amino acid CDRH3 + CDRL3 seuqence, identical germline usage, or seqeunce homology requirements. This function works both for V3 of platypus.

global.clonotype If set to TRUE, clonotypes will be generated across all samples of the vgm. This may be useful if more than one sample has been taken from the same animal (e.g. spleen and bone marrow)

A note for cells with aberrant numbers of VDJ and VJ chains: the functions implements a hierarchical clonotyping strategy for such cells. In brief, clonotypes are determined based on all cells with exactly 1VDJ and 1VJ chain. Other cells are then “joined in” in post. Two examples for illustration of the hierarchical mode “single.chains”: 1. A cell with 0VDJ and 1VJ chain is joined to an existing clonotype containing the same or a homologous VJ chain, depending of the clonotyping strategy. 2. A cell with 1VDJ and 2VJ chain will be matched to a clonotype that contains either of the two possible combinations of the chains of that aberrant cell. In case that no existing clonotype matches the abberant cell in question, a new clonotype will be assigned In case of multiple matching clonotypes, abberant cells will be assigned to the most frequent of them. Cells with 2VDJ and 2VJ chains are excluded as these most likely correspond to doublets. Further the function implements the hierarchical mode “double.and.single.chains” Here the function will proceed as if set to “single.chains” but include two more steps 3. Check the frequency of each cell 1 VDJ 2 VJ chain exact clone (by exact nucleotide CDR3 matching). Only if this count exceeds the triple.chain.count.threshold, the clone is used as a “hub clone”. This protects from merging clonotypes on the basis of rare doublets. 4. Merge existing clonotypes into the 1 VDJ 2 VJ clonotypes as they match with the assumption that e.g. a cell with 1 VDJ 1 VJ is part of that same clonotype, but missing a VJ chain due to stochastical sampling

In the example that follows, we use VDJ_clonotype to group cells into clones based on identical CDRH3 + CDRL3 amino acid sequence. We will compare this to the case in which we group B cells by using the same germline genes (both heavy chain and light chain).

The function adds new columns containing clonotype ids, features and frequencies, named based on the clonotyping strategy. This allows multiple clonotyping strategy outputs to be stored in the same VGM object

The columns clonotype_frequency and clonotype_id are considered the “active” slot, as functions downstream pull from them. These two columns are updated every time the VDJ_clonotype function is run

vgm[[1]] <- VDJ_clonotype(VDJ = vgm[[1]], clone.strategy = "cdr3.aa", global.clonotype = F, VDJ.VJ.1chain = F, hierarchical = "single.chains") #Not filtering cells with counts other than 1VDJ 1VJ chain and integrating these cells hierarchically into clonotypes

cat(" Nr and distribution of clonotypes using exact CDR3.aa matching \n")
print(length(unique(vgm[[1]]$clonotype_id_cdr3.aa)))
print(table(vgm[[1]]$clonotype_frequency_cdr3.aa)) #Check distribution of clonotypes with identical CDR3 aa sequences

vgm[[1]] <- VDJ_clonotype(VDJ = vgm[[1]], clone.strategy = "hvj.lvj", global.clonotype = F, output.format = "vgm", VDJ.VJ.1chain = F, hierarchical = "single.chains")

cat("\n Nr and distribution of clonotypes using germline gene matching \n ")
print(length(unique(vgm[[1]]$clonotype_id_hvj.lvj)))
print(table(vgm[[1]]$clonotype_frequency_hvj.lvj)) #Check distribution of clonotypes with identical germline genes

6.3 Extracting full-length sequences from the VDJRegion

At the start of this vignette the VGM function was run with trim.and.align = F. Here we run the function again, but this time with trim.and.align = T. This does take significantly longer

To accelerate runtime, parallel computing options are available. If running a Windows machine parallel.processing should be set to “parlapply” This does require the package Parallel and its dependencies. The additional parameter numcores allows to set the numbers of cores to use. This is important when running the function on a cluster. If on a MAC or Linux machine, set parallel.processing to “mclapply”

Reference sequences are optained from the cellranger output and thereby the 10x Genomics reference. Further trimming is made based on annotations by cellranger (feature$region_type). Alignment is done via Biostrings::pairwiseAlignment and the alignment with maximum score is returned. If alignments are not complete, the vgm parameters gap.opening.cost and gap.extension.cost can be modified

vgm <- VDJ_GEX_matrix(VDJ.out.directory.list = VDJ.out.directory.list,
                               GEX.out.directory.list = GEX.out.directory.list,
                               GEX.integrate = T,
                               VDJ.combine = T,
                               integrate.GEX.to.VDJ = T,
                               integrate.VDJ.to.GEX = T,
                               exclude.GEX.not.in.VDJ = F,
                               filter.overlapping.barcodes.GEX = T,
                               filter.overlapping.barcodes.VDJ = T,
                               exclude.on.cell.state.markers = c("CD3E"),
                               get.VDJ.stats = T,
                               parallel.processing = "none",
                               trim.and.align = T, #Set this to TRUE for full length sequence recovery
                               group.id = c(1,2),
                               gap.opening.cost = 10,
                               gap.extension.cost = 4) #Tweak to optimize alignments if neccessary

#saving this for later
#saveRDS(vgm, "VDJ_GEX_matrix_agedCNS.rds")

No trimmed and reference sequence columns should be filled.

print(vgm[[1]][1,])

A different way to get germline reference sequences is by using MIXCR. MIXCR can be called directly from R on both Windows and MAC machines. Given that output files have to be read in, it is important to set the working directory correctly and (Essential for Windows users) have a version of MIXCR available in that working directory.

Moreover to quantify the number of somatic variants or to extract full-length sequences for expression, it is often useful to have the nucleotide sequence from framework region 1 (FR1) to framework region 4 (FR4). Using the call_MIXCR function, the full-length VDJRegion sequences can be added to the clonal information and easily extracted thereafter. This function works on UNIX/Mac and furthermore requires that mixcr is already downloaded locally (with license agreement).

FOR MAC/UNIX users

One needs to supply the directory to the executable in the call_MIXCR function as below. Either “mmu” or “hsa” for mouse and human, respectively. Again, the format is similar to the input, in that the outer list corresponds to the individual repertoire and the inner list is a dataframe with various information, including the full-length VDJ sequences (e.g., VDJ.AA.LC and VDJ.AA.HC for the light and heavy chain amino acid sequence). One will notice that the germline sequence is still very long (e.g., in the example below the “full_HC_germline” length is over 600 nucleotides). This will be filled in using the separate function VDJ_extract_germline.

### WARNING: You will need to download MiXCR and change the mixcr.directory to the location of MiXCR
#VDJ_mixcr_out <- VDJ_call_MIXCR(VDJ.matrix = vgm[[1]], mixcr.directory = "~/Downloads/mixcr.jar",species = "mmu", platypus.version = "v3", operating.system = "Darwin", simplify = T)
#set simplify to T to append only a selected columns of the MIXCR output to the vgm matrix

FOR Windows users The mixcr.jar executable needs to be in the current working directory!

#check working directory
#getwd()

### WARNING: You will need to download MiXCR and have it in your current working directory
VDJ_mixcr_out <- VDJ_call_MIXCR(VDJ = vgm[[1]], mixcr.directory = "Is set automatically to current working directory",species = "mmu", platypus.version = "v3", operating.system = "Windows", simplify = T)
#set simplify to T to append only a selected columns of the MIXCR output to the vgm matrix
#set to False to save as separate object with the complete MIXCR output as in this case

6.5 Plotting SHMs

Directly from the MIXCR output we can plot the frequency of somatic hypermutations and run basic tests of significance between groups

VDJ_plot_SHM(VDJ = VDJ_mixcr_out, group.by = "sample_id", quantile.label = 0.95)

7.1 Visualizing clonal frequencies (V3 only)

For a basic view of clonal expansion the VDJ_clonal_donut produces circular plots per sample. Label position and size should only be adjusted after deciding on a plot export format. This function uses the clonotype_id column as input. If VDJ_clonotyping function was used, its result are stored in that column. If not, the default 10x clonotyping is stored there. To retrieve default 10x clonotyping, use the column clonotyping_id_10x

donuts <- VDJ_clonal_donut(VDJ = vgm[[1]], expanded.colors = c("grey50", "grey65", "grey80"), non.expanded.color = "black", counts.to.use = "freq_column")
#Counts to use = "freq_column" uses the counts in the clonotype_frequency column. Counts to use = "vgm" simply counts the rows of a given clonotype in the VGM table (these counts may differ if cells have been filtered out due to overlapping barcodes or if another clonotyping strategy was used)

donuts[[1]]

The black section shows cells in unexpanded clones. The rest shows expanded clones by their frequency. The middle label shows the total number of clones (cells)

7.2 Calculating common repertoire diversity metrics

This requires the Vegan package. We can calculate the diversity for any column(s) in the vgm. If more than one column is provided, the content will be pasted together before calculating diversity metrics

#Shannon Evenness for the VDJ chain CDR3
diversity_plot <- VDJ_diversity(VDJ = vgm[[1]],feature.columns = c("VDJ_cdr3s_aa"),grouping.column = "sample_id",metric = c("shannonevenness"), platypus.version = "v3", subsample.to.same.n = T)
diversity_plot

#Gini-Simpson index for pasted VDJ and VJ chain CDR3s
diversity_plot <- VDJ_diversity(VDJ = vgm[[1]],feature.columns = c("VDJ_cdr3s_aa", "VJ_cdr3s_aa"),grouping.column = "sample_id",metric = c("ginisimpson"), platypus.version = "v3", subsample.to.same.n = T)
diversity_plot

#exact values can be retrived as
print(head(diversity_plot$data))

#Jaccard index between repertoires of the two samples
diversity_plot <- VDJ_diversity(VDJ = vgm[[1]],feature.columns = c("VDJ_cgene"),grouping.column = "sample_id",metric = c("jaccard"), platypus.version = "v3", subsample.to.same.n = T)
diversity_plot

7.3 Examining isotypes

Additional VDJ functions include plotting the clonal expansion for each clone. This can be performed by the VDJ_clonal_expansion function and setting the number of clones to show.

The function has two operating modes: Set color.by to “isotype” to exctract and color bars by isotype. Set subtypes to TRUE to view IgG subtypes. Set color.by to any other column in vgm[[1]] to examine the distribution of other parameters (This will be used during VDJ - GEX integration)

The following code extracts and plots the isotype distribution of the top 30 clones. We can see a clear IgA isotype majority of the top four clones in the first patient when using the default clonotyping strategy. We can additionally use the new clonotyping strategies to compare how changing the clonal defintion impacts the clonal expansion profiles. We simply supply the output from VDJ_clonotype

clonal_out <- VDJ_clonal_expansion(VDJ = vgm[[1]],celltype = "Bcells",clones = "30", platypus.version = "v3", group.by = "sample_id", color.by = "isotype", isotypes.to.plot = "all", treat.incomplete.clones = "exclude", treat.incomplete.cells = "proportional")
#group by specifies how many separate plots should be generated. If vgm contains global clonotype information this can be set to "none"
print(clonal_out[[1]])

7.4 Sequence similarity networks

Other functions are specifically tailored to repertoire analysis - such as VDJ_network, which creates a sequence similarity network between repertoires or within a repertoire by connecting those clones with sequence similarity. This function relies upon igraph to visually display and construct the graph - which means that networks with high number of sequences will not display easily. In the following example we are using a small dataset. In case of a bigger dataset one may subsample by using the sample.n function on vgm[[1]]. Setting the per.mouse argument to false indicates that one network for multiple repertoires should be produced.

#subsampled_VGM <- dplyr::sample_n(vgm[[1]], 60)

agedCNS_igraph <- VDJ_network(VDJ = vgm[[1]],per.sample = T,distance.cutoff = 8, platypus.version = "v3")

for(i in 1:2){
igraph::plot.igraph(agedCNS_igraph[[1]][[i]],vertex.label=NA,vertex.size=7)
}

#For a plot including clonotype frequencies
agedCNS_igraph <- VDJ_network(VDJ = vgm[[1]],per.sample = F,distance.cutoff = 8, platypus.version = "v3")

igraph::plot.igraph(agedCNS_igraph[[4]],vertex.label=NA,vertex.size=3+(.03*as.numeric(agedCNS_igraph[[2]]$clonotype_frequency)),vertex.color= as.factor(vgm[[1]]$sample_id))

For more details see the documentation of the VDJ_network function, but essentially information such as clonal frequency and which sample (here still indicated by the “sample_id” column) are stored in the second element of the output list. Here we can see only a few clones that are showing connections (produced by edges between those with a distance of 8 amino acids or less between heavy and light chain paired CDR3 sequence homology)

7.5 Germline gene usage heatmaps

It is also possible to produce heatmaps of the germline gene usage in the context of heavy chain V gene and light chain V gene. The output of the VDJ_Vgene_usage function is a matrix for each repertoire corresponding to the order specified by VDJ_analyze. The outer list corresponds to the sample and the inner list corresponds to a matrix, where the rows correspond to the heavy chain V genes and the columns correspond to the light chains of the V genes. Therefore the output[[1]][i,j] corresponds to the number of clones using the combination of IGH-Vgene[i] and IGK/L-Vgene[j].

#First calculate adjacency matrix for V gene usage
agedCNS_Vgene_usage <- VDJ_Vgene_usage(VDJ = vgm[[1]], platypus.version = "v3")

#library(pheatmap)
pheatmap::pheatmap(agedCNS_Vgene_usage[[1]],show_rownames = F,show_colnames = F)

print(class(agedCNS_Vgene_usage[[1]]))
print(head(rownames(agedCNS_Vgene_usage[[1]])))
print(head(colnames(agedCNS_Vgene_usage[[1]])))

This can then be easily plotted as a heatmap to observe patterns between repertoires or can be used to calculate V gene correlation using the “pheatmap” package.

Platypus also allows a separate analysis of V gene usage for HC and LC. The VDJ_Vgene_usage_barplot allows the user to plot most frequently used IgH or IgK/L V genes. By default, this function only provides visualizations for the HC V genes, but can also provide for the LC if LC. Vgene is set to TRUE. The User can also select the number of most used genes to be depicted.

agedCNS_Vgene_usage_barplot <- VDJ_Vgene_usage_barplot(VDJ = vgm[[1]], HC.gene.number = 10, LC.Vgene = T, LC.gene.number = 10, platypus.version = "v3")
agedCNS_Vgene_usage_barplot[[1]]

#VDJ chains
agedCNS_Vgene_usage_stackedbarplot <- VDJ_Vgene_usage_stacked_barplot(VDJ = vgm[[1]], LC.Vgene = F,HC.gene.number = 10, Fraction.HC = 1, platypus.version = "v3")
agedCNS_Vgene_usage_stackedbarplot

Furthermore, we can also produce a circular visualization of how V and J genes are combined throughout the repertoire. In the example that follows we use VDJ_VJ_usage_circos to look at the V gene with the corresponding J gene for each expanded clonotype.

#VDJ and VJ V and J gene pairing
vj_circos_bcells <- VDJ_VJ_usage_circos(VDJ = vgm[[1]], c.threshold = 1,label.threshold=2,cell.level = T, A.or.B = "both", platypus.version = "v3")

#VDJ and VJ pairing
HL_circos_bcells <- VDJ_alpha_beta_Vgene_circos(VDJ = vgm[[1]], c.threshold = 1,label.threshold=2,cell.level = T, V.or.J= "both", platypus.version = "v3")

7.6. Assessing Nr of variants per clone

Depending on the clonotype strategy, within a clone several different variants of CDR3s or full length sequences may be contained. To count variants, the function VDJ_variants_per_clone is used. It returns a table per sample_id with stats about variants

variants_agedCNS <- VDJ_variants_per_clone(VDJ = vgm[[1]], variants.of = c("VDJ_cdr3s_aa", "VJ_cdr3s_aa"), clonotypes.col = "clonotype_id_10x", split.by = "sample_id", stringDist.method = "levenshtein")

head(variants_agedCNS[[1]])

#set split.by to "none" if clonotyping was conducted across all samples

7.7. Searching for overlap between repertoires

Public clones or sequences are sequences that are shared between individuals. The VDJ_overlap_heatmap function can quantify this overlap and return a list of overlapping items

#overlap of VDJ V genes
VDJv_overlap <- VDJ_overlap_heatmap(VDJ = vgm[[1]],feature.columns = c("VDJ_vgene") ,grouping.column = "sample_id", axis.label.size = 20, pvalues.label.size = 12, platypus.version = "v3", add.barcode.table = T, plot.type = "ggplot")

VDJv_overlap[[2]] #summary of overlap
VDJv_overlap[[2]] #Table of overlapping items

#overlap of clones
VDJv_overlap <- VDJ_overlap_heatmap(VDJ = vgm[[1]],feature.columns = c("VDJ_cdr3s_aa","VJ_cdr3s_aa") ,grouping.column = "sample_id", axis.label.size = 20, pvalues.label.size = 12, platypus.version = "v3", add.barcode.table = T, plot.type = "ggplot")
#Pheatmap plots will function only with length(unique(grouping.column)) > 3

Given the small size of the two repertoires, it is unsurprising that no clones are shared

8.1. Integrating transcriptional clusters to the VDJ objects

One thing we may ask is how similar the B or T cells in a given clonal family are on the transcriptional level. The basic work of matching barcodes from V(D)J sequencing to those from GEX sequencing was already done in the VGM function. This part is therefore dedicated to functions that help explore and exploite this data

First we get basic stats on the overlap between GEX and VDJ

#Nr of cells for which VDJ info is available
nrow(vgm[[1]])

#Nr of cells for which GEX info is available
ncol(vgm[[2]])

#VDJ sequences for which GEX is available
sum(vgm[[2]]$VDJ_available)

#We can also plot this
Seurat::DimPlot(vgm[[2]],reduction = "umap", group.by = "VDJ_available", shuffle = T)

8.2. Relating clonal expansion to transcriptional cluster membership

The VDJ_clonal_expansion function has been used above to plot isotype information. Here its exdented functionality is used to color clones by their transcriptional cluster

clonal_out <- VDJ_clonal_expansion(VDJ = vgm[[1]],celltype = "Bcells",clones = "30", platypus.version = "v3", group.by = "sample_id", color.by = "seurat_clusters")
#group by specifies how many separate plots should be generated. If vgm contains global clonotype information this can be set to "none"
print(clonal_out[[1]][[2]])

Similar plots can be generated using the GEX_phenotype_per_clone function

clonal_out <- GEX_phenotype_per_clone(GEX = vgm[[2]], GEX.clonotypes = "topclones", GEX.group.by = "seurat_clusters", platypus.version = "v3")
#If vgm contains global clonotype information this can be set global.clonotypes to TRUE
clonal_out[[1]]

8.3. Visualizing clones on the 2 dimensional landscape

To better visualize this, the cells of any clone can also be highlighted on a dimensional reduction

Because the legend of these plots can get quite large, we can split the plot and the legend and draw them separately. For this the gridExtra and cowplot package is required

#here we overlay the top 5 clones
plot_out <- VDJ_GEX_overlay_clones(GEX = vgm[[2]], reduction = "umap", n.clones = 5, by.sample = F, ncol.facet = 1, split.plot.and.legend = F, pt.size = 0.5)

plot_out[[1]] # the plot

#We can also plot any clone of interest by adding a column were TRUE values select which clones are to be plotted
interesting_clones <- c("clonotype7", "clonotype42")
vgm[[2]]@meta.data$clones_to_plot <- FALSE
vgm[[2]]@meta.data$clones_to_plot[which(vgm[[2]]$clonotype_id_10x %in% interesting_clones)] <- TRUE

plot_out <- VDJ_GEX_overlay_clones(GEX = vgm[[2]], reduction = "umap", clones.to.plot = "clones_to_plot", by.sample = F, ncol.facet = 1, split.plot.and.legend = T, pt.size = 0.5)

plot_out[[1]] # the plot
plot(plot_out[[2]]) # the legend. Alternatively use gridExtra::grid.arrange(plot_out[[2]])

8.4 Specific gene expression information on the clonal level

We previously integrated the GEX information into the format of the VDJ output. However, we may want to ask how the gene expression looks for certain clonotypes (e.g., how many of the cells in the top clone are expressing markers of activated B cells). This is done by adding an extra column to the expanded clones and examining transcritional differences to cells of non expanded clones

vgm[[2]]$Expanded <- FALSE
vgm[[2]]$Expanded[which(vgm[[2]]$clonotype_frequency > 1)] <- TRUE

table(vgm[[2]]$Expanded)

#We can use the dottile function to look at selected genes
GEX_dottile_plot(GEX = vgm[[2]], genes = c("CD19", "CD74","IL2RA", "CD27","CD80"), group.by = "Expanded", threshold.to.plot = 1) + ggplot2::theme(legend.position = "bottom")

For a more unbiased view, DEGs can be calculated between expanded and non-expanded clones

DE_byexpansion <- GEX_DEgenes(GEX = vgm[[2]], min.pct = 0.25, group1 = "TRUE", group2 = "FALSE", grouping.column = "Expanded", return.plot = "volcano", label.n.top.genes = 10, platypus.version = "v3")

head(DE_byexpansion[[1]])
DE_byexpansion[[2]]

Indeed inflammation associated markers as well as MHC clomplex components are upregulated in expanded clones.