TCR Analysis-7-11-2025

1. Introduction

Single-cell sequencing is an emerging technology in the field of immunology and oncology that allows researchers to couple RNA quantification and other modalities, like immune cell receptor profiling at the level of an individual cell. A number of workflows and software packages have been created to process and analyze single-cell transcriptomic data. These packages allow users to take the vast dimensionality of the data generated in single-cell-based experiments and distill the data into novel insights. Unlike the transcriptomic field, there is a lack of options for software that allow for single-cell immune receptor profiling. Enabling users to easily combine RNA and immune profiling, the scRepertoire framework supports use of 10x, single-cell clonal formats and interaction with popular R-based single-cell data pipelines.

scRepertoire is designed to take filter contig outputs from the 10x Genomics Cell Ranger pipeline, process that data to assign clonotype based on two TCR or Ig chains and analyze the clonotype dynamics. The latter can be separated into 1) clonotype-only analysis functions, such as unique clonotypes or clonal space quantification, and 2) interaction with mRNA expression data using Seurat or SingleCellExperiment packages.

References:

• scRepertoire Vignette
  
• Borch.dev scRepertoire

Prerequisite: Ensure the All_samples_Merged Seurat object is loaded into your R environment before running the chunks below.

1.1 Load libraries

library(Seurat)
library(SeuratObject)
library(SeuratObject)
library(SeuratData)
library(patchwork)
library(dplyr)
library(ggplot2)
library(tidyverse)
library(rmarkdown)
library(tinytex)
library(grid)
library(cowplot)
library(presto)

#TCR Analysis
library(scRepertoire)
library(SingleCellExperiment)
library(circlize)
library(scales)

1.2 Load Seurat Object

#Load Seurat Object merged from cell lines and a control(PBMC) after filtration
All_samples_Merged <- readRDS("../0-Seurat_RDS_OBJECT_FINAL/All_samples_Merged_with_Renamed_Clusters_final-26-10-2025.rds")

All_samples_Merged

An object of class Seurat 
62900 features across 49305 samples within 6 assays 
Active assay: RNA (36601 features, 0 variable features)
 2 layers present: data, counts
 5 other assays present: ADT, prediction.score.celltype.l1, prediction.score.celltype.l2, prediction.score.celltype.l3, SCT
 5 dimensional reductions calculated: integrated_dr, ref.umap, pca, umap, harmony

1.3. Load contigs

scRepertoire functions using the filtered_contig_annotations.csv output from the 10x Genomics Cell Ranger. This file is located in the ./outs/ directory of the VDJ alignment folder. To generate a list of contigs to use for scRepertoire:

- load the filtered_contig_annotations.csv for each of the samples.

- make a list in the R environment.

Create display_sample column for uniform naming

# Create a new metadata column for plotting with simplified sample names
All_samples_Merged$display_sample <- All_samples_Merged$orig.ident

# Rename L3_B and L4_B to L3 and L4
All_samples_Merged$display_sample <- gsub("^L3_B$", "L3", All_samples_Merged$display_sample)
All_samples_Merged$display_sample <- gsub("^L4_B$", "L4", All_samples_Merged$display_sample)

# Rename PBMC to CD4T_lab
All_samples_Merged$display_sample <- gsub("^PBMC$", "CD4T_lab", All_samples_Merged$display_sample)

2. Combining Contigs into Clones

There are varying definitions of clones in the literature. For the purposes of scRepertoire, we define a clone as cells with shared/trackable complementarity-determining region 3 (CDR3) sequences. Within this definition, one might use amino acid (aa) sequences of one or both chains to define a clone. Alternatively, we could use nucleotide (nt) or the V(D)JC genes (genes) to define a clone. The latter, genes, would be a more permissive definition of “clones,” as multiple amino acid or nucleotide sequences can result from the same gene combination. Another option to define a clone is the use of the V(D)JC and nucleotide sequence (strict). scRepertoire allows for the use of all these definitions of clones and enables users to select both or individual chains to examine.

2.1 Combining Contigs into Clones

contig_list <- list(L1, L2, L3_B, L4_B,L5,L6, L7)
contig_list_with_PBMC <- list(L1, L2, L3_B, L4_B,L5,L6, L7, PBMC)

combined.TCR <- combineTCR(contig_list,
                           samples = c("L1", "L2", "L3", "L4", "L5", "L6", "L7"),
                           removeNA = FALSE,
                           removeMulti = FALSE,
                           filterMulti = FALSE)

combined.TCR_with_PBMC <- combineTCR(contig_list_with_PBMC, 
                           samples = c("L1", "L2", "L3", "L4", "L5", "L6", "L7", 
                                       "CD4T_lab"),
                           removeNA = FALSE, 
                           removeMulti = FALSE, 
                           filterMulti = FALSE)

 exportClones(combined.TCR,
              write.file = TRUE,
              dir = "TCR_analysis-7-11-2025/",
             file.name = "clones.csv")

exportClones(combined.TCR_with_PBMC,
             write.file = TRUE,
             dir = "TCR_analysis-7-11-2025/",
             file.name = "clones_with_PBMC.csv")

head(combined.TCR[[1]])

head(combined.TCR_with_PBMC[[1]])

# Combine all data frames in the list into a single data frame
combined_TCR_df <- do.call(rbind, combined.TCR)

# Write the combined data frame to a CSV file
write.csv(combined_TCR_df, file = "TCR_analysis-7-11-2025/combined_TCR.csv", row.names = FALSE)


# Combine all data frames in the list into a single data frame
combined_TCR_with_PBMC_df <- do.call(rbind, combined.TCR_with_PBMC)

# Write the combined data frame to a CSV file
write.csv(combined_TCR_with_PBMC_df, file = "TCR_analysis-7-11-2025/combined_TCR_with_PBMC.csv", row.names = FALSE)

2.2 Write the fasta file for phylogeny

# Define the FASTA file for alpha chains
alpha_fasta_file <- "TCR_analysis-7-11-2025/Phylogeny/cdr3_alpha_sequences.fasta"

# Open a connection to the file
f <- file(alpha_fasta_file, open = "w")

# Loop through the data frame
for (i in seq_len(nrow(combined_TCR_with_PBMC_df))) {
  if (!is.na(combined_TCR_with_PBMC_df$cdr3_nt1[i])) {
    # Write the header and sequence
    writeLines(paste0(">", combined_TCR_df$sample[i], "_alpha"), f)
    writeLines(combined_TCR_df$cdr3_nt1[i], f)
  }
}

# Close the connection
close(f)

cat("Alpha chain FASTA file created at:", alpha_fasta_file, "\n")

Alpha chain FASTA file created at: TCR_analysis-7-11-2025/Phylogeny/cdr3_alpha_sequences.fasta 

# Define the FASTA file for beta chains
beta_fasta_file <- "TCR_analysis-7-11-2025/Phylogeny/cdr3_beta_sequences.fasta"

# Open a connection to the file
f <- file(beta_fasta_file, open = "w")

# Loop through the data frame
for (i in seq_len(nrow(combined_TCR_with_PBMC_df))) {
  if (!is.na(combined_TCR_with_PBMC_df$cdr3_nt2[i])) {
    # Write the header and sequence
    writeLines(paste0(">", combined_TCR_with_PBMC_df$sample[i], "_beta"), f)
    writeLines(combined_TCR_df$cdr3_nt2[i], f)
  }
}

# Close the connection
close(f)

cat("Beta chain FASTA file created at:", beta_fasta_file, "\n")

Beta chain FASTA file created at: TCR_analysis-7-11-2025/Phylogeny/cdr3_beta_sequences.fasta 

3. Basic Clonal Visualizations

3.1. clonalQuant: Quantifying Unique Clones

#clonalQuant
 clonalQuant(combined.TCR, 
            cloneCall="strict", 
            chain = "both", 
            scale = TRUE)

clonalQuant(combined.TCR_with_PBMC, 
            cloneCall="strict", 
            chain = "both", 
            scale = TRUE)

#clonalQuant
 clonalQuant(combined.TCR, 
            cloneCall="strict", 
            chain = "both", 
            scale = FALSE)

clonalQuant(combined.TCR_with_PBMC, 
            cloneCall="strict", 
            chain = "both", 
            scale = FALSE)

3.2 clonalAbundance: Distribution of Clones by Size

#clonalAbundance

clonalAbundance(combined.TCR_with_PBMC, 
                cloneCall = "gene",  palette = "Zissou 1",
                scale = FALSE)

  
clonalAbundance(combined.TCR_with_PBMC, 
                cloneCall = "gene", palette = "Zissou 1",
                scale = TRUE)

clonalAbundance(combined.TCR, 
                cloneCall = "gene",  
                scale = FALSE)

 

clonalAbundance(combined.TCR, 
                cloneCall = "gene", 
                scale = TRUE)

3.3 clonalLength: Distribution of Sequence Lengths

#clonalLength
clonalLength(combined.TCR_with_PBMC, 
             cloneCall="nt", 
             chain = "both") 

clonalLength(combined.TCR_with_PBMC, 
             cloneCall="aa", 
             chain = "both") 

clonalLength(combined.TCR, 
             cloneCall="nt", 
             chain = "both") 

clonalLength(combined.TCR, 
             cloneCall="aa", 
             chain = "both") 

#TRA
clonalLength(combined.TCR_with_PBMC, 
             cloneCall="nt", 
             chain = "TRA", 
             scale = TRUE) 

clonalLength(combined.TCR_with_PBMC, 
             cloneCall="aa", 
             chain = "TRA", 
             scale = TRUE) 

clonalLength(combined.TCR, 
             cloneCall="nt", 
             chain = "TRA",  
             scale = TRUE) 

clonalLength(combined.TCR, 
             cloneCall="aa", 
             chain = "TRA", 
             scale = TRUE) 

#TRB
clonalLength(combined.TCR_with_PBMC, 
             cloneCall="nt", 
             chain = "TRB", 
             scale = TRUE) 

clonalLength(combined.TCR_with_PBMC, 
             cloneCall="aa", 
             chain = "TRB", 
             scale = TRUE) 

clonalLength(combined.TCR, 
             cloneCall="nt", 
             chain = "TRB",  
             scale = TRUE) 

clonalLength(combined.TCR, 
             cloneCall="aa", 
             chain = "TRB", 
             scale = TRUE) 

3.5 clonalCompare: Clonal Dynamics Between Categorical Variables

# clonalCompare
clonalCompare(combined.TCR, 
                  top.clones = 10, 
                  samples = c("L1", "L2"), 
                  cloneCall="aa", 
                  graph = "alluvial")

clonalCompare(combined.TCR, 
                  top.clones = 10, 
                  samples = c("L3", "L4"), 
                  cloneCall="aa", 
                  graph = "alluvial")

clonalCompare(combined.TCR, 
                  top.clones = 10, 
                  samples = c("L5", "L6","L7"), 
                  cloneCall="aa", 
                  graph = "alluvial")

clonalCompare(combined.TCR, 
                  top.clones = 10, 
                  samples = c("L5","L6"), 
                  cloneCall="aa", 
                  graph = "alluvial")

clonalCompare(combined.TCR, 
                  top.clones = 10, 
                  samples = c("L5","L7"), 
                  cloneCall="aa", 
                  graph = "alluvial")

clonalCompare(combined.TCR, 
                  top.clones = 10, 
                  samples = c("L6","L7"), 
                  cloneCall="aa", 
                  graph = "alluvial")

# clonalCompare


clonalCompare(combined.TCR_with_PBMC, 
                  top.clones = 10, 
                  samples = c("L1", "L2", "L3", "L4", "L5", "L6", "L7", "CD4T_lab"), 
                  cloneCall="aa", 
                  graph = "alluvial")

clonalCompare(combined.TCR, 
                  top.clones = 10, 
                  samples = c("L1", "L2", "L3", "L4", "L5", "L6", "L7"), 
                  cloneCall="aa", 
                  graph = "alluvial")

3.6 clonalScatter: Scatterplot of Two Variables

#clonalScatter

clonalScatter(combined.TCR, 
              cloneCall ="gene", 
              x.axis = "L1", 
              y.axis = "L2",
              dot.size = "total",
              graph = "proportion")

clonalScatter(combined.TCR, 
              cloneCall ="gene", 
              x.axis = "L3", 
              y.axis = "L4",
              dot.size = "total",
              graph = "proportion")

clonalScatter(combined.TCR, 
              cloneCall ="gene", 
              x.axis = "L5", 
              y.axis = "L6",
              dot.size = "total",
              graph = "proportion")

clonalScatter(combined.TCR, 
              cloneCall ="gene", 
              x.axis = "L5", 
              y.axis = "L7",
              dot.size = "total",
              graph = "proportion")

clonalScatter(combined.TCR_with_PBMC, 
              cloneCall ="gene", 
              x.axis = "L6", 
              y.axis = "L7",
              dot.size = "total",
              graph = "proportion")

4. Visualizing Clonal Dynamics

4.1 clonalHomeostasis: Examining Clonal Space

# Visualizing Clonal Dynamics

# clonalHomeostasis

clonalHomeostasis(combined.TCR, 
                  cloneCall = "gene")

clonalHomeostasis(combined.TCR, 
                  cloneCall = "gene",
                  cloneSize = c(Rare = 0.001, Small = 0.01, Medium = 0.1, Large = 0.3, Hyperexpanded =1))

clonalHomeostasis(combined.TCR, 
                  cloneCall = "gene",
                  cloneSize = c(Rare = 0.1, Small = 1, Medium = 10, Large = 30, Hyperexpanded =100))

4.2 clonalProportion: Examining Space Occupied by Ranks of Clones

# clonalProportion

clonalProportion(combined.TCR, 
                 cloneCall = "gene") 

clonalProportion(combined.TCR, 
                 cloneCall = "nt",
                 clonalSplit = c(1, 5, 10, 100, 1000, 10000)) 

clonalProportion(combined.TCR, 
                 cloneCall = "nt",
                 clonalSplit = c(10, 100, 1000, 10000, 30000, 100000)) 

NA
NA

5. Summarizing Repertoires

# Summarizing Repertoires
# percentAA

percentAA(combined.TCR, 
          chain = "TRB", 
          aa.length = 20)

# positionalEntropy

positionalEntropy(combined.TCR, 
                  chain = "TRB", 
                  aa.length = 20)

positionalEntropy(combined.TCR_with_PBMC, 
                  chain = "TRB", 
                  aa.length = 20)

# positionalProperty

positionalProperty(combined.TCR[c(1,2)], 
                  chain = "TRB", 
                  aa.length = 20, 
                  method = "atchleyFactors") + 
  scale_color_manual(values = hcl.colors(5, "inferno")[c(2,4)])

NA
NA

# percentGeneUsage
vizGenes(combined.TCR, 
         x.axis = "TRBV",
         y.axis = NULL,
         plot = "barplot",  
         summary.fun = "proportion")

vizGenes(combined.TCR, 
        x.axis = "TRAV",
        y.axis = NULL,
        plot = "barplot",  
        summary.fun = "proportion")

vizGenes(combined.TCR_with_PBMC, 
         x.axis = "TRBV",
         y.axis = NULL,
         plot = "barplot",  
         summary.fun = "proportion")

vizGenes(combined.TCR_with_PBMC, 
        x.axis = "TRAV",
        y.axis = NULL,
        plot = "barplot",  
        summary.fun = "proportion")

vizGenes(combined.TCR[c(1,2,3,4,5,6,7)], 
         x.axis = "TRBV",
         y.axis = "TRBJ",
         plot = "heatmap",  
          summary.fun = "percent")

vizGenes(combined.TCR[c(1,2,3,4,5,6,7)], 
         x.axis = "TRAV",
         y.axis = "TRAJ",
         plot = "heatmap",  
          summary.fun = "percent")

vizGenes(combined.TCR[c(1,2)], 
         x.axis = "TRBV",
         y.axis = "TRAV",
         plot = "heatmap", 
          summary.fun = "percent")

vizGenes(combined.TCR[c(3,4)], 
         x.axis = "TRBV",
         y.axis = "TRAV",
         plot = "heatmap", 
          summary.fun = "percent")

vizGenes(combined.TCR[c(5,6,7)], 
         x.axis = "TRBV",
         y.axis = "TRAV",
         plot = "heatmap", 
          summary.fun = "percent")

# percentGenes: Quantifying Single Gene Usage
percentGenes(combined.TCR, 
             chain = "TRB", 
             gene = "Vgene",
              summary.fun = "percent")

percentGenes(combined.TCR, 
             chain = "TRA", 
             gene = "Vgene",
              summary.fun = "percent")

percentGenes(combined.TCR, 
             chain = "TRB", 
             gene = "Jgene")

percentGenes(combined.TCR, 
             chain = "TRA", 
             gene = "Jgene",
              summary.fun = "percent")

# percentGenes: Quantifying Single Gene Usage
percentGenes(combined.TCR_with_PBMC, 
             chain = "TRB", 
             gene = "Vgene",
              summary.fun = "percent")

percentGenes(combined.TCR_with_PBMC, 
             chain = "TRA", 
             gene = "Vgene",
              summary.fun = "percent")

percentGenes(combined.TCR_with_PBMC, 
             chain = "TRB", 
             gene = "Jgene")

percentGenes(combined.TCR_with_PBMC, 
             chain = "TRA", 
             gene = "Jgene",
              summary.fun = "percent")

# percentVJ
# Quantify the proportion of V and J gene usage with percentVJ(). Like percentGenes(), this function will quantify the percentage of V and J paired together across individual repertoires. The output can be visualized using a heatmap or as input for further dimensional reduction.

percentVJ(combined.TCR_with_PBMC, 
          chain = "TRB",
           summary.fun = "percent")

percentVJ(combined.TCR_with_PBMC, 
          chain = "TRA",
           summary.fun = "percent")

# percentKmer

percentKmer(combined.TCR_with_PBMC, 
            cloneCall = "aa",
            chain = "TRB", 
            motif.length = 3, 
            top.motifs = 25)

percentKmer(combined.TCR_with_PBMC, 
            cloneCall = "aa",
            chain = "TRA", 
            motif.length = 3, 
            top.motifs = 25)

percentKmer(combined.TCR_with_PBMC, 
            cloneCall = "nt",
            chain = "TRB", 
            motif.length = 3, 
            top.motifs = 25)

percentKmer(combined.TCR_with_PBMC, 
            cloneCall = "nt",
            chain = "TRA", 
            motif.length = 3, 
            top.motifs = 25)

6. Comparing Clonal Diversity and Overlap

# clonalDiversity
clonalDiversity(combined.TCR_with_PBMC, 
                cloneCall = "gene")

clonalDiversity(combined.TCR, 
                cloneCall = "gene")

#clonalSizeDistribution
clonalSizeDistribution(combined.TCR, 
                       cloneCall = "aa", 
                       method= "ward.D2")

clonalSizeDistribution(combined.TCR, 
                       cloneCall = "nt", 
                       method= "ward.D2")

#clonalSizeDistribution
clonalSizeDistribution(combined.TCR_with_PBMC, 
                       cloneCall = "aa", 
                       method= "ward.D2")

clonalSizeDistribution(combined.TCR_with_PBMC, 
                       cloneCall = "nt", 
                       method= "ward.D2")

#clonalOverlap

clonalOverlap(combined.TCR_with_PBMC, 
              cloneCall = "strict", 
              method = "morisita")

clonalOverlap(combined.TCR_with_PBMC, 
              cloneCall = "strict", 
              method = "raw")

#clonalOverlap

clonalOverlap(combined.TCR, 
              cloneCall = "strict", 
              method = "morisita")

clonalOverlap(combined.TCR, 
              cloneCall = "strict", 
              method = "raw")

clonalRarefaction(combined.TCR,
                  plot.type = 1,
                  hill.numbers = 0,
                  n.boots = 2)

clonalRarefaction(combined.TCR_with_PBMC,
                  plot.type = 1,
                  hill.numbers = 0,
                  n.boots = 2)

clonalRarefaction(combined.TCR,
                  plot.type = 2,
                  hill.numbers = 0,
                  n.boots = 2)

clonalRarefaction(combined.TCR_with_PBMC,
                  plot.type = 2,
                  hill.numbers = 0,
                  n.boots = 2)

clonalRarefaction(combined.TCR,
                  plot.type = 3,
                  hill.numbers = 0,
                  n.boots = 2)

clonalRarefaction(combined.TCR_with_PBMC,
                  plot.type = 3,
                  hill.numbers = 0,
                  n.boots = 2)

# Rarefaction using Shannon Diversity (q = 1)
clonalRarefaction(combined.TCR,
                  plot.type = 1,
                  hill.numbers = 1,
                  n.boots = 2)

clonalRarefaction(combined.TCR_with_PBMC,
                  plot.type = 1,
                  hill.numbers = 1,
                  n.boots = 2)