AntibodyForests vignette: building and analyzing B-cell lineage trees from 10x sc-V(D)J seq data

The VDJ_build() function imports Cell Ranger output into an R dataframe, priming the data for downstream analyses. It is optimized for Cell Ranger v7 and older versions and is suitable for B and T cell repertoires. Upon execution, the function generates a VDJ dataframe where VDJ and VJ transcripts are paired (if present) in each row, with each row representing a single cell. Post VDJ dataframe construction, a table is displayed on the console, delineating the number of filtered cells, contingent upon the remove.divergent.cells and complete.cells.only parameters.

For seamless integration with functions such as the Af_build() function (discussed next), it is strongly recommended to set the remove.divergent.cells and complete.cells.only parameters to TRUE. Cells with more than one VDJ or VJ transcript, referred to as ‘divergent cells’, typically represent doublets that may arise during the single-cell sorting, or show dual expression of kappa and lambda light chains. Such occurrences are rare, accounting for only 0.2-0.5% of peripheral blood B cells, as reported by Giachino et al. (1995). Cells harboring both a single VDJ and a VJ transcript may be included for comprehensive downstream analyses.

The output of 10x Genomics contains one output folder per run/sample. In each sample folder, the following files should be present in order to run the VDJ_build() function successfully.

• filtered_contig_annotations.csv
  This CSV file contains the amino acid and nucleotide framework (FWRs) and complementarity determining regions (CDRs) of the contig sequences stored in the filtered_contig.fasta file. These annotations are required to obtain the trimmed V(D)J sequences.
  
• filtered_contig.fasta
  This FASTA file contains the contig sequences from barcodes that passed the filtering steps. These contigs are annotated in the filtered_contig_annotations.csv.
  
• consensus_annotations.csv
  This CSV file contains the amino acid and nucleotide framework (FWRs) and complementarity determining regions (CDRs) of the consensus sequences stored in the consensus.fasta file.
  
• consensus.fasta
  This FASTA file contains the consensus sequence of each assembled contig, which is identical to the sequence of the most frequent exact subclonotype, in which an exact subclonotype is defined as a group of cells with identical V(D)J transcripts.
  
• concat_ref.fasta
  This FASTA file contains the concatenated V(D)J reference segments for the segments detected on each consensus sequence. Note, these reference sequences serve as an approximate reference for each consensus sequence, and need to be aligned and trimmed in order to have the same length as the consensus sequence (see trim.and.align parameter).
  
  NOTE: In older versions of Cell Ranger (older than v6), the FWR1-4 and CDR1-3 regions are not annotated. In this scenario, the all_contig_annotations.json file will be employed to obtain the ‘L-REGION+V-REGION’ and ‘J-REGION’ annotations. The raw sequences can now be trimmed using the ‘L-REGION+V-REGION’ start and ‘J-REGION’ end indices. If the all_contig_annotations.json file is not found in the sample/run directory, the columns with the trimmed sequences will remain empty and a warning message is printed.
  
  

VDJ.directory

• Description: Path to the parent directory containing the output folders (one folder for each run/sample) of Cell Ranger (note: when specified, the VDJ.sample.list parameter should not be specified).
  
• Example: "/path/to/VDJ/directory"
  
• Default: None
  
  

VDJ.sample.list

• Description: Character vector of the paths to the output folders (one folder for each run/sample) of Cell Ranger (note: when specified, the VDJ.directory parameter should not be specified).
  
• Example: c("/path/to/sample1", "/path/to/sample2")
  
• Default: None
  
  

remove.divergent.cells

• Description: If TRUE, cells with more than one VDJ transcript or more than one VJ transcript will be excluded from the VDJ dataframe.
  
• Default: FALSE
  
  

complete.cells.only

• Description: If TRUE, only cells with both a VDJ transcript and a VJ transcript are included in the VDJ dataframe.
  
• Default: FALSE
  
  

trim.germlines

• Description: If TRUE, the raw germline sequence of each clone will be aligned with the consensus sequence of that clone serving as the reference sequence, utilizing the Biostrings::pairwiseAlignment() function. The alignment conducted will be of the ‘global-local’ type, therebyy trimming the raw germline sequence.
  
• Default: FALSE
  
  

fill.germline.CDR3

• Description: If TRUE, the trimmed germline sequence of each clone will be aligned with the most frequently observed sequence for the CDR3 region of that clone serving as the reference sequence, utilizing the Biostrings::pairwiseAlignment() function. The alignment conducted will be of the ‘global-local’ type, thereby selecting the CDR3 region from the trimmed germline sequence. After the alignment, the germline CDR3 sequence is replaced by the most frequently observed CDR3 sequence, in order to obtain germline sequences that are more likely to encode producible and productive antibodies.
  
• Default: FALSE
  
  

gap.opening.cost

• Description: The cost for opening a gap when aligning and trimming germline sequences using the consensus sequences. A cost of Inf will result in a gapless alignment.
  
• Default: 10
  
  

gap.extension.cost

• Description: The cost for extending a gap in when aligning and trimming germline sequences. A cost of Inf will result in a gapless alignment.
  
• Default: 4
  
  

In this example, divergent cells are filtered out, and only cells with both a VDJ and VJ transcript are included in the VDJ dataframe. The resulting VDJ dataframe contains 17,891 cells with both a single VDJ and VJ transcript. The number of cells filtered out due to the lack of VDJ and VJ transcripts, referred to as ‘incomplete cells’, is shown in the right column of the printed table. Since a cell can have a double VDJ or VJ transcript or may lack a VDJ or VJ transcript, such a cell can be included in both the counts of ‘divergent cells’ and ‘incomplete cells’. Note that, as a result, the sum of these counts may not exactly match the difference in the number of rows between the dataframes. Additionally, for each clone, the germline sequences are trimmed, and an additional germline sequence column is added in which the CDR3 region is replaced by the most frequently observed CDR3 sequence in that clone.

# Read in data, keep only complete, non-divergent cells and trim germline sequences
VDJ_OVA <- Platypus::VDJ_build(VDJ.directory = "10x_output/VDJ/",
                     remove.divergent.cells = TRUE,
                     complete.cells.only = TRUE,
                     trim.germlines = TRUE,
                     fill.germline.CDR3 = TRUE)

For illustration, the dummy VDJ dataframe below contains a few columns of the first 10 rows of S1:

The Af_build() function infers B/T cell evolutionary networks for all clonotypes in a VDJ dataframe, thereby providing insights into the evolutionary relationships between BCR/TCR sequences from each clonotype. Upon execution, the function generates an AntibodyForests object containing lineage trees, sequences, and other specified features for downstream analyses. The resulting object can be used for visualization, comparison, and further analysis of B cell repertoires.

For all clonotypes, the unique (combination of) sequences are selected, which will be nodes in the lineage tree. Until now, five algorithms have been implemented to construct B cell lineage trees from this list of sequences. Three construction algorithms are based on a pairwise distance matrix: ‘phylo.network.default’, ‘phylo.network.mst’, and ‘phylo.tree.nj’. If the string.dist.metric is specified, the distance matrices are calculated using the stringdist::stringdistmatrix() function; if the dna.model or aa.model is specified, the distance matrices are calculated using the ape::dist.dna() or phangorn::dist.ml() function, respectively. Two construction algorithms are based on a multiple sequence alignment (msa): ‘phylo.tree.mp’ and ‘phylo.tree.ml’. The msa is created with the msa::msa() function. The five algorithms are explained in the description of the construction.method parameter below.

In addition, there is also the option to import pre-constructed lineage trees from an IgPhyML output file into an AntibodyForests object, which enables to comparison of trees created with different construction methods on the same sequencing data in a repertoire-wide manner. IgPhyML is a command-line tool that infers maximum likelihood trees using a phylogenetic codon substitution model specific for antibody lineages, as explained by Hoehn et al. (2017). Some functions have been included to create the appropriate input files from IgPhyML from 10x Genomics outputs, by running the IgBLAST tool and integrating the output with the VDJ dataframe. These function are discussed in the supplementary of this vignette.

VDJ

• Description: VDJ dataframe as obtained from the VDJ_build() function or a dataframe that contain the columns specified in the sequence.columns, germline.columns, and node.features column.
  
• Example: VDJ_df
  
• Default: None
  
  

sequence.columns

• Description: Sequence column(s) in the VDJ dataframe that contain the sequences that will be used to infer B cell lineage trees. The sequences should be of the same type (DNA or protein), otherwise, an error message will be returned.
  
• Example: c("VDJ_sequence_aa_trimmed", "VJ_sequence_aa_trimmed")
  
• Default: c("VDJ_sequence_nt_trimmed", "VJ_sequence_nt_trimmed")
  
  

germline.columns

• Description: Germline column(s) in the VDJ dataframe that contain the sequences that will be the starting points (roots) of the trees. The germline columns should be in the corresponding order with the columns specified by the `sequence.columns’ parameter, and should be of the same type. Note: if the germline sequences are not trimmed using the consensus sequences, the distance between the nodes and the germline node will be overestimated.
  
• Example: c("VDJ_germline_aa_trimmed", "VJ_germline_aa_trimmed")
  
• Default: c("VDJ_germline_nt_trimmed", "VJ_germline_nt_trimmed")
  
  

concatenate.sequences

• Description: If TRUE, sequences from multiple sequence columns are concatenated into one sequence for single distance matrix calculations/multiple sequence alignments. If FALSE, a distance matrix is calculated/multiple sequence alignment is performed for each sequence column separately.
  
• Default: FALSE
  
  

node.features

• Description: Column name(s) in the VDJ dataframe which should be imported into the AntibodyForests object, i.e. for plotting of lineage trees later on.
  
• Example: c("VDJ_vgene", "VDJ_dgene", "VDJ_jgene", "isotype")
  
• Default: "isotype" (if present)
  
  

string.dist.metric

• Description: Specifies the metric used to calculate pairwise string distances between sequences with the stringdist::stringdistmatrix() function. This parameter is applicable only when a pairwise DNA or protein distance matrix will be calculated and when the dna.model and aa.model parameters are not specified. Options include:
  
  
  • "lv" (Levenshtein distance / edit distance)
    
  • "dl" (Damerau-Levenshtein distance)
    
  • "osa" (Optimal String Alignment distance)
    
  • "hamming" (Hamming distance)
    
  • "qgram" (Q-gram distance)
    
  • "cosine" (Cosine distance)
    
  • "jaccard" (Jaccard distance)
    
  • "jw" (Jaro-Winkler distance)
    
• Default: "lv"
  
  

dna.model

• Description: When a distance-based method is specified in the construction.method parameter, an evolutionary DNA model can be specified that is to be used during the pairewise distance calculation. Currently, the following models are included: ‘raw’, ‘N’, ‘TS’, ‘TV’, ‘JC69’, ‘K80’, ‘F81’, ‘K81’, ‘F84’, ‘BH87’, ‘T92’, ‘TN93’, ‘GG95’, ‘logdet’, ‘paralin’, ‘indel’, and ‘indelblock’. The pairwise-distance matrix is calculated using the ape::dist.dna() function. When the construction.method parameter is set to phylo.tree.ml, a nucleotide subsitution model can be specified that is to be use during the maximum likelihood tree-inference. Currently, the following nucleotide substitution models are included: ‘JC’, ‘F81’, ‘K80’, ‘HKY’, ‘TrNe’, ‘TrN’, ‘TPM1’, ‘K81’, ‘TPM1u’, ‘TPM2’, ‘TPM2u’, ‘TPM3’, ‘TPM3u’, ‘TIM1e’, ‘TIM1’, ‘TIM2e’, ‘TIM2’, ‘TIM3e’, ‘TIM3’, ‘TVMe’, ‘TVM’, ‘SYM’, and ‘GTR’. When set to "all", all available nucleotide substitution models will be tested using the phangorn::modelTest() function, after which the result is given to the phangorn::pml_bb() function for the maximum likelihood tree inference.
  
• Example: "all"
  
• Default: None (when a distance-based method is specified with the construction.methodparameter) and "all" (when the construction.method parameter is set to “phylo.tree.ml”).
  
  

aa.model

• Description: When a distance-based method is specified in the construction.method parameter, an evolutionary protein model can be specified that is to be used during the pairwise distance calculation. Currently, the following models are included: “WAG”, “JTT”, “LG”, “Dayhoff”, “VT”, “Dayhoff_DCMut” and “JTT-DCMut”. The pairwise-distance matrix is calculated using the phangorn::dist.ml() function. When the construction.method parameter is set to ‘phylo.tree.ml’, an amino acid subsitution model can be specified that is to be use during the maximum likelihood tree-inference. Currently, the following amino acid substitution models are included: “WAG”, “JTT”, “LG”, “Dayhoff”, “VT”, “Dayhoff_DCMut” and “JTT-DCMut”. When set to "all", all available nucleotide substitution models will be tested using the phangorn::modelTest() function, after which the result is given to the phangorn::pml_bb() function for the maximum likelihood tree inference.
  
• Example: "all"
  
• Default: None (when a distance-based method is specified with the construction.methodparameter) and "all" (when the construction.method parameter is set to “phylo.tree.ml”).
  
  

codon.model

• Description: Specifies the codon substitution model to be used during maximum likelihood tree infrence. This parameter can only be specified when the construction.method is set to “phylo.tree.ml” and when columns with nucleotide sequences are specified in the sequence.columns and germline.columns parameters. Currently, three codon substution models are available: ‘M0’ (which assumes a non-distinct non-synonymous/synonymous ratio (ω), ‘M1a’ (that estimates two different ω classes: ω = 1 & ω < 1), and ‘M2a’ (that estimates three different ω classes: ω < 1, ω = 1, positive selection ω > 1).
• Example: "M0"
  
• Default: None
  
  

construction.method

• Description: Denotes the approach and algorithm used to convert the distance matrix or multiple sequence alignment into a lineage tree. Currently, there are six options for this paramater:
  
  
  • "phylo.network.default"
    This approach employs a mst-like algorithm to construct a tree evolutionary network. In this method, the germline node is positioned at the top of the tree, and nodes with the minimum distance to any existing node in the tree are linked iteratively.
    
  • "phylo.network.mst"
    This method utilizes the minimum spanning tree (MST) algorithm from the ape::mst() function. It constructs networks with the minimum sum of edge lengths/weights by iteratively adding edges to the network in ascending order of edge weights, while ensuring that no cycles are formed. The network is then reorganized into a germline-rooted lineage tree.
    
  • "phylo.tree.nj"
    This approach employs the neighbor-joining (NJ) algorithm from the ape::nj() function. It constructs phylogenetic trees by joining pairs of nodes with the minimum distance, resulting in a bifurcating tree consisting of internal nodes (representing unrecovered sequences) and terminal nodes (representing the recovered sequences).
    
  • "phylo.tree.mp"
    This method utilizes the maximum-parsimony (MP) algorithm from the phangorn::pratchet() function. It constructs phylogenetic trees by minimizing the total number of edits required to explain the observed differences among sequences.
    
  • "phylo.tree.ml"
    This approach utilizes the maximum-likelihood (ML) algorithm from the phangorn::pml_bb() function. It constructs phylogenetic trees by estimating the tree topology and branch lengths that maximize the likelihood of observing the given sequence data under a specified evolutionary model.
    
  • "phylo.tree.IgPhyML"
    With this option, no trees/networks are inferred directly. Instead, trees are imported from an IgPhyML output file specified by the IgPhyML.output.file parameter.
    
• Default: "phylo.network.default" (if the IgPhyML.output.file parameter is not specified) and "phylo.tree.IgPhyML" (if the IgPhyML.output.file parameter is specified)
  
  

IgPhyML.output.file

• Description: Specifies the path to the IgPhyML output file, from which the trees will be imported if construction.method is not specified or set to "phylo.tree.IgPhyML".
• Example: "path/to/IgPhyML/output/file.tab"
  
• Default: None
  
  

resolve.ties

• Description: Denotes how ties are handled during the conversion of the distance matrix into lineage trees when construction.method is set to “phylo.network.default”. Ties occur when an unlinked node, which is to be linked to the tree next, shares identical distances with multiple previously linked nodes in the lineage tree. If a vector is provided, ties will be resolved in a hierarchical manner, following the order specified in the vector. If a tie could not be resolved, the node is connected to all nodes, thereby creating a cyclic structure, and a warning message is printed after creation of the AntibodyForests object. There are seven ways to handle ties:
  
  
  • "min.expansion"
    Selects the node(s) with the smallest size.
    
  • "max.expansion" Selects the node(s) with the biggest size.
    
  • "min.germline.dist"
    Selects the node(s) with the smallest string distance to the germline node.
    
  • "max.germline.dist" Selects the node(s) with the biggest string distance to the germline node.
    
  • "min.germline.edges"
    Selects the node(s) with the lowest possible number of edges to the germline node.
    
  • "max.germline.edges"
    Selects the node(s) with the highest possible number of edges to the germline node.
    
  • "min.descendants"
    Selects the node(s) with the lowest possible number of descendants.
    
  • "max.descendants"
    Selects the node(s) with the highest possible number of descendants.
    
  • "random" Selects a random node.
    
• Default: c("max.expansion", "close.germline.dist", "close.germline.edges", "random") (if the construction.method parameter is set to "phylo.network.default")
  
  

remove.internal.nodes

• Description: Denotes if and how internal nodes should be removed after constructing phylogenetic trees. There are four algorithms included to remove internal nodes:
  
  
  • "zero.length.edges.only"
    This option removes internal nodes that only have zero-length edges to terminal nodes.
    
  • "connect.to.parent"
    This option first removes internal nodes that have zero-length edges to terminal nodes. Subsequently, it connects all terminal nodes directly to the first parental sequence-recovered node found higher in the tree.
    
  • "minimum.length" This algorithm iteratively removes internal nodes by prioritizing edges with the minimum length for deletion.
    
  • "minimum.cost"
    Similar to the "minimum.length" option, this algorithm iteratively removes internal nodes. However, instead of prioritizing edges based solely on their length, it prioritizes edges that result in the minimum increase in the sum of all edge lengths.
• Default: "minimum.cost" (when the construction.method parameter is set to "phylo.tree.nj") and "connect.to.parent" (when the construction.method parameter is set to "phylo.tree.mp", "phylo.tree.ml", or "phylo.tree.IgPhyML")

The Af_plot_tree() function retrieves the igraph object from the provided AntibodyForests object for the specified clone within the specified sample and plots the lineage tree using the specified plotting parameters. This visualization offers insights into sequence evolution, enabling the observation of evolutionary processes in relation to sequence distances, clonal expansion, and other attributes such as isotype or evolutionary likelihood, as well as additional features appended to the VDJ dataframe.

AntibodyForests_object

• Description: A list of AntibodyForests objects as obtained from the Af_build() function.
  
• Example: list("Default" = AntibodyForests_OVA_default, "MST" = AntibodyForests_OVA_mst, "IgPhyML" = AntibodyForests_OVA_IgPhyML)
  
• Default: None
  
  

sample

• Description: The sample name for which the lineage tree should be plotted.
• Example: "S1"
• Default: None
  

clonotype

• Description: The clonotype name for which the lineage tree should be plotted.
• Example: "clonotype1"
• Default: None
  

label.by

• Description: Specifies the label to plot on the nodes. Options include:
  • "name"
    The node names are displayed on the nodes.
    
  • "size" The number of cells per node are displayed.
  • "none"
    No labels are displayed on the nodes.
• Default: "name"
  

color.by

• Description: Specifies the node feature to color the nodes.
• Example: "celltype"
• Default: "isotype" or NULL (if isotype is not present)
  

edge.label

• Description: Specifies the label to be plotted on the edges. Options include:
  • "original"
    The distance that is stored in the igraph object.
  • "none"
    No labels are displayed on the edges.

First, node labels need to be synchronized between the different AntibodyForests objects. This can be done using the Af_sync_nodes() function. Here, the nodes of the AntibodyForests objects created with the IgPhyML tool, MST, NJ, and ML algorithms are synchronized with the nodes of the AntibodyForests object created with the default algorithm.

# Synchronize the nodes of the different AntibodyForests objects
AntibodyForests_OVA_IgPhyML <- Af_sync_nodes(reference = AntibodyForests_OVA_default,
                                                    subject = AntibodyForests_OVA_IgPhyML)
AntibodyForests_OVA_mst <- Af_sync_nodes(reference = AntibodyForests_OVA_default,
                                subject = AntibodyForests_OVA_mst)
AntibodyForests_OVA_nj <- Af_sync_nodes(reference = AntibodyForests_OVA_default,
                                subject = AntibodyForests_OVA_nj)
AntibodyForests_OVA_mp <- Af_sync_nodes(reference = AntibodyForests_OVA_default,
                                subject = AntibodyForests_OVA_mp)
AntibodyForests_OVA_ml <- Af_sync_nodes(reference = AntibodyForests_OVA_default,
                                subject = AntibodyForests_OVA_ml)

Now we can calculate the distance between the trees of the different AntibodyForests objects using the Af_compare_methods() function using the Generalized Branch Length Distance.

# Calculate the euclidean distance between different trees based on the GBLD
out <- Af_compare_methods(input = list("Default" = AntibodyForests_OVA_default,
                                       "MST" = AntibodyForests_OVA_mst,
                                       "NJ" = AntibodyForests_OVA_nj,
                                       "MP" = AntibodyForests_OVA_mp,
                                       "ML" = AntibodyForests_OVA_ml, 
                                       "IgPhyML" = AntibodyForests_OVA_IgPhyML),
                          min.nodes = 10,
                          distance.method = "GBLD",
                          visualization.methods = "heatmap",
                          include.average = T)
# Plot the average heatmap
out$average$Heatmap

As an example, we calculate the metric dataframe of an AntibodyForests object for trees with at least 10 nodes. The rownames include the sample ID, followed by the clonotype ID.

# Calculate the metric dataframe for the AntibodyForests object
metric_df <- Af_metrics(input = AntibodyForests_OVA_default,
                        min.nodes = 10,
                        metrics = c("mean.depth", "sackin.index", "spectral.density"))

load("objects/metric_df.RData")

We also investigate if nodes from the IgA isotype are closer to the germline in these trees than nodes from the IgG isotype. We visualize the results in a boxplot and observe no significant difference.

# Cluster the trees in the repertoire based on the euclidean distance between various metrics of each tree
Af_distance_boxplot(AntibodyForests_object = AntibodyForests_OVA_default,
                    node.feature = "isotype",
                    distance = "node.depth",
                    min.nodes = 8,
                    groups = c("IgG1", "IgA"),
                    text.size = 16,
                    unconnected = T,
                    significance = T)

# Cluster the trees in the repertoire based on the euclidean distance between various metrics of each tree
out <- Af_compare_within_repertoires(input = AntibodyForests_OVA_default,
                                     min.nodes = 8,
                                     distance.method = "euclidean",
                                     distance.metrics = c("mean.depth", "sackin.index", "spectral.density"),
                                     clustering.method = "mediods",
                                     visualization.methods = "PCA")
# Plot the PCA colored on the assigned clusters
out$plots$PCA_clusters

# Analyze the difference between the clusters based on the calculated sackin.index
plots <- Af_cluster_metrics(input = AntibodyForests_OVA_default,
                   clusters = out$clustering,
                   metrics = c("sackin.index", "spectral.density"),
                   min.nodes = 8,
                   significance = T)

plots$sackin.index
plots$modalities

# Analyze the difference between the clusters based on the calculated sackin.index
plots <- Af_cluster_node_features(input = AntibodyForests_OVA_default,
                   clusters = out$clustering,
                   features = "isotype",
                   fill = "max",
                   significance = T)

plots$isotype

# Cluster the trees in the repertoire based on the Jensen-Shannon divergence between the Spectral Density profiles
out <- Af_compare_within_repertoires(input = AntibodyForests_OVA_default,
                                     min.nodes = 8,
                                     distance.method = "jensen-shannon",
                                     clustering.method = "mediods",
                                     visualization.methods = "heatmap")
# Plot the heatmap
out$plots$heatmap_clusters

# Analyze the difference between the clusters
plots <- Af_cluster_metrics(input = AntibodyForests_OVA_default,
                   clusters = out$clustering,
                   metrics = "spectral.density",
                   min.nodes = 8,
                   significance = T)

plots$spectral.asymmetry
plots$modalities

Pseudolikelihoods were calculated with the PLM-pipeline for the sequences in the VDJ dataframe. These pseudolikelihoods first have to be added to the AntibodyForests object. Then we plot the pseudolikelihood against the distance to the germline based on the edge.length.

# Add pseudolikelihood as node feature to the AntibodyForests object
AntibodyForests_OVA_default <- Af_add_node_feature(AntibodyForests_object = AntibodyForests_OVA_default, 
                                                                feature.df = VDJ_OVA,
                                                                feature.names = "pseudolikelihood")

# Plot pseudolikelihood against distance to the germline
Af_distance_scatterplot(AntibodyForests_object = AntibodyForests_OVA_default, 
                                     node.features = "pseudolikelihood",
                                     distance = "edge.length",
                                     min.nodes = 5,
                                     correlation = "pearson",
                                     color.by = "sample")

To investigate the PLM likelihoods of somatic hypermutation, we first need to extract sequences from the AntibodyForests object. These sequences can be used as input CSV file for the PLM-pipeline.

#Create a dataframe of the sequences with a sequence ID containing the sample, clonotype, and node number
df <- Af_get_sequences(AntibodyForests_object = AntibodyForests_OVA_default, 
                                               sequence.name = "VDJ_sequence_aa_trimmed", 
                                               min.edges = 1)
#Save the data frame as CSV to serve as input for the PLM-pipeline
write.csv(df, file = "/directory/PLM_input.csv", row.names = FALSE)

The generated probability matrices can be used to create a dataframe of the ranks and probabilities of the mutations along the edges of the lineage trees.

#Create a PLM dataframe
PLM_dataframe <- Af_PLM_dataframe(AntibodyForests_object = AntibodyForests_OVA_default, 
                                               sequence.name = "VDJ_sequence_aa_trimmed", 
                                               path_to_probabilities = "/directory/ProbabilityMatrix")

Below is an example of the first 20 rows of a PLM dataframe. Each row in this data frame represents an edge between node1 and node2 in the lineage tree of sample and clonotype. The number of mutations along this edge is given by n_subs. When there are multiple mutations, the average of the ranks and probabilities are reported.

The ranks and probabilities of the original and substitution redisues can be visualized in distribution plots.

# Plot the substitution rank
Af_plot_PLM(PLM_dataframe = PLM_dataframe, 
                         values = "substitution_rank", 
                         group_by = "sample_id")

Protein structures of the antibodies of one of the trees from the OVA dataset were determined with AlphaFold3, together with the Ovalbumin sequence. These PDB files were superimposed and used as input for AntibodyForests.

#Biophysical properties of the heavy and light chain sequence of the antibodies are calculated
vdj_structure_AbAg <- VDJ_3d_properties(VDJ = vdj,
                                   pdb.dir = "~/path/PDBS_superimposed/",
                                   file.df = files,
                                   properties = c("charge", 
                                                  "3di_germline", 
                                                  "hydrophobicity", 
                                                  "RMSD_germline", 
                                                  "pKa_shift", 
                                                  "free_energy", 
                                                  "pLDDT"),
                                   chain = "whole.complex",
                                   sequence.region = "binding.residues",
                                   propka.dir = "~/path/Propka_output/",
                                   germline.pdb = "~/path/PDBS_superimposed/germline_5_model_0.pdb",
                                   foldseek.dir = "~/path/3di_sequences/")

#Add the biophysical properties to the AntibodyForests object
af_structure_AbAg <- Af_add_node_feature(AntibodyForests_object = af,
                                                 feature.df = vdj_structure_AbAg,
                                                 feature.names = c("average_charge", 
                                                                   "3di_germline", 
                                                                   "average_hydrophobicity", 
                                                                   "RMSD_germline", 
                                                                   "average_pKa_shift", 
                                                                   "free_energy", 
                                                                   "average_pLDDT"))

Af_distance_scatterplot(af_structure_AbAg,
                        node.features = "free_energy",
                        correlation = "pearson",
                        color.by = "average_pLDDT",
                        color.by.numeric = T,
                        geom_smooth.method = "lm")

We investigate the change in protein structure during somatic hypermutation for the heavy and light chain sequences of the antibody.

rmsd <- Af_edge_RMSD(af,
             VDJ = vdj,
             pdb.dir = "~/path/PDBS_superimposed/",
             file.df = files,
             sequence.region = "full.sequence",
             chain = "HC+LC")
rmsd$plot