PlatypusV3 Structure Vignette

2.2.1 Selection of top clonotypes with the VDJ_select_clonotypes function.

Prior to the structure prediction, desired antibody sequences of the most expanded clonotypes can be selected with the VDJ_select_clonotypes function. Furthermore, this function will call MIXCR to determine the germline reference sequences, which will be used to annotate the predicted sequences. Therefore a copy of MIXCR must be installed locally on your computer and the path has to be specified in the mixcr.direcotry argument. In order to integrate the UMI counts to the VGM data, this function needs the VGM and the raw data from the PlatypusDB_fetch function as an input. The function selects the most expanded clonotypes per sample by default (overall most expanded clonotypes can be selected by setting clonotypes.per.sample to FALSE). Among the clonotypes, the entries are ranked according to UMI counts, where the most prominent unique VDJ/VJ sequences are selected per clonotype. The number of clonotypes can be specified in the top.clontypes argument and the number of sequences per clonotype in the seq.per.clonotype argument.

Here the top three most expanded clonotypes per sample with the two most prominent sequences per clonotype are selected.

VDJ_top_clonotypes_mixcr_out <- VDJ_select_clonotypes(VGM = vgm, 
                                                      raw.data = neumeier2021_raw, 
                                                      clone.strategy = "10x.default", 
                                                      top.clonotypes = 3, 
                                                      seq.per.clonotype = 2, 
                                                      donut.plot = F, 
                                                      mixcr.directory = "/usr/local/Cellar/mixcr/3.0.13-2", 
                                                      species = "mmu")

## MAC system detected

## Clonotyping strategy: 10x.default

## Processing sample 1 of 5

## Processing sample 2 of 5

## Processing sample 3 of 5

## Processing sample 4 of 5

## Processing sample 5 of 5

## Backing up 10x default clonotyping in columns clonotype_id_10x and clonotype_frequency_10x before updating clonotype_id and clonotype_frequency columns

kable(VDJ_top_clonotypes_mixcr_out %>% head(),caption = "MIXCR output")

2.2.2 Custom selection of VGM entries

Instead of using the VDJ_select_clonotypes function, sequences can also be selected based on custom filters from the vgm object. The germline reference sequences can then be added by calling the VDJ_call_MIXCR_full function. Here only unique reads from sample s1 with the clonotype ID “clonotype4” are considered.

vgm_selection <- vgm[[1]] %>%  
                 filter(sample_id == "s1" & clonotype_id_10x == "clonotype4") %>%  
                 distinct(VDJ_raw_ref, .keep_all = T)

VDJ_mixcr_out <- VDJ_call_MIXCR_full(VDJ = vgm_selection, 
                                     mixcr.directory = "/usr/local/Cellar/mixcr/3.0.13-2", 
                                     species = "mmu")

## MAC system detected

2.3.1 Structure prediction from the vgm MIXCR output on Euler

The MIXCR output can be used as an input for the AlphaFold_prediction function. For ETH members with access to the Euler cluster and GPU access, the function will connect to the Euler cluster and send the AlphaFold structure predictions to the queue. The function will therefore ask for your password, which is handled in a safe manner by the shh package of R.

In the VDJ.mixcr.out argument the MIXCR output data frame either from the VDJ_select_clonotypes function or from calling the VDJ_call_MIXCR_full function can be specified as an input. In the cells.to.predict argument one can additionally specify a list of barcodes from which cells the structure shall be predicted. If it is set to “ALL”, all the structures from of the input data frame are predicted.

In the euler.user.name argument the ETH username for the Euler cluster has to be specified, which is the regular ETH username. In order to connect to server one has either to be connected to the ETH network or use a VPN connection. The function will then create a folder on your computer with all the FASTA files, which are then copied to the Euler Cluster in your scratch directory. By default the directory will be named AlphaFold_Fasta, however by specifying the argument dir.name the directory can be given any name.

Here we use the output from the VDJ_select_clonotypes function and predict the structure of three cells specified in cells.to.predict

#AlphaFold structure prediction
AlphaFold_prediction(VDJ.mixcr.out = VDJ_top_clonotypes_mixcr_out, 
                     cells.to.predict = c("s1_AAACCTGCATTGGCGC-1", 
                                          "s2_AAAGATGAGCATCATC-1",
                                          "s3_AAAGATGTCTCTAGGA-1"),
                     euler.user.name = "lucstalder", 
                     dir.name = "top_clonotypes_OVA")

This generates a single AlphaFold prediction job for each structure. Depending on the size and the number of structures this will take several hours to complete. Once the prediction is finished, the structures can be imported using the same function.
In the VDJ.mixcr.out argument the same input is specified as used previously for the prediction. In this case the VDJ_top_clonotypes_mixcr_out data frame. In order to import data, the import argument must be specified. Here the structures are imported directly from the Euler cluster and therefore import = “euler” is used.
As the directory was named “top_clonotypes_OVA” previously, the euler.dirname argument must be specified to tell the function in which directory the predicted outputs can be found.
AlphaFold produces multiple outputs for each structure which are then ranked according to the overall confidence, where the prediction with the highest confidence is named “ranked_0”. In the import function the number of ranked predictions per structure can be specified in the n.ranked argument, where n.ranked = 5 will import the top five ranked predictions per structure. The structures will be copied from the server to a directory on your local computer before they are loaded to the environment. By default, the local copy will be deleted afterwards. Here we want to keep them and therefore the rm.local.output is set to FALSE.

#AlphaFold import
VDJ_top_clonotypes_structure <- AlphaFold_prediction(VDJ.mixcr.out = VDJ_top_clonotypes_mixcr_out, 
                                                     import = "euler",
                                                     euler.user.name = "lucstalder",
                                                     euler.dirname = "top_clonotypes_OVA",
                                                     n.ranked = 5, 
                                                     rm.local.output = F)

The functions output is a list object which has in its first element the MIXCR output data frame that was used as input for structure prediction and in its second element a list of the predicted structures. As there are multiple ranked predictions per structure, the second list element of the output contains a list for each structure with the ranked predictions.

The second element contains a list for every structure with the ranked models. So here the structure from cell “s1_AAACCTGCATTGGCGC-1” is selected and from all predictions the one with the highest confidence “ranked_0.pdb” is chosen. This is a pdb object with the following attributes: atom, xyz, calpha, call

VDJ_top_clonotypes_structure[[2]]$`s1_AAACCTGCATTGGCGC-1`$ranked_0.pdb

## 
##  Call:  bio3d::read.pdb(file = paste0(CurDir, "/", OutDir, "/", out[i], 
##     "_ranked/", "ranked_", n - 1, ".pdb"))
## 
##    Total Models#: 1
##      Total Atoms#: 3419,  XYZs#: 10257  Chains#: 2  (values: A B)
## 
##      Protein Atoms#: 3419  (residues/Calpha atoms#: 225)
##      Nucleic acid Atoms#: 0  (residues/phosphate atoms#: 0)
## 
##      Non-protein/nucleic Atoms#: 0  (residues: 0)
##      Non-protein/nucleic resid values: [ none ]
## 
##    Protein sequence:
##       EVQLLETGGGLVQPGGSRGLSCEGSGFTFSGFWMSWVRQTPGKTLEWIGDINSDGSAINY
##       APSIKDRFTIFRDNDKSTLYLQMSNVRSEDTATYFCMRYSAYWYFDVWGAGTTVTVSSDI
##       KMTQSPSSMYASLGERVTITCKASQDIKSYLSWYQQKPWKSPKTLIYYATSLADGVPSRF
##       SGSGSGQDYSLTISSLESDDTATYYCLQHGESPYTFGGGTKLEIK
## 
## + attr: atom, xyz, calpha, call

2.3.2 Structure prediction from local FASTA files

The AlphaFold_prediction function can also be used to predict a protein of interest form a local FASTA file. For this purpose the argument fasta.directory.path can be specified with a path to a directory containing the FASTA files. The function is written in a way that the structure is predicted for all files in the specified directory with the ending .fasta. Like this a lot of structures can be predicted with calling the function only once and without adding a separate path for every file. Here we are only interested in the structure of ovalbumin, so the fasta.directory.path points to a directory which only contains the FASTA file of ovalbumin.

#AlphaFold structure prediction from local fasta file
AlphaFold_prediction(fasta.directory.path = "OVA", 
                     euler.user.name = "lucstalder", 
                     dir.name = "OVA_only")

The structure can then be imported with the same function as for prediction. If there is no vgm object specified, the output will just be a list of structures with a sublist for each prediction containing the top ranked pdb files.

OVA_structure <- AlphaFold_prediction(import = "euler", 
                                      euler.user.name = "lucstalder", 
                                      euler.dirname = "OVA", 
                                      n.ranked = 5)

OVA_structure[[1]]$OVA$ranked_0.pdb

## 
##  Call:  bio3d::read.pdb(file = paste0(CurDir, "/", OutDir, "/", out[i], 
##     "_ranked/", "ranked_", n - 1, ".pdb"))
## 
##    Total Models#: 1
##      Total Atoms#: 5999,  XYZs#: 17997  Chains#: 1  (values: A)
## 
##      Protein Atoms#: 5999  (residues/Calpha atoms#: 386)
##      Nucleic acid Atoms#: 0  (residues/phosphate atoms#: 0)
## 
##      Non-protein/nucleic Atoms#: 0  (residues: 0)
##      Non-protein/nucleic resid values: [ none ]
## 
##    Protein sequence:
##       MGSIGAASMEFCFDVFKELKVHHANENIFYCPIAIMSALAMVYLGAKDSTRTQINKVVRF
##       DKLPGFGDSIEAQCGTSVNVHSSLRDILNQITKPNDVYSFSLASRLYAEERYPILPEYLQ
##       CVKELYRGGLEPINFQTAADQARELINSWVESQTNGIIRNVLQPSSVDSQTAMVLVNAIV
##       FKGLWEKAFKDEDTQAMPFRVTEQESKPVQMMYQIGLFRVASMAS...<cut>...CVSP
## 
## + attr: atom, xyz, calpha, call

2.3.3 Structure prediction with Antigen interaction

It might be of interest not only to get the structure of the antibody but also to see AlphaFold’s prediction of antigen interaction. Therefore the antigen can be added to the prediction as an additional domain.
This can be easily done by specifying the path to the antigen FASTA file in the antigen.fasta.path argument of the function.

Here we are interested to see the prediction of some experimentally identified binders to the OVA protein. So we first filter the vgm object for the given samples

VDJ_OVA_binder <- vgm[[1]] %>% 
                  filter(barcode %in% c("s1_AAACCTGCACGTCTCT-1", 
                                        "s1_AACCGCGGTTCAGTAC-1", 
                                        "s2_TTTGCGCTCGGCTACG-1", 
                                        "s2_AACTGGTAGTGCAAGC-1", 
                                        "s5_ATTACTCAGGCTCTTA-1"))

For prediction we need to get the full length amino acid HC and LC sequence with MIXCR

VDJ_OVA_binder_mixcr_out <- VDJ_call_MIXCR_full(VDJ = VDJ_OVA_binder, 
                                                mixcr.directory = "/usr/local/Cellar/mixcr/3.0.13-2", 
                                                species = "mmu")

## MAC system detected

Now the structure of these antibodies can be predicted with the OVA antigen as an additional domain. The only thing to be added is the path to the antigen amino acid FASTA file in the antigen.fasta.path argument of the function.

#AlphaFold structure prediction
AlphaFold_prediction(VDJ.mixcr.out = VDJ_OVA_binder_mixcr_out, 
                     cells.to.predict = "ALL", 
                     euler.user.name = "lucstalder", 
                     dir.name = "OVA_binder_with_OVA", 
                     antigen.fasta.path = "OVA.fasta")

After prediction is finished the data is imported from the server:

#AlphaFold  import
OVA_binder_with_OVA_structure <- AlphaFold_prediction(VDJ.mixcr.out = VDJ_OVA_binder_mixcr_out, 
                                                      import = "euler", 
                                                      euler.user.name = "lucstalder", 
                                                      euler.dirname = "OVA_binder_with_OVA", 
                                                      n.ranked = 5, 
                                                      rm.local.output = F)

2.3.1 Antibody visualisation

The predicted structures can be visualized and analysed in R using a variety of packages for protein structure analysis. For an easy visualization, annotation and analysis of antibody structures, the Platypus VDJ_structure_analysis function can be used. It is manly based on the bio3d and r3dmol packages. The VDJ_structure_analysis function visualizes the antibody structure and automatically annotates the framework and CDR regions.

Lets have a look at the three structures of the top most expanded clonotypes that were previously predicted in 2.3.1.

out <- VDJ_structure_analysis(VDJ.structure = VDJ_top_clonotypes_structure, 
                              cells.to.vis = "ALL")

The function automatically annotates the structure based on the VDJ MIXCR sequences. The annotation can be disabled by setting the VDJ.anno argument to FALSE.

out <- VDJ_structure_analysis(VDJ.structure = VDJ_top_clonotypes_structure, 
                       cells.to.vis = "s1_AAACCTGCATTGGCGC-1", 
                       VDJ.anno = F)

The function has a variety of options which can be seen in the description. The annotation label can for example be disabled, or the label size can be changed, the color of the different regions can be set individually and so on.

The framework color can be set to gray, for a better visualization of the CDR’s.

out <- VDJ_structure_analysis(VDJ.structure = VDJ_top_clonotypes_structure, 
                       cells.to.vis = "s1_AAACCTGCATTGGCGC-1", 
                       label = F, 
                       color.frameworks = "#C0C0C0", 
                       angle.x = 0, 
                       angle.z = -10, 
                       angle.y = 170)

The function can also be used to visualize structures directly from pdb files. Previously the structure of ovalbumin was predicted with AlphaFold and can be visualized now by the PDB.file option

out <- VDJ_structure_analysis(PDB.file = OVA_structure[[1]]$OVA$ranked_0.pdb)

If you would like to visualize a local .pdb file then you have to use the function bio3d::read.pdb() function to read in the PDB file as follows:

VDJ_structure_analysis(PDB.file = bio3d::read.pdb("Path/to/file.pdb"))

2.3.2 Binding Interaction

For the OVA binders the structure was predicted together with the antigen. The Platypus VDJ_structure_analysis function can be used to visualize the binding interaction and determine some binding site metrics.

First let’s have a glimpse at the structures

out <- VDJ_structure_analysis(VDJ.structure = OVA_binder_with_OVA_structure, 
                       cells.to.vis = "s1_AACCGCGGTTCAGTAC-1", 
                       angle.x = 40, 
                       angle.y = 200, 
                       angle.z = -50, 
                       label = F)

Experimentally there were some distinct epitops determined on the ovalbumin antigen. With the anno.seq argument of the VDJ_plot_structure function there is the possibility to annotate custom regions on the sequence. It can either be a list of multiple regions to annotate or a single region. For every annotation there are four mandatory elements: The residue index of the annotation’s start, the residue index of the annotation’s end, the protein chain (this is alphabetic A:Z and in the same order as the sequences in the FASTA file) and the color of the annotation. There is an optional fifth element if you would like to have a label for the annotation.

For the OVA there are three bins to be annotated so a list of annotations is specified. First the OVA structure is set to gray, then annotations for Bin1 to Bin3 are specified with the respective color and the label. Bin3 consists of two regions so there is a second annotation with the same color but no label.

OVA_Bins <- list(
  c(1,392,"C","#b3b3b3"),
  c(186,201,"C","#9900ff","Bin1"),
  c(133,148,"C","#ffcc00","Bin2"),
  c(265,280,"C","#0000ff","Bin3"),
  c(214,229,"C","#0000ff")
)

Now the created custom annotation can be added by the anno.seq argument.

out <- VDJ_structure_analysis(VDJ.structure = OVA_binder_with_OVA_structure, 
                       cells.to.vis = "s1_AACCGCGGTTCAGTAC-1", 
                       angle.x = 40, 
                       angle.y = 200, 
                       angle.z = -50, 
                       label = F, 
                       anno.seq = OVA_Bins)

The three bins are now visible and the antibody seems to interact with the antigen somewhere close to Bin3. To get some more information about the binding site there is the antigen.interaction option. This will determine the binding site residues based on the bio3d::binding.site function and then calculate some metrics of the binding site residues which are summarized in a data frame.

By setting the BindingResidues.plot argument to TRUE, a plot is shown where the binding residues are colored on the antigen and on the antibody.

out <- VDJ_structure_analysis(VDJ.structure = OVA_binder_with_OVA_structure, 
                       cells.to.vis = "s1_AACCGCGGTTCAGTAC-1", 
                       structure.plot = F, 
                       anno.seq = OVA_Bins, 
                       antigen.interaction = T, 
                       BindingResidues.plot = T, 
                       angle.x = 40, 
                       angle.y = 200, 
                       angle.z = -50)

Alpha Fold has a confidence score for every residue in form a pLDDT value. This is a scale from 1 to 100 where a higher score means more confidence in the model prediction. The per residue pLDDT scores can also be visualized with the VDJ_structure_analysis function by setting the plddt.plot argument to TRUE. In this plot blue means high confidence and red low confidence. It can be seen that the binding residues are normally of lower confidence compared to the rest of the structure.

out <- VDJ_structure_analysis(VDJ.structure = OVA_binder_with_OVA_structure, 
                       cells.to.vis = "s1_AACCGCGGTTCAGTAC-1", 
                       structure.plot = F, 
                       anno.seq = OVA_Bins, 
                       antigen.interaction = T, 
                       plddt.plot = T, 
                       angle.x = 40, 
                       angle.y = 200, 
                       angle.z = -50)

Now let’s have a look at the binding site metrics data frame that is produced when the antigen.interaction label is set to TRUE. The first two columns are just the cellular barcode and the rank of the model that was analyzed. Then the mean confidence (pLDDT) is summarized for the non binding residues over the whole structure as well as for the binding residues separately for heavy chain, light chain and antigen.
Then the mean minimal distance is shown in the Mean_dist_bind_resi column. For every binding residue on the antibody the distance to its closest partner on the antigen is calculated and then the mean over all binding residues is calculated. So it is the mean of all the minimal distances.
In the last three columns the residue indexes for all the binding sites are summarized as a semi colon separated list.

out <- VDJ_structure_analysis(VDJ.structure = OVA_binder_with_OVA_structure, 
                       cells.to.vis = "s1_AACCGCGGTTCAGTAC-1", 
                       structure.plot = F, 
                       anno.seq = OVA_Bins, 
                       antigen.interaction = T)

It might be of interest to see to which region of the antibody the binding site residues belong to. It would be expected that manly the CDR’s are involved in binding and most prominently the CDR3. The Platypus VDJ_structure_analysis function has an argument binding.residue.barplot, which will plot the distribution of the biding site residues to the regions of the antigen in a bar plot.

out <- VDJ_structure_analysis(VDJ.structure = OVA_binder_with_OVA_structure, 
                       cells.to.vis = "s1_AACCGCGGTTCAGTAC-1", 
                       structure.plot = F, 
                       anno.seq = OVA_Bins, 
                       antigen.interaction = T, 
                       binding.residue.barplot = T)

A bar plot is produced for every structure that is analyzed with the absolute counts of binding residues for each region of the antibody. Furthermore, a summary bar plot is produced where a percentage of binding site residues is given for each region. This plot is the average over all the analyzed structures. In this example, only one structure was analyzed, so the summary plot is only based on this one structure. The data of for this summary baplot is also returned as a dataframe.

For this vignette we might be interested in the summary plot of all the predicted structures

out <- VDJ_structure_analysis(VDJ.structure = OVA_binder_with_OVA_structure, 
                              cells.to.vis = "ALL", 
                              structure.plot = F, 
                              anno.seq = OVA_Bins, 
                              antigen.interaction = T, 
                              binding.residue.barplot = T)

It can be nicely seen that for the heavy chain most of the binding residues are actually part of the CDR3 and in the light chain they are part of the CDR’s as well.

The barplot can also be shown in a simplyfied version where only the distribution between frame work vs CDR is sown. For this the binding.residue.barplot.style argument can be set to “FR_CDR”

out <- VDJ_structure_analysis(VDJ.structure = OVA_binder_with_OVA_structure, 
                              cells.to.vis = "ALL", 
                              structure.plot = F, 
                              anno.seq = OVA_Bins, 
                              antigen.interaction = T, 
                              binding.residue.barplot = T, 
                              binding.residue.barplot.style = "FR_CDR")

## Not for every entry a structure is defined! All the defined structures are analysed.

Here it can be seen that for the analyzed models in this vignette 80% of all the binding residues are actually part of the CDR regions. Which is quite astonishing to see AlphaFold modeling the binding of an antibody with the CDRs as main binding interaction partners.