inst/doc/AntibodyForests_vignette.R

In AntibodyForests: Delineating Inter- And Intra-Antibody Repertoire Evolution

## ----setup, echo = FALSE, eval = TRUE, include = FALSE------------------------
# Set global options
knitr::opts_chunk$set(eval = TRUE, include = TRUE, echo = TRUE, message = FALSE, warning = FALSE, error = FALSE, dpi = 72)
# Define custom hooks for handling errors, warnings, and messages during knitting
knitr::knit_hooks$set(
   error = function(x, options){
     paste('\n\n<div class="alert alert-danger">',
           gsub('##', '\n', gsub('^##\ Error', '**Error**', x)),
           '</div>', sep = '\n')},
   warning = function(x, options){
     paste('\n\n<div class="alert alert-warning">',
           gsub('##', '\n', gsub('^##\ Warning:', '**Warning**', x)),
           '</div>', sep = '\n')},
   message = function(x, options){
     paste('\n\n<div class="alert alert-info">',
           gsub('##', '\n', x),
           '</div>', sep = '\n')})

#Load AntibodyForests functions
devtools::load_all(".")

## ----eval = FALSE-------------------------------------------------------------
#  #Install from CRAN
#  install.packages("Platypus")
#  install.packages("AntibodyForests")

## ----eval = FALSE-------------------------------------------------------------
#  #Load the libraries
#  library(Platypus)
#  library(AntibodyForests)

## ----installation, echo = FALSE, eval = FALSE---------------------------------
#  # List of CRAN packages required to run all functions
#  CRAN_packages <- c("Platypus", "alakazam", "ape", "base",
#                     "base64enc", "BiocManager", "dplyr",
#                     "grDevices", "gtools", "igraph",
#                     "jsonlite", "knitr", "parallel",
#                     "phangorn", "scales", "stats",
#                     "stringdist", "stringr", "utils")
#  
#  # List of Bioconductor packages required to run all functions
#  Bioconductor_packages <- c("Biostrings", "msa")
#  
#  # Iterate through packages
#  for(package in c(CRAN_packages, Bioconductor_packages)){
#    # Check if the current packages is already installed
#    if(!requireNamespace(package, quietly = TRUE)){
#      # Install CRAN packages if not installed
#      if(package %in% CRAN_packages){install.packages(package)}
#      # Install Bioconductor packages if not installed
#      if(package %in% Bioconductor_packages){BiocManager::install(package)}
#    }
#  }
#  

## ----quick-start1, eval = FALSE-----------------------------------------------
#  # Import 10x Genomics output files into VDJ dataframe and only keep cells with one VDJ and one VJ transcript
#  # For each clone, trim germline sequence and replace CDR3 region in germline with most-frequently observed CDR3 sequence
#  VDJ_OVA <- Platypus::VDJ_build(VDJ.directory = "10x_output/VDJ/",
#                       remove.divergent.cells = TRUE,
#                       complete.cells.only = TRUE,
#                       trim.germlines = TRUE)
#  
#  # Build lineage trees for all clones present in the VDJ dataframe with the default algorithm
#  AntibodyForests_OVA_default <- Af_build(VDJ = VDJ_OVA,
#                                            construction.method = "phylo.network.default")

## ----quick-start2, eval = F---------------------------------------------------
#  # Plot one of the lineage trees, constructed with the `phylo.network.default' method
#  Af_plot_tree(AntibodyForests_object = AntibodyForests_OVA_default,
#                       sample = "S1",
#                       clonotype = "clonotype2")

## ----quick-start3, eval = F---------------------------------------------------
#  # 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 = 5,
#                                       distance.method = "euclidean",
#                                       distance.metrics = c("mean.depth", "sackin.index", "spectral.density"),
#                                       clustering.method = "mediods",
#                                       visualization.methods = "PCA",
#                                       plot.label = T)
#  
#  # 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("spectral.density"),
#                     min.nodes = 5,
#                     significance = T)
#  

## ----quick-start4, eval=FALSE-------------------------------------------------
#  # Plot the PCA colored on the assigned clusters
#  out$plots$PCA_clusters

## ----quick-start5, eval=FALSE-------------------------------------------------
#  # Plot the different number of modalities between the clusters
#  plots$modalities

## ----vdj_build-example-3-code, eval = FALSE-----------------------------------
#  # 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)

## ----vdj_build-example-3-message, echo = FALSE, results = "asis"--------------
# Print number of filtered cells/barcodes 
warning("Please be aware of the number of cells excluded, when 'remove.divergent.cells' or 'complete.cells.only' is set to TRUE:\n")
print(knitr::kable(data.frame(divergent.cells = c(360, 295, 79, 183 ,526), 
                              incomplete.cells = c(708, 503, 387, 380, 1113),
                              row.names = c("S1", "S2", "S3", "S4", "S5"))))

## ----vdj_build-example-3-df, echo = FALSE-------------------------------------
# Load the example VDJdataframe
load("objects/VDJ_OVA_dummy_3.RData")
knitr::kable(head(VDJ_OVA_dummy_3[,c("sample_id", "barcode", "isotype", "clonotype_id", "clonotype_frequency", "VDJ_vgene", "VDJ_jgene")], 10))

## ----AntibodyForests-example-default, eval = FALSE----------------------------
#  # Infer lineage trees for all clones using the 'phylo.network.default' construction method and using all default parameter settings
#  AntibodyForests_OVA_default <- Af_build(VDJ = VDJ_OVA,
#                                                 sequence.columns = c("VDJ_sequence_nt_trimmed", "VJ_sequence_nt_trimmed"),
#                                                 germline.columns = c("VDJ_germline_nt_trimmed", "VJ_germline_nt_trimmed"),
#                                                 concatenate.sequences = FALSE,
#                                                 node.features = "isotype",
#                                                 construction.method = "phylo.network.default",
#                                                 string.dist.metric = "lv",
#                                                 resolve.ties = c("max.expansion", "close.germline.dist", "close.germline.edges", "random"),
#                                                 parallel = TRUE)

## ----AntibodyForests-example-ML, eval = FALSE---------------------------------
#  # Infer lineage trees for all clones using the 'phylo.tree.ml' construction method and using all default parameter settings
#  AntibodyForests_OVA_ml <- Af_build(VDJ = VDJ_OVA,
#                                            sequence.columns = c("VDJ_sequence_nt_trimmed", "VJ_sequence_nt_trimmed"),
#                                            germline.columns = c("VDJ_germline_nt_trimmed", "VJ_germline_nt_trimmed"),
#                                            concatenate.sequences = FALSE,
#                                            construction.method = "phylo.tree.ml",
#                                            dna.model = "all",
#                                            remove.internal.nodes = "connect.to.parent",
#                                            parallel = TRUE)

## ----AntibodyForests-example-IgPhyML, eval = FALSE----------------------------
#  # Infer lineage trees for all clones using the 'phylo.tree.IgPhyML' construction method and using all default parameter settings
#  AntibodyForests_OVA_IgPhyML <- Af_build(VDJ = VDJ_OVA,
#                                                 sequence.columns = "VDJ_sequence_nt_trimmed",
#                                                 germline.columns = "VDJ_germline_nt_trimmed",
#                                                 construction.method = "phylo.tree.IgPhyML",
#                                                 IgPhyML.output.file = "data/OVA_scRNA-seq_data/VDJ/airr_rearrangement_igphyml-pass.tab",
#                                                 remove.internal.nodes = "connect.to.parent",
#                                                 parallel = TRUE)

## ----AntibodyForests_plot_ML_with_IN_example, eval = FALSE--------------------
#  # Plot lineage tree, constructed with the `phylo.tree.ml' method, without internal nodes
#  Af_plot_tree(AntibodyForests_object = AntibodyForests_OVA_ml,
#                       sample = "S1",
#                       clonotype = "clonotype8",
#                       show.inner.nodes = TRUE,
#                       label.by = "name",
#                       color.by = "isotype",
#                       node.size = "expansion",
#                       main.title = "Lineage tree - ML",
#                       sub.title = "OVA - S1 - clonotype8")

## ----AntibodyForests_plot_ML_without_IN_example, eval = FALSE-----------------
#  # Plot lineage tree, constructed with the `phylo.tree.ml' method, without internal nodes
#  Af_plot_tree(AntibodyForests_object = AntibodyForests_OVA_ml,
#                       sample = "S1",
#                       clonotype = "clonotype8",
#                       show.inner.nodes = FALSE,
#                       label.by = "name",
#                       color.by = "isotype",
#                       node.size = "expansion",
#                       main.title = "Lineage tree - ML",
#                       sub.title = "OVA - S1 - clonotype8")

## ----Af_sync_nodes, eval= FALSE-----------------------------------------------
#  # 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)

## ----compare_methods_example2, eval = FALSE-----------------------------------
#  # 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

## ----metric_dataframe1, eval = FALSE------------------------------------------
#  # 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"))
#  

## ----metric_dataframe2--------------------------------------------------------
load("objects/metric_df.RData")

## ----metric_dataframe3, echo=FALSE--------------------------------------------
# Print the dataframe
DT::datatable(data = metric_df,
              rownames = TRUE,
              options = list(scrollX = TRUE))

## ----compare_distance_isotype, eval=FALSE-------------------------------------
#  # 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)

## ----compare_within_repertoire_example1, eval = FALSE-------------------------
#  # 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

## ----cluster_metrics_example1, eval = FALSE-----------------------------------
#  # 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

## ----cluster_node_features_example1, eval = FALSE-----------------------------
#  # 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

## ----compare_within_repertoire_example2, eval = FALSE-------------------------
#  # 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

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

## ----compare_across_repertoire_example1, eval=FALSE---------------------------
#  #Calculate the difference between the two groups of repertoires
#  boxplots <- Af_compare_across_repertoires(list("Blood" = AntibodyForests_example_blood,
#                                                 "Lymph Nodes" = AntibodyForests_example_LN),
#                                            metrics = c("betweenness", "mean.depth"), significance = T)
#  
#  #Plot the boxplots
#  boxplots$betweenness
#  boxplots$mean.depth

## ----add_pseudolikelihood_example, eval=FALSE---------------------------------
#  # 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_example1, eval=FALSE-------------------------------
#  # 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")

## ----get_sequences_probability_example, eval=FALSE----------------------------
#  #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)

## ----probability_matrix_example, eval=FALSE-----------------------------------
#  #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")

## ----load_plm_dataframe, echo=FALSE-------------------------------------------
load("objects/PLM_dataframe.RData")

## ----vdj_plm-df, echo = FALSE-------------------------------------------------
knitr::kable(PLM_dataframe)

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

## ----load_structure_data, echo=FALSE------------------------------------------
load("objects/VDJ_ova_structure.RData") #vdj
load("objects/VDJ_structure_AbAg.RData") #vdj_structure_AbAg
load("objects/AntibodyForests_OVA_structure_AbAg.RData") #af_structure_AbAg

## ----structure_example2, eval = F---------------------------------------------
#  #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"))

## ----structure_example4_plot, eval = T, fig.width=5, fig.height=5-------------
Af_distance_scatterplot(af_structure_AbAg,
                        node.features = "free_energy",
                        correlation = "pearson",
                        color.by = "average_pLDDT",
                        color.by.numeric = T,
                        geom_smooth.method = "lm")

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

## ----write-AIRR-rearrangement, eval = FALSE-----------------------------------
#  # Append the IgBLAST annotations and alignment to the VDJ dataframe
#  VDJ_OVA <- VDJ_import_IgBLAST_annotations(VDJ = VDJ_OVA,
#                                            VDJ.directory = "../data/OVA_scRNA-seq_data/VDJ/",
#                                            method = "append")
#  
#  # Write the VDJ dataframe into an AIRR-formatted TSV file
#  VDJ_to_AIRR(VDJ = VDJ_OVA,
#              output.file = "../data/OVA_scRNA-seq_data/VDJ/airr_rearrangement.tsv")

## ----bulk_example, eval = FALSE-----------------------------------------------
#  # Integrate bulk RNA-seq data into the VDJ dataframe, if multiple clonotypes match the bulk sequence, the bulk transcript is added to all clonotypes (tie.resolvement = "all)
#  VDJ_OVA <- VDJ_integrate_bulk(sc.VDJ = VDJ_OVA,
#                     bulk.tsv = "path/to/bulk/OVA_bulkRNA.tsv",
#                     bulk.tsv.sequence.column="VDJ_sequence_nt_raw",
#                     bulk.tsv.sample.column="sample_id",
#                     bulk.tsv.barcode.column="barcode",
#                     bulk.tsv.isotype.column="isotype",
#                     organism="mouse",
#                     igblast.dir="/path/to/anaconda3/envs/igphyml/share/igblast",
#                     tie.resolvement = "all")
#  
#  #Create the lineage tree of the integrated VDJ dataframe
#  af <- Af_build(VDJ = VDJ_OVA,
#                 sequence.columns = "VDJ_sequence_aa_trimmed",
#                 germline.columns = "VDJ_germline_aa_trimmed",
#                 construction.method = "phylo.network.default",
#                 node.features = c("dataset"))
#  
#  #Plot lineage tree
#  Af_plot_tree(af,
#               sample = "S1",
#               clonotype = "clonotype3",
#               color.by = "dataset",
#               label.by = "none",
#               main.title = "OVA S1",
#               sub.title = "clonotype3")