Immune repertoires are comprised of a diverse collection of BCRs and TCRs, which enable molecular recognition to a vast number of pathogen and disease antigens. Immune repertoire diversity is initially generated as a result of lymphocyte V(D)J recombination and, in the case of B cells, can undergo further sequence diversification in the form of somatic hypermutation. Targeted deep sequencing of BCRs and TCRs from bulk populations of lymphocytes has paved the way to quantify the diversity, distribution and evolution of immune repertoires [1] [2] [3] [4] . However, one major challenge in immune repertoire sequencing is acquiring information on correct receptor chain pairing [variable light (VL) and variable heavy (VH) for BCR and variable alpha (Va) and variable beta (Vb) for TCR], which greatly complicates identification of clonal groups and antigen-specificity 5, 6 . Until only recently it was not possible to directly integrate a lymphocyte's phenotypic gene expression information (i.e., activation, exhaustion, antibody secretion) with its immune receptor sequence.
Recent developments in microfluidic and scSeq technologies have now made it possible to obtain information on immune repertoires or transcriptional profiles at high-throughput [7] [8] [9] . Several of these methods have been tailored specifically for lymphocytes, thus making it possible to perform parallel sequencing of immune repertoires and whole transcriptomes 10, 11 . Furthermore, commercially available instruments and protocols [10X Genomics Chromium and VDJ and gene expression (GEX) libraries] are further accelerating progress in this field. This simultaneous sequencing of immune repertoires and transcriptomes produces single-cell datasets with features ranging from quantitative gene profiles, cellular phenotypes, transcriptional clustering, clonal diversity and expansion, germline gene usage, somatic hypermutation, among many others. These high-dimensional datasets can be mined to discover novel insight on lymphocyte immunobiology, function, and specificity. For example, one recent study leveraged scSeq to discover the distinct transcriptional profiles and specificity of B cells following influenza vaccination 8 . In another study, clonal expansion and activation signatures of tumor-infiltrating T cells were profiled 12 . Other studies have leveraged this technology to answer fundamental questions across a variety of areas in immunology such as bacterial infection responses 13 , tumor-immune microenvironment 14 , clonal expansion in Alzheimer's disease 9 , and B cell differentiation 15 .
While multiple bioinformatic tools exist to facilitate rapid analysis of gene expression from scSeq [16] [17] [18] [19] , they do not allow the incorporation of immune repertoire information. Analogously, existing software packages to analyze immune repertoires do not allow the user to supply accompanying gene expression and transcriptome data [20] [21] [22] . Taken together, these considerations complicate the analysis of BCR and TCR repertoires for those with little bioinformatics experience and who are unfamiliar with the output data from scSeq experiments. To address the lack of software specifically tailored to single-cell lymphocyte sequencing data, we developed Platypus, an open-source R package that contains an automated pipeline to analyze and integrate single-cell immune repertoires with transcriptome data. With only a few lines of code, Platypus allows users to easily extract immune repertoire features such as clonal expansion, somatic hypermutation, isotype switching and integrate it with transcriptome features such as differential gene expression. We subsequently demonstrate the value of this package using scSeq data from convalescent coronavirus infectious disease 2019 (COVID-19) patients. Our analysis revealed clonal expansion in B and T cells, and within these we could identify distinct patterns of somatic hypermutation, amino acid usage, clonal convergence and transcriptional heterogeneity. Taken together, Platypus helps facilitate the analysis of single cell immune repertoires and transcriptomes and reveal novel insights such as the transcriptional profile of clonal expanded and potentially pathogen-reactive lymphocytes.
Platypus allows the user to integrate single-cell immune repertoire and transcriptome sequencing data, which includes automation of pre-processing, filtering, and data visualization ( Figure 1 ). While Platypus is optimized for data generated by the 10X Genomics System, it is also adaptable to other cell-barcode based scSeq data (e.g. RAGE-seq, Split-Seq 10,23 ). Users can simply supply the path to the output directory from the 10X cellranger analysis as input to Platypus, which then extracts and annotates key immune repertoire metrics such as clonal diversity, clonal expansion, somatic hypermutation, reference germline gene usage, and sequence motifs ( Figure 1 ). Platypus can perform additional clonotyping, either increasing or relaxing the pre-determined stringency of upstream alignment tools by incorporating information regarding germline gene usage or sequence homology thresholds. In addition to clonal sequence information [based on complementarity determining region 3 (CDR3)], it also extracts full-length sequences of both immune receptor variable chains (VH and VL for BCRs and Va and Vb for TCRs). Furthermore, Platypus enables the quantification, comparison, and visualization of more advanced repertoire features such as sequence similarity networks 1 , phylogenetic tree construction 24, 25 , isotype quantification, and diversity metrics 26 .
To demonstrate the potential of Platypus, we performed single-cell immune repertoire and transcriptome sequencing on B and T cells isolated from peripheral blood mononuclear cells (PBMCs) of two convalescent COVID-19 patients (Figure 2A ). 10X Genomics' basic alignment tool is cellranger, which has a default clonotyping strategy that groups identical CDRH3 + CDRL3 nucleotide sequences into clonal families; this approach would be too restrictive in identifying clonotypes of B cells that have undergone somatic hypermutation in the CDR3s. We thereby demonstrated the ability and impact of changing the clonotyping strategy to include germline genes, CDR3 length restrictions, and sequence homology requirements for the B cell repertoires of the two COVID-19 patients, which resulted in a decrease in the number of unique clones when additional repertoire features were included in the clonotyping definition ( Figure 2B ). Next, using Platypus, we were able to detect and visualize clonal expansion for both B and T cells ( Figures 2C, 2D ). Furthermore, we were able to relate isotype information with clonal expansion at single-cell resolution, thereby observing that the majority of most expanded B cell clones were of the IgA isotype ( Figure 2D ). We next used built-in functions of Platypus to extract other common immune repertoire statistics and features, such as CDR3 length distribution and common sequence space motifs (sequence logo plots). This revealed tremendous diversity in the B cell response at the most common paired CDRH3 + CDRL3 amino acid sequence length in both COVID-19 patients ( Figures 2E, 2F ). While we focused on the most frequent CDRH3 + CDRL3 sequence length, such an analysis using Platypus could theoretically be applied to other single-cell subsets, such as those cells with known antigen-specificity or those B cells that underwent extensive somatic hypermutation.
A critical feature of Platypus is that it can seamlessly integrate single-cell immune repertoire data with transcriptome sequencing data. It allows users to directly interact with the commonly used scSeq transcriptome analysis program Seurat 16 , while tuning parameters specifically relevant for immune repertoires. Therefore, we next investigated additional repertoire and transcriptional data of highly expanded B and T cell clonal groups, which allows us to relate repertoire information (e.g. expansion, CDR3 sequence, isotype) to phenotypic cellular behavior (e.g. whether a cell is proliferating, differentiated, activated, exhausted, etc.). We first integrated the transcriptome sequencing data from both COVID-19 patients by normalizing and scaling the data using default parameters in Seurat (although Platypus supports other normalization methods such as SCTransform and Harmony) 16, 17 . We could thereby compute clusters based on gene expression and subsequently visualize the cells from each patient on the same two-dimensional uniform manifold approximation projection (UMAP) plot ( Figure 3A ). Overlaying the gene expression of CD4, CD8, CD3, and CD19 revealed a separation between B and T cells, in addition to demonstrating similarly distributed lymphocytes populations arising from both patients ( Figure S1A ). Of note is that we detected CD3 expression in the B cell clusters on the UMAP, in addition to minor CD19 expression in the T cell clusters ( Figures 3A, S1A) , together suggesting the possibility off doublets in which B and T cells were present in the same microfluidic droplet. We next used Platypus to investigate which transcriptional-based cluster contained the most expanded T cell clonal families. We could demonstrate in one of the COVID-19 patients that some of the most expanded T cell clones were located across multiple transcriptional clusters, with 6 of the most expanded clones located in the CD8+ T cell clusters and the remaining 4 were located in the CD4+ T cell clusters ( Figures 3B, S1A) . Furthermore, we observed that the majority of expanded CD8 T cells were located in cluster 6, which corresponded to high expression of granzymes and perforin ( Figure S1B ). Another option to visualize transcriptional heterogeneity using Platypus is to overlay clonal information directly onto the UMAP plot. To investigate the transcriptional heterogeneity of clonally expanded T cells, we supplied the clonal index of the top 10 most expanded clones from a single patient and visualized where these cells lie within the UMAP plot, again revealing that many of the expanded clones were located in cluster 6 ( Figure 3C ). This approach can be leveraged to profile transcriptome signatures of well-defined clones, for example ones in which antigen-specificity is known.
Next we determined whether we could use Platypus to identify potentially virus-reactive B and T cell clones by integrating repertoire metrics such as somatic hypermutation and clonal similarity with phenotypic markers of activation and differentiation. We first noticed that 30 most expanded B cell clones had undergone somatic hypermutation in their VH segments, and since these patients were convalescent but still actively infected, highly mutated antibodies represent potential specific clones to SARS-CoV-2 ( Figure   4A ). Next, we used Platypus to infer the phylogenetic tree for the most expanded B cell clone while also annotating information about clonal expansion (based on identified number of cell barcodes) ( Figure 4B) .
Surprisingly, we uncovered that 59 cells produced the exact same, full-length nucleotide antibody sequence and was actually the least mutated relative to the unmutated germline ancestral sequence ( Figure 4B ). The potential specificity of an antibody that has few somatic hypermutation is compatible with the recent discovery that B cells from COVID-19 patients that have germline-like antibodies with specificity to SARS-CoV-2 22 .
We next questioned whether we could detect similar CDR3 sequences shared between the two COVID-19 patient repertoires by inferring sequence similarity networks for the top 60 most expanded B and T cell clones ( Figures 4C, 4D) 1 . This involved first calculating the edit distance between the CDR3s of clones from both patients and then drawing edges between those clones that were separated by less than 7 amino acid mutations in either the CDRH3 or the CDRL3 (or CDRb3 + CDRa3 in T cells). This revealed that within B cell repertoires, but not T cell repertoires, there were B cell clones that were connected across patients, thereby suggesting convergence of antibody sequences which may be a result of specificity to SARS-CoV-2. Highly similar antibody sequences from the memory B cells of different COVID-19 patients have recently been shown to occur 27 . The lack of connected T cell clones may be expected since it is likely these patients were not HLA-matched, in addition to the small portion of the blood repertoire sampled when performing single-cell sequencing.
We lastly investigated whether signatures of activation or differentiation would reveal potential B and T cell clones that may have recently interacted with viral antigen. Using B and T cell-specific gene sets contained internally in Platypus, we could investigate phenotypic and differentiation (ontogeny) markers for expanded clones ( Figure 5 ). In the case of B cells, this included the following markers: B220 (PTPRC), MS4A1, CD27, CD38, and CD138 (SDC1), which distinguish between naïve, memory, and plasma cells subsets ( Figure   5A ). This analysis further supported a transcriptional heterogeneity within the most expanded clones, as Figure 5B ). Together, these genes allowed us to distinguish between naïve (CD44-, TCF1+, SELL+), memory (IL7R+), effector (KLRG1+), and exhausted (PDCD1+, LAG3+) subtypes that coexist within clonotypes ( Figure 5B ). Of interest was that the cells arising from expanded clones (e.g. clonotypes 1, 2, 3) expressed CCL5 and granzymes, whereas several clones corresponding to either one or two cells (clonotype 64, 69, 71, 72, 73) did not express these markers ( Figure 5B ). Platypus enables the user to relate clonotype information to user-defined transcriptional signatures, thereby connecting phenotypes involving antibody secretion, T cell exhaustion, amongst others to immune receptor sequence.
Importantly, the user can customize the genes of interest to explore in the context of expanded clones.
Taken together, Platypus enabled us to investigate and visualize repertoire and transcriptome parameters that can help facilitate the rapid discovery of relevant B and T cell clones.
As in many areas of biology, researchers in adaptive immunity are poised to heavily utilize scSeq as an essential method to address basic and translational questions, therefore creating a substantial need for suitable bioinformatics tools. We have demonstrated here that Platypus enables rapid extraction, integration, and analysis of scSeq of lymphocyte repertoires and transcriptomes to uncover meaningful insight on lymphocyte immunobiology and function. By using convalescent COVID-19 patients as an example, we were able to detect clonally expanded B and T cells at the single-cell resolution from less than 10,000 cells, suggesting these clones are present at high-levels in the remainder of the blood repertoire. Interestingly, we observed many expanded clones were of the IgA isotype and often did not coexist with other isotypes, highlighting that IgA antibodies may play an important role in the response to SARS-CoV-2. We could relate these expanded clones to diverse cellular phenotypes based on district transcriptional cluster groups, which revealed that even when cells belonged to the same BCR or TCR clonal group, surprisingly they could still be found in different transcriptional-clusters. This highlights that lymphocytes sharing the same immune receptor specificity may still undergo very different cell fates and function in the course of an immune response. We could furthermore analyze various repertoire parameters including somatic hypermutation, clonal evolution, and clonal convergence at the single-cell resolution.
Platypus allowed us to analyze and extract full-length, paired, heavy-light chain information for clonally expanded and heavily mutated clones, which offers a possible approach to identify SARS-CoV-2-specific lymphocytes. Leveraging Platypus, we were able to visualize this information and have full-length BCR and TCR sequences ready for cloning in less than 10 lines of code -something that can greatly accelerate the discovery of adaptive immune therapeutics. In conclusion, Platypus enables a broad range of immunologists and bioinformaticians alike to gain quantitative insight at a single-cell resolution of how immune repertoire parameters relate to heterogeneous transcriptome information. 
The R package, accompanying code, and processed sequencing data used in this study are publicly available at github.com/alexyermanos/Platypus and doi.org/10.5281/zenodo.4140161. Briefly, clonotyping information was extracted directly from the output directory of cellranger using the function analyze_VDJ in Platypus. Quantifying the number of unique clones was performed using the VDJ_clonotype function in Platypus, with clone.strategy set to either "cdr3.aa", "hvj.lvj", hvj.lvj.cdr3length.cdr3homology", or "hvj.lvj.cdr3lengths". The isotype distributions for the top thirty B cell clones were calculated using the VDJ_isotypes_per_clone function in Platypus. CDR3 length distribution and sequence logo plots were calculated on the output of the VDJ_per_clone function. Sequence logos were calculated based on the function ggseqlogo from the R package ggseqlogo, with method set to "prob" and seq_type set to "aa".
The output directory from cellranger count was supplied as input to the function automate_GEX, which analyzes and integrates transcription data using functions from the R package Seurat. Briefly, the GEX libraries were integrated using the SCTransform function from Seurat. Cells containing more than 20% mitochondrial genes were removed. TCR and BCR genes were filtered prior to integration and gene expression analysis. The number of variable features selected was 2000 for the RunPCA function using the first 10 dimensions and Cluster resolution was set to the default 0.5. Transcriptional cluster and clonotype information were integrated using the VDJ_GEX_integrate function in Platypus under default parameters.
Quantification of the transcriptional cluster distribution for the 10 most expanded clones from patient 1 were calculated and visualised using the VDJ_GEX_expansion function in Platypus. Those cells from the two most expanded clones containing barcodes in both VDJ and GEX datasets were highlighted on the UMAP plot using the function visualize_clones_GEX in Platypus under default parameters. Somatic hypermutation was calculated in the VDJRegion as defined by MiXCR via the function call_MIXCR in Platypus, which utilized MiXCR (v3.0.12) 22 . The number of nucleotide alignment mismatches between the reference germline and the full-length VDJRegion for both heavy and light chain nucleotide sequence was then computed based on the best alignment determined by MIXCR. The phylogenetic tree was inferred by appending the full-length VDJRegion of the heavy and light chain for each unique sequence after appending VH VDJRegion and VL VDJRegion, as determined from the output function of call_MIXCR in Platypus. The reference germline sequence was first extracted from the initial cellranger alignment using the function VDJ_extract_germline in Platypus and added to the set of input sequences which were supplied to VDJ_tree. Sequence similarity networks were calculated using the function VDJ_networks in Platypus by calculating the edit distance separately for CDRH3 and CDRL3 amino acid sequences and then summing the two matrices. Edges were then drawn between those clones separated by less than 7 amino acid mutations. The heatmap integrating clonotype membership with user-defined gene lists was performed using the GEX_heatmap function in Platypus. Normalized expression of common T cell genes distributed by cluster. Each column represents a single cell and column identity refers to transcriptional cluster as previously determined ( Figure 3A ). The width of the column identity corresponds to the number of cells within each cluster.