TCRView.Rmd

---
title: "TCRView: Visualizing relationships in a comprehensive reference dataset of paired-chain T cell receptor sequences and epitope specificities"
author: "Maia Bennett-Boehm"
date: "`r Sys.Date()`"
output: html_document
---

# TCRView: Visualizing relationships in a comprehensive reference dataset of paired-chain T cell receptor (TCR) sequences and epitope specificities

## Introduction

This document facilitates the visualization of paired-chain T cell receptor (TCR) datasets prepared using TCRPaired. TCRView provides a complementary set of tools to explore the distribution of, and relationships between, TCR:pMHC complex components in a benchmark dataset. These components include epitopes, major histocompatibility complexes (MHCs), and TCR features such as VJ genes and sequence segments. 


```{r setup, include=FALSE}
source("./util/library.R")
source("./util/plotting.R")
source("./util/clustering.R")
source("./util/filtering.R")
library(rstudioapi)
library(plotly)

knitr::opts_chunk$set(
  fig.width = 10  # Set figure width
)

```

```{r data-import}

tcr.view.file <- "./paired-data/all-distinct-cdr-paired-sequences.csv"
tcr.view.file <- "./paired-data/high-confidence-paired-sequences.csv"
tcr.view.file <- "./paired-data/high-confidence-cdr3-similarity-90-paired-sequences.csv"
tcr.view.set <- read.csv(tcr.view.file)

# To make plots more legible, string wrap Epitope species values to 12 characters
tcr.view.set$Epitope.species <- str_wrap(tcr.view.set$Epitope.species, 12)

```

TCRView is intended to run on entire datasets as well as data subsets for both 1) investigation of data composition and bias in TCR:pMHC complex datasets and 2) exploration of relationships between TCR:pMHC complex components. As such, a number of parameters can be set to filter and facet the data, as well as to save clustering data for further analysis. 

The following parameters can be set to customize the data visualization:

- `plotly`: Set to `TRUE` to generate interactive plotly plots, or `FALSE` to generate basic ggplot2 plots.

- `filter.by`: A filtering expression to subset the data. For example, `filter.by = "Epitope.species == 'Influenza A'"` will filter the data to include only receptors with epitopes from Influenza A in all plots. This can be useful for exploring relationships between TCR:pMHC complex components for specific epitopes, epitope species, MHC loci, or V/J gene profiles.

- `group.by`: The name of a column to group the data by. For example, `group.by = "Epitope.species"` will create separate distribution lines for each epitope species in count distribution plots and limit composition plots to those epitopes. This can be useful for comparing the distribution of TCR:pMHC complex components across different subgroups. 

- `top`: The number of top group.by values to display in count distribution and composition plots. For example, `top = 10` will display the top 10 epitopes in the dataset. Additionally, this has use in the clustering section, which is specified later.

- `option`: The name of the viridis palette to use for coloring the plots. For example, `option = "turbo"` will color the plots using the turbo palette. This can be useful for distinguishing between different group.by values in the plots.

- `begin`: The beginning of the viridis palette to use for coloring the plots. For example, `begin = 0.2` will color the plots using the turbo palette starting at the 20th percentile. This can be useful for distinguishing between different group.by values in the plots.

- `end`: The end of the viridis palette to use for coloring the plots. For example, `end = 0.8` will color the plots using the turbo palette ending at the 80th percentile. This can be useful for distinguishing between different group.by values in the plots.

- `out.path`: The path to save the clustering data. This allows you to save plots and data for further analysis, especially when running multiple subsets of the same data.


Due to the reuse of some of these parameters in differing plot types across components, there are multiple `setup` sections throughout TCRView to allow ease of analysis while specifying intended plot orientations. These also introduce new parameters as needed.


```{r setup-epitope}

# Do you want the plots to be a basic ggplot or interactive plotly plot?
plotly = FALSE

# Do you want to filter the data using a filtering expression (e.g., str_detect(Epitope.species, "SARS-CoV-2"))?
filter.by = NULL

# Do you want to group the data by a factor (e.g., "Epitope.species")?
group.by = "Epitope.species"

# Do you want to limit the grouped data to the top n values?
top = 5

# How do you want to color the plots (viridis palettes)?
option = "turbo"
begin = 0
end = 1

# Where do you want to save any plot and clustering data?
out.path = "./view-data/all-data/"
out.path = "./view-data/high-confidence-data/"
out.path = "./view-data/high-confidence-similarity-90-data/"

dir.create(out.path, showWarnings = FALSE, recursive = TRUE)

```


## Components of the TCR:pMHC complex

### Epitopes

Epitopes are the antigenic peptides that are presented by MHC molecules and recognized by TCRs. The distribution of epitopes in a dataset can provide insight into the diversity of the dataset and the potential biases in the data. The following tables and plots provide an overview of the epitopes in the dataset.

```{r epitope-table, echo = FALSE}

if (!is.null(filter.by)) {
    filter.set <- tcr.view.set %>% filter(!!rlang::parse_expr(filter.by))
} else {
    filter.set <- tcr.view.set
}

# Epitope counts table (highest to lowest abundance)
epitope.table <- filter.set %>%
    select(Epitope, Epitope.gene, Epitope.species) %>%
    add_count(Epitope, name = "Epitope Count") %>%
    add_count(Epitope.gene, name = "Epitope Gene Count") %>%
    add_count(Epitope.species, name = "Epitope Species Count") %>%
    distinct(.keep_all = TRUE) %>%
    arrange(desc("Epitope Count"))

write.csv(epitope.table, file = paste0(out.path, "epitope-table.csv"), row.names = FALSE)

datatable((epitope.table %>%
        mutate(`Epitope Species` = paste0(Epitope.species, ": ", `Epitope Species Count`, " receptors")) %>%
        arrange(`Epitope Species`)), 
    extensions=c('Scroller', 'RowGroup'),
    options = list(dom='t', scrollX=TRUE, deferRender=TRUE, scrollY=200, scroller=TRUE, 
    rowGroup = list(dataSrc = 7), columnDefs = list(list(visible=FALSE, targets=c(7)))),
    caption="Epitope distribution in the TCR dataset")


```

```{r epitope-count-dist}

# Epitope count distribution
plotEpitopeCountDistribution(tcr.view.set, plotly = plotly, filter.by = filter.by, option = option) 
ggsave(paste0(out.path, "epitope-count-distribution.png"), width=10, height=5)

# Epitope count distribution per epitope species (top species, all other as "Other")
plotEpitopeCountDistribution(tcr.view.set, group.by = group.by, top = top, plotly = plotly, filter.by = filter.by, option = option, begin = begin, end = end)

# All epitopes
plotEpitopeComposition(tcr.view.set, plotly = plotly, filter.by = filter.by, option = option, begin = begin, end = end)

# Epitopes for top epitope species 
plotEpitopeComposition(tcr.view.set, group.by = group.by, top = top, plotly = plotly, filter.by = filter.by, option = option, begin = begin, end = end)

# All epitope species
plotEpitopeSpeciesComposition(tcr.view.set %>% 
    mutate(Epitope.species = ifelse(str_detect(Epitope.species, "SARS-CoV-2"), 
        "Severe acute\nrespiratory syndrome\ncoronavirus 2 (SARS-CoV-2)", 
        Epitope.species), 
    Epitope.species = ifelse(str_detect(Epitope.species, "SARS-CoV-1"), 
        "Severe acute\nrespiratory syndrome\ncoronavirus 1 (SARS-CoV-1)", 
        Epitope.species),
    Epitope.species = ifelse(str_detect(Epitope.species, "Salmonella"), 
        "Salmonella enterica", 
        Epitope.species),
    Epitope.species = ifelse(str_detect(Epitope.species, "\\("), 
        str_extract(Epitope.species, "(?<=\\().+?(?=\\))"), 
        Epitope.species)), 
    plotly = plotly, filter.by = filter.by, option = option, begin = begin, end = end) + 
    theme(title = element_text(size = 14), axis.text.x = element_text(size=12))
ggsave(paste0(out.path, "epitope-species-composition.png"), width=12, height=6)

```

### Major histocompatibility complexes (MHCs)

MHC molecules are responsible for presenting epitopes to TCRs. Although not directly present involved in the TCR:epitope binding interface, MHC identity plays an important role in epitope processing and assist in determining protein complex binding conformations. The following tables and plots provide an overview of the MHCs in the dataset.


```{r epitope-mhc-table}

# MHC + epitope counts table
mhc.table <- filter.set %>%
    select(Epitope, Epitope.gene, Epitope.species, MHC.allele) %>%
    count(Epitope, Epitope.gene, Epitope.species, MHC.allele) %>% 
    arrange(desc(n))

write.csv(mhc.table, file = paste0(out.path, "epitope-mhc-table.csv"), row.names = FALSE)

datatable((mhc.table %>% arrange(Epitope.species)), extensions=c('Scroller', 'RowGroup'),
    options = list(dom='t', scrollX=TRUE, deferRender=TRUE, scrollY=200, scroller=TRUE, 
    rowGroup = list(dataSrc = 3)),
    caption="Epitope + MHC distribution in the TCR dataset")



```


```{r mhc-table, echo = FALSE}

# MHC counts table
mhc.table <- filter.set %>%
    select(MHC.locus, MHC.allele, MHC) %>%
    add_count(MHC.locus, name = "MHC Locus Count") %>%
    add_count(MHC.allele, name = "MHC Allele Count") %>%
    add_count(MHC, name = "MHC Group Count") %>%
    distinct(.keep_all = TRUE) %>%
    na.omit() %>%
    filter(MHC.locus != "") %>%
    mutate(`MHC Locus` = paste0(MHC.locus, ": ", `MHC Locus Count`, " receptors")) %>%
    arrange(MHC.allele)

write.csv(mhc.table, file = paste0(out.path, "mhc-table.csv"), row.names = FALSE)

datatable(mhc.table, extensions=c('Scroller', 'RowGroup'),
    options = list(dom='t', scrollX=TRUE, deferRender=TRUE, scrollY=200, scroller=TRUE, 
    rowGroup = list(dataSrc = 1)),
    caption="MHC distribution in the TCR dataset")

mhc.locus.table <- filter.set %>% 
    count(MHC.locus) %>%
    arrange(desc(n))

write.csv(mhc.locus.table, file = paste0(out.path, "mhc-locus-table.csv"), row.names = FALSE)

rm(filter.set)


```




```{r setup-mhc}

# Do you want the plots to be a basic ggplot or interactive plotly plot?
plotly = FALSE

# Do you want to filter the data using a filtering expression (e.g., str_detect(Epitope.species, "SARS-CoV-2"))?
filter.by = NULL

# Do you want to group the data by a factor (e.g., "Epitope.species")?
group.by = "MHC.locus"

# Do you want to limit the grouped data to the top n values?
top = 5

# How do you want to color the plots (viridis palettes)?
option = "plasma"
begin = 0.9
end = 0.1

```


```{r mhc}

# MHC distribution
# MHC locus
plotMHCLocusCountDistribution(tcr.view.set, plotly = plotly, filter.by = filter.by, option = option)

# MHC allele
plotMHCAlleleCountDistribution(tcr.view.set, plotly = plotly, filter.by = filter.by, option = option)

# MHC 
plotMHCCountDistribution(tcr.view.set, plotly = plotly, filter.by = filter.by, option = option)

# MHC allele distribution, by loci
plotMHCAlleleCountDistribution(tcr.view.set, group.by = group.by, plotly = plotly, filter.by = filter.by, option = option, begin = begin, end = end)

# MHC composition
# MHC locus
plotMHCLocusComposition(tcr.view.set, plotly = plotly, filter.by = filter.by, option = option, begin = begin, end = end) + 
    theme(title = element_text(size = 14), axis.text.x = element_text(size=14))
ggsave(paste0(out.path, "mhc-locus-composition.png"), width=12, height=6)

# MHC allele
plotMHCAlleleComposition(tcr.view.set, plotly = plotly, filter.by = filter.by, option = option, begin = begin, end = end) + 
    theme(title = element_text(size = 14), axis.text.x = element_text(size=14))
ggsave(paste0(out.path, "mhc-allele-composition.png"), width=12, height=6)

# MHC, by loci
plotMHCComposition(tcr.view.set, plotly = plotly, filter.by = filter.by, group.by = group.by, option = option, begin = begin, end = end)


```

### T cell receptors (TCRs)

#### Variable and joining (VJ) genes

VJ genes are the variable and joining gene segments that make up the TCR. These genes contribute directly to TCR sequence through the process of variable (diversity) joining recombination, and as such play an important role in understanding and predicting TCR specificity. The following tables and plots provide an overview of the VJ genes in the dataset. 

```{r vj-gene-table, echo = FALSE}

if (!is.null(filter.by)) {
    vj.set <- tcr.view.set %>% filter(!!rlang::parse_expr(filter.by))
} else {
    vj.set <- tcr.view.set
}

# Datatable of VJ alleles, genes, families, with counts
vj.table <- vj.set %>%
    pivot_longer(cols = c(AV, AJ, BV, BJ), names_to = "Gene Segment", values_to = "Allele") %>%
    mutate(keep = case_when(
        `Gene Segment` == "AV" ~ "AV.gene",
        `Gene Segment` == "AJ" ~ "AJ.gene",
        `Gene Segment` == "BV" ~ "BV.gene",
        `Gene Segment` == "BJ" ~ "BJ.gene"
    )) %>%
    pivot_longer(cols = c(AV.gene, AJ.gene, BV.gene, BJ.gene), names_to = "Gene Fam. Segment", values_to = "Gene") %>%
    filter(keep == `Gene Fam. Segment`) %>%
    mutate(keep = case_when(
        `Gene Fam. Segment` == "AV.gene" ~ "AV.family",
        `Gene Fam. Segment` == "AJ.gene" ~ "AJ.family",
        `Gene Fam. Segment` == "BV.gene" ~ "BV.family",
        `Gene Fam. Segment` == "BJ.gene" ~ "BJ.family"
    )) %>%
    pivot_longer(cols = c(AV.family, AJ.family, BV.family, BJ.family), names_to = "Gene Chain", values_to = "Gene Family") %>%
    filter(keep == `Gene Chain`) %>%
    select(Allele, Gene, `Gene Family`, `Gene Segment`) %>%
    add_count(Allele, name = "Allele Count") %>%
    add_count(Gene, name = "Gene Count") %>%
    add_count(`Gene Family`, name = "Family Count") %>%
    distinct(.keep_all = TRUE) %>%
    group_by(`Gene Segment`) %>%
    add_count(`Gene Segment`, name = "Gene Segment Count") %>%
    ungroup() %>%
    mutate(`Gene Segment` = paste0(`Gene Segment`, ": ", `Gene Segment Count`, " alleles")) %>%
    select(-`Gene Segment Count`) %>%
    arrange(`Gene Segment`) %>%
    distinct(Allele, Gene, `Gene Family`, .keep_all = TRUE) 

write.csv(vj.table %>% filter(str_detect(Allele, "AV")) %>% select(-"Gene Segment") %>% ungroup() %>% arrange(desc(`Allele Count`)), file = paste0(out.path, "av-gene-table.csv"), row.names = FALSE)
write.csv(vj.table %>% filter(str_detect(Allele, "AJ")) %>% select(-"Gene Segment") %>% arrange(desc(`Allele Count`)), file = paste0(out.path, "aj-gene-table.csv"), row.names = FALSE)
write.csv(vj.table %>% filter(str_detect(Allele, "BV")) %>% select(-"Gene Segment") %>% arrange(desc(`Allele Count`)), file = paste0(out.path, "bv-gene-table.csv"), row.names = FALSE)
write.csv(vj.table %>% filter(str_detect(Allele, "BJ")) %>% select(-"Gene Segment") %>% arrange(desc(`Allele Count`)), file = paste0(out.path, "bj-gene-table.csv"), row.names = FALSE)

datatable(vj.table, extensions=c('Scroller', 'RowGroup'), 
    options = list(dom='t', scrollX=TRUE, deferRender=TRUE, scrollY=200, scroller=TRUE, 
    rowGroup = list(dataSrc = 4), columnDefs = list(list(visible=FALSE, targets=c(4)))),
    caption="V and J gene usage in the TCR dataset")


co.gene.table <- vj.set %>% 
    count(AV.gene, AJ.gene, BV.gene, BJ.gene) %>%
    arrange(desc(n))

write.csv(co.gene.table, file = paste0(out.path, "vj-gene-cooccurrence-table.csv"), row.names = FALSE)


rm(vj.set)

```



```{r setup-vj-gene}

# Do you want the plots to be a basic ggplot or interactive plotly plot?
plotly = FALSE

# Do you want to filter the data using a filtering expression (e.g., str_detect(Epitope.species, "SARS-CoV-2"))?
filter.by = NULL

# Do you want to facet the data by a factor (e.g., "Epitope.species")?
group.by = "`Gene Family`"

# Do you want to limit the filtered and/or faceted data to the top n values?
top = 10

# How do you want to color the plots (viridis palettes)?
option = "viridis"
begin = 0
end = 1

```


```{r vj-genes}

# VJ composition
# V J alleles 
plotVJAllele(tcr.view.set, plotly = plotly, filter.by = filter.by, option = option, begin = begin, end = end)

# VJ genes 
plotVJGene(tcr.view.set, plotly = plotly, filter.by = filter.by, option = option, begin = begin, end = end)

# VJ gene families
plotVJFamily(tcr.view.set, plotly = plotly, filter.by = filter.by, option = option, begin = begin, end = end) + 
    theme(title = element_text(size = 14), axis.text.x = element_text(size=12), strip.text = element_text(size = 14))
ggsave(paste0(out.path, "vj-family-composition.png"), width=20, height=8)

```


#### Sequence

TCR sequences are highly variable, with the complementarity determining regions (CDRs) being the most variable regions. This high degree of variability directly informs the vast diversity of TCR specificity and binding mechanisms. As such, multiple facets of the TCR sequence are explored in the following sections. 


The following parameters are changed or added for the analysis of TCR sequences, particularly clustering of complementarity determining regions (CDRs) and full sequences:

- `top.col`: If provided, the name of a column to limit the coloring of the data to. For example, `top.col = "Epitope.species"` will limit the coloring of the data to the epitope species in clustering plots. All other values will be colored as "Other". Cannot be used concurrently with 'threshold.col'.

- `top`: The number of top facet.by values to display in count distribution and composition plots. For example, `top = 10` will display the top 10 epitopes in the dataset. Additionally, this is the number of top values to display in clustering plots, *e.g.*, the top 10 epitope species in the dataset.

- `threshold.col`: If provided, the name of a column to limit the coloring of the data to some feature with > n receptors. For example, `threshold.col = "Epitope"` and `threshold = 50` will limit the coloring of the data to receptors with epitopes that have > 50 receptors in clustering plots. All other values will be colored as "Other". Cannot be used concurrently with 'top.col'.

- `threshold`: The number of receptors to limit the coloring to. For example, `threshold = 50` will limit the coloring of the data to receptors with epitopes that have > 50 receptors in clustering plots. All other values will be colored as "Other".

- `highlight`: A vector of values to highlight in clustering plots **without** discarding non-matching data, as is accomplished using `filter.by`. For example, `highlight = c("Epitope == 'GILGFVFTL'")` will highlight receptors with the epitope "GILGFVFTL" in clustering plots, coloring all other receptors as "Other". It can also be used to highlight specific TCRs using the clone.id parameter for use in benchmark case studies.

- `color`: If used in conjunction with highlight, the name of a column to label the highlighted values by. For example, `color = "Epitope"` will label the highlighted values by epitope in clustering plots. If not used in conjunction with highlight, this parameter will be used to color the data in clustering plots. The default value of `color` is "Epitope". 

- `save`: Set to `TRUE` to save clustering data in the initial TCRView dataset, or `FALSE` to not save clustering data. If not run prior to the clustering functions, plots will include the clustering step, but UMAP coordinates cannot be guaranteed to be the same across plots. 


As before, these parameters may be changed on a per-plot basis. This is shown in the cluster plots below, which aim to highlight the utility of clustering functions. 


```{r setup-tcr}

# Do you want the plots to be a basic ggplot or interactive plotly plot?
plotly = FALSE

# Do you want to filter the data using a filtering expression (e.g., str_detect(Epitope.species, "SARS-CoV-2"))?
filter.by = NULL

# Do you want to limit the coloring of the data to a specific factor (e.g., "Epitope.species")?
top.col = NULL

# How many top values do you want to display?
top = 10

# Do you want to limit the coloring of the data to some feature with > n receptors?
threshold.col = NULL

# How many receptors do you want to limit the coloring to?
threshold = NULL

# Do you want to highlight the position of specific values (e.g., clone.id %in% ..., Epitope == "GILGFVFTL") without hiding the rest of the data?
highlight = NULL

# What column do you want to label these values by?
color = NULL

# How do you want to color the plots (viridis palettes)?
option = "turbo"
begin = 0
end = 1

# Do you want to save clustering data?
save = FALSE


```


##### Complementarity determining region 3 (CDR3)

CDR3 is the most variable region of the TCR, and is formed by random recombination rather than V(D)J recombination. Structurally, CDR3 sits directly at the binding interface of the TCR:pMHC complex, and as such is a key determinant of TCR specificity. Many binding prediction methods focus exclusively on the CDR3 region, and as such, it is important to understand the distribution of CDR3 lengths and sequences in a dataset.

```{r cdr3-length}

# CDR3 alpha and beta lengths
# By top 10 epitopes
plotCDR3SeqLength(tcr.view.set, color.by = "Epitope", top = top, plotly = plotly, filter.by = filter.by, option = "turbo")

# By epitope species
plotCDR3SeqLength(tcr.view.set, color.by = "Epitope.species", plotly = plotly, filter.by = filter.by, option = "turbo")

# By MHC locus
plotCDR3SeqLength(tcr.view.set, color.by = "MHC.locus", plotly = plotly, filter.by = filter.by, option = "plasma", begin = 0.9, end = 0.1)


```


```{r cdr3-seq-sim}

# CDR3a sequence similarity
tcr.view.set <- computeSimilarity(tcr.view.set, tcr.view.set$CDR3a, tcr.view.set$CDR3a) %>%
    select(-dist.matrix, -similarity.level) %>%
    dplyr::rename(CDR3a.max.similarity = max.similarity.level)

# CDR3b sequence similarity
tcr.view.set <- computeSimilarity(tcr.view.set, tcr.view.set$CDR3b, tcr.view.set$CDR3b) %>%
    select(-dist.matrix, -similarity.level) %>%
    dplyr::rename(CDR3b.max.similarity = max.similarity.level)

# Density plot of max similarity in CDR3a and CDR3b
plotCDR3SimilarityDistribution(tcr.view.set, plotly = plotly, filter.by = filter.by, option = "mako", begin=0.2, end=0.8) +
    theme(title = element_text(size = 14), axis.text.x = element_text(size=12))
ggsave(paste0(out.path, "max-similarity-density.png"), width=12, height=6)

# Box plot per top epitopes
plotCDR3SimilarityByFactor(tcr.view.set, group.by = "Epitope", top = top, plotly = plotly, filter.by = filter.by, option = "mako", begin=0.2, end=0.8) +
    theme(title = element_text(size = 14), axis.text.x = element_text(size=14))
ggsave(paste0(out.path, "max-similarity-boxplot.png"), width=12, height=6)

```


```{r cdr3-clust}

set.seed(77482951)

# CDR3 
# Combined chains
# Run clustering 
# tcr.view.set <- getCDR3SeqCluster(tcr.view.set, combine.chains = TRUE, filter.by = filter.by)

# All, colored by cluster (default)
plotCDR3SeqClusterResults(tcr.view.set, combine.chains = TRUE, plotly = plotly, filter.by = filter.by, option  = option, begin = begin, end = end)

# Top 10 epitope species
plotCDR3SeqClusterResults(tcr.view.set, combine.chains = TRUE, top = top, top.col = "Epitope.species", plotly = plotly, filter.by = filter.by, option  = option, begin = begin, end = end)

# Receptors from epitopes with > 50 receptors
plotCDR3SeqClusterResults(tcr.view.set, combine.chains = TRUE, threshold = 50, threshold.col = "Epitope", plotly = plotly, filter.by = filter.by, option  = option, begin = begin, end = end)

# Receptors limited to SARS-CoV-2 epitopes, colored by epitope: this will not save clusters to the original data file; if you wanted to do that, you'd have to re-run the clustering with the filter.by argument
plotCDR3SeqClusterResults(tcr.view.set, combine.chains = TRUE, plotly = plotly, filter.by = "Epitope.species == 'Influenza A'", color = "Epitope", option  = option, begin = begin, end = end)

# SARS-CoV-2 epitope, overlaid on all receptors
plotCDR3SeqClusterResults(tcr.view.set, combine.chains = TRUE, plotly = plotly, highlight = "Epitope.species == 'Influenza A'", color = "Epitope", option  = option, begin = begin, end = end)

# Per chain
# Run clustering 
# tcr.view.set <- getCDR3SeqCluster(tcr.view.set, combine.chains = FALSE, filter.by = filter.by)

# All, colored by cluster (default)
plotCDR3SeqClusterResults(tcr.view.set, combine.chains = FALSE, plotly = plotly, filter.by = filter.by, option  = option, begin = begin, end = end)

# Top 10 epitope species
plotCDR3SeqClusterResults(tcr.view.set, combine.chains = FALSE, top = top, top.col = "Epitope.species", plotly = plotly, filter.by = filter.by, option  = option, begin = begin, end = end)

# Receptors from epitopes with > 50 receptors
plotCDR3SeqClusterResults(tcr.view.set, combine.chains = FALSE, threshold = 50, threshold.col = "Epitope", plotly = plotly, filter.by = filter.by, option  = option, begin = begin, end = end)

# Receptors limited to SARS-CoV-2 epitopes
plotCDR3SeqClusterResults(tcr.view.set, combine.chains = FALSE, plotly = plotly, filter.by = "Epitope.species == 'Influenza A'", color = "Epitope", option  = option, begin = begin, end = end)

# SARS-CoV-2 epitope, overlaid on all receptors
plotCDR3SeqClusterResults(tcr.view.set, combine.chains = FALSE, plotly = plotly, highlight = "Epitope.species == 'Influenza A'", color = "Epitope", option  = option, begin = begin, end = end)

```

##### All complementarity determining regions (CDRs)

Although not as immensely diverse as the CDR3 region, the remaining CDRs (CDR1, CDR2, and CDR2.5) are still highly diverse and contribute to TCR binding conformations. While commonly assumed that only the CDR3 makes direct contact with epitopes during binding, multiple studies have shown other CDRs to be involved in binding. 

```{r cdr-clust}

set.seed(77482951)

# CDR clustering
# Combined chains
# Run clustering 
# tcr.view.set <- getCDRSeqCluster(tcr.view.set, combine.chains = TRUE, filter.by = filter.by)

# All, colored by cluster (default)
plotCDRSeqClusterResults(tcr.view.set, combine.chains = TRUE, plotly = plotly, filter.by = filter.by, option  = option, begin = begin, end = end)

# Top 10 epitope species
plotCDRSeqClusterResults(tcr.view.set, combine.chains = TRUE, top = top, top.col = "Epitope.species", plotly = plotly, filter.by = filter.by, option  = option, begin = begin, end = end)

# Receptors from epitopes with > 50 receptors
plotCDRSeqClusterResults(tcr.view.set, combine.chains = TRUE, threshold = 50, threshold.col = "Epitope", plotly = plotly, filter.by = filter.by, option  = option, begin = begin, end = end)

# Receptors limited to SARS-CoV-2 epitopes, colored by epitope
plotCDRSeqClusterResults(tcr.view.set, combine.chains = TRUE, plotly = plotly, filter.by = "Epitope.species == 'Influenza A'", color = "Epitope", option  = option, begin = begin, end = end)

# SARS-CoV-2 epitope, overlaid on all receptors
plotCDRSeqClusterResults(tcr.view.set, combine.chains = TRUE, plotly = plotly, highlight = "Epitope.species == 'Influenza A'", color = "Epitope", option  = option, begin = begin, end = end)

# Per chain
# Run clustering 
# tcr.view.set <- getCDRSeqCluster(tcr.view.set, combine.chains = FALSE, filter.by = filter.by)

# All, colored by cluster (default)
plotCDRSeqClusterResults(tcr.view.set, combine.chains = FALSE, plotly = plotly, filter.by = filter.by, option  = option, begin = begin, end = end)

# Top 10 epitope species
plotCDRSeqClusterResults(tcr.view.set, combine.chains = FALSE, top = top, top.col = "Epitope.species", plotly = plotly, filter.by = filter.by, option  = option, begin = begin, end = end)

# Receptors from epitopes with > 50 receptors
plotCDRSeqClusterResults(tcr.view.set, combine.chains = FALSE, threshold = 50, threshold.col = "Epitope", plotly = plotly, filter.by = filter.by, option  = option, begin = begin, end = end)

# Receptors limited to SARS-CoV-2 epitopes
plotCDRSeqClusterResults(tcr.view.set, combine.chains = FALSE, plotly = plotly, filter.by = "Epitope.species == 'Influenza A'", color = "Epitope", option  = option, begin = begin, end = end)

# SARS-CoV-2 epitope, overlaid on all receptors
plotCDRSeqClusterResults(tcr.view.set, combine.chains = FALSE, plotly = plotly, highlight = "Epitope.species == 'Influenza A'", color = "Epitope", option  = option, begin = begin, end = end)


```

##### Full sequence

The full TCR sequence includes the V and J gene segments, the CDRs, and the constant region. The full sequence is important for understanding the complete TCR structure and function, and can provide insight into the diversity of the dataset.

```{r full-seq-clust}

# Full sequence clustering
# Run clustering 
# tcr.view.set <- getFullSeqCluster(tcr.view.set, combine.chains = TRUE, filter.by = filter.by)

# All, colored by cluster (default)
plotFullSeqClusterResults(tcr.view.set, plotly = plotly, filter.by = filter.by, option  = option, begin = begin, end = end)

# Top 10 epitope species
plotFullSeqClusterResults(tcr.view.set, top = top, top.col = "Epitope.species", plotly = plotly, filter.by = filter.by, option  = option, begin = begin, end = end)

# Receptors from epitopes with > 50 receptors
plotFullSeqClusterResults(tcr.view.set, threshold = 50, threshold.col = "Epitope", plotly = plotly, filter.by = filter.by, option  = option, begin = begin, end = end)

# Receptors limited to SARS-CoV-2 epitopes
plotFullSeqClusterResults(tcr.view.set, plotly = plotly, filter.by = "Epitope.species == 'Influenza A'", color = "Epitope", option  = option, begin = begin, end = end)

# SARS-CoV-2 epitope, overlaid on all receptors
plotFullSeqClusterResults(tcr.view.set, plotly = plotly, highlight = "Epitope.species == 'Influenza A'", color = "Epitope", option  = option, begin = begin, end = end)

```


```{r save-data}

if (save) {
    write.csv(tcr.view.set, file = paste0(out.path, "TCRView-data.csv"), row.names = FALSE)
}


```