Starting work with scRepertoire.
Compiled: June 29, 2026
Package
scRepertoire 2.9.1
1 Introduction
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.
1.1 Loading Libraries
suppressMessages(library(scRepertoire))1.2 Loading Data
1.2.1 What data to load into scRepertoire?
scRepertoire primarily functions using the filtered_contig_annotations.csv output generated by 10x Genomics Cell Ranger. This file is typically found in the ./outs/ directory of your VDJ alignment folder.
1.2.2 Demonstrating Manual Data Loading (10x Genomics)
To prepare data for scRepertoire from 10x Genomics outputs, you would:
- Load the
filtered_contig_annotations.csvfile for each of your samples. - Combine these loaded data frames into a single list in your R environment.
S1 <- read.csv(".../Sample1/outs/filtered_contig_annotations.csv")
S2 <- read.csv(".../Sample2/outs/filtered_contig_annotations.csv")
S3 <- read.csv(".../Sample3/outs/filtered_contig_annotations.csv")
S4 <- read.csv(".../Sample4/outs/filtered_contig_annotations.csv")
contig_list <- list(S1, S2, S3, S4)1.3 Other alignment workflows
Beyond the default 10x Genomics Cell Ranger pipeline outputs, scRepertoire supports various other single-cell immune receptor sequencing formats through the loadContigs() function.
1.3.1 Supported Formats and Expected File Names
10X: “filtered_contig_annotations.csv”AIRR: “airr_rearrangement.tsv”BD: “Contigs_AIRR.tsv”Dandelion: “all_contig_dandelion.tsv”Immcantation: "_data.tsv" (or similar)JSON: “.json”MiXCR: “clones.tsv”ParseBio: “barcode_report.tsv”TRUST4: “barcode_report.tsv”WAT3R: “barcode_results.csv”
Key Parameter(s) for loadContigs()
input: A directory path containing your contig files (the function will recursively search) or a list/data frame of pre-loaded contig data.format: A string specifying the data format (e.g.,10X,TRUST4,WAT3R). If set to “auto”, the function will attempt to automatically detect the format based on file names or data structure.
Regardless of format, loadContigs() standardizes a full-length sequence (nucleotide) and sequence_aa (amino acid) column. AIRR-family inputs carry these natively, 10X and MiXCR reconstruct them from the framework and CDR regions, and CDR3-only formats leave them NA. These columns feed the optional retain.sequences argument of combineTCR() and combineBCR() described below, which preserves them through to the combined object for downstream tools like dowser.
You can provide loadContigs() with a directory where your sequencing experiments are located, and it will recursively load and process the contig data based on the file names:
# Directory example
contig.output <- c("~/Documents/MyExperiment")
contig.list <- loadContigs(input = contig.output,
format = "TRUST4")Alternatively, loadContigs() can be given a list of pre-loaded data frames and process the contig data based on the specified format:
# List of data frames example
S1 <- read.csv("~/Documents/MyExperiment/Sample1/outs/barcode_results.csv")
S2 <- read.csv("~/Documents/MyExperiment/Sample2/outs/barcode_results.csv")
S3 <- read.csv("~/Documents/MyExperiment/Sample3/outs/barcode_results.csv")
S4 <- read.csv("~/Documents/MyExperiment/Sample4/outs/barcode_results.csv")
contig.list <- list(S1, S2, S3, S4)
contig.list <- loadContigs(input = contig.list,
format = "WAT3R")1.4 Multiplexed Experiment
It is now easy to create the contig list from a multiplexed experiment by first generating a single-cell RNA object (either Seurat or Single Cell Experiment), loading the filtered contig file, and then using createHTOContigList(). This function will return a list separated by the group.by variable(s).
1.4.1 Important Considerations for createHTOContigList()
- This function depends on the match of barcodes between the single-cell object and contigs. If there is a prefix or different suffix added to the barcode, this will result in no contigs being recovered.
- It is currently recommended to perform this step before integration workflows, as integration commonly alters the barcodes.
- There is a multi.run variable that can be used on an integrated object. However, it assumes you have modified the barcodes with the Seurat pipeline (automatic addition of _# to the end), and your contig list is in the same order.
To create a contig list separated by HTO (Hash Tag Oligo) IDs from a single-cell object:
contigs <- read.csv(".../outs/filtered_contig_annotations.csv")
contig.list <- createHTOContigList(contigs,
Seurat.Obj,
group.by = "HTO_maxID")1.5 Example Data in scRepertoire
scRepertoire includes a built-in example dataset to demonstrate the functionality of the R package. This dataset consists of T cells derived from four patients with acute respiratory distress with paired peripheral-blood (B) and bronchoalveolar lavage (L), effectively creating 8 distinct runs for T cell receptor (TCR) enrichment. More information on the data set can be found in the corresponding manuscript.
The built-in example data is derived from the 10x Cell Ranger pipeline, so it is ready to go for downstream processing and analysis.
To load and preview the example data built into scRepertoire:
data("contig_list") #the data built into scRepertoire
head(contig_list[[1]])## barcode is_cell contig_id high_confidence length
## 1 AAACCTGAGTACGACG-1 True AAACCTGAGTACGACG-1_contig_1 True 500
## 2 AAACCTGAGTACGACG-1 True AAACCTGAGTACGACG-1_contig_2 True 478
## 4 AAACCTGCAACACGCC-1 True AAACCTGCAACACGCC-1_contig_1 True 506
## 5 AAACCTGCAACACGCC-1 True AAACCTGCAACACGCC-1_contig_2 True 470
## 6 AAACCTGCAGGCGATA-1 True AAACCTGCAGGCGATA-1_contig_1 True 558
## 7 AAACCTGCAGGCGATA-1 True AAACCTGCAGGCGATA-1_contig_2 True 505
## chain v_gene d_gene j_gene c_gene full_length productive
## 1 TRA TRAV25 None TRAJ20 TRAC True True
## 2 TRB TRBV5-1 None TRBJ2-7 TRBC2 True True
## 4 TRA TRAV38-2/DV8 None TRAJ52 TRAC True True
## 5 TRB TRBV10-3 None TRBJ2-2 TRBC2 True True
## 6 TRA TRAV12-1 None TRAJ9 TRAC True True
## 7 TRB TRBV9 None TRBJ2-2 TRBC2 True True
## cdr3 cdr3_nt
## 1 CGCSNDYKLSF TGTGGGTGTTCTAACGACTACAAGCTCAGCTTT
## 2 CASSLTDRTYEQYF TGCGCCAGCAGCTTGACCGACAGGACCTACGAGCAGTACTTC
## 4 CAYRSAQAGGTSYGKLTF TGTGCTTATAGGAGCGCGCAGGCTGGTGGTACTAGCTATGGAAAGCTGACATTT
## 5 CAISEQGKGELFF TGTGCCATCAGTGAACAGGGGAAAGGGGAGCTGTTTTTT
## 6 CVVSDNTGGFKTIF TGTGTGGTCTCCGATAATACTGGAGGCTTCAAAACTATCTTT
## 7 CASSVRRERANTGELFF TGTGCCAGCAGCGTAAGGAGGGAAAGGGCGAACACCGGGGAGCTGTTTTTT
## reads umis raw_clonotype_id raw_consensus_id
## 1 8344 4 clonotype123 clonotype123_consensus_2
## 2 65390 38 clonotype123 clonotype123_consensus_1
## 4 18372 8 clonotype124 clonotype124_consensus_1
## 5 34054 9 clonotype124 clonotype124_consensus_2
## 6 5018 2 clonotype1 clonotype1_consensus_2
## 7 25110 11 clonotype1 clonotype1_consensus_12 Combining Contigs into Clones
The combineTCR() function processes a list of TCR sequencing results, consolidating them to the level of individual cell barcodes. It handles potential issues with repeated barcodes by adding prefixes from samples and ID parameters. The output includes combined reads into clone calls by nucleotide sequence (CTnt), amino acid sequence (CTaa), VDJC gene sequence (CTgene), or a combination of nucleotide and gene sequence (CTstrict).
Key Parameter(s) for combineTCR()
input.data: A list of filtered contig annotations (e.g., filtered_contig_annotations.csv from 10x Cell Ranger) or outputs fromloadContigs().samples: Labels for your samples (recommended).ID: Additional sample labels (optional).removeNA: IfTRUE, removes any cell barcode with an NA value in at least one chain (default isFALSE).removeMulti: IfTRUE, removes any cell barcode with more than two immune receptor chains (default isFALSE).filterMulti: IfTRUE, isolates the top two expressed chains in cell barcodes with multiple chains (default isFALSE).filterNonproductive: IfTRUE, removes non-productive chains if the variable exists in the contig data (default isTRUE).retain.sequences: Retains full-length sequences as additional per-chain columns without altering the clone calls (default isFALSE). See Retaining Full-Length Sequences.
To combine TCR contigs from contig_list and apply sample prefixes:
combined.TCR <- combineTCR(contig_list,
samples = c("P17B", "P17L", "P18B", "P18L",
"P19B","P19L", "P20B", "P20L"),
removeNA = FALSE,
removeMulti = FALSE,
filterMulti = FALSE)
head(combined.TCR[[1]])## barcode sample TCR1 cdr3_aa1
## 1 P17B_AAACCTGAGTACGACG-1 P17B TRAV25.TRAJ20.TRAC CGCSNDYKLSF
## 3 P17B_AAACCTGCAACACGCC-1 P17B TRAV38-2/DV8.TRAJ52.TRAC CAYRSAQAGGTSYGKLTF
## 5 P17B_AAACCTGCAGGCGATA-1 P17B TRAV12-1.TRAJ9.TRAC CVVSDNTGGFKTIF
## 7 P17B_AAACCTGCATGAGCGA-1 P17B TRAV12-1.TRAJ9.TRAC CVVSDNTGGFKTIF
## 9 P17B_AAACGGGAGAGCCCAA-1 P17B TRAV20.TRAJ8.TRAC CAVRGEGFQKLVF
## 10 P17B_AAACGGGAGCGTTTAC-1 P17B TRAV12-1.TRAJ9.TRAC CVVSDNTGGFKTIF
## cdr3_nt1
## 1 TGTGGGTGTTCTAACGACTACAAGCTCAGCTTT
## 3 TGTGCTTATAGGAGCGCGCAGGCTGGTGGTACTAGCTATGGAAAGCTGACATTT
## 5 TGTGTGGTCTCCGATAATACTGGAGGCTTCAAAACTATCTTT
## 7 TGTGTGGTCTCCGATAATACTGGAGGCTTCAAAACTATCTTT
## 9 TGTGCTGTGCGAGGAGAAGGCTTTCAGAAACTTGTATTT
## 10 TGTGTGGTCTCCGATAATACTGGAGGCTTCAAAACTATCTTT
## TCR2 cdr3_aa2
## 1 TRBV5-1.None.TRBJ2-7.TRBC2 CASSLTDRTYEQYF
## 3 TRBV10-3.None.TRBJ2-2.TRBC2 CAISEQGKGELFF
## 5 TRBV9.None.TRBJ2-2.TRBC2 CASSVRRERANTGELFF
## 7 TRBV9.None.TRBJ2-2.TRBC2 CASSVRRERANTGELFF
## 9 <NA> <NA>
## 10 TRBV9.None.TRBJ2-2.TRBC2 CASSVRRERANTGELFF
## cdr3_nt2
## 1 TGCGCCAGCAGCTTGACCGACAGGACCTACGAGCAGTACTTC
## 3 TGTGCCATCAGTGAACAGGGGAAAGGGGAGCTGTTTTTT
## 5 TGTGCCAGCAGCGTAAGGAGGGAAAGGGCGAACACCGGGGAGCTGTTTTTT
## 7 TGTGCCAGCAGCGTAAGGAGGGAAAGGGCGAACACCGGGGAGCTGTTTTTT
## 9 <NA>
## 10 TGTGCCAGCAGCGTAAGGAGGGAAAGGGCGAACACCGGGGAGCTGTTTTTT
## CTgene
## 1 TRAV25.TRAJ20.TRAC_TRBV5-1.None.TRBJ2-7.TRBC2
## 3 TRAV38-2/DV8.TRAJ52.TRAC_TRBV10-3.None.TRBJ2-2.TRBC2
## 5 TRAV12-1.TRAJ9.TRAC_TRBV9.None.TRBJ2-2.TRBC2
## 7 TRAV12-1.TRAJ9.TRAC_TRBV9.None.TRBJ2-2.TRBC2
## 9 TRAV20.TRAJ8.TRAC_NA
## 10 TRAV12-1.TRAJ9.TRAC_TRBV9.None.TRBJ2-2.TRBC2
## CTnt
## 1 TGTGGGTGTTCTAACGACTACAAGCTCAGCTTT_TGCGCCAGCAGCTTGACCGACAGGACCTACGAGCAGTACTTC
## 3 TGTGCTTATAGGAGCGCGCAGGCTGGTGGTACTAGCTATGGAAAGCTGACATTT_TGTGCCATCAGTGAACAGGGGAAAGGGGAGCTGTTTTTT
## 5 TGTGTGGTCTCCGATAATACTGGAGGCTTCAAAACTATCTTT_TGTGCCAGCAGCGTAAGGAGGGAAAGGGCGAACACCGGGGAGCTGTTTTTT
## 7 TGTGTGGTCTCCGATAATACTGGAGGCTTCAAAACTATCTTT_TGTGCCAGCAGCGTAAGGAGGGAAAGGGCGAACACCGGGGAGCTGTTTTTT
## 9 TGTGCTGTGCGAGGAGAAGGCTTTCAGAAACTTGTATTT_NA
## 10 TGTGTGGTCTCCGATAATACTGGAGGCTTCAAAACTATCTTT_TGTGCCAGCAGCGTAAGGAGGGAAAGGGCGAACACCGGGGAGCTGTTTTTT
## CTaa
## 1 CGCSNDYKLSF_CASSLTDRTYEQYF
## 3 CAYRSAQAGGTSYGKLTF_CAISEQGKGELFF
## 5 CVVSDNTGGFKTIF_CASSVRRERANTGELFF
## 7 CVVSDNTGGFKTIF_CASSVRRERANTGELFF
## 9 CAVRGEGFQKLVF_NA
## 10 CVVSDNTGGFKTIF_CASSVRRERANTGELFF
## CTstrict
## 1 TRAV25.TRAJ20.TRAC;TGTGGGTGTTCTAACGACTACAAGCTCAGCTTT_TRBV5-1.None.TRBJ2-7.TRBC2;TGCGCCAGCAGCTTGACCGACAGGACCTACGAGCAGTACTTC
## 3 TRAV38-2/DV8.TRAJ52.TRAC;TGTGCTTATAGGAGCGCGCAGGCTGGTGGTACTAGCTATGGAAAGCTGACATTT_TRBV10-3.None.TRBJ2-2.TRBC2;TGTGCCATCAGTGAACAGGGGAAAGGGGAGCTGTTTTTT
## 5 TRAV12-1.TRAJ9.TRAC;TGTGTGGTCTCCGATAATACTGGAGGCTTCAAAACTATCTTT_TRBV9.None.TRBJ2-2.TRBC2;TGTGCCAGCAGCGTAAGGAGGGAAAGGGCGAACACCGGGGAGCTGTTTTTT
## 7 TRAV12-1.TRAJ9.TRAC;TGTGTGGTCTCCGATAATACTGGAGGCTTCAAAACTATCTTT_TRBV9.None.TRBJ2-2.TRBC2;TGTGCCAGCAGCGTAAGGAGGGAAAGGGCGAACACCGGGGAGCTGTTTTTT
## 9 TRAV20.TRAJ8.TRAC;TGTGCTGTGCGAGGAGAAGGCTTTCAGAAACTTGTATTT_NA;NA
## 10 TRAV12-1.TRAJ9.TRAC;TGTGTGGTCTCCGATAATACTGGAGGCTTCAAAACTATCTTT_TRBV9.None.TRBJ2-2.TRBC2;TGTGCCAGCAGCGTAAGGAGGGAAAGGGCGAACACCGGGGAGCTGTTTTTTcombineTCR() is the essential first step for organizing raw TCR sequencing data into a structured format for scRepertoire analyses. It allows for flexible handling of single and paired chains, barcode disambiguation, and initial filtering, producing a list of data frames where each row represents a single cell and its associated TCR clonotypes.
2.1 combineBCR
The combineBCR() function is the primary tool for processing raw B cell receptor contig data into a format ready for analysis. It is analogous to combineTCR() but includes specialized logic for handling the complexities of BCRs, such as somatic hypermutation. The function consolidates contigs into a single data frame per sample, with each row representing a unique cell.
By default (call.related.clones = TRUE), combineBCR() groups B cells into clones based on the similarity of their CDR3 sequences.
2.1.2 Advanced: Grouping by Alignment
For more biological accuracy—specifically when analyzing amino acid sequences—you can use alignment metrics. This allows the clustering to penalize conservative amino acid changes less than radical changes.
Here we use Needleman-Wunsch alignment with the BLOSUM80 substitution matrix:
combined.BCR.aligned <- combineBCR(BCR.contigs,
samples = "Patient1",
sequence = "aa",
dist_type = "nw",
dist_mat = "BLOSUM80",
threshold = 0.85)
head(combined.BCR.aligned[[1]][, c("barcode", "CTstrict", "IGH", "cdr3_aa1")])## barcode
## 1 Patient1_AAACCTGAGGGCACTA-1
## 2 Patient1_AAACCTGAGTACGTTC-1
## 3 Patient1_AAACCTGAGTCCGGTC-1
## 4 Patient1_AAACCTGCACCAGGTC-1
## 5 Patient1_AAACCTGGTAAATGAC-1
## 6 Patient1_AAACGGGTCACCTCGT-1
## CTstrict
## 1 NA_IGLV1-51.CGTWDSSLSAGVF
## 2 IGHV3-64.CAKSYSRDLPRYFGSW_IGKV3-15.CQQYSNWPLTF
## 3 NA_IGKV3-20.CQQYGDSLPF
## 4 IGHV1-69-2.CVRDGAGHYGSGIGYYGMDVW_IGLV2-14.CSSYISTSTLEVLF
## 5 IGHV1-46.CSRGRLPMPAAGSGLVSW_IGKV3-15.CQHYHNWPPYTF
## 6 NA_IGKV3-11.CQQRSNWPLTF
## IGH cdr3_aa1
## 1 <NA> <NA>
## 2 IGHV3-64.IGHD6-6.IGHJ4.IGHG1 CAKSYSRDLPRYFGSW
## 3 <NA> <NA>
## 4 IGHV1-69-2.IGHD3-10.IGHJ6.IGHM CVRDGAGHYGSGIGYYGMDVW
## 5 IGHV1-46.IGHD6-13.IGHJ5.IGHA2 CSRGRLPMPAAGSGLVSW
## 6 <NA> <NA>2.1.3 Filtering and Cleaning Data
combineBCR() includes several arguments to filter and clean the contig data during processing.
* filterNonproductive = TRUE (Default): Removes any contigs that are not classified as productive, ensuring that only functional receptor chains are included in the analysis.
* filterMulti = TRUE (Default): For cells with more than one heavy or light chain detected, this automatically selects the chain with the highest UMI count (read abundance) and discards the others. This helps resolve cellular multiplets or technical artifacts.
Here is an example applying these filters (though they are on by default):
cleaned.BCR <- combineBCR(BCR.contigs,
samples = "Patient1",
filterNonproductive = TRUE,
filterMulti = TRUE)
head(cleaned.BCR[[1]])## barcode sample IGH
## 1 Patient1_AAACCTGAGGGCACTA-1 Patient1 <NA>
## 2 Patient1_AAACCTGAGTACGTTC-1 Patient1 IGHV3-64.IGHD6-6.IGHJ4.IGHG1
## 3 Patient1_AAACCTGAGTCCGGTC-1 Patient1 <NA>
## 4 Patient1_AAACCTGCACCAGGTC-1 Patient1 IGHV1-69-2.IGHD3-10.IGHJ6.IGHM
## 5 Patient1_AAACCTGGTAAATGAC-1 Patient1 IGHV1-46.IGHD6-13.IGHJ5.IGHA2
## 6 Patient1_AAACGGGTCACCTCGT-1 Patient1 <NA>
## cdr3_aa1
## 1 <NA>
## 2 CAKSYSRDLPRYFGSW
## 3 <NA>
## 4 CVRDGAGHYGSGIGYYGMDVW
## 5 CSRGRLPMPAAGSGLVSW
## 6 <NA>
## cdr3_nt1
## 1 <NA>
## 2 TGTGCGAAATCGTATAGCAGAGACCTGCCGCGGTACTTTGGCTCCTGG
## 3 <NA>
## 4 TGTGTGAGAGATGGGGCGGGTCACTATGGTTCGGGGATAGGCTACTACGGTATGGACGTCTGG
## 5 TGTTCGAGGGGTCGTCTCCCGATGCCAGCAGCTGGCAGTGGTTTGGTCTCCTGG
## 6 <NA>
## IGLC cdr3_aa2
## 1 IGLV1-51.IGLJ3.IGLC2 CGTWDSSLSAGVF
## 2 IGKV3-15.IGKJ4.IGKC CQQYSNWPLTF
## 3 IGKV3-20.IGKJ4.IGKC CQQYGDSLPF
## 4 IGLV2-14.IGLJ2.IGLC2 CSSYISTSTLEVLF
## 5 IGKV3-15.IGKJ2.IGKC CQHYHNWPPYTF
## 6 IGKV3-11.IGKJ1.IGKC CQQRSNWPLTF
## cdr3_nt2
## 1 TGCGGAACATGGGATAGCAGCCTGAGTGCTGGCGTGTTC
## 2 TGTCAGCAGTATAGTAACTGGCCGCTCACTTTC
## 3 TGTCAACAGTATGGTGACTCTCTCCCTTTC
## 4 TGCAGTTCATATATAAGTACCAGCACTCTCGAGGTCCTATTC
## 5 TGTCAGCACTATCATAACTGGCCTCCGTACACTTTT
## 6 TGTCAGCAGCGTAGCAACTGGCCTCTGACGTTC
## CTgene
## 1 NA_IGLV1-51.IGLJ3.IGLC2
## 2 IGHV3-64.IGHD6-6.IGHJ4.IGHG1_IGKV3-15.IGKJ4.IGKC
## 3 NA_IGKV3-20.IGKJ4.IGKC
## 4 IGHV1-69-2.IGHD3-10.IGHJ6.IGHM_IGLV2-14.IGLJ2.IGLC2
## 5 IGHV1-46.IGHD6-13.IGHJ5.IGHA2_IGKV3-15.IGKJ2.IGKC
## 6 NA_IGKV3-11.IGKJ1.IGKC
## CTnt
## 1 NA_TGCGGAACATGGGATAGCAGCCTGAGTGCTGGCGTGTTC
## 2 TGTGCGAAATCGTATAGCAGAGACCTGCCGCGGTACTTTGGCTCCTGG_TGTCAGCAGTATAGTAACTGGCCGCTCACTTTC
## 3 NA_TGTCAACAGTATGGTGACTCTCTCCCTTTC
## 4 TGTGTGAGAGATGGGGCGGGTCACTATGGTTCGGGGATAGGCTACTACGGTATGGACGTCTGG_TGCAGTTCATATATAAGTACCAGCACTCTCGAGGTCCTATTC
## 5 TGTTCGAGGGGTCGTCTCCCGATGCCAGCAGCTGGCAGTGGTTTGGTCTCCTGG_TGTCAGCACTATCATAACTGGCCTCCGTACACTTTT
## 6 NA_TGTCAGCAGCGTAGCAACTGGCCTCTGACGTTC
## CTaa
## 1 NA_CGTWDSSLSAGVF
## 2 CAKSYSRDLPRYFGSW_CQQYSNWPLTF
## 3 NA_CQQYGDSLPF
## 4 CVRDGAGHYGSGIGYYGMDVW_CSSYISTSTLEVLF
## 5 CSRGRLPMPAAGSGLVSW_CQHYHNWPPYTF
## 6 NA_CQQRSNWPLTF
## CTstrict
## 1 NA_IGLV1-51.TGCGGAACATGGGATAGCAGCCTGAGTGCTGGCGTGTTC
## 2 IGHV3-64.TGTGCGAAATCGTATAGCAGAGACCTGCCGCGGTACTTTGGCTCCTGG_IGKV3-15.TGTCAGCAGTATAGTAACTGGCCGCTCACTTTC
## 3 NA_IGKV3-20.TGTCAACAGTATGGTGACTCTCTCCCTTTC
## 4 IGHV1-69-2.TGTGTGAGAGATGGGGCGGGTCACTATGGTTCGGGGATAGGCTACTACGGTATGGACGTCTGG_IGLV2-14.TGCAGTTCATATATAAGTACCAGCACTCTCGAGGTCCTATTC
## 5 IGHV1-46.TGTTCGAGGGGTCGTCTCCCGATGCCAGCAGCTGGCAGTGGTTTGGTCTCCTGG_IGKV3-15.TGTCAGCACTATCATAACTGGCCTCCGTACACTTTT
## 6 NA_IGKV3-11.TGTCAGCAGCGTAGCAACTGGCCTCTGACGTTCcombineBCR() is designed for processing B cell repertoire data, going beyond simple contig aggregation to incorporate advanced clustering based on CDR3 sequence similarity. This enables the identification of clonally related B cells, crucial for studying B cell development, affinity maturation, and humoral immune responses. Its filtering options further ensure the quality and interpretability of the processed data.
2.2 Retaining Full-Length Sequences
By default, combineTCR() and combineBCR() keep only the CDR3 in the CTnt and CTaa columns. Some downstream tools need the full-length receptor sequence. B cell lineage tracing with dowser, for example, reconstructs ancestral sequences from IMGT-gapped alignments rather than the CDR3 alone.
The retain.sequences argument carries the standardized full-length columns produced by loadContigs() (sequence/sequence_aa, and sequence_alignment/germline_alignment when present) into the combined object. It is additive and never touches CTaa, CTnt, CTgene, or CTstrict.
retain.sequences = FALSE(default): legacy output with no extra columns.retain.sequences = TRUE: retainssequence_nt1/sequence_nt2andsequence_aa1/sequence_aa2per chain.- A character vector retains specific source columns, for example
c("sequence_alignment", "germline_alignment")fordowser.
# Retain the standardized full-length sequences for both chains
combined.BCR.seq <- combineBCR(BCR.contigs,
samples = "Patient1",
retain.sequences = TRUE)
head(combined.BCR.seq[[1]][, c("barcode", "sequence_nt1", "sequence_aa1")])These columns flow through combineExpression() into single-cell metadata, and exportClones(format = "airr") writes them into the matching AIRR rearrangement fields. For B cell lineage analysis, exportDowser() (or exportClones(format = "dowser")) reshapes a combined object built with retain.sequences = c("sequence_alignment", "germline_alignment") into a data frame ready for dowser::formatClones().
3 Additional Processing
3.1 addVariable: Adding Variables for Plotting
What if there are more variables to add than just sample and ID? We can add them by using the addVariable() function. For each element, the function will add a column (labeled by variable.name) with the variable. The length of the variables parameter needs to match the length of the combined object.
Key Parameter(s) for addVariable()
variable.name: A character string that defines the new variable to add (e.g., “Type”, “Treatment”).variables: A character vector defining the desired column value for each list element. Its length must match the number of elements in the input.data list.
As an example, here we add the Type in which the samples were processed and sequenced to the combined.TCR object:
combined.TCR <- addVariable(combined.TCR,
variable.name = "Type",
variables = rep(c("B", "L"), 4))
head(combined.TCR[[1]])## barcode sample TCR1 cdr3_aa1
## 1 P17B_AAACCTGAGTACGACG-1 P17B TRAV25.TRAJ20.TRAC CGCSNDYKLSF
## 3 P17B_AAACCTGCAACACGCC-1 P17B TRAV38-2/DV8.TRAJ52.TRAC CAYRSAQAGGTSYGKLTF
## 5 P17B_AAACCTGCAGGCGATA-1 P17B TRAV12-1.TRAJ9.TRAC CVVSDNTGGFKTIF
## 7 P17B_AAACCTGCATGAGCGA-1 P17B TRAV12-1.TRAJ9.TRAC CVVSDNTGGFKTIF
## 9 P17B_AAACGGGAGAGCCCAA-1 P17B TRAV20.TRAJ8.TRAC CAVRGEGFQKLVF
## 10 P17B_AAACGGGAGCGTTTAC-1 P17B TRAV12-1.TRAJ9.TRAC CVVSDNTGGFKTIF
## cdr3_nt1
## 1 TGTGGGTGTTCTAACGACTACAAGCTCAGCTTT
## 3 TGTGCTTATAGGAGCGCGCAGGCTGGTGGTACTAGCTATGGAAAGCTGACATTT
## 5 TGTGTGGTCTCCGATAATACTGGAGGCTTCAAAACTATCTTT
## 7 TGTGTGGTCTCCGATAATACTGGAGGCTTCAAAACTATCTTT
## 9 TGTGCTGTGCGAGGAGAAGGCTTTCAGAAACTTGTATTT
## 10 TGTGTGGTCTCCGATAATACTGGAGGCTTCAAAACTATCTTT
## TCR2 cdr3_aa2
## 1 TRBV5-1.None.TRBJ2-7.TRBC2 CASSLTDRTYEQYF
## 3 TRBV10-3.None.TRBJ2-2.TRBC2 CAISEQGKGELFF
## 5 TRBV9.None.TRBJ2-2.TRBC2 CASSVRRERANTGELFF
## 7 TRBV9.None.TRBJ2-2.TRBC2 CASSVRRERANTGELFF
## 9 <NA> <NA>
## 10 TRBV9.None.TRBJ2-2.TRBC2 CASSVRRERANTGELFF
## cdr3_nt2
## 1 TGCGCCAGCAGCTTGACCGACAGGACCTACGAGCAGTACTTC
## 3 TGTGCCATCAGTGAACAGGGGAAAGGGGAGCTGTTTTTT
## 5 TGTGCCAGCAGCGTAAGGAGGGAAAGGGCGAACACCGGGGAGCTGTTTTTT
## 7 TGTGCCAGCAGCGTAAGGAGGGAAAGGGCGAACACCGGGGAGCTGTTTTTT
## 9 <NA>
## 10 TGTGCCAGCAGCGTAAGGAGGGAAAGGGCGAACACCGGGGAGCTGTTTTTT
## CTgene
## 1 TRAV25.TRAJ20.TRAC_TRBV5-1.None.TRBJ2-7.TRBC2
## 3 TRAV38-2/DV8.TRAJ52.TRAC_TRBV10-3.None.TRBJ2-2.TRBC2
## 5 TRAV12-1.TRAJ9.TRAC_TRBV9.None.TRBJ2-2.TRBC2
## 7 TRAV12-1.TRAJ9.TRAC_TRBV9.None.TRBJ2-2.TRBC2
## 9 TRAV20.TRAJ8.TRAC_NA
## 10 TRAV12-1.TRAJ9.TRAC_TRBV9.None.TRBJ2-2.TRBC2
## CTnt
## 1 TGTGGGTGTTCTAACGACTACAAGCTCAGCTTT_TGCGCCAGCAGCTTGACCGACAGGACCTACGAGCAGTACTTC
## 3 TGTGCTTATAGGAGCGCGCAGGCTGGTGGTACTAGCTATGGAAAGCTGACATTT_TGTGCCATCAGTGAACAGGGGAAAGGGGAGCTGTTTTTT
## 5 TGTGTGGTCTCCGATAATACTGGAGGCTTCAAAACTATCTTT_TGTGCCAGCAGCGTAAGGAGGGAAAGGGCGAACACCGGGGAGCTGTTTTTT
## 7 TGTGTGGTCTCCGATAATACTGGAGGCTTCAAAACTATCTTT_TGTGCCAGCAGCGTAAGGAGGGAAAGGGCGAACACCGGGGAGCTGTTTTTT
## 9 TGTGCTGTGCGAGGAGAAGGCTTTCAGAAACTTGTATTT_NA
## 10 TGTGTGGTCTCCGATAATACTGGAGGCTTCAAAACTATCTTT_TGTGCCAGCAGCGTAAGGAGGGAAAGGGCGAACACCGGGGAGCTGTTTTTT
## CTaa
## 1 CGCSNDYKLSF_CASSLTDRTYEQYF
## 3 CAYRSAQAGGTSYGKLTF_CAISEQGKGELFF
## 5 CVVSDNTGGFKTIF_CASSVRRERANTGELFF
## 7 CVVSDNTGGFKTIF_CASSVRRERANTGELFF
## 9 CAVRGEGFQKLVF_NA
## 10 CVVSDNTGGFKTIF_CASSVRRERANTGELFF
## CTstrict
## 1 TRAV25.TRAJ20.TRAC;TGTGGGTGTTCTAACGACTACAAGCTCAGCTTT_TRBV5-1.None.TRBJ2-7.TRBC2;TGCGCCAGCAGCTTGACCGACAGGACCTACGAGCAGTACTTC
## 3 TRAV38-2/DV8.TRAJ52.TRAC;TGTGCTTATAGGAGCGCGCAGGCTGGTGGTACTAGCTATGGAAAGCTGACATTT_TRBV10-3.None.TRBJ2-2.TRBC2;TGTGCCATCAGTGAACAGGGGAAAGGGGAGCTGTTTTTT
## 5 TRAV12-1.TRAJ9.TRAC;TGTGTGGTCTCCGATAATACTGGAGGCTTCAAAACTATCTTT_TRBV9.None.TRBJ2-2.TRBC2;TGTGCCAGCAGCGTAAGGAGGGAAAGGGCGAACACCGGGGAGCTGTTTTTT
## 7 TRAV12-1.TRAJ9.TRAC;TGTGTGGTCTCCGATAATACTGGAGGCTTCAAAACTATCTTT_TRBV9.None.TRBJ2-2.TRBC2;TGTGCCAGCAGCGTAAGGAGGGAAAGGGCGAACACCGGGGAGCTGTTTTTT
## 9 TRAV20.TRAJ8.TRAC;TGTGCTGTGCGAGGAGAAGGCTTTCAGAAACTTGTATTT_NA;NA
## 10 TRAV12-1.TRAJ9.TRAC;TGTGTGGTCTCCGATAATACTGGAGGCTTCAAAACTATCTTT_TRBV9.None.TRBJ2-2.TRBC2;TGTGCCAGCAGCGTAAGGAGGGAAAGGGCGAACACCGGGGAGCTGTTTTTT
## Type
## 1 B
## 3 B
## 5 B
## 7 B
## 9 B
## 10 B3.2 subsetClones: Filter Out Clonal Information
Likewise, we can remove specific list elements after combineTCR() or combineBCR() using the subsetClones() function. In order to subset, we need to identify the column header we would like to use for subsetting (name) and the specific values to include (variables).
Key Parameter(s) for subsetClones()
name: The column header/name in the metadata of input.data to use for subsetting (e.g., “sample”, “Type”).variables: A character vector of the specific values within the chosen name column to retain in the subsetted data.
Below, we isolate just the two sequencing results from “P18L” and “P18B” samples:
subset1 <- subsetClones(combined.TCR,
name = "sample",
variables = c("P18L", "P18B"))
head(subset1[[1]][,1:4])## barcode sample TCR1 cdr3_aa1
## 1 P18B_AAACCTGAGGCTCAGA-1 P18B TRAV26-1.TRAJ37.TRAC CIVRGGSSNTGKLIF
## 3 P18B_AAACCTGCATGACATC-1 P18B TRAV3.TRAJ20.TRAC CAVQRSNDYKLSF
## 5 P18B_AAACCTGGTATGCTTG-1 P18B TRAV26-1.TRAJ53.TRAC CIGSSGGSNYKLTF
## 8 P18B_AAACGGGCAGATGGGT-1 P18B <NA> <NA>
## 9 P18B_AAACGGGTCTTACCGC-1 P18B TRAV20.TRAJ9.TRAC CAVQAKRYTGGFKTIF
## 12 P18B_AAAGATGAGTTACGGG-1 P18B TRAV8-3.TRAJ8.TRAC CAVGGDTGFQKLVFAlternatively, we can also just select the list elements after combineTCR() or combineBCR().
subset2 <- combined.TCR[c(3,4)]
head(subset2[[1]][,1:4])## barcode sample TCR1 cdr3_aa1
## 1 P18B_AAACCTGAGGCTCAGA-1 P18B TRAV26-1.TRAJ37.TRAC CIVRGGSSNTGKLIF
## 3 P18B_AAACCTGCATGACATC-1 P18B TRAV3.TRAJ20.TRAC CAVQRSNDYKLSF
## 5 P18B_AAACCTGGTATGCTTG-1 P18B TRAV26-1.TRAJ53.TRAC CIGSSGGSNYKLTF
## 8 P18B_AAACGGGCAGATGGGT-1 P18B <NA> <NA>
## 9 P18B_AAACGGGTCTTACCGC-1 P18B TRAV20.TRAJ9.TRAC CAVQAKRYTGGFKTIF
## 12 P18B_AAAGATGAGTTACGGG-1 P18B TRAV8-3.TRAJ8.TRAC CAVGGDTGFQKLVF3.3 exportClones: Save Clonal Data
After assigning the clone by barcode, we can export the clonal information using exportClones() to save for later use or to integrate with other bioinformatics pipelines. This function supports various output formats tailored for different analytical needs.
Key Parameter(s) for exportClones()
* format: The desired output format for the clonal data.
* airr: Exports data in an Adaptive Immune Receptor Repertoire (AIRR) Community-compliant format, with each row representing a single receptor chain.
* immunarch: Exports a list containing a data frame and metadata formatted for use with the immunarch package.
* paired: Exports a data frame with paired chain information (amino acid, nucleotide, genes) per barcode. This is the default.
* TCRMatch: Exports a data frame specifically for the TCRMatch algorithm, containing TRB chain amino acid sequence and clonal frequency.
* tcrpheno: Exports a data frame compatible with the tcrpheno pipeline, with TRA and TRB chains in separate columns.
* dowser: Exports an AIRR data frame ready for dowser::formatClones() for B cell lineage analysis. Requires full-length IMGT-gapped alignments retained via retain.sequences. See exportDowser() and the Retaining Full-Length Sequences section.
* write.file: If TRUE (default), saves the output to a CSV file. If FALSE, returns the data frame or list to the R environment.
* dir: The directory where the output file will be saved. Defaults to the current working directory.
* file.name: The name of the CSV file to be saved.
To export the combined clonotypes as a paired data frame and save it to a specified directory:
exportClones(combined,
write.file = TRUE,
dir = "~/Documents/MyExperiment/Sample1/"
file.name = "clones.csv")To return an immunarch-formatted data frame directly to your R environment without saving a file:
immunarch <- exportClones(combined.TCR,
format = "immunarch",
write.file = FALSE)
head(immunarch[[1]][[1]])## Clones Proportion CDR3.nt
## 1 1 0.0003565062 NA;TGCGCCAGCAGTCGGGGACTAGCGGGATACAATGAGCAGTTCTTC
## 2 1 0.0003565062 NA;TGTGCCATCAGCGCGGACCCCCGCTACAATGAGCAGTTCTTC
## 3 1 0.0003565062 NA;TGTGCCAGCAGCTTGAGGGACAGCTATCGGTACTATGGCTACACCTTC
## 4 2 0.0007130125 NA;TGTGCCAGCAGCCGGCAGGGCGCAGATACGCAGTATTTT
## 5 1 0.0003565062 NA;TGTGCCAGCAGTCCCTTTACAGGGTTCTATGGCTACACCTTC
## 6 1 0.0003565062 NA;TGTGCCAGCTCATCCGGGATCAATCAGCCCCAGCATTTT
## CDR3.aa V.name D.name J.name C.name
## 1 NA;CASSRGLAGYNEQFF NA;TRBV10-2 NA;None NA;TRBJ2-1 NA;TRBC2
## 2 NA;CAISADPRYNEQFF NA;TRBV10-3 NA;None NA;TRBJ2-1 NA;TRBC2
## 3 NA;CASSLRDSYRYYGYTF NA;TRBV11-3 NA;None NA;TRBJ1-2 NA;TRBC1
## 4 NA;CASSRQGADTQYF NA;TRBV11-3 NA;None NA;TRBJ2-3 NA;TRBC2
## 5 NA;CASSPFTGFYGYTF NA;TRBV12-4 NA;None NA;TRBJ1-2 NA;TRBC1
## 6 NA;CASSSGINQPQHF NA;TRBV18 NA;None NA;TRBJ1-5 NA;TRBC1
## Barcode
## 1 P17B_AGCGGTCCAAAGGAAG-1
## 2 P17B_GGCTCGAGTCGCGGTT-1
## 3 P17B_CGCGTTTTCGGCTACG-1
## 4 P17B_CACCAGGGTTCCTCCA-1;P17B_TCTATTGCAGGTGCCT-1
## 5 P17B_GCTGGGTGTACGAAAT-1
## 6 P17B_AGGGTGACATTGGTAC-13.4 clonalBin: Bin Clones by Frequency or Proportion
The clonalBin() function adds a clonal grouping variable (cloneSize) to the output of combineTCR(), combineBCR(), or combineExpression(). This function calculates the clonal frequency and proportion, then bins clones into categories based on customizable thresholds. This is useful for categorizing clones prior to downstream analysis or visualization, without needing to attach the data to a single-cell object.
Key Parameter(s) for clonalBin()
clone.call: Defines the clonal sequence grouping (gene,nt,aa, orstrict).chain: The TCR/BCR chain to use (both,TRA,TRB,TRG,TRD,IGH,IGL).group.by: A column header in the metadata to group the analysis by (e.g., “sample”, “Type”).proportion: Whether to use proportion (TRUE) or total frequency (FALSE) for binning.clone.size: The bins for grouping based on proportion or frequency.
3.4.1 Basic Usage with Proportion
By default, clonalBin() uses proportion-based binning with the following default bins: c(Rare = 1e-4, Small = 0.001, Medium = 0.01, Large = 0.1, Hyperexpanded = 1).
combined.TCR <- clonalBin(combined.TCR, clone.call = "strict")
# Check the cloneSize column
head(combined.TCR[[1]][, c("barcode", "CTstrict", "clonalProportion", "clonalFrequency", "cloneSize")])## barcode
## 1 P17B_AGCGGTCCAAAGGAAG-1
## 2 P17B_GGCTCGAGTCGCGGTT-1
## 3 P17B_CGCGTTTTCGGCTACG-1
## 4 P17B_CACCAGGGTTCCTCCA-1
## 5 P17B_TCTATTGCAGGTGCCT-1
## 6 P17B_GCTGGGTGTACGAAAT-1
## CTstrict
## 1 NA;NA_TRBV10-2.None.TRBJ2-1.TRBC2;TGCGCCAGCAGTCGGGGACTAGCGGGATACAATGAGCAGTTCTTC
## 2 NA;NA_TRBV10-3.None.TRBJ2-1.TRBC2;TGTGCCATCAGCGCGGACCCCCGCTACAATGAGCAGTTCTTC
## 3 NA;NA_TRBV11-3.None.TRBJ1-2.TRBC1;TGTGCCAGCAGCTTGAGGGACAGCTATCGGTACTATGGCTACACCTTC
## 4 NA;NA_TRBV11-3.None.TRBJ2-3.TRBC2;TGTGCCAGCAGCCGGCAGGGCGCAGATACGCAGTATTTT
## 5 NA;NA_TRBV11-3.None.TRBJ2-3.TRBC2;TGTGCCAGCAGCCGGCAGGGCGCAGATACGCAGTATTTT
## 6 NA;NA_TRBV12-4.None.TRBJ1-2.TRBC1;TGTGCCAGCAGTCCCTTTACAGGGTTCTATGGCTACACCTTC
## clonalProportion clonalFrequency cloneSize
## 1 0.0003565062 1 Small (1e-04 < X <= 0.001)
## 2 0.0003565062 1 Small (1e-04 < X <= 0.001)
## 3 0.0003565062 1 Small (1e-04 < X <= 0.001)
## 4 0.0007130125 2 Small (1e-04 < X <= 0.001)
## 5 0.0007130125 2 Small (1e-04 < X <= 0.001)
## 6 0.0003565062 1 Small (1e-04 < X <= 0.001)# View the distribution of clone sizes
table(combined.TCR[[1]]$cloneSize)##
## Hyperexpanded (0.1 < X <= 1) Large (0.01 < X <= 0.1)
## 952 784
## Medium (0.001 < X <= 0.01) Small (1e-04 < X <= 0.001)
## 340 729
## Rare (0 < X <= 1e-04) None ( < X <= 0)
## 0 03.4.2 Using Frequency-Based Binning
When using frequency-based binning (proportion = FALSE), the clone.size values must be integers. If the maximum frequency exceeds the upper bin limit, it will be automatically adjusted.
combined.TCR.freq <- combineTCR(contig_list,
samples = c("P17B", "P17L", "P18B", "P18L",
"P19B","P19L", "P20B", "P20L"))
combined.TCR.freq <- clonalBin(combined.TCR.freq,
clone.call = "strict",
proportion = FALSE,
clone.size = c(Rare = 1, Small = 5, Medium = 20,
Large = 100, Hyperexpanded = 500))
# Check the frequency-based binning
head(combined.TCR.freq[[1]][, c("barcode", "clonalFrequency", "cloneSize")])## barcode clonalFrequency cloneSize
## 1 P17B_AGCGGTCCAAAGGAAG-1 1 Rare (0 < X <= 1)
## 2 P17B_GGCTCGAGTCGCGGTT-1 1 Rare (0 < X <= 1)
## 3 P17B_CGCGTTTTCGGCTACG-1 1 Rare (0 < X <= 1)
## 4 P17B_CACCAGGGTTCCTCCA-1 2 Small (1 < X <= 5)
## 5 P17B_TCTATTGCAGGTGCCT-1 2 Small (1 < X <= 5)
## 6 P17B_GCTGGGTGTACGAAAT-1 1 Rare (0 < X <= 1)3.4.3 Grouping by a Variable
You can calculate clonal frequency and proportion across samples grouped by a variable (e.g., “Type” or “sample”). This is particularly useful when you want consistent binning across experimental conditions.
combined.TCR.grouped <- combineTCR(contig_list,
samples = c("P17B", "P17L", "P18B", "P18L",
"P19B","P19L", "P20B", "P20L"))
combined.TCR.grouped <- addVariable(combined.TCR.grouped,
variable.name = "Type",
variables = rep(c("B", "L"), 4))
combined.TCR.grouped <- clonalBin(combined.TCR.grouped,
clone.call = "strict",
group.by = "Type")
# Check grouped binning
head(combined.TCR.grouped[[1]][, c("barcode", "Type", "clonalProportion", "cloneSize")])## barcode Type clonalProportion cloneSize
## 1 P20B_TACTCGCAGAGGTTGC-1 B 4.983554e-05 Rare (0 < X <= 1e-04)
## 2 P18L_CGTGTCTGTGTATGGG-1 L 1.404692e-04 Small (1e-04 < X <= 0.001)
## 3 P18L_AGTGTCAGTAGCGCAA-1 L 1.404692e-04 Small (1e-04 < X <= 0.001)
## 4 P18L_TTTACTGCAAGTTCTG-1 L 1.404692e-04 Small (1e-04 < X <= 0.001)
## 5 P19L_GTGTGCGCAAAGGAAG-1 L 1.404692e-04 Small (1e-04 < X <= 0.001)
## 6 P20L_CCCTCCTAGTAGCGGT-1 L 1.404692e-04 Small (1e-04 < X <= 0.001)3.5 annotateInvariant
The annotateInvariant() function enables the identification of mucosal-associated invariant T (MAIT) cells and invariant natural killer T (iNKT) cells in single-cell sequencing datasets. These specialized T-cell subsets are defined by their characteristic TCR usage, making them distinguishable within single-cell immune profiling data. The function extracts TCR chain information from the provided single-cell dataset and evaluates it against known invariant TCR criteria for either MAIT or iNKT cells. Each cell is assigned a score indicating the presence (1) or absence (0) of the specified invariant T-cell population.
Key Parameter(s) for annotateInavriant()
type: Character string specifying the type of invariant T cell to annotate (MAIToriNKT).species: Character string specifying the species (mouseorhuman).
combined <- annotateInvariant(combined,
type = "MAIT",
species = "human")
combined <- annotateInvariant(combined,
type = "iNKT",
species = "human")4 Basic Clonal Visualizations
4.1 cloneCall: How to call clones.
gene- use the VDJC genes comprising the TCR/Ig
nt- use the nucleotide sequence of the CDR3 region
aa- use the amino acid sequence of the CDR3 region
strict- use the VDJC genes comprising the TCR + the nucleotide sequence of the CDR3 region. This is the proper definition of clonotype. ForcombineBCR()strict refers to the edit distance clusters + Vgene of the Ig.
It is important to note that the clonotype is called using essentially the combination of genes or nt/aa CDR3 sequences for both loci. As of this implementation of scRepertoire, clonotype calling is not incorporating small variations within the CDR3 sequences. As such the gene approach will be the most sensitive, while the use of nt or aa is moderately so, and the most specific for clonotypes being strict. Additionally, the clonotype call is trying to incorporate both loci, i.e., both TCRA and TCRB chains and if a single cell barcode has multiple sequences identified (i.e., 2 TCRA chains expressed in one cell). Using the 10x approach, there is a subset of barcodes that only return one of the immune receptor chains. The unreturned chain is assigned an NA value.
4.2 clonalQuant: Quantifying Unique Clones
The clonalQuant() function is used to explore the clones by returning the total or relative numbers of unique clones.
Key Parameter(s) for clonalQuant()
scale: IfTRUE, converts the output to the relative percentage of unique clones scaled by the total repertoire size; ifFALSE(default), reports the total number of unique clones.
To visualize the relative percent of unique clones across all chains ("both") using the strict clone definition:
clonalQuant(combined.TCR,
cloneCall="strict",
chain = "both",
scale = TRUE)
Another option is to define the visualization by data classes using the group.by parameter. Here, we’ll use the "Type" variable, which was previously added to the combined.TCR list.
clonalQuant(combined.TCR,
cloneCall="gene",
group.by = "Type",
scale = TRUE)
4.3 clonalAbundance: Distribution of Clones by Size
clonalAbundance() allows for the examination of the relative distribution of clones by abundance. It produces a line graph showing the total number of clones at specific frequencies within a sample or group.
Key Parameter(s) for clonalAbundance()
scale: IfTRUE, converts the graphs into density plots to show relative distributions; ifFALSE(default), displays raw counts.
To visualize the raw clonal abundance using the gene clone definition:
clonalAbundance(combined.TCR,
cloneCall = "gene",
scale = FALSE)
clonalAbundance() output can also be converted into a density plot, which may allow for better comparisons between different repertoire sizes, by setting scale = TRUE.
clonalAbundance(combined.TCR,
cloneCall = "gene",
scale = TRUE)
4.4 clonalLength: Distribution of Sequence Lengths
clonalLength() allows you to look at the length distribution of the CDR3 sequences. Importantly, unlike the other basic visualizations, the cloneCall can only be nt (nucleotide) or aa (amino acid). Due to the method of calling clones as outlined previously (e.g., using NA for unreturned chain sequences or multiple chains within a single barcode), the length distribution may reveal a multimodal curve.
To visualize the amino acid length distribution for both chains (“both”):
clonalLength(combined.TCR,
cloneCall="aa",
chain = "both") 
To visualize the amino acid length distribution for the TRA chain, scaled as a density plot:
clonalLength(combined.TCR,
cloneCall="aa",
chain = "TRA",
scale = TRUE) 
4.5 clonalCompare: Clonal Dynamics Between Categorical Variables
clonalCompare() allows you to look at clones between samples and changes in dynamics. It is useful for tracking how the proportions of top clones change between conditions.
Key Parameters for clonalCompare()
samples: A character vector to isolate specific samples by their list element name.clones: A character vector of specific clonal sequences to visualize. If used,top.cloneswill be ignored. *top.clones: The top n number of clones to graph, calculated based on the frequency within individual samples.highlight.clones: A character vector of specific clonal sequences to color; all other clones will be greyed out.relabel.clones: IfTRUE, simplifies the legend by labeling isolated clones numerically (e.g., “Clone: 1”).graph: The type of plot to generate;alluvial(default) orarea.proportion: IfTRUE(default), the y-axis represents proportional abundance; ifFALSE, it represents raw clone counts.
To compare the top 10 clones between samples “P17B” and “P17L” using amino acid sequences as an alluvial plot:
clonalCompare(combined.TCR,
top.clones = 10,
samples = c("P17B", "P17L"),
cloneCall="aa",
graph = "alluvial")
We can also choose to highlight specific clones, such as in the case of “CVVSDNTGGFKTIF_CASSVRRERANTGELFF” and “NA_CASSVRRERANTGELFF” using the highlight.clones parameter. In addition, we can simplify the plot to label the clones as “Clone: 1”, “Clone: 2”, etc., by setting relabel.clones = TRUE.
clonalCompare(combined.TCR,
top.clones = 10,
highlight.clones = c("CVVSDNTGGFKTIF_CASSVRRERANTGELFF", "NA_CASSVRRERANTGELFF"),
relabel.clones = TRUE,
samples = c("P17B", "P17L"),
cloneCall="aa",
graph = "alluvial")
Alternatively, if we only want to show specific clones, we can use the clones parameter.
clonalCompare(combined.TCR,
clones = c("CVVSDNTGGFKTIF_CASSVRRERANTGELFF", "NA_CASSVRRERANTGELFF"),
relabel.clones = TRUE,
samples = c("P17B", "P17L"),
cloneCall="aa",
graph = "alluvial")
4.6 clonalScatter: Scatterplot of Two Variables
clonalScatter() organizes two repertoires, quantifies their relative clone sizes, and produces a scatter plot comparing the two samples. Clones are categorized by counts into singlets or expanded, either exclusively present or shared between the selected samples.
Key Parameter(s) for clonalScatter()
x.axis,y.axis: Names of the list elements or meta data variable to place on the x-axis and y-axis.dot.size: Specifies how dot size is determined;total(default) displays the total number of clones between the x- and y-axis, or a specific list element name for size calculation.graph: The type of graph to display;proportionfor the relative proportion of clones (default) orcountfor the total count of clones by sample.
To compare samples “P18B” and “P18L” based on gene clone calls, with dot size representing the total number of clones, and plotting clone proportions:
clonalScatter(combined.TCR,
cloneCall ="gene",
x.axis = "P18B",
y.axis = "P18L",
dot.size = "total",
graph = "proportion")
5 Visualizing Clonal Dynamics
5.1 clonalHomeostasis: Examining Clonal Space
By examining the clonal space, we effectively look at the relative space occupied by clones at specific proportions. Another way to think about this would be to consider the total immune receptor sequencing run as a measuring cup. In this cup, we will fill liquids of different viscosity—or different numbers of clonal proportions. Clonal space homeostasis asks what percentage of the cup is filled by clones in distinct proportions (or liquids of different viscosity, to extend the analogy). The proportional cut points are set under the cloneSize variable in the function and can be adjusted.
Default cloneSize Bins
- Rare: 0 to 0.0001
- Small: 0.0001 to 0.001
- Medium: 0.001 to 0.01
- Large: 0.01 to 0.1
- Hyperexpanded: 0.1 to 1
To visualize clonal homeostasis using gene clone calls with default cloneSize bins:
clonalHomeostasis(combined.TCR,
cloneCall = "gene")
We can reassign the proportional cut points for cloneSize, which can drastically alter the visualization and analysis
clonalHomeostasis(combined.TCR,
cloneCall = "gene",
cloneSize = c(Rare = 0.001, Small = 0.01, Medium = 0.1,
Large = 0.3, Hyperexpanded = 1))
In addition, we can use the group.by parameter to look at the relative proportion of clones between groups, such as by tissue type. First, ensure the “Type” variable is added to combined.TCR:
combined.TCR <- addVariable(combined.TCR,
variable.name = "Type",
variables = rep(c("B", "L"), 4))Now, visualize clonal homeostasis grouped by “Type”:
clonalHomeostasis(combined.TCR,
group.by = "Type",
cloneCall = "gene")
clonalHomeostasis() provides an assessment of how different “sizes” of clones (based on their proportional abundance) contribute to the overall repertoire. This visualization helps to identify shifts in repertoire structure, such as expansion of large clones in response to infection or contraction in chronic conditions, offering insights into immune system activity and health.
5.2 clonalProportion: Examining Space Occupied by Ranks of Clones
Like clonal space homeostasis, clonalProportion() also categorizes clones into separate bins. The key difference is that instead of looking at the relative proportion of the clone to the total, clonalProportion() ranks the clones by their total count or frequency of occurrence and then places them into predefined bins.
The clonalSplit parameter represents the ranking of clonotypes. For example, 1:10 refers to the top 10 clonotypes in each sample. The default bins are set under the clonalSplit variable and can be adjusted.
Default clonalSplit Bins
- 10 (top 1-10 clones)
- 100: (top 11-100 clones)
- 1000: (top 101-1000 clones)
- 10000: (top 1001-10000 clones)
- 30000: (top 10001-30000 clones)
- 100000: (top 30001-100000 clones)
To visualize the clonal proportion using gene clone calls with default clonalSplit bins:
clonalProportion(combined.TCR,
cloneCall = "gene") 
To visualize clonal proportion using nt (nucleotide) clone calls with custom clonalSplit bins:
clonalProportion(combined.TCR,
cloneCall = "nt",
clonalSplit = c(1, 5, 10, 100, 1000, 10000)) 
clonalProportion() complements clonalHomeostasis() by providing a perspective on how the “richest” (most abundant) clones contribute to the overall repertoire space. By segmenting clones into rank-based bins, it helps identify whether a few highly expanded clones or a larger number of moderately expanded clones dominate the immune response, offering distinct insights into repertoire structure and dynamics.
6 Summarizing Repertoires
6.1 percentGeneUsage
Gene quantification and visualization has been redesigned to offer a more robust and translatable function under percentGeneUsage(). We have maintained the functionality of the previous functions for gene-level quantification under the aliases vizGenes(), percentGenes(), and percentVJ().
6.1.1 vizGenes: Flexible Gene Usage Visualization
The vizGenes() function offers a highly adaptable approach to visualizing the relative usage of TCR or BCR genes. It acts as a versatile alias for percentGeneUsage(), allowing for comparisons across different chains, scaling of values, and selection between bar charts and heatmaps.
x.axis: Specifies the gene segment to display along the x-axis (e.g., “TRBV”, “TRBD”, “IGKJ”).y.axis- Another gene segment (e.g., “TRAV”, “TRBJ”) for paired gene analysis. When x.axis and y.axis are both gene segments, vizGenes() internally calls percentGeneUsage() with genes = c(x.axis, y.axis), resulting in a heatmap.
- A categorical variable (e.g., “sample”, “orig.ident”) to visualize gene usage across different groups. When y.axis is a categorical variable, vizGenes() maps it to the group.by parameter of percentGeneUsage(), creating facets for each category.
plot: Determines the visualization type:- “barplot”: Ideal for visualizing the distribution of a single gene segment.
- “heatmap”: Suitable for single gene usage (with a group.by or y.axis categorical variable) or for paired gene analysis.
summary.fun: (Inherited frompercentGeneUsage) Defines the statistic to display: “percent” (default), “proportion”, or “count”. This implicitly handles the scaling of values.
vizGenes(combined.TCR,
x.axis = "TRBV",
y.axis = NULL, # No specific y-axis variable, will group all samples
plot = "barplot",
summary.fun = "proportion") 
This plot shows the proportion of each TRBV gene segment observed across the entire combined.TCR dataset. Since y.axis is NULL, samples are grouped by the list element of the combined.TCR data.
vizGenes() is particularly useful for examining gene pairings. Let’s look at the differences in TRBV and TRBJ usage between peripheral blood and lung samples from your dataset. We’ll subset combined.TCR for this analysis.
# Peripheral Blood Samples
vizGenes(combined.TCR[c("P17B", "P18B", "P19B", "P20B")],
x.axis = "TRBV",
y.axis = "TRBJ",
plot = "heatmap",
summary.fun = "percent") # Display percentages
# Lung Samples
vizGenes(combined.TCR[c("P17L", "P18L", "P19L", "P20L")],
x.axis = "TRBV",
y.axis = "TRBJ",
plot = "heatmap",
summary.fun = "percent") # Display percentages
In these examples, by providing both x.axis and y.axis as gene segments (“TRBV” and “TRBJ”), vizGenes() automatically performs a paired gene analysis, generating a heatmap where the intensity reflects the percentage of each V-J pairing.
Beyond V-J pairings within a single chain, vizGenes() can also visualize gene usage across different chains. For instance, to examine TRBV and TRAV pairings for patient P17’s samples:
vizGenes(combined.TCR[c("P17B", "P17L")],
x.axis = "TRBV",
y.axis = "TRAV",
plot = "heatmap",
summary.fun = "count") 
6.1.2 percentGenes: Quantifying Single Gene Usage
The percentGenes() function is a specialized alias for percentGeneUsage() designed to quantify the usage of a single V, D, or J gene locus for a specified immune receptor chain. By default, it returns a heatmap visualization.
Key Parameters for percentGenes():
chain: Specifies the immune receptor chain (e.g., “TRB”, “TRA”, “IGH”, “IGL”).gene: Indicates the gene locus to quantify: “Vgene”, “Dgene”, or “Jgene”.group.by,order.by,summary.fun,exportTable,palette: These parameters are directly passed topercentGeneUsage().
To quantify and visualize the percentage of TRBV gene usage across your samples:
percentGenes(combined.TCR,
chain = "TRB",
gene = "Vgene",
summary.fun = "percent")
This generates a heatmap showing the percentage of each TRBV gene segment within each sample, allowing for easy visual comparison of gene usage profiles across your samples.
The raw data returned by percentGenes() (when exportTable = TRUE) can be a powerful input for further downstream analysis, such as dimensionality reduction techniques like Principal Component Analysis. This allows you to summarize the complex gene usage patterns and identify samples with similar or distinct repertoires.
df.genes <- percentGenes(combined.TCR,
chain = "TRB",
gene = "Vgene",
exportTable = TRUE,
summary.fun = "proportion")
# Performing PCA on the gene usage matrix
pc <- prcomp(t(df.genes) )
# Getting data frame to plot from
df_plot <- as.data.frame(cbind(pc$x[,1:2], colnames(df.genes)))
colnames(df_plot) <- c("PC1", "PC2", "Sample")
df_plot$PC1 <- as.numeric(df_plot$PC1)
df_plot$PC2 <- as.numeric(df_plot$PC2)
ggplot(df_plot, aes(x = PC1, y = PC2)) +
geom_point(aes(fill = Sample), shape = 21, size = 5) +
guides(fill=guide_legend(title="Samples")) +
scale_fill_manual(values = hcl.colors(nrow(df_plot), "inferno")) +
theme_classic() +
labs(title = "PCA of TRBV Gene Usage")
This PCA plot visually clusters samples based on their TRBV gene usage profiles, helping to identify underlying patterns or relationships between samples.
6.1.3 percentVJ: Quantifying V-J Gene Pairings
The percentVJ() function is another specialized alias for percentGeneUsage(), specifically designed to quantify the proportion or percentage of V and J gene segments paired together within individual repertoires. It always produces a heatmap visualization.
Key Parameters for percenVJ():
chain: Specifies the immune receptor chain (e.g., “TRB”, “TRA”, “IGH”, “IGL”). This dictates which V and J gene segments are analyzed (e.g., TRBV and TRBJ forchain= “TRB”)
percentVJ(combined.TCR[1:2],
chain = "TRB",
summary.fun = "percent")
This heatmap displays the percentage of each TRBV-TRBJ pairing for each sample, providing a detailed view of the V-J recombination landscape.
Similar to percentGenes(), the quantitative output of percentVJ() can be used for dimensionality reduction to summarize V-J pairing patterns across samples.
df.vj <- percentVJ(combined.TCR,
chain = "TRB",
exportTable = TRUE,
summary.fun = "proportion") # Export proportions for PCA
# Performing PCA on the V-J pairing matrix
pc.vj <- prcomp(t(df.vj))
# Getting data frame to plot from
df_plot_vj <- as.data.frame(cbind(pc.vj$x[,1:2], colnames(df.vj)))
colnames(df_plot_vj) <- c("PC1", "PC2", "Sample")
df_plot_vj$PC1 <- as.numeric(df_plot_vj$PC1)
df_plot_vj$PC2 <- as.numeric(df_plot_vj$PC2)
# Plotting the PCA results
ggplot(df_plot_vj, aes(x = PC1, y = PC2)) +
geom_point(aes(fill = Sample), shape = 21, size = 5) +
guides(fill=guide_legend(title="Samples")) +
scale_fill_manual(values = hcl.colors(nrow(df_plot_vj), "inferno")) +
theme_classic() +
labs(title = "PCA of TRBV-TRBJ Gene Pairings")
6.2 percentAA: Amino Acid Composition by Residue
Quantify the proportion of amino acids along the CDR3 sequence with percentAA(). By default, the function will pad the sequences with NAs up to the maximum of aa.length. Sequences longer than aa.length will be removed before visualization (default aa.length = 20).
Key Parameter(s) for percentAA()
aa.length: The maximum length of the CDR3 amino acid sequence to consider.
To visualize the relative percentage of amino acids at each position of the CDR3 sequences for the TRB chain, up to a length of 20 amino acids:
percentAA(combined.TCR,
chain = "TRB",
aa.length = 20)
This plot displays the relative proportion of each amino acid at different positions across the CDR3 sequence for the specified chain. It provides insights into the amino acid composition and variability at each residue, which can be indicative of functional constraints or selection pressures on the CDR3
6.3 positionalEntropy: Entropy across CDR3 Sequences
We can also quantify the level of entropy/diversity across amino acid residues along the CDR3 sequence. positionalEntropy() combines the quantification by residue of percentAA() with diversity calculations. Positions without variance will have a value reported as 0 for the purposes of comparison.
Key Parameter(s) for positionalEntropy()
methodshannon- Shannon Indexinv.simpson- Inverse Simpson Indexgini.simpson- Gini-Simpson Indexnorm.entropy- Normalized Entropypielou- Pielou’s Evennesshill1,hill2,hill3- Hill Numbers
To visualize the normalized entropy across amino acid residues for the TRB chain, up to a length of 20 amino acids:
positionalEntropy(combined.TCR,
chain = "TRB",
aa.length = 20)
The plot generated by positionalEntropy() illustrates the diversity or entropy at each amino acid position within the CDR3 sequence. Higher entropy values indicate greater variability in amino acid usage at that position, suggesting less selective pressure or more promiscuous binding, while lower values suggest conserved positions critical for structural integrity or antigen recognition.
6.4 positionalProperty: Amino Acid Properties across CDR3 Sequence
Like positionalEntropy(), we can also examine a series of amino acid properties along the CDR3 sequences using positionalProperty(). Important differences for positionalProperty() is dropping NA values as they would void the mean calculation and displaying a ribbon with the 95% confidence interval surrounding the mean value for the selected properties.
Key Parameter(s) for positionalEntropy()
method
To examine the Atchley factors of amino acids across the CDR3 sequence for the first two samples:
positionalProperty(combined.TCR[c(1,2)],
chain = "TRB",
aa.length = 20,
method = "atchleyFactors") +
scale_color_manual(values = hcl.colors(5, "inferno")[c(2,4)])
This function provides a detailed view of how physicochemical properties of amino acids change along the CDR3 sequence. By visualizing the mean property value and its confidence interval, one can identify positions with distinct characteristics or significant variations between groups, offering insights into structural and functional aspects of the CDR3.
6.5 percentKmer: Motif Quantification
Another quantification of the composition of the CDR3 sequence is to define motifs by sliding across the amino acid or nucleotide sequences at set intervals, resulting in substrings or kmers.
Key Parameter(s) for percentKmer():
cloneCall: Defines the clonal sequence grouping; only acceptsaa(amino acids) ornt(nucleotides)motif.length: The length of the kmer to analyze.top.motifs: Displays the n most variable motifs as a function of median absolute deviation.
To visualize the percentage of the top 25 most variable 3-mer amino acid motifs for the TRB chain:
percentKmer(combined.TCR,
cloneCall = "aa",
chain = "TRB",
motif.length = 3,
top.motifs = 25)
To perform the same analysis but for nucleotide motifs:
percentKmer(combined.TCR,
cloneCall = "nt",
chain = "TRB",
motif.length = 3,
top.motifs = 25)