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

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.

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_B", "L4_B", "L5", "L6", "L7"),
                           removeNA = FALSE,
                           removeMulti = FALSE,
                           filterMulti = FALSE)

combined.TCR_with_PBMC <- combineTCR(contig_list_with_PBMC, 
                           samples = c("L1", "L2", "L3_B", "L4_B", "L5", "L6", "L7", 
                                       "PBMC"),
                           removeNA = FALSE, 
                           removeMulti = FALSE, 
                           filterMulti = FALSE)

 exportClones(combined.TCR,
              write.file = TRUE,
              dir = "TCR_analysis_Harmony/",
             file.name = "clones.csv")

exportClones(combined.TCR_with_PBMC,
             write.file = TRUE,
             dir = "TCR_analysis_Harmony/",
             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_Harmony/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_Harmony/combined_TCR_with_PBMC.csv", row.names = FALSE)

Remove TCR VDJ genes

DefaultAssay(All_samples_Merged) <- "SCT"

TCR <- All_samples_Merged



# Combining Clones and Single-Cell Objects

#Getting a sample of a Seurat object
TCR <- All_samples_Merged

# Check the first 10 variable features before removal
VariableFeatures(TCR)[1:10]

 [1] "CCL17"   "CCL1"    "CCL4"    "TRBV7-2" "PPBP"    "CCL3"    "OASL"    "IFIT2"   "GZMB"    "XCL1"   

# Remove TCR VDJ genes
TCR <- quietTCRgenes(TCR)

# Check the first 10 variable features after removal
VariableFeatures(TCR)[1:10]

 [1] "CCL17"    "CCL1"     "CCL4"     "PPBP"     "CCL3"     "OASL"     "IFIT2"    "GZMB"     "XCL1"     "HLA-DRB1"

3. Combining Clones and Single-Cell Objects

DefaultAssay(TCR) <- "SCT"

# If All_samples_Merged is already loaded in the environment:
TCR <- All_samples_Merged

# Define color palette
colorblind_vector <- hcl.colors(n=7, palette = "inferno", fixup = TRUE)

# Combine expression with filtering of NA clonotype cells
TCR <- combineExpression(
    combined.TCR_with_PBMC,
    TCR,
    cloneCall = "gene",
    group.by = "sample",
    proportion = TRUE,
    filterNA = TRUE  # This will exclude cells without clonotype info
)

# You no longer need the manual barcode matching or NA replacements



# Plot UMAP colored by cloneSize
DimPlot(TCR, group.by = "cloneSize", reduction = "umap") +
    scale_color_manual(values = rev(colorblind_vector[c(1, 3, 4, 5, 7)]))

DimPlot(TCR, group.by = "cloneSize", reduction = "umap")

DimPlot(TCR, group.by = "cloneSize", reduction = "umap") +
    scale_color_manual(values=rev(colorblind_vector[c(1,3,4,5,6)]))

#Define color palette 
colorblind_vector <- hcl.colors(n=9, palette = "inferno", fixup = TRUE)

Seurat::DimPlot(TCR, group.by = "cloneSize", reduction = "umap") +
    scale_color_manual(values=rev(colorblind_vector[c(1,3,4,5,7)]))

TCR <- combineExpression(combined.TCR_with_PBMC, 
                                   TCR, 
                                   cloneCall="gene", 
                                   group.by = "sample", 
                                   proportion = FALSE, 
                                   cloneSize=c(Single=1, Small=5, Medium=20, Large=100, Hyperexpanded=500))

Seurat::DimPlot(TCR, group.by = "cloneSize") +
    scale_color_manual(values=rev(colorblind_vector[c(1,3,4,5,7)]))

4. Visualizations for Single-Cell Objects

clonalOverlay(TCR, 
              reduction = "umap", 
              cutpoint = 1, 
              bins = 10, 
              facet.by = "orig.ident") + 
              guides(color = "none")

#clonalNetwork
#ggraph needs to be loaded due to issues with ggplot
library(ggraph)

clonalNetwork(TCR, 
              reduction = "umap", 
              group.by = "seurat_clusters",
              filter.clones = NULL,
              filter.identity = NULL,
              cloneCall = "aa")

#Examining Cluster 3 only
clonalNetwork(TCR, 
              reduction = "umap", 
              group.by = "seurat_clusters",
              filter.identity = 8,
              cloneCall = "aa")

shared.clones <- clonalNetwork(TCR, 
                               reduction = "umap", 
                               group.by = "seurat_clusters",
                               cloneCall = "aa", 
                               exportClones = TRUE)
head(shared.clones)

#ggraph needs to be loaded due to issues with ggplot
library(ggraph)

#No Identity filter
clonalNetwork(TCR, 
              reduction = "umap", 
              group.by = "seurat_clusters",
              filter.clones = NULL,
              filter.identity = NULL,
              cloneCall = "aa")

# clonalOccupy
#clonalOccupy
clonalOccupy(TCR, 
              x.axis = "seurat_clusters")

NA
NA
NA

# clonalOccupy
clonalOccupy(TCR, 
              x.axis = "orig.ident")

clonalOccupy(TCR, 
                     x.axis = "orig.ident", 
                     proportion = TRUE, 
                     label = FALSE)

# getCirclize

library(circlize)
library(scales)

circles <- getCirclize(TCR, 
                       group.by = "seurat_clusters")

#Just assigning the normal colors to each cluster
grid.cols <- hue_pal()(length(unique(TCR$seurat_clusters)))
names(grid.cols) <- unique(TCR$seurat_clusters)

#Graphing the chord diagram
chordDiagram(circles, self.link = 1, grid.col = grid.cols)

circles <- getCirclize(TCR, group.by = "orig.ident")

grid.cols <- scales::hue_pal()(length(unique(TCR@active.ident)))
names(grid.cols) <- levels(TCR@active.ident)

chordDiagram(circles, 
             self.link = 1, 
             grid.col = grid.cols)

5. Quantifying Clonal Bias

# # StartracDiversity # From the excellent work by Lei Zhang, et al., the authors introduce new methods for looking at clones by cellular origins and cluster identification. Their STARTRAC software has been adapted to work with scRepertoire and please read and cite their excellent work. # # In order to use the StartracDiversity() function, you will need to include the product of the combinedExpression() function. The second requirement is a column header in the meta data of the Seurat object that has tissue of origin. In the example data, type corresponds to the column “Type”, which includes the “P” and “T” classifiers. The indices can be subsetted for a specific patient or examined overall using the by variable. Importantly, the function uses only the strict definition of a clone of the VDJC genes and the CDR3 nucleotide sequence. # # The indices output includes: # # expa - Clonal Expansion # migr - Cross-tissue Migration # tran - State Transition

Idents(TCR) <- "seurat_clusters"

StartracDiversity(TCR, 
                  type = "orig.ident",
                  group.by = "orig.ident")

NA
NA

6. clonalBias

clonalBias(TCR, 
           cloneCall = "aa", 
           split.by = "orig.ident", 
           group.by = "seurat_clusters",
           n.boots = 10, 
           min.expand =5)

7. save the TCR object for future Use

 # saveRDS(TCR, file = "TCR_Seurat_with_cloneSize_12_11_2025_finalized.rds")