2. Loading the VGM
library(Platypus)##
## Attaching package: 'Platypus'## The following object is masked _by_ '.GlobalEnv':
##
## VDJ_call_MIXCRlibrary(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──## ✔ ggplot2 3.3.5 ✔ purrr 0.3.4
## ✔ tibble 3.1.6 ✔ dplyr 1.0.7
## ✔ tidyr 1.1.4 ✔ stringr 1.4.0
## ✔ readr 2.1.1 ✔ forcats 0.5.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()library(foreach)##
## Attaching package: 'foreach'## The following objects are masked from 'package:purrr':
##
## accumulate, whenlibrary(reticulate)
#Fetch only samples 1 (bone marrow, 3 months-old mice), 6 (bone marrow, 12 months-old), 9 (bone marrow, 18 months-old), 12 (spleen, 3 months-old mice), 17 (spleen, 12 months-old mice), and 20 (spleen, 18 months-old mice) - one sample from every group.
PlatypusDB_fetch(PlatypusDB.links = c("agrafiotis2021a/TNFR2.BM.3m.S1/ALL",
"agrafiotis2021a/TNFR2.BM.12m.S6/ALL",
"agrafiotis2021a/TNFR2.BM.18m.S9/ALL",
"agrafiotis2021a/TNFR2.SP.3m.S12/ALL",
"agrafiotis2021a/TNFR2.SP.12m.S17/ALL",
"agrafiotis2021a/TNFR2.SP.18m.S20/ALL"
),
load.to.enviroment = T, combine.objects = T)
#Integrate samples into a single VGM object
VGM <- VDJ_GEX_matrix(Data.in = list(agrafiotis2021a_TNFR2.BM.3m.S1_VDJGEXdata,
agrafiotis2021a_TNFR2.BM.12m.S6_VDJGEXdata,
agrafiotis2021a_TNFR2.BM.18m.S9_VDJGEXdata,
agrafiotis2021a_TNFR2.SP.3m.S12_VDJGEXdata,
agrafiotis2021a_TNFR2.SP.12m.S17_VDJGEXdata,
agrafiotis2021a_TNFR2.SP.18m.S20_VDJGEXdata
),
parallel.processing = "mclapply", trim.and.align = T)The output of the VDJ_GEX_matrix is the main object for all downstream functions in Platypus, as well as Steropodon. We can create this directly into the R session by using the public data available on PlatypusDB. Using PlatypusDB_fetch, we were able to download only a subset of our dataset.
The dataset we will be focusing on throughout this vignette comprises 22 total mouse single-cell repertoires (with both sequenced receptors and transcriptome information) following immunizations with the tumor necrosis factor 2 (TNFR2) - https://www.biorxiv.org/content/10.1101/2021.11.09.467876v2: 5 bone-marrow samples from 3-month old mice (s1-s5), 3 bone-marrow samples from 12-month old mice (s6-s8), 3 bone-marrow samples from 18-month old mice (s9-s11), 5 spleen samples from 3-month old mice (s12-s16), 3 spleen samples from 12-month old mice (s17-s19), and 3 spleen samples from 18-month old mice (s17-s19). For a faster computation time and to potentially observe a greater structural diversity, we have decided to select a single sample from each group (s1, s6, s9, s12, s7, and s20 were selected) when creating our VGM object.
#Annotate the VDJ and GEX objects with the site and age information
VGM[[1]]$site <- NA
VGM[[1]]$age <- NA
VGM[[1]]$site[VGM[[1]]$sample_id %in% c('s1', 's2', 's3')] <- 'BM'
VGM[[1]]$site[VGM[[1]]$sample_id %in% c('s4', 's5', 's6')] <- 'SP'
VGM[[1]]$age[VGM[[1]]$sample_id %in% c('s1', 's4')] <- '3 months'
VGM[[1]]$age[VGM[[1]]$sample_id %in% c('s2', 's5')] <- '12 months'
VGM[[1]]$age[VGM[[1]]$sample_id %in% c('s3', 's6')] <- '18 months'
VGM[[2]]@meta.data$site <- NA
VGM[[2]]@meta.data$age <- NA
VGM[[2]]@meta.data$site[VGM[[2]]@meta.data$sample_id %in% c('s1', 's2', 's3')] <- 'BM'
VGM[[2]]@meta.data$site[VGM[[2]]@meta.data$sample_id %in% c('s4', 's5', 's6')] <- 'SP'
VGM[[2]]@meta.data$age[VGM[[2]]@meta.data$sample_id %in% c('s1', 's4')] <- '3 months'
VGM[[2]]@meta.data$age[VGM[[2]]@meta.data$sample_id %in% c('s2', 's5')] <- '12 months'
VGM[[2]]@meta.data$age[VGM[[2]]@meta.data$sample_id %in% c('s3', 's6')] <- '18 months'We have annotated the VDJ and GEX objects (parts of the VGM object, see the Platypus introductory vignette!) with the site and age information, which will be used in the downstream plots and analyses.
Lastly, we will check whether our annotations are viable in both the VDJ and GEX, by plotting the richness of CDR3 sequences per age and the UMAP projection per site, for each object respectively.
#Check the VDJ annotations
Platypus::VDJ_diversity(VGM[[1]], feature.columns = c('VDJ_cdr3s_aa', 'VJ_cdr3s_aa'), grouping.column = 'age', metric = 'richness')## Used group 3 months as a reference for subsampling with 5193 entries## Subsampled group 12 months to 5193 entries## Subsampled group 18 months to 5193 entries
#Check the GEX annotations
Seurat::DimPlot(VGM[[2]], group.by = 'site')
Similarly, the TM-score also requires the common residues between two
structures, therefore seq.struct.superpose must be set to TRUE.
The
weighted RMSD sum by region offers a greater flexibility and can
incorporate some meaningful B-cell biology priors, such as the greater
structural and sequence diversity of the CDR regions. The
custom.specific.values defines which regions will be used from the
entire structure, with the column names defining those regions in our
PDB object being specified in custom.features. The custom.weights
parameter is a vector of weights for each region in the
custom.specific.values - in this case, we will assign a greater weight
to the heavy chain CDR3 region. Lastly, each region’s RMSD distance can
be normalized by its number of residues if
custom.normalize.sequence.length is set to T.
We
can also trim our structures beforehand (using Steropodon_trim) and then
calculate the RMSD - here we select only the CDR regions of both the
heavy and light chains! In this case, we have kept only the CDR regions
of both heavy and light chains - this increases the overall RMSD as
opposed to the previous RMSD heatmap on whole structure.
If
plot.results is set to F, Steropodon_distances will output the symmetric
distance matrix instead of a heatmap. This distance matrix can be used
for multidimensional scaling in Steropodon_distances_MDS, or for
hierarchical and k-medoids clustering in Steropodon_cluster.
When
considering MDS distances across the entire structure, we can see how
performing MDS is not very informative, with no specific clusters
forming across our structures.
We have
also done the same analysis as above, but for the aligned coordinate
vectors of the CDRH3 structures. However, these vectors are less
informative, mostly due to many coordinates being removed when
performing sequence alignment on the highly variable CDRH3 regions (as
the pdbs object keeps only the C-alpha coordinates of the aligned
residues).
However,
when using UMAP for dimensionality reduction and 2D projection, we can
notice some distinct clusters forming, with some novel clusters also for
each sampling site (first plot).
We have
done the same for the CDRH3 structures - some clusters appear, yet they
are not correlated with our features of interested (sample site and age
specifically).
Again, we
can see that the obtained clusters are mostly correlated with the V gene
of each structure, as expected. We can see no specific clustering across
sample sites or age for our TNFR2-immunized mice.
This plot
illustrates the best clustering - with two distinct clusters forming
depending on the V gene used. However, our other 2 features of interest
are still not correlated with the overall receptor structure, especially
with the structures of only the CDR regions.
We can see
how the multimodal integration refines the overall cell embeddings when
considering all the features highlighted (site, age, sample,
transcriptomic clusters, and V genes). Surprisingly, it appears that the
Seurat/transcriptomic clusters are more clearly defined in the
multimodal embedding, and that structural clustering on the coordinate
vectors is also correlated with transcriptome identity (e.g., clusters
10 and 11 clearly separating from the largest mass of points when
comparing the right side plot to the left side, colored by the Seurat
cluster identity obtained for the gene expression of all 22 TNFR2
samples). This illustrates a tight correlation between the structural
features of a receptor and the gene expression patterns of its
corresponding immune cell.
The
‘cluster’ slots contains the cluster ID for each structure.
We
might also be interested in investigating the CDR3 sequence diversity
across structural clusters, using the VDJ_diversity function from the
Platpyus package.
##
10. Calculating physicochemical properties via Steropodon_properties The
Steropodon_properties function can be used to calculate physicochemical
and developability parameters such as the solvent-accessible surface
area (SASA), charge, hydrophobicity, pKa, as well as obtain secondary
structure annotations. The resulting Steropodon object will have a
modified bio3d pdb object, with new columns corresponding to the
aforementioned properties. More calculators/properties will be
added!
We
can get proportions for the residue counts above by setting proportions
to T (proportions of variable regions across each group = across each
chain).
We
can also get the proportions of residues per region across other
grouping factors obtained via Steropodon, such as the core column from
Steropodon_find_core.
If
the feature parameter is given a column/feature name for continuous
features, then it will calculate the mean of that feature across the
values determined by the grouping parameter. For example, we can see the
B-value averages across the variable regions, pooled for all structures
in the nested list of Steropodon objects (plot.level = ‘global’).
To
get the mean of B-values across regions pooled across each sample, we
can set plot.level to ‘sample’.
If
plot.level is not set to ‘global’, we can obtain a list of plots per
sample/ samples and clonotypes.
Lastly, if
combine.lineplot is set to T, the line plots from above are combined
into a single plot and colored by the factor in plot.level (sample in
this case).
Same as above, but for B-values.
##
12. Docking structures via Steropodon_dock We can obtain receptor-ligand
complexes via docking (with ZDOCK) using Steropodon_dock. This requires
a single Steropodon object or a nested list as obtained from the
Steropodon_model function, as well as the PDB ID or PDB file of the
ligand we wish to dock. In this case, we will try docking the
experimental structure of the tumor necrosis factor used in immunizing
our mice - PDB ID = 2tnf.