The use of non-functional clonotypes as a natural calibrator for quantitative bias correction in adaptive immune receptor repertoire profiling

Version of Record

Accepted for publication after peer review and revision.

Download
Cite
Share
CommentOpen annotations (there are currently 0 annotations on this page).

731 views

Version of Record updated: February 9, 2023 (This version)
Version of Record published: February 6, 2023 (Go to version)
Accepted Manuscript updated: January 26, 2023 (Go to version)
Accepted Manuscript published: January 24, 2023 (Go to version)
Accepted: January 22, 2023
Received: April 6, 2021
Preprint posted: March 25, 2021 (Go to version)

1. Of interest
Gut Health: The value of connections

Vanessa Rossetto Marcelino

Insight Sep 19, 2023
Further reading

Abstract
Editor's evaluation
Introduction
Results
Discussion
Methods
Data availability
References
Decision letter
Author response
Article and author information
Metrics

Abstract

High-throughput sequencing of adaptive immune receptor repertoires is a valuable tool for receiving insights in adaptive immunity studies. Several powerful TCR/BCR repertoire reconstruction and analysis methods have been developed in the past decade. However, detecting and correcting the discrepancy between real and experimentally observed lymphocyte clone frequencies are still challenging. Here, we discovered a hallmark anomaly in the ratio between read count and clone count-based frequencies of non-functional clonotypes in multiplex PCR-based immune repertoires. Calculating this anomaly, we formulated a quantitative measure of V- and J-genes frequency bias driven by multiplex PCR during library preparation called Over Amplification Rate (OAR). Based on the OAR concept, we developed an original software for multiplex PCR-specific bias evaluation and correction named iROAR: immune Repertoire Over Amplification Removal (https://github.com/smiranast/iROAR). The iROAR algorithm was successfully tested on previously published TCR repertoires obtained using both 5’ RACE (Rapid Amplification of cDNA Ends)-based and multiplex PCR-based approaches and compared with a biological spike-in-based method for PCR bias evaluation. The developed approach can increase the accuracy and consistency of repertoires reconstructed by different methods making them more applicable for comparative analysis.

Editor's evaluation

This paper describes a newly developed, publicly available algorithm (iROAR) that was tested on pre-exisiting datasets and is of interest to T and B cell immunologists who perform repertoire analysis via multiplex PCR based techniques. iROAR utilises naturally occurring non-functional sequences to improve and partially correct the amplification bias inherent in multiplex PCR based sequencing technologies.

https://doi.org/10.7554/eLife.69157.sa0

Introduction

Adaptive immune receptor (TCR – T-cell receptor and BCR – B-cell receptor) repertoire is usually defined as a set of TCR or BCR sequences obtained from an individual’s blood, bone marrow, or specific lymphocyte population. Reflecting the T/B cell’s clonal composition, the repertoire is characterized by a high degree of specificity for each individual and substantial variation in clone frequencies. The accuracy of both sequences and frequencies of TCR/BCR genes in the obtained repertoire is essential to receiving the correct biological information from immune repertoire analysis.

High-throughput sequencing (HTS) of adaptive immune receptor repertoires is widely used in immunological studies (reviewed in Minervina et al., 2019) for the investigation of immune response to vaccines (Minervina et al., 2021; Pogorelyy et al., 2018; Sycheva et al., 2022), tumor-infiltrating lymphocytes (Gee et al., 2018; Goncharov et al., 2022; Oliveira et al., 2021), new therapeutic agents (Huang et al., 2019; Wang et al., 2018; Wilson et al., 2022), leukemia clonality, and minimal residual disease monitoring (Brüggemann et al., 2019; Komkov et al., 2020; Nazarov et al., 2016; Tirtakusuma et al., 2022; Wood et al., 2018). HTS-based methods for immune repertoire profiling use either RNA or DNA as a starting material and, in most cases, use PCR for the selective enrichment of receptor sequences. DNA-based methods generally use two-sided multiplex PCR with primers annealing to multiple V- and J-genes of the rearranged receptor (Brüggemann et al., 2019; Komkov et al., 2020; Robins et al., 2009). RNA-based methods start with cDNA synthesis, usually with TCR/BCR C(constant)-genes specific oligonucleotides, followed by one-side multiplex amplification with a set of V-gene specific primers and a universal C-gene specific primer (Wang et al., 2010). Alternatively, two universal primers are used for amplification if an artificial sequence is added to the 5’ end during synthesis using a template-switch (5’-Rapid Amplification of cDNA Ends [RACE]; Mamedov et al., 2013) or ligation (Oakes et al., 2017). DNA-based methods protect the repertoire from gene transcription bias and provide more comprehensive results (Barennes et al., 2020) which include most non-functional (out-of-frame) as well as functional (in-frame) rearrangements but produce high-amplification bias in the course of multiplex PCR. Additionally, each T/B cell contains a single DNA copy (i.e. two target strands) of the receptor molecule in contrast to tens of single-stranded RNA copies. RNA-based methods using 5’-RACE or ligation are characterized by the lowest PCR bias as they need a single primer pair for the amplification. However, the low efficiency of adding a universal oligo to the 5’ end makes its sensitivity comparable to or even lower than DNA-based methods. The compromise between these two approaches is the RNA-based method with a one-side multiplex that has moderate amplification bias yet sufficient sensitivity (Ma et al., 2018). Most bias in one-side multiplex RNA-based approaches could be removed by using unique molecular identifiers (UMIs; Ma et al., 2018). Unfortunately, for DNA-based methods, efficient incorporation of UMIs into the initial molecule before PCR is still challenging. The only method for DNA multiplex bias correction (Carlson et al., 2013) is undirected and cost-ineffective due to the utilization of an expensive synthetic spike-in control repertoire. Here, we propose an orthogonal solution for this challenge: the first fully computational algorithm for amplification bias detection and correction in adaptive immune receptor repertoires named iROAR (immune Repertoire Over Amplification Removal).

Results

The rationale for the Over Amplification Rate measure

Since out-of-frame TCR/BCR rearrangements do not form a functional receptor, they are not subjected to any specific clonal expansions and selection (Murugan et al., 2012). Being a passenger genomic variation, they change their initial (recombinational) clonal frequencies just randomly following the frequency changes of the second functional (in-frame) TCR/BCR allele present in the same T/B cell clone. According to the TCR/BCR loci rearrangement mechanism, the formation of in-frame and out-of-frame allele combinations in the same cell is also a stochastic and independent process in terms of V- and J-genes frequency. It leads to the conclusion that V- and J-gene frequencies among out-of-frame rearrangements must be sufficiently stable and must be equal to the initial recombination frequencies despite repertoire changes caused by various immune challenges (Figure 1). Thus, reproducible deviation of out-of-frame V- and J-gene frequencies (for the same multiplex PCR primer set) from the initial recombinational frequencies observed in the sequenced repertoire dataset is a result of artificial aberration caused by PCR amplification rather than immune repertoire evolution. Thus out-of-frame clonotypes can be considered a natural calibrator that can be used to measure amplification bias and quantitatively correct immune repertoire data.

Figure 1

Download asset Open asset

Stability of TRBV and TRBJ genes frequencies calculated based on unique out-of-frame rearrangements after Yellow fever vaccination (model of acute viral infection).

Out-of-frame clonotypes for frequencies calculation were extracted from low-biased 5’ Rapid Amplification of cDNA Ends (RACE) TRB repertoires of PBMC (peripheral blood mononuclear cells) samples obtained in two replicates for six time points: 0, 5, 10, 15, 21, and 45 d after YFV (Yellow Fever Vaccine) injection (donor M1, SRA accession number PRJNA577794 Minervina et al., 2020).

Figure 1—source data 1 XLSX table. TRB V and J frequencies for plots in Figure 1.: https://cdn.elifesciences.org/articles/69157/elife-69157-fig1-data1-v4.xlsx
Download elife-69157-fig1-data1-v4.xlsx

Formulating this observation, we developed the Over Amplification Rate (OAR) measure, which we define as a ratio of the observed and expected frequency of a V- (OAR[Vi]) or a J-gene (OAR[Ji]) among identified out-of-frame rearrangements. Observed frequency represents a value calculated as read counts (RCs) for each V- and J-gene (related to out-of-frames) divided by the sum of all out-of-frame clones RC in the obtained repertoire sequencing dataset. The expected frequency is a value before amplification calculated as a number of unique out-of-frame clones (UCN) having each V- or J-gene divided by the total number of unique out-of-frame clones in the repertoire. At the final stage, each OAR is normalized by dividing by the average OAR.

O A R (V i) = \frac{\frac{R C (V i)}{\sum_{1}^{N} R C (V i)}}{\frac{U C N (V i)}{\sum_{1}^{N} U C N (V i)}}

O A R (J i) = \frac{\frac{R C (J i)}{\sum_{1}^{N} R C (J i)}}{\frac{U C N (J i)}{\sum_{1}^{N} U C N (J i)}}

OAR value tends to be equal to 1 under ideal conditions (low or no amplification bias). It deviates from 1 as amplification bias increases in line: 5’-RACE with a single universal primer pair, one-side multiplex PCR (VMPlex), and two-side multiplex PCR (VJMPlex; Figure 2).

Figure 2 with 1 supplement see all

Download asset Open asset

Dependence of the Over Amplification Rate on the TCR repertoire profiling method.

(a) Comparison of Over Amplification Rate (OAR) values variances for TRB repertoires obtained with 5’-Rapid Amplification of cDNA Ends (RACE), one-side multiplex (VMPlex), and two-side multiplex (VJMPlex) PCR. The Levene’s test was performed to compare OAR variances: ****p<0.0001 and ***p<0.001. The bar and whiskers indicate a mean and SD. (b) Average (bold lines) OAR values for TRBV and TRBJ genes in repertoires obtained with 5’-RACE, one-side multiplex (VMPlex), and two-side multiplex (VJMPlex) PCR. Pale lines illustrate OARs of individual repertoires. Datasets: six repertoires for RACE from PRJNA847436 (Sycheva et al., 2022), six repertoires for VMPlex from PRJNA427746 (Ma et al., 2018), and six repertoires for VJMPlex from 27483#.XpCuQ1MzZQI (zenodo.org; Weinberger et al., 2015).

Figure 2—source data 1 XLSX table. Over Amplification Rate (OAR) values used for the plots in Figure 2.: https://cdn.elifesciences.org/articles/69157/elife-69157-fig2-data1-v4.xlsx
Download elife-69157-fig2-data1-v4.xlsx

The versatility of OAR measure

OAR measurement is a universal approach and can be applied to different types of immune repertoire data. To demonstrate this versatility, we calculated OAR values for low-biased (5’ RACE) repertoires of different adaptive immune receptor chains obtained from bulk human PBMC: TCR alpha (TRA), TCR beta (TRB), and BCR heavy chains (Figure 3a). The results show that OARs for both TCR and BCR repertoires obtained by 5’ RACE are close to 1 and stay within the range of 0.5–2, which is much narrower than OAR for multiplex PCR-based repertoires (see main text Figure 2).

Figure 3

Download asset Open asset

Over Amplification Rate of low-biased TCR and BCR repertoires in different lymphocyte subsets.

(a) Distribution of Over Amplification Rates of V- and J-genes in Rapid Amplification of cDNA Ends (RACE)-based repertoires of TCR and BCR (the empty dots represent average OARs among TCR repertoires: SRA accession numbers: PRJNA577794, PRJNA316572, PRJEB27352, and BCR repertoires: SRA accession number: PRJNA297771 and PRJNA494572). (b) Over Amplification Rates of V- and J-genes of TCR alpha chains in RACE-based repertoires of different types of T-cells (donors M1 and P30, 45 d after booster vaccination, SRA accession number PRJNA577794; Minervina et al., 2020). The bar and whiskers indicate a mean and SD.

Figure 3—source data 1 XLSX table. Over Amplification Rate (OAR) values used for the plots in Figure 3.: https://cdn.elifesciences.org/articles/69157/elife-69157-fig3-data1-v4.xlsx
Download elife-69157-fig3-data1-v4.xlsx

We also analyzed OARs for low-biased (5’ RACE) TCR repertoires of different T cell subpopulations, including T-helper, cytotoxic, central memory, effector memory, and naïve T cells. As shown in Figure 3, the OAR values demonstrate much less differences between analyzed T cell types then between RACE and multiplex PCR and are close enough to 1 similarly to the repertoire of bulk T cell mix obtained from PBMC.

Herewith, the variance of IGHV’s OARs compared TCRs’ and the variance of TCR subpopulations OARs compared PBMCs’ is slightly higher. This phenomenon may be linked to well-known differences in clonal expansion intensities of B/T-cell subsets which can affect indirectly the OAR values. However, the proof of this hypothesis demands separate deep analysis which is beyond the main focus of this research.

Despite it, our results demonstrate that OAR is a sufficiently universal measure of repertoires and can be applied to most adaptive immune receptors and cell types.

Factors affecting OAR measure accuracy

In the case of insufficient sequencing coverage, high-PCR bias can lead to the dramatic loss of clones and thus an incorrect measurement of V- and J-genes frequencies. In this instance, for the majority of V- and J-genes, the population frequencies can approximate the real frequencies better than multiplex repertoire-based ones (Figure 2—figure supplement 1). If upon comparison samples’ UCN-based frequencies significantly differ from the average frequencies calculated for the population (i.e. exceeds 99% CI), OAR calculation should be based on the latter.

Also, the balance of V- and J-genes frequencies can be disrupted by accidentally arisen abnormally large non-functional clonotypes generated in the course of abnormal clonal expansion in various lymphoproliferative disorders or stochastic spike in normal lymphocyte population. To reduce the impact of this anomaly on OAR value, the top clone of each V- and J-gene containing subgroups must be excluded from OAR calculation. Since V- and J-specific bias affects all clones non-selectively, the remaining large part of clones after top clones exclusion should be still representative for PCR bias calculation. As shown in Figure 4a, the exclusion of one top clonotype from OAR calculation for RACE-based TRB repertoire is enough to restore OAR calculation accuracy for TRBV2, TRBV5-6, TRBV7-9, and TRBV11-3. The further top clones exclusion has no significant effect on OAR values.

Figure 4

Download asset Open asset

Factors impacting Over Amplification Rate (OAR) calculation accuracy.

(a) Impact of highly proliferated top non-functional clonotype on OAR calculation accuracy in low-biased Rapid Amplification of cDNA Ends (RACE)-based TRB repertoire (Data: SRR19594184). (b) Impact of sequencing depth on OARs calculation error. (c) PCR bias independent changes of TRB V-genes OARs as a function of sequencing depth. Data: two-sided multiplex-based TRB repertoire (Data: RACE - SRR3129976, VJMPlex – SRR3129972).

Figure 4—source data 1 XLSX table. Over Amplification Rate (OAR) values used for the plots in Figure 4.: https://cdn.elifesciences.org/articles/69157/elife-69157-fig4-data1-v4.xlsx
Download elife-69157-fig4-data1-v4.xlsx

Another aspect impacting the accuracy of OAR calculation is the low-sequencing coverage of the TCR/BCR repertoire. The ratio of total RCs and the sum of unique clone counts can affect OAR value despite PCR bias solely because of the mathematical properties of the OAR formula. In the extreme case, the OAR value (OAR = 1) for V- and J-genes represented in a single out-of-frame clone with only one read will not reflect the real amplification bias. To address this issue, we analyzed the OAR calculation error as a function of the number of reads per clone used for the OARs evaluation (Figure 4b). For this purpose, we performed a serial down-sampling of TCR datasets generated by RACE and two-side multiplex PCR and calculated OAR measurement error for each dataset portion. OAR calculated for the entire dataset was taken as a benchmark. The result shows that 1.8 (for MPlex) and 2.5 (for RACE) reads per out-of-frame clonotype are a minimal sufficient sequencing coverage to get adequate OAR values with an acceptable error rate of ~10%.

It is also important to note that errors in nucleotide sequences occurring during library preparation and sequencing could lead to an artificial increase in both in-frame and out-of-frame clone diversity. Single nucleotide substitutions generate artificial clones as a branch of real most abundant clones inside of each in-frame and out-of-frame group independently. Single nucleotide indels lead to cross-generation artificial clones between groups: real in-frame clones generate false out-of-frame clones and vice versa. Artificial clones compromise the accuracy of both repertoire itself and OAR value. To eliminate such clones generated by single-nucleotide substitutions, we filtered them out by the VDJTOOLS software (see Methods section). To eliminate artificial clones produced by indels, we searched for in-frame and out-of-frame clone pairs which differ by one indel (Levenshtein distance = 1). If their ratio is less than 1:500, the smaller clone in pair is discarded, and its count is added to the count of the larger clone (this procedure guarantees to discard most sequencing errors present in 1 per 1000 nucleotides average).

OAR-index

To estimate the value of immune repertoire structure disruption by amplification bias, we proposed the OAR-index, which represents the mean square deviation of OARs for each V- and J-gene from the value characteristic for repertoire with no bias (OAR = 1). OAR-index is directly proportional to the amplification bias and thus can be used for rapid estimation and comparison of immune repertoire bias. The less OAR-index is, the less PCR bias is with an ideally unbiased repertoire having OAR-index=0.

O A R - i n d e x = \sqrt[2]{\frac{\sum_{0}^{n} {(O A R i - 1)}^{2}}{n}}

Using OAR for the removal of amplification bias

Normalization coefficients for each VJ combination are estimated by multiplication of corresponding V- and J-gene OARs for two-side multiplex and V-gene OAR for one-side multiplex (Figure 2a). The corrected RC for each clonotype with the particular V-J gene combination is obtained simply by dividing the observed RC by the corresponding normalization coefficient. OAR of V- and J-genes could be co-dependent, which can be a reason for overcorrection. To avoid this issue, the procedure can be recursively repeated with a modified normalization coefficient defined as described coefficient raised to the power of a number in the range from 0 to 1 (parameter ‘mt’). The corrected RCs are used to estimate the real percentage of each clonotype in the repertoire. However, the all multiplex-based repertoires analyzed in actual study required just one iteration with mt = 1. A detailed flowchart of the OAR-based amplification bias correction algorithm named iROAR is shown in Figure 5a.

Figure 5

Download asset Open asset

Implementation of OAR-based correction of quantitative bias in TCR repertoire.

(a) Flowchart of iROAR (immune Repertoire Over Amplification Removal) algorithm. UCF – Frequency calculated using unique clones counts (denominator of Over Amplification Rate [OAR]), pUCF – population UCF, RC – read count, RCn – normalized RC, RCo – observed RC, and mt – the number in the range from 0 to 1 for the iterative procedure. (b) Clone frequencies in the low biased 5’-Rapid Amplification of cDNA Ends (RACE)-based repertoire (ENA database, accession number ERR2869430) vs. the same repertoire with introduced artificial bias: before and after iROAR processing. (c) Out-of-frame (non-functional) clone frequencies in low biased 5’-RACE-based repertoire vs. two-side multiplex (VJMPlex)-based repertoire obtained for the same RNA sample (SRA database, accession numbers SRR3129976 and SRR3129972): before and after iROAR processing. (d) In-frame (functional) clone frequencies in low biased 5’-RACE-based repertoire vs. two-side multiplex (VJMPlex)-based repertoire obtained for the same RNA sample: before and after iROAR processing. SRA database, accession numbers SRR3129976 and SRR3129970. R² is the squared Pearson correlation coefficient. iROAR was applied only for biased repertoires: artificially biased RACE and VJMPlex. (e) OARV and OARJ of test 5’-RACE-based TRB repertoire (b) before artificial bias introduction (green dots and line), biased one (red dots and line) and corrected one by iROAR (blue dots and line).

Figure 5—source data 1 XLSX table. Clonal frequencies used for the plots in Figure 5.: https://cdn.elifesciences.org/articles/69157/elife-69157-fig5-data1-v4.xlsx
Download elife-69157-fig5-data1-v4.xlsx

OAR-based approach validation

The validation of OAR-based amplification bias correction was performed on the TRB dataset with in silico introduced bias generated from real (experimental) low-biased (5’-RACE) repertoire (Figure 5b). After correction, the OAR-index indicating general repertoire bias expectedly decreased from 1.81 to 0.76. Interestingly, the OAR independent measure R² value of in silico biased and original repertoire correlation raised from 0.5350 to 0.8915, confirming the substantial reduction of in silico introduced quantitative bias. Afterward, we tested our approach on real paired experimental datasets obtained from the same RNA sample by two different method types: 5’-RACE and multiplex PCR (Barennes et al., 2020; Liu et al., 2016; Figure 5c–d, Figure 6).

Figure 6

Download asset Open asset

Effect of iROAR (immune Repertoire Over Amplification Removal)-based PCR bias correction in MPlex repertoire on similarity with low-biased Rapid Amplification of cDNA Ends (RACE)-based repertoire obtained from the same RNA sample.

(a) Pearson correlation coefficient; (b) R² measure; (c) Morisita-Horn similarity index. **** p<0.0001 (two-tailed Wilcoxon matched-pairs signed rank test, CI = 0.95). Dataset: PRJNA548335 three different RACE (RACE-2, RACE-3, and RACE-4 in six replicates each) protocols vs. RNA-based MPlex (Multiplex-3) protocol for Donor1 and Donor2 (100 ng RNA input): total 36 points; PRJNA309577 (one RACE protocol vs. one MPlex protocol) for Donors S01 (four MPlex replicates vs. two RACE replicates), S02 (two MPlex replicates vs. four RACE replicates), and donor S03 (one MPlex replicate vs. one RACE replicate): total 17 points.

Figure 6—source data 1 XLSX table. Pearson R, R², and Morisita-Horn index values used for the plots in Figure 6.: https://cdn.elifesciences.org/articles/69157/elife-69157-fig6-data1-v4.xlsx
Download elife-69157-fig6-data1-v4.xlsx

As a result of amplification bias correction, OAR-index for multiplex-based repertoire decreased 1.5-fold average. At the exact time, the correlation of clonal frequencies obtained with RACE and multiplex significantly increased (Pearson correlation measure and R² value increased 1.5-fold average each) with a significant rise of repertoires similarity (Morisita-Horn index increased 1.7-fold average; Figure 6). Importantly, amplification bias decreased in both out-of-frame and in-frame clone sets, although normalization coefficients were calculated using out-of-frame ones only.

Comparison of iROAR and spike-in-based approach for amplification bias detection

Biological spike-in is considered a classical technique for multiplex PCR bias evaluation. Several options for this technique including synthetic repertoire (Carlson et al., 2013; Wu et al., 2020), lymphoid cell lines DNA mix, and DNA from human blood, tonsil, and thymus (Kallemeijn et al., 2018; Knecht et al., 2019) were established to measure V- and J-segment specific primers performance during TCR/BCR rearrangements amplification in multiplex PCR. In this study, we compared iROAR-based amplification bias evaluation with a spike-in-based approach. Similarly to Kallemeijn et al., 2018; Knecht et al., 2019, we were using natural thymic cell-derived spike-ins rather than synthetic ones. Human CD8 T-cells derived DNA was used as a target input for the libraries’ preparation. TRA rearrangements library of thymocytes were used as a source of spike-ins. Two different random mixes of TRAV- and TRAJ-specific primers (0.18–4.7 μM each) were used for multiplex PCR amplification of target DNA with spike-in added. Each test library was prepared in two replicas (four test libraries total). The obtained libraries were sequenced with an average coverage of 9.88 reads per clonotype and contained 35,818–40,209 target and 2298–3571 spike-in clonotypes after pseudogenes removal (Table 1).

Table 1

The number of spike-in and target clonotypes in test TRA libraries.

Sample	Spike-in clonotypes			Target clonotypes
Sample	Number	Read count	Coverage	Number	Read count	Coverage
Primer mix 1 Replica 1	3571	30,474	8.53	39,911	348,792	8.74
Primer mix 1 Replica 2	2698	19,420	7.20	35,818	303,494	8.47
Primer mix 2 Replica 1	3439	34,717	10.10	40,209	425,508	10.58
Primer mix 2 Replica 2	2298	24,823	10.80	33,406	383,615	11.48

Multiplex PCR bias of each separate V- and J-gene was calculated using both iROAR and biological spike-in approaches demonstrating the high-correlation level (Figure 7a and b) for the matched OAR/spike-in pairs (Pearson’s r=0.78 average) in contrast to mismatched ones (Pearson’s r=0.46 average). VJ combination bias for both approaches was calculated by multiplying V- and J-segment biases and compared using correlation analysis (Figure 7c–e). iROAR and spike-in detected VJ biases showed a strong positive correlation (Pearson’s r=0.7182 average) for all four test TRA libraries (Figure 7c). Based on replicas comparison, the reproducibility of iROAR detected VJ bias was higher than one detected using spike-in control (Figure 7d and e).

Figure 7

Download asset Open asset

Comparison of Over Amplification Rate (OAR)-based and biological spike-in-based approaches for multiplex PCR bias detection.

Pearson’s correlation coefficient for V-segments bias measure (a) and J-segments bias measure (b). Column and row titles: PM = Primer mix, R=replica. (c) Correlation of VJ combination bias calculated by iROAR (immune Repertoire Over Amplification Removal) and biological spike-ins. (d) Reproducibility of iROAR-based VJ combination bias detection. (e) Reproducibility of spike-in-based VJ combination bias detection. Data: PRJNA825832.

Figure 7—source data 1 XLSX table. iROAR (immune Repertoire Over Amplification Removal)-based and spike-in-based PCR bias measures used for the plots in Figure 7.: https://cdn.elifesciences.org/articles/69157/elife-69157-fig7-data1-v4.xlsx
Download elife-69157-fig7-data1-v4.xlsx

Impact of iROAR on a similarity of repertoires prepared by different multiplex PCR systems

To further test the iROAR approach’s ability to raise the uniformity of repertoires by reducing multiplex PCR-specific bias, we analyzed changes in the similarity of repertoires prepared for the same individual but using different multiplex methods. For this purpose, we compared OARs, V/J, and clonotype frequencies before and after bias correction using iROAR in test TRA libraries prepared with Primer mix 1 and Primer mix 2 (after spike-in removal). As a result of iROAR-based bias correction, the difference between OARs for these two library types significantly decreases, and OARs themselves approach a value of one. By default, iROAR does not affect the diversity of repertoires and does not remove any clonotypes. Meanwhile, V and J frequencies are subject to substantial changes (Figure 8b) depending on the initial bias level. These changes occur in both biased repertoires (Primer mix 1 and Primer mix 2) and lead to an increase its convergence (Figure 8d). Herewith, R² measure increased 1.31-fold and 1.5-fold for V- and J-gene frequencies, respectively. Moreover, bias correction using iROAR also increases similarities of clone frequencies (Figure 8c). In this case, both the Morisita-Horn index and Pearson correlation coefficient increase twofold and R² measure increases 1.5-fold.

Figure 8

Download asset Open asset

Convergence of Over Amplification Rate (OAR), clonotype, and V/J frequencies between two TRA repertoire before and after iROAR (immune Repertoire Over Amplification Removal) based bias correction.

(a) OAR values changes in four test TRA libraries after PCR bias correction using iROAR. (b) TRAV and TRAJ frequency changes after PCR bias correction using iROAR (Sample: Primer mix 1 Replica 1). (c) Correlation of clonal frequencies of two different types of test TRA repertoires before and after iROAR-based PCR bias correction. (d) Correlation of V- and J-gene frequencies of test TRA repertoires before and after iROAR-based PCR bias correction.

Figure 8—source data 1 XLSX table. Clonal frequencies, TRB V and J frequencies and Over Amplification Rate (OAR) values used for the plots in Figure 8.: https://cdn.elifesciences.org/articles/69157/elife-69157-fig8-data1-v4.xlsx
Download elife-69157-fig8-data1-v4.xlsx

It is important to note that OARs calculation and bias correction for each of the analyzed test TRA repertoires was performed entirely independently without the involvement of any common normalization coefficients or spike-in controls. Therefore, each repertoire contains enough information to correct it adequately, increasing the consistencies of interrogated repertoires obtained even by different multiplex PCR protocols.

All observed results can be considered evidence of the actual capacity of iROAR approach to accurately detect and reduce multiplex-specific quantitative bias in adaptive immune receptor repertoires.

Discussion

Even a small difference in amplification efficiencies can lead to a massive bias after multiple amplification cycles due to the exponential nature of PCR. Thus, most of the existing immune repertoire library preparation methods are subjected to amplification bias. The effect of distinct PCR bias-generating factors can be reduced experimentally by varying reaction mixture content and introducing special protocols (UMI, crafty primer structures, and spike-in controls). However, the criteria for estimating and removing the residual bias after applying these optimization approaches are lacking. Here, we close this gap by introducing the OAR value and OAR-index, which score PCR bias for both V- and J-genes separately (OAR values) and the whole repertoire dataset (OAR-index). Based on OAR values, we developed the first fully computational approach to decipher and correct amplification bias in adaptive immune receptor repertoires produced by one- or two-side multiplex PCR-based methods, using RNA or DNA as a template. Due to the inability to use UMI-based correction for DNA-based multiplex, the developed approach is the only currently available technique allowing direct measuring and correcting PCR bias in such repertoires without additional experiments.

In contrast to cell-line mix spike-in (Knecht et al., 2019) or synthetic repertoire-based (Carlson et al., 2013; Wu et al., 2020) PCR bias correction, the proposed approach operates with hundreds and even thousands of natural calibrators (out-of-frame clones) for each V-J gene pair. It makes this method potentially more reliable due to the ability to minimize the impact of CDR3 structure on PCR bias calculation since out-of-frame captures significantly higher CDR3 diversity than biological spike-ins. Moreover, similarly to a previously described method (Carlson et al., 2013), the OAR-based approach can also be used for primer efficacy evaluation to optimize their structures and concentrations, which in turn will straighten the coverage of various V- and J-genes and minimize the number of experimentally lost clones. Being fully computational, the developed PCR bias correction algorithm can be easily implemented in any TCR/BCR repertoire analysis pipeline, noticeably improving the quantitative parts of the analysis. Even though it’s not possible to fully substitute the low-biased RACE methods, iROAR is capable to make multiplex-PCR-based repertoires more consistent with RACE-based ones. Therefore, the developed approach can provide the opportunity to compare the immune repertoire datasets generated using different library preparation methods.

Methods

Raw-data processing and immune repertoires reconstruction

All sequencing data used in this study represent human TCR and BCR repertoires. The repertoires (see Supplementary file 1) were reconstructed from fastq data using MiXCR v2 software (Bolotin et al., 2017; Bolotin et al., 2015) after primers and adapters trimming using FASTP software (Chen et al., 2018). All obtained repertoires were converted to VDJTOOLS (Shugay et al., 2015) format for unification. Erroneous clones generated by single nucleotide substitution were removed from the repertoire using the ‘Correct’ function from the VDJTOOLS software package. Erroneous clones generated by single nucleotide indels were removed from repertoires using the ‘Filter’ function from developed iROAR software. V and J pseudogenes were removed from repertoires using the ‘FilterBySegment’ function of VDJTOOLS.

TRA repertoires preparation

The peripheral blood was collected from a healthy volunteer from the article’s co-authors with informed consent in a certified clinical lab. PBMC was separated from whole blood using the Ficoll-Paque approach. CD8+ T cells were isolated using Dynabeads CD8 Positive Isolation Kit (Invitrogen). DNA for library preparation was extracted from CD8+ T cells using FlexiGene DNA Kit (Qiagen). 150 ng aliquots of obtained DNA were used as input to prepare each out of four TRA libraries. Each DNA aliquot was premixed with 0.1 pg of serial diluted low-biased TRA library (prepared using MiLaboratories Human TCR kit) of thymic cells (spike-in matrix) as biological spike-ins. Two pools of previously designed (Komkov et al., 2020) TRAV- and TRAJ-specific primers (MiLaboratories LLC) with randomly selected concentrations (0.18–4.7 μM each) were generated to produce two types of TRA libraries with different quantitative bias status simulating libraries produced by different multiplex PCR methods. Library preparation was performed according to the protocol from Komkov et al., 2020. Both types of TRA libraries were prepared in two replicas and sequenced along with a spike-in matrix library on the MiSeq Illumina instrument (SE 150 nt) with moderate coverage 480,000 reads per library.

Biological spike-in detection and analysis

TRA repertoires were extracted from FASTQ files using MIXCR software. All obtained MIXCR output files were converted to VDJTOOLS format as described above. Extraction of the spike-in sequences and spike-in-free repertoires from sequenced libraries was performed using the VDJTOOLS function ‘ApplySampleAsFilter’ and the sequenced spike-in library as a filter. Spike-in-based amplification bias was calculated as the quotient of V- and J-frequency in spike-ins extracted from target libraries and corresponding V- and J-frequency in the spike-in matrix, which was not subjected to multiplex amplification. OARs for obtained TRA libraries were calculated using iROAR software and spike-in-free repertoires as input. VJ bias values were calculated by multiplying V- to J-segment-specific biases. Correlation analysis of iROAR and spike-in VJ bias values was performed using GraphPad Prism9 software.

Step-by-step pipeline for the OAR evaluation used in this study

Single nucleotide error correction in read1/read2 intersected sequences and Illumina adapters removal (optional): fastp -c -i input_R1.fastq.gz -I input_R2.fastq.gz -o fp.input_R1.fastq.gz -O fp.input_R2.fastq.gz.
Raw reads alignment (essential):
1. For TCR beta chainsmixcr align -c TRB fp.input_R1.fastq.gz fp.input_R2.fastq.gz output1.vdjca.
2. For TCR beta chainsmixcr align -c TRA fp.input_R1.fastq.gz fp.input_R2.fastq.gz output1.vdjca.
Clonotypes assemble (essential): mixcr assemble output1.vdjca output2.clns
TCR repertoire export in a human-readable format (essential): mixcr export Clones output2.clns clones.txt.
Convert repertoire into VDJtools format (essential):java -jar vdjtools.jar Convert -S mixcr clones.txt vdjtools.
Artificial diversity removal by single nucleotide substitutions correction (optional):java -jar vdjtools.jar Correct vdjtools.clones.txt correct.
Pseudogenes removal (optional):
1. For TCR beta chains.java -jar vdjtools.jar FilterBySegment –-j-segment TRBJ2-2P --v-segment TRBV1,TRBV12-1,TRBV12-2,TRBV17,TRBV21-1,TRBV22-1,TRBV23-1,TRBV26,TRBV5-2,TRBV5-3,TRBV5-7,TRBV6-7,TRBV7-1,TRBV7-5,TRBV8-1,TRBV8-2 --negative correct.vdjtools.clones.txt filter
2. For TCR alpha chainsjava -jar vdjtools.jar FilterBySegment -–j-segment TRAJ1,TRAJ19,TRAJ2,TRAJ25,TRAJ51,TRAJ55,TRAJ58,TRAJ59,TRAJ60,TRAJ61 --v-segment TRAV11,TRAV11-1,TRAV14-1,TRAV15,TRAV28,TRAV31,TRAV32,TRAV33,TRAV37,TRAV46,TRAV8-5,TRAV8-6-1,TRAV8-7 --negative correct.vdjtools.clones.txt filter
Artificial diversity removal by single nucleotide indels correction (optional):iroar Filter -se 0.01 filter.correct.vdjtools.clones.txt filter2.txt
OARs calculation and quantitative bias correction (essential):iroar Count -min_outframe 15 r -z 1 -iter 1 -mt 1 input_folder output_folder (input_folder must contain filter2.txt file).

OAR evaluation and statistical analysis

For the OAR and OAR-index calculation and amplification bias removal, we used the command-line-based iROAR software designed in this study and freely available for non-profit use at GitHub (https://github.com/smiranast/iROAR; Komkov, 2023; copy archived at swh:1:rev:2362c4f41d40519154e1c2dc6ce7af619f15fb4b). For the OAR comparison between 5’-RACE, one-side, and two-side multiplex PCRs, an equal number of out-of-frame clones (50,000) was randomly selected from TCR repertoires of 15 healthy individuals (for each approach). Average population V- and J-gene frequencies (unweighted) were calculated based on out-of-frame clones from 105 TRB repertoires obtained by two methods: 5’-RACE (95 repertoires) and single-cell TCR profiling (10× Genomics; 10 repertoires; Supplementary file 1) using the ‘CalcSegmentUsage’ function with ‘-u’ parameter of VDJTOOLS. All statistical tests were performed using Prism9 GraphPad software (https://www.graphpad.com/).

iROAR software requirement

Recommended system configuration for iROAR running: Linux or MacOS, 2 CPU, 8 GB RAM, programming language: python = 3.7.3, required Python packages: matplotlib = 3.0.3, numpy = 1.16.2, pandas = 0.24.2, and requests = 2.21.0. Starting iROAR package includes list of average populational frequencies with SDs of TRB and TRA V- and J-genes related to the European population. iROAR run command: iroar Count (optional parameters) <input> <output>. Recommended parameters for most tasks: -min_outframe 15r -z 1 -iter 1 -mt 1. Full list of available parameters is deposited in project directory at github (https://github.com/smiranast/iROAR).

Data access

All analyzed datasets were downloaded from open-source databases: NCBI SRA (https://www.ncbi.nlm.nih.gov/sra), ENA (https://www.ebi.ac.uk/ena), and Zenodo project (https://zenodo.org/). A complete list of web links and accession numbers is summarized in Supplementary file 1. TRA repertoire dataset generated in this study for iROAR validation is available under access number PRJNA825832.

The iROAR software and its documentation are available at the link: https://github.com/smiranast/iROAR. The additional software used in this study is available in the GitHub repository: MiXCR v2 (Bolotin et al., 2015) (https://github.com/milaboratory/mixcr; MiLaboratories, 2023), VDJTOOLS (Shugay et al., 2015) ( https://github.com/mikessh/vdjtools; Shugay, 2022), and FASTP (Chen et al., 2018) (https://github.com/OpenGene/fastp; OpenGene - Open Source Genomics Toolbox, 2022).

Data availability

Sequencing data have been deposited in SRA under accession code PRJNA825832. All other sequencing data analyzed during this study are previously published and fully available under links or access numbers included in the manuscript and supporting files.

The following data sets were generated

1. Smirnova A
(2022) NCBI BioProject
ID PRJNA825832. TRA repertoire.

https://www.ncbi.nlm.nih.gov/bioproject/PRJNA825832/

The following previously published data sets were used

1. Warren RL
2. Freeman JD
3. Zeng T
4. Choe G
5. Munro S
(2011) NCBI Sequence Read Archive
ID SRA020989. Exhaustive T-cell repertoire sequencing of human peripheral blood samples reveals signatures of antigen selection and a directly measured repertoire size of at least 1 million clonotypes.

https://www.ncbi.nlm.nih.gov/sra/?term=SRA020989
1. Zvyagin IV
(2014) NCBI BioProject
ID PRJNA214848. Homo sapiens T-cell repertoire - MZ twins.

https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA214848
1. Rosati E
(2020) NCBI BioProject
ID PRJEB27352. TCR repertoire in IBD twins.

https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJEB27352
1. de Greef PC
(2020) NCBI BioProject
ID PRJNA390125. TCR repertoire sequencing of T cell subsets from healthy individuals.

https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA390125
1. Pogorelyy MV
(2018) NCBI BioProject
ID PRJNA493983. Precise tracking of vaccine-responding T-cell clones reveals convergent and personalized response in identical twins.

https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA493983
1. Turchaninova et al
(2016) NCBI BioProject
ID PRJNA297771. Protocol for full length profiling of IG repertoires.

https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA297771
1. Minervina et al
(2020) NCBI BioProject
ID PRJNA577794. Comprehensive analysis of antiviral adaptive immunity formation and reactivation down to single cell level.

https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA577794
1. Simon et al
(2018) NCBI BioProject
ID PRJNA494572. Sequencing the Peripheral Blood B and T cell Repertoire - Quantifying robustness and limitations.

https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA494572
1. Pan et al
(2019) NCBI BioProject
ID PRJNA550004. Identification of drug-specific public TCR driving severe cutaneous adverse reactions.

https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA550004
1. Ma et al
(2018) NCBI BioProject
ID PRJNA427746. Homo sapiens Raw sequence reads.

https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA427746
1. Truong et al
(2019) NCBI BioProject
ID PRJEB31283. TCR diversity and clonality of human CD4+ memory T cells.

https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJEB31283
1. Simnica et al
(2019) NCBI BioProject
ID PRJEB33490. Immunoaging.

https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJEB33490
(2015) Zenodo
Immune repertoire profiling reveals that clonally expanded B and T cells infiltrating diseased human kidneys can also be tracked in the blood.
https://doi.org/10.5281/zenodo.27483
1. Tanno et al
(2019) NCBI BioProject
ID PRJNA593622. Paired TCR alpha:TCR beta sequencing at the single-cell level.

https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA593622
1. Liu et al
(2016) NCBI BioProject
ID PRJNA309577. TRB and IGH are captured from peripheral blood using Multiplex PCR and 5'RACE.

https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA309577
1. Barennes T
(2020) NCBI BioProject
ID PRJNA548335. Benchmarking of T cell receptor repertoire profiling methods reveals large systematic biases.

https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA548335

References

1. Barennes P
2. Quiniou V
3. Shugay M
4. Egorov ES
5. Davydov AN
6. Chudakov DM
7. Uddin I
8. Ismail M
9. Oakes T
10. Chain B
11. Eugster A
12. Kashofer K
13. Rainer PP
14. Darko S
15. Ransier A
16. Douek DC
17. Klatzmann D
18. Mariotti-Ferrandiz E
(2020) Benchmarking of T cell receptor repertoire profiling methods reveals large systematic biases
Nature Biotechnology 39:236–245.
https://doi.org/10.1038/s41587-020-0656-3
- PubMed
- Google Scholar
(2015) MiXCR: software for comprehensive adaptive immunity profiling
Nature Methods 12:380–381.
https://doi.org/10.1038/nmeth.3364
- PubMed
- Google Scholar
(2017) Antigen receptor repertoire profiling from RNA-Seq data
Nature Biotechnology 35:908–911.
https://doi.org/10.1038/nbt.3979
- PubMed
- Google Scholar
1. Brüggemann M
2. Kotrová M
3. Knecht H
4. Bartram J
5. Boudjogrha M
6. Bystry V
7. Fazio G
8. Froňková E
9. Giraud M
10. Grioni A
11. Hancock J
12. Herrmann D
13. Jiménez C
14. Krejci A
15. Moppett J
16. Reigl T
17. Salson M
18. Scheijen B
19. Schwarz M
20. Songia S
21. Svaton M
22. van Dongen JJM
23. Villarese P
24. Wakeman S
25. Wright G
26. Cazzaniga G
27. Davi F
28. García-Sanz R
29. Gonzalez D
30. Groenen PJTA
31. Hummel M
32. Macintyre EA
33. Stamatopoulos K
34. Pott C
35. Trka J
36. Darzentas N
37. Langerak AW
38. EuroClonality-NGS working group
(2019) Standardized next-generation sequencing of immunoglobulin and T-cell receptor gene recombinations for MRD marker identification in acute lymphoblastic leukaemia; a euroclonality-NGS validation study
Leukemia 33:2241–2253.
https://doi.org/10.1038/s41375-019-0496-7
- PubMed
- Google Scholar
(2013) Using synthetic templates to design an unbiased multiplex PCR assay
Nature Communications 4:2680.
https://doi.org/10.1038/ncomms3680
- PubMed
- Google Scholar
1. Chen S
2. Zhou Y
3. Chen Y
4. Gu J
(2018) Fastp: an ultra-fast all-in-one FASTQ preprocessor
Bioinformatics 34:i884–i890.
https://doi.org/10.1093/bioinformatics/bty560
- PubMed
- Google Scholar
1. Gee MH
2. Han A
3. Lofgren SM
4. Beausang JF
5. Mendoza JL
6. Birnbaum ME
7. Bethune MT
8. Fischer S
9. Yang X
10. Gomez-Eerland R
11. Bingham DB
12. Sibener LV
13. Fernandes RA
14. Velasco A
15. Baltimore D
16. Schumacher TN
17. Khatri P
18. Quake SR
19. Davis MM
20. Garcia KC
(2018) Antigen identification for orphan T cell receptors expressed on tumor-infiltrating lymphocytes
Cell 172:549–563.
https://doi.org/10.1016/j.cell.2017.11.043
- PubMed
- Google Scholar
(2022) Pinpointing the tumor-specific T cells via TCR clusters
eLife 11:e77274.
https://doi.org/10.7554/eLife.77274
- PubMed
- Google Scholar
1. Huang J
2. Alexey S
3. Li J
4. Jones T
5. Grande G
6. Douthit L
7. Xie J
8. Chen D
9. Wu X
10. Michael M
11. Xiao C
12. Zhao J
13. Xie X
14. Xie J
15. Chen XL
16. Fu G
17. Alexander G
18. Tzeng C-M
(2019) Unique CDR3 epitope targeting by CAR-T cells is a viable approach for treating T-cell malignancies
Leukemia 33:2315–2319.
https://doi.org/10.1038/s41375-019-0455-3
- PubMed
- Google Scholar
(2018) Next-generation sequencing analysis of the human TCRγδ+ T-cell repertoire reveals shifts in Vγ- and Vδ-usage in memory populations upon aging
Frontiers in Immunology 9:448.
https://doi.org/10.3389/fimmu.2018.00448
- PubMed
- Google Scholar
1. Knecht H
2. Reigl T
3. Kotrová M
4. Appelt F
5. Stewart P
6. Bystry V
7. Krejci A
8. Grioni A
9. Pal K
10. Stranska K
11. Plevova K
12. Rijntjes J
13. Songia S
14. Svatoň M
15. Froňková E
16. Bartram J
17. Scheijen B
18. Herrmann D
19. García-Sanz R
20. Hancock J
21. Moppett J
22. van Dongen JJM
23. Cazzaniga G
24. Davi F
25. Groenen PJTA
26. Hummel M
27. Macintyre EA
28. Stamatopoulos K
29. Trka J
30. Langerak AW
31. Gonzalez D
32. Pott C
33. Brüggemann M
34. Darzentas N
35. EuroClonality-NGS Working Group
(2019) Quality control and quantification in IG/TR next-generation sequencing marker identification: protocols and bioinformatic functionalities by euroclonality-NGS
Leukemia 33:2254–2265.
https://doi.org/10.1038/s41375-019-0499-4
- PubMed
- Google Scholar
(2020) High-Throughput sequencing of T-cell receptor alpha chain clonal rearrangements at the DNA level in lymphoid malignancies
British Journal of Haematology 188:723–731.
https://doi.org/10.1111/bjh.16230
- PubMed
- Google Scholar
Software
1. Komkov A
(2023) IROAR, version swh:1:rev:2362c4f41d40519154e1c2dc6ce7af619f15fb4b
Software Heritage.

https://archive.softwareheritage.org/swh:1:dir:6d7201302952a15597cfb50c2d4917c22efd966f;origin=https://github.com/smiranast/iROAR;visit=swh:1:snp:e92f096a8db6efcdc50263dba10a70e759b225f7;anchor=swh:1:rev:2362c4f41d40519154e1c2dc6ce7af619f15fb4b
1. Liu X
2. Zhang W
3. Zeng X
4. Zhang R
5. Du Y
6. Hong X
7. Cao H
8. Su Z
9. Wang C
10. Wu J
11. Nie C
12. Xu X
13. Kristiansen K
(2016) Systematic comparative evaluation of methods for investigating the TCRβ repertoire
PLOS ONE 11:e0152464.
https://doi.org/10.1371/journal.pone.0152464
- PubMed
- Google Scholar
1. Ma KY
2. He C
3. Wendel BS
4. Williams CM
5. Xiao J
6. Yang H
7. Jiang N
(2018) Immune repertoire sequencing using molecular identifiers enables accurate clonality discovery and clone size quantification
Frontiers in Immunology 9:33.
https://doi.org/10.3389/fimmu.2018.00033
- PubMed
- Google Scholar
(2013) Preparing unbiased T-cell receptor and antibody cdna libraries for the deep next generation sequencing profiling
Frontiers in Immunology 4:456.
https://doi.org/10.3389/fimmu.2013.00456
- PubMed
- Google Scholar
Software
1. MiLaboratories
(2023) Mixcr, version a2a8ed8
GitHub.

https://github.com/milaboratory/mixcr
(2019) T-Cell receptor and B-cell receptor repertoire profiling in adaptive immunity
Transplant International 32:1111–1123.
https://doi.org/10.1111/tri.13475
- PubMed
- Google Scholar
1. Minervina AA
2. Pogorelyy MV
3. Komech EA
4. Karnaukhov VK
5. Bacher P
6. Rosati E
7. Franke A
8. Chudakov DM
9. Mamedov IZ
10. Lebedev YB
11. Mora T
12. Walczak AM
(2020) Primary and secondary anti-viral response captured by the dynamics and phenotype of individual T cell clones
eLife 9:e53704.
https://doi.org/10.7554/eLife.53704
- PubMed
- Google Scholar
1. Minervina AA
2. Komech EA
3. Titov A
4. Bensouda Koraichi M
5. Rosati E
6. Mamedov IZ
7. Franke A
8. Efimov GA
9. Chudakov DM
10. Mora T
11. Walczak AM
12. Lebedev YB
13. Pogorelyy MV
(2021) Longitudinal high-throughput TCR repertoire profiling reveals the dynamics of T-cell memory formation after mild COVID-19 infection
eLife 10:e63502.
https://doi.org/10.7554/eLife.63502
- PubMed
- Google Scholar
1. Murugan A
2. Mora T
3. Walczak AM
4. Callan CG
(2012) Statistical inference of the generation probability of T-cell receptors from sequence repertoires
PNAS 109:16161–16166.
https://doi.org/10.1073/pnas.1212755109
- PubMed
- Google Scholar
(2016) Reliability of immune receptor rearrangements as genetic markers for minimal residual disease monitoring
Bone Marrow Transplantation 51:1408–1410.
https://doi.org/10.1038/bmt.2016.148
- PubMed
- Google Scholar
1. Oakes T
2. Heather JM
3. Best K
4. Byng-Maddick R
5. Husovsky C
6. Ismail M
7. Joshi K
8. Maxwell G
9. Noursadeghi M
10. Riddell N
11. Ruehl T
12. Turner CT
13. Uddin I
14. Chain B
(2017) Quantitative characterization of the T cell receptor repertoire of naïve and memory subsets using an integrated experimental and computational pipeline which is robust, economical, and versatile
Frontiers in Immunology 8:1267.
https://doi.org/10.3389/fimmu.2017.01267
- PubMed
- Google Scholar
1. Oliveira G
2. Stromhaug K
3. Klaeger S
4. Kula T
5. Frederick DT
6. Le PM
7. Forman J
8. Huang T
9. Li S
10. Zhang W
11. Xu Q
12. Cieri N
13. Clauser KR
14. Shukla SA
15. Neuberg D
16. Justesen S
17. MacBeath G
18. Carr SA
19. Fritsch EF
20. Hacohen N
21. Sade-Feldman M
22. Livak KJ
23. Boland GM
24. Ott PA
25. Keskin DB
26. Wu CJ
(2021) Phenotype, specificity and avidity of antitumour CD8+ T cells in melanoma
Nature 596:119–125.
https://doi.org/10.1038/s41586-021-03704-y
- PubMed
- Google Scholar
Software
1. OpenGene - Open Source Genomics Toolbox
(2022) Fastp, version 1bc2c53
GitHub.

https://github.com/OpenGene/fastp
1. Pogorelyy MV
2. Minervina AA
3. Touzel MP
4. Sycheva AL
5. Komech EA
6. Kovalenko EI
7. Karganova GG
8. Egorov ES
9. Komkov AY
10. Chudakov DM
11. Mamedov IZ
12. Mora T
13. Walczak AM
14. Lebedev YB
(2018) Precise tracking of vaccine-responding T cell clones reveals convergent and personalized response in identical twins
PNAS 115:12704–12709.
https://doi.org/10.1073/pnas.1809642115
- PubMed
- Google Scholar
1. Robins HS
2. Campregher PV
3. Srivastava SK
4. Wacher A
5. Turtle CJ
6. Kahsai O
7. Riddell SR
8. Warren EH
9. Carlson CS
(2009) Comprehensive assessment of T-cell receptor beta-chain diversity in alphabeta T cells
Blood 114:4099–4107.
https://doi.org/10.1182/blood-2009-04-217604
- PubMed
- Google Scholar
(2015) VDJtools: unifying post-analysis of T cell receptor repertoires
PLOS Computational Biology 11:e1004503.
https://doi.org/10.1371/journal.pcbi.1004503
- PubMed
- Google Scholar
Software
1. Shugay M
(2022) Vdjtools, version 1655b10
GitHub.

https://github.com/mikessh/vdjtools
(2022) Inactivated tick-borne encephalitis vaccine elicits several overlapping waves of T cell response
Frontiers in Immunology 13:970285.
https://doi.org/10.3389/fimmu.2022.970285
- PubMed
- Google Scholar
1. Tirtakusuma R
2. Szoltysek K
3. Milne P
4. Grinev VV
5. Ptasinska A
6. Chin PS
7. Meyer C
8. Nakjang S
9. Hehir-Kwa JY
10. Williamson D
11. Cauchy P
12. Keane P
13. Assi SA
14. Ashtiani M
15. Kellaway SG
16. Imperato MR
17. Vogiatzi F
18. Schweighart EK
19. Lin S
20. Wunderlich M
21. Stutterheim J
22. Komkov A
23. Zerkalenkova E
24. Evans P
25. McNeill H
26. Elder A
27. Martinez-Soria N
28. Fordham SE
29. Shi Y
30. Russell LJ
31. Pal D
32. Smith A
33. Kingsbury Z
34. Becq J
35. Eckert C
36. Haas OA
37. Carey P
38. Bailey S
39. Skinner R
40. Miakova N
41. Collin M
42. Bigley V
43. Haniffa M
44. Marschalek R
45. Harrison CJ
46. Cargo CA
47. Schewe D
48. Olshanskaya Y
49. Thirman MJ
50. Cockerill PN
51. Mulloy JC
52. Blair HJ
53. Vormoor J
54. Allan JM
55. Bonifer C
56. Heidenreich O
57. Bomken S
(2022) Epigenetic regulator genes direct lineage switching in MLL/AF4 leukemia
Blood 140:1875–1890.
https://doi.org/10.1182/blood.2021015036
- PubMed
- Google Scholar
1. Wang C
2. Sanders CM
3. Yang Q
4. Schroeder HW
5. Wang E
6. Babrzadeh F
7. Gharizadeh B
8. Myers RM
9. Hudson JR
10. Davis RW
11. Han J
(2010) High throughput sequencing reveals a complex pattern of dynamic interrelationships among human T cell subsets
PNAS 107:1518–1523.
https://doi.org/10.1073/pnas.0913939107
- PubMed
- Google Scholar
1. Wang B
2. DeKosky BJ
3. Timm MR
4. Lee J
5. Normandin E
6. Misasi J
7. Kong R
8. McDaniel JR
9. Delidakis G
10. Leigh KE
11. Niezold T
12. Choi CW
13. Viox EG
14. Fahad A
15. Cagigi A
16. Ploquin A
17. Leung K
18. Yang ES
19. Kong W-P
20. Voss WN
21. Schmidt AG
22. Moody MA
23. Ambrozak DR
24. Henry AR
25. Laboune F
26. Ledgerwood JE
27. Graham BS
28. Connors M
29. Douek DC
30. Sullivan NJ
31. Ellington AD
32. Mascola JR
33. Georgiou G
(2018) Functional interrogation and mining of natively paired human VH: VL antibody repertoires
Nature Biotechnology 36:152–155.
https://doi.org/10.1038/nbt.4052
- Google Scholar
(2015) Immune repertoire profiling reveals that clonally expanded B and T cells infiltrating diseased human kidneys can also be tracked in blood
PLOS ONE 10:e0143125.
https://doi.org/10.1371/journal.pone.0143125
- PubMed
- Google Scholar
1. Wilson TL
2. Kim H
3. Chou C-H
4. Langfitt D
5. Mettelman RC
6. Minervina AA
7. Allen EK
8. Métais J-Y
9. Pogorelyy MV
10. Riberdy JM
11. Velasquez MP
12. Kottapalli P
13. Trivedi S
14. Olsen SR
15. Lockey T
16. Willis C
17. Meagher MM
18. Triplett BM
19. Talleur AC
20. Gottschalk S
21. Crawford JC
22. Thomas PG
(2022) Common trajectories of highly effective CD19-specific CAR T cells identified by endogenous T-cell receptor lineages
Cancer Discovery 12:2098–2119.
https://doi.org/10.1158/2159-8290.CD-21-1508
- PubMed
- Google Scholar
1. Wood B
2. Wu D
3. Crossley B
4. Dai Y
5. Williamson D
6. Gawad C
7. Borowitz MJ
8. Devidas M
9. Maloney KW
10. Larsen E
11. Winick N
12. Raetz E
13. Carroll WL
14. Hunger SP
15. Loh ML
16. Robins H
17. Kirsch I
(2018) Measurable residual disease detection by high-throughput sequencing improves risk stratification for pediatric B-ALL
Blood 131:1350–1359.
https://doi.org/10.1182/blood-2017-09-806521
- PubMed
- Google Scholar
1. Wu J
2. Wang X
3. Lin L
4. Li X
5. Liu S
6. Zhang W
7. Luo L
8. Wan Z
9. Fang M
10. Zhao Y
11. Wang X
12. Mai H
13. Yuan X
14. Wen F
15. Li C
16. Liu X
(2020) Developing an unbiased multiplex PCR system to enrich the TRB repertoire toward accurate detection in leukemia
Frontiers in Immunology 11:1631.
https://doi.org/10.3389/fimmu.2020.01631
- PubMed
- Google Scholar

Decision letter

Gabrielle T Belz

Reviewing Editor; University of Queensland, Australia
Tadatsugu Taniguchi

Senior Editor; University of Tokyo, Japan
Lindsay Cowell

Reviewer

Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "The use of non-functional clonotypes as a natural spike-in for multiplex PCR bias correction in immune receptor repertoire profiling" for consideration by eLife. Your article has been reviewed by 2 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Tadatsugu Taniguchi as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Lindsay Cowell (Reviewer #2).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

The reviewers are enthusiastic about this manuscript. They have identified a number of issues that require clarification and are generally thought to be important to reach the broad immunology audience.

Reviewer #1 (Recommendations for the authors):

The conclusions are mostly well supported. However, there are some concerns that should be addressed before this paper is accepted for publication.

The authors have kindly provided the iROAR software free for non-profit use on Github. As this entire study was performed to validate and share this software to other scientists, I would suggest that the iROAR documentation be improved such that it can be used by a wide audience. More detailed instructions for use and perhaps even a step by step example on how to replicate the results in this study included in the documentation would be very helpful.

The authors utilised pre-existing datasets that performed both multiplex PCR (high bias) and 5'-RACE (low bias) on the same sample to test the functionality of iROAR. The iROAR algorithm improved the PCR bias in repertoires that were obtained by multiplex PCR which better correlated with the same repertoire determined by 5'-RACE. It is important to note that although the iROAR algorithm improved PCR bias, it is by no means perfect and cannot substitute for a low bias approach. Even for a 5'-RACE repertoire in which artificial bias was introduced in silico, iROAR improved the correlation of the in silico biased 5'-RACE from R2 = 0.4052 to R2 = 0.6286 which is an improvement but certainly cannot substitute for a repertoire that was determined using a low bias approach.

The authors very nicely showed the extent of amplification bias when using multiplex PCR based technologies by plotting the OAR distributions in Figure 2A and Figure 2B. It would be very informative to plot the iROAR corrected OAR and OAR distributions as well in order to visualise the improvement in bias correction in the same way..

The authors state that an F-test was performed to compare OAR distributions in Figure 2A. Please state exactly which F-test was used to calculate statistical significance.

The authors state that they use a z-test to exclude outliers from the OAR calculation that may have resulted from large clonotypes introduced by PCR bias or naturally generated from clonal expansions. The authors should describe further how this affects analysis of repertoires in which large clonal expansions are expected such as during an antiviral immune response.

The authors nicely compared clonal frequency between iROAR corrected and 5'-RACE repertoires in Figure 4B-D. I think it would also be important to know how the iROAR algorithm would affect other measurements commonly used when analysing TCR/BCR repertoires. For example, does iROAR also affect measurements of diversity? This would be especially useful when comparing corrected repertoires vs. 5'-RACE.

As per eLife policy, "Regardless of whether authors use original data or are reusing data available from public repositories, they must provide program code, scripts for statistical packages, and other documentation sufficient to allow an informed researcher to precisely reproduce all published results." I would recommend that the authors very much improve the documentation provided in the Github repository for iROAR such that anyone can reproduce the published results before publication. Additionally, eLife suggests that the authors should license their code using an open source license.

Line 158: What do you mean by z-test? Can you explain under which circumstances you can decide to exclude abnormally large clonotypes? How has this threshold been calculated? Would this negatively affect measurements of OAR during an immune response such as viral infections where large clonotypes might be observed?

Line 167, 2.5 reads is a minimal sufficient sequencing coverage. Does the package provide this quality control check? This would be useful information that could be included in the manuscript as well as within iROAR documentation.

Line 203 Figure 4: Can the authors also generate the OAR distributions (a per Figure 2B) for these datasets before and after iROAR? Would be very useful as a comparison between 5' RACE and other multiplex PCR methods as well as to assess the efficacy of iROAR in correcting bias.

Line 219, How does the in silico introduced bias look with regards to the OAR? Can you generate the OAR distributions of in silico biased and unbiased samples using the same graphing method in Figure2B?

The Discussion section needs much more discussion on its applications as well as limitations. Ie., how does the software deal with instances where clonal expansions occur in cell clones that are also bearing a non-functional receptor? How does this affect the OAR calculation?

Reviewer #2 (Recommendations for the authors):

Enthusiasm is dampened by the fact that the proposed method is not directly compared to the gold standard of biological spike-ins. The results would be stronger if the authors could use data generated with a biological spike-in and compare correction using data from the spike-in versus their proposed algorithm. If such a data set is publicly available or can be obtained by request, this should be done. If such a data set cannot be obtained, then the Discussion section should directly address the fact that this gold standard validation remains to be done.

Figure 2:

– More interpretation of the results is needed in the Discussion. For example, Figure 2a appears to show that for VMPlex and VJMPlex, most genes suffer from under- rather than over-amplification. This seems unexpected. I would have expected "zero-sum-game" behavior. Also, it is surprising that VJMPlex shows less overall bias than VMPlex. Does this make sense? Why or why not? I think a plot like Figure 1 is needed for VMPlex and VJMPlex to ensure that the base assumption underlying OAR calculation holds in the same way and to the same degree across the three technologies.

– Regarding the input data: "for each method type 15 PBMC TCR repertoires were chosen randomly from" and this is followed by a list of 11 study/project identifiers. More information is needed. Fifteen seems like a small number, but it depends on the details. How do the 11 studies distribute across the three methods? How many studies were included in each method? Were 15 repertoires total used for each method? Or 15 from each study? Did the different studies within a method use the same primer set? The same depth coverage?

– How much within-study between-repertoire variability is there for a single V or J gene? And how much variability is there across studies? It seems important to understand this basic behavior of the metric before pooling into a single figure, as the pooling can average out important behavior.

Figure 3: More interpretation of the results is needed.

– In Figure 3a, the range of values is greater for TRBV and IGHV. Could this be attributed to a larger number of genes or primers, or to more sequence variability (and presumably therefore greater variability in primer hybridization efficiency)?

– In Figure 3b, while there may be no difference in the mean or median OAR value across cell types, there does appear to be a difference in the range, with PBMCs and THs showing much less variability than the other cell types. What could explain this? Is this an artifact of the way samples were pooled for the figure? Are both donors represented in all six cell types? Or does this point to some interesting biology? If the latter, then exploration of this is understandably beyond the scope of this paper, but it should be mentioned and possible explanations should be put forward.

Lines 197 – "the procedure can be recursively repeated with a modified normalization coefficient defined as described coefficient raised to the power of a number in the range from 0 to 1.": How is the value of this number determined? This sounds like an optimization procedure, in which case much more information is needed. In particular: what algorithm was used? what is the objective function? What are the stopping criteria? In which of the presented results was this procedure applied?

Line 87 "constant during a lifetime": this needs to be substantiated with either a reference or data. The authors show relevant data in Supplementary Figure 1, but there is no description of how the data shown in the figure were derived. For example: are all data points from a single study? do all data points used in a single regression (i.e., all those corresponding to the same gene segment) derive from the same patient? do all data in the figure derive from one single patient? what sequencing protocol was used? Has the proposed correction been applied? The strength of the claim should reflect the strength of the data. Note: I don't think the truth of this claim is a requirement for the proposed method to be valid. So if the data are such that the statement cannot be generalized to "all" repertoires, I don't see this as a problem.

Line 89 "reproducible": reproducible in what way? Over multiple aliquots of a sample? Over multiple samples for the same person? Over multiple people?

Figure 1: Is the same stability observed for the sequencing protocols with even more amplification bias (VMPlex and VJMPlex)?

Lines 102 – 109: this paragraph is not clearly written, but it is the heart of the manuscript. The paragraph reads as if only the two terms with summations are based on out-of-frame rearrangements. I assume all four terms are based on out-of-frame rearrangements? Because otherwise the equation does not make sense …. Also, I recommend removing "a percentage of" in line 105, as this reads as if there is a percentage in the numerator of the numerator (i.e., the RC(Vi) term). Finally, in line 107, I assume you mean PCR amplification instead of clonal expansion?

Figure 2: More information about the data is needed. Specific details are below.

Lines 149 – 153 starting with "the average population frequencies": It isn't clear what is meant by "average population frequencies", and it isn't clear what the rest of the paragraph implies about calculations described in the paper and results displayed in the figures. It sounds as if different repertoires may have been subjected to different calculations.

Lines 157 – 158: The last sentence of this paragraph suggests that V- or J-genes with large or small relative frequencies may have been excluded. More details are needed. Were only large clones excluded? How were they identified for exclusion? By what method and what exclusion threshold was used? How many were excluded from each repertoire? If more than "a few", then before and after data need to be shown for some repertoires to show how this impacts the OAR distributions and averages.

Lines 167 – 169: how is "adequate" defined? Within the 10% error range discussed in the context of Supplemental Figure 3?

Lines 257 – 262: These sentences suggest that the proposed approach is more reliable than biological spike-ins. This should be substantiated before being claimed. The authors state that this is because "the impact of CDR3 structure" is minimized. Why is this? Because the number of out-of-frame rearrangements is much much larger than the number of cell-line or synthetic spike-ins and so captures more CDR3 diversity? Or some other reason? More explanation of why this is being suggested should be given and the wording should make clear that this is an untested hypothesis.

https://doi.org/10.7554/eLife.69157.sa1

Author response

Reviewer #1 (Recommendations for the authors):

The conclusions are mostly well supported. However, there are some concerns that should be addressed before this paper is accepted for publication.

The authors have kindly provided the iROAR software free for non-profit use on Github. As this entire study was performed to validate and share this software to other scientists, I would suggest that the iROAR documentation be improved such that it can be used by a wide audience. More detailed instructions for use and perhaps even a step by step example on how to replicate the results in this study included in the documentation would be very helpful.

We edited the documentation at GitHub and provided step-by-step instructions for the general usage in the Methods section “Step-by-step pipeline for the OAR evaluation used in this study”:

Single nucleotide error correction in read1/read2 intersected sequences and Illumina adapters removal (optional):

fastp -c -i input_R1.fastq.gz -I input_R2.fastq.gz -o fp.input_R1.fastq.gz -O fp.input_R2.fastq.gz

Raw reads alignment (essential):

For TCR β chains

mixcr align -c TRB fp.input_R1.fastq.gz fp.input_R2.fastq.gz output1.vdjca

For TCR β chains

mixcr align -c TRA fp.input_R1.fastq.gz fp.input_R2.fastq.gz output1.vdjca

Clonotypes assemble (essential):

mixcr assemble output1.vdjca output2.clns

TCR repertoire export in a human-readable format (essential):

mixcr exportClones output2.clns clones.txt

Convert repertoire into VDJtools format (essential):

java -jar vdjtools.jar Convert -S mixcr clones.txt vdjtools

Artificial diversity removal by single nucleotide substitutions correction (optional):

java -jar vdjtools.jar Correct vdjtools.clones.txt correct

Pseudogenes removal (optional):

For TCR β chains

java -jar vdjtools.jar FilterBySegment –-j-segment TRBJ2-2P --v-segment TRBV1,TRBV12-1,TRBV12-2,TRBV17,TRBV21-1,TRBV22-1,TRBV23-1,TRBV26,TRBV5-2,TRBV5-3,TRBV5-7,TRBV6-7,TRBV7-1,TRBV7-5,TRBV8-1,TRBV8-2 –-negative correct.vdjtools.clones.txt filter

For TCR α chains

java -jar vdjtools.jar FilterBySegment -–j-segment TRAJ1,TRAJ19,TRAJ2,TRAJ25,TRAJ51,TRAJ55,TRAJ58,TRAJ59,TRAJ60,TRAJ61 --v-segment TRAV11,TRAV11-1,TRAV14-1,TRAV15,TRAV28,TRAV31,TRAV32,TRAV33,TRAV37,TRAV46,TRAV8-5,TRAV8-6-1,TRAV8-7 –-negative correct.vdjtools.clones.txt filter

Artificial diversity removal by single nucleotide indels correction (optional):

iroar Filter -se 0.01 filter.correct.vdjtools.clones.txt filter2.txt

OARs calculation and quantitative bias correction (essential):

iroar Count -min_outframe 15 -z 1 -mt 1 -iter 1 input_folder output_folder

(input_folder must contain filter2.txt file)

The authors utilised pre-existing datasets that performed both multiplex PCR (high bias) and 5'-RACE (low bias) on the same sample to test the functionality of iROAR. The iROAR algorithm improved the PCR bias in repertoires that were obtained by multiplex PCR which better correlated with the same repertoire determined by 5'-RACE. It is important to note that although the iROAR algorithm improved PCR bias, it is by no means perfect and cannot substitute for a low bias approach. Even for a 5'-RACE repertoire in which artificial bias was introduced in silico, iROAR improved the correlation of the in silico biased 5'-RACE from R2 = 0.4052 to R2 = 0.6286 which is an improvement but certainly cannot substitute for a repertoire that was determined using a low bias approach.

We added a suggested sentence: “Even though it’s not possible to fully substitute the low-biased RACE methods, iROAR is capable to make multiplex-PCR-based repertoires more consistent with RACE based ones.”

The authors very nicely showed the extent of amplification bias when using multiplex PCR based technologies by plotting the OAR distributions in Figure 2A and Figure 2B. It would be very informative to plot the iROAR corrected OAR and OAR distributions as well in order to visualise the improvement in bias correction in the same way.

Done. Figures 5 and 8.

The authors state that an F-test was performed to compare OAR distributions in Figure 2A. Please state exactly which F-test was used to calculate statistical significance.

We used nonparametric Levene's test. We specified it in the Figure 2 legend.

The authors state that they use a z-test to exclude outliers from the OAR calculation that may have resulted from large clonotypes introduced by PCR bias or naturally generated from clonal expansions. The authors should describe further how this affects analysis of repertoires in which large clonal expansions are expected such as during an antiviral immune response.

It’s just a simple mathematical issue due to the formula which we used for OAR calculation. The dramatic error in OAR calculation can occur only in the case of enormous lymphoproliferation far beyond normal antiviral immune response.

We corrected the related statement in the manuscript:

Also, the balance of V- and J-genes frequencies can be disrupted by accidentally arisen abnormally large non-functional clonotypes either artificially amplified due to the extreme PCR bias or naturally generated in the course of abnormal clonal expansion in various lymphoproliferative disorders or stochastic spike in normal lymphocyte population. To reduce the impact of this anomaly on OAR value, the top clone of each V and J-gene containing subgroups must be excluded from OAR calculation. Since V- and J-specific bias affects all clones non-selectively, the remaining large part of clones after top clones exclusion is still representative for PCR bias calculation. As shown in Figure 4A, the exclusion of one top clonotype from OAR calculation for RACE-based TRB repertoire is enough to restore OAR calculation accuracy for TRBV2, TRBV5-6, TRBV7-9, TRBV11-3. The further top clones exclusion has no significant effect on OAR values.

The authors nicely compared clonal frequency between iROAR corrected and 5'-RACE repertoires in Figure 4B-D. I think it would also be important to know how the iROAR algorithm would affect other measurements commonly used when analysing TCR/BCR repertoires. For example, does iROAR also affect measurements of diversity? This would be especially useful when comparing corrected repertoires vs. 5'-RACE.

By default, iROAR does not affect diversity since it does not remove any clones. This statement was added to the manuscript.

As per eLife policy, "Regardless of whether authors use original data or are reusing data available from public repositories, they must provide program code, scripts for statistical packages, and other documentation sufficient to allow an informed researcher to precisely reproduce all published results." I would recommend that the authors very much improve the documentation provided in the Github repository for iROAR such that anyone can reproduce the published results before publication. Additionally, eLife suggests that the authors should license their code using an open source license.

We provided more detailed documentation of iROAR on GitHub and added step-by-step instructions in the Methods section. iROAR is totally free to the scientific community for non-profit usage.

Line 158: What do you mean by z-test? Can you explain under which circumstances you can decide to exclude abnormally large clonotypes? How has this threshold been calculated? Would this negatively affect measurements of OAR during an immune response such as viral infections where large clonotypes might be observed?

In the current version of iROAR we deleted z-test to improve calculation speed. However, we kept the option to discount the user-specified number of top non-functional clones (for each V and J gene) from the OAR calculation. It could significantly improve the accuracy of OAR calculation in case of strong clonal expansion, especially in lymphoproliferative disorders. To illustrate the impact of one top clone exclusion on OAR calculation improvement we added Figure 4a and the related text.

Line 167, 2.5 reads is a minimal sufficient sequencing coverage. Does the package provide this quality control check? This would be useful information that could be included in the manuscript as well as within iROAR documentation.

Done. This information is added to the iROAR report file.

Line 203 Figure 4: Can the authors also generate the OAR distributions (a per Figure 2B) for these datasets before and after iROAR? Would be very useful as a comparison between 5' RACE and other multiplex PCR methods as well as to assess the efficacy of iROAR in correcting bias.

Done. Figures 5 and 8.

Line 219, How does the in silico introduced bias look with regards to the OAR? Can you generate the OAR distributions of in silico biased and unbiased samples using the same graphing method in Figure2B?

Done. Figure 5.

The Discussion section needs much more discussion on its applications as well as limitations. Ie., how does the software deal with instances where clonal expansions occur in cell clones that are also bearing a non-functional receptor? How does this affect the OAR calculation?

To deal with this situation we recommend excluding the strongly expanded (top) non-functional clones from OAR calculation, as we mention above in responses as well as in the manuscript.

Reviewer #2 (Recommendations for the authors):

Enthusiasm is dampened by the fact that the proposed method is not directly compared to the gold standard of biological spike-ins. The results would be stronger if the authors could use data generated with a biological spike-in and compare correction using data from the spike-in versus their proposed algorithm. If such a data set is publicly available or can be obtained by request, this should be done. If such a data set cannot be obtained, then the Discussion section should directly address the fact that this gold standard validation remains to be done.

During manuscript revision, we designed and performed a wet-lab experiment to directly compare the iROAR approach with biological spike-ins.

Figure 2:

– More interpretation of the results is needed in the Discussion. For example, Figure 2a appears to show that for VMPlex and VJMPlex, most genes suffer from under- rather than over-amplification. This seems unexpected. I would have expected "zero-sum-game" behavior. Also, it is surprising that VJMPlex shows less overall bias than VMPlex. Does this make sense? Why or why not? I think a plot like Figure 1 is needed for VMPlex and VJMPlex to ensure that the base assumption underlying OAR calculation holds in the same way and to the same degree across the three technologies.

Thanks a lot for this observation. OAR coefficient is actually a relative measure. In the initial manuscript version, we used unnormalized OARs which can lead to misinterpretation of OAR values compared to other PCR bias measures. We corrected this inconvenience by using OAR normalized by average OAR value. This change can’t affect the bias removal process and at the same time will help to present bias coefficients in a more convenient way.

– Regarding the input data: "for each method type 15 PBMC TCR repertoires were chosen randomly from" and this is followed by a list of 11 study/project identifiers. More information is needed. Fifteen seems like a small number, but it depends on the details. How do the 11 studies distribute across the three methods? How many studies were included in each method? Were 15 repertoires total used for each method? Or 15 from each study? Did the different studies within a method use the same primer set? The same depth coverage?

Despite the number of published studies which used multiplex PCR-based methods for TCR profiling the number of publicly available raw read datasets is dramatically limited especially deep-sequenced ones. For this reason, the number of repertoires used for the analysis presented in Figure 2 was limited by available multiplex-based datasets (15 repertoires total for each method, i.e. 45 repertoires in sum, one repertoire for one individual). However, to address the concerns pointed out below we had to reduce this number to 6 repertoires per method. We corrected the figure legend according to the recommendation.

To equalize the sequencing coverage, we performed the down sampling as described in the Methods section.

– How much within-study between-repertoire variability is there for a single V or J gene? And how much variability is there across studies? It seems important to understand this basic behavior of the metric before pooling into a single figure, as the pooling can average out important behavior.

We corrected Figure 2 in order to show within-study between-repertoire variabilities.

Figure 3: More interpretation of the results is needed.

– In Figure 3a, the range of values is greater for TRBV and IGHV. Could this be attributed to a larger number of genes or primers, or to more sequence variability (and presumably therefore greater variability in primer hybridization efficiency)?

Thanks for the questions. We truly appreciate such attention to even small details of our research. The variability of IGHV’s OARs is smaller compared to multiplex PCR but a little bit greater than TCR’s OARs. Probably, it’s a consequence of differences between clonal dynamics of particular subsets of Т-cell and B-cell. We agreed, that it’s interesting phenomenon, it’s indeed needed to be investigated separately. But we can only speculate about it for now.

Added to manuscript:

“Herewith, the variance of IGHV’s OARs compared TCRs’ and the variance of TCR subpopulations' OARs compared PBMCs’ is slightly higher. This phenomenon may be linked to well-known differences in clonal expansion intensities of B/T-cell subsets. But the proof of this hypothesis demands separate deep analysis which is beyond the main focus of this research.”

– In Figure 3b, while there may be no difference in the mean or median OAR value across cell types, there does appear to be a difference in the range, with PBMCs and THs showing much less variability than the other cell types. What could explain this? Is this an artifact of the way samples were pooled for the figure? Are both donors represented in all six cell types? Or does this point to some interesting biology? If the latter, then exploration of this is understandably beyond the scope of this paper, but it should be mentioned and possible explanations should be put forward.

Added to manuscript:

Lines 197 – "the procedure can be recursively repeated with a modified normalization coefficient defined as described coefficient raised to the power of a number in the range from 0 to 1.": How is the value of this number determined? This sounds like an optimization procedure, in which case much more information is needed. In particular: what algorithm was used? what is the objective function? What are the stopping criteria? In which of the presented results was this procedure applied?

The most quantitative bias can be corrected in the first iteration. (We added the corresponding sentence to the main text). All bias correction procedures in the manuscript were performed using a single iteration. But in terms of the flexibility of iROAR, we added an opportunity for the potential user to optimize the process of bias correction for their own multiplex system.

Line 87 "constant during a lifetime": this needs to be substantiated with either a reference or data. The authors show relevant data in Supplementary Figure 1, but there is no description of how the data shown in the figure were derived. For example: are all data points from a single study? do all data points used in a single regression (i.e., all those corresponding to the same gene segment) derive from the same patient? do all data in the figure derive from one single patient? what sequencing protocol was used? Has the proposed correction been applied? The strength of the claim should reflect the strength of the data. Note: I don't think the truth of this claim is a requirement for the proposed method to be valid. So if the data are such that the statement cannot be generalized to "all" repertoires, I don't see this as a problem.

We agreed with Reviewer 2, that this part is quite irrelevant to this study and has to be analyzed in more detail, so we deleted the statement "constant during a lifetime" along with Supplementary Figure 1.

Line 89 "reproducible": reproducible in what way? Over multiple aliquots of a sample? Over multiple samples for the same person? Over multiple people?

Reproducible in all these ways for each particular multiplex PCR primer system. The phrase “for the same multiplex PCR primer set” was added to the sentence.

Figure 1: Is the same stability observed for the sequencing protocols with even more amplification bias (VMPlex and VJMPlex)?

The point of the figure is that we show this stability on low-biased repertoires as a benchmark.

Lines 102 – 109: this paragraph is not clearly written, but it is the heart of the manuscript. The paragraph reads as if only the two terms with summations are based on out-of-frame rearrangements. I assume all four terms are based on out-of-frame rearrangements? Because otherwise the equation does not make sense …. Also, I recommend removing "a percentage of" in line 105, as this reads as if there is a percentage in the numerator of the numerator (i.e., the RC(Vi) term). Finally, in line 107, I assume you mean PCR amplification instead of clonal expansion?

Corrected

Figure 2: More information about the data is needed. Specific details are below.

Lines 149 – 153 starting with "the average population frequencies": It isn't clear what is meant by "average population frequencies", and it isn't clear what the rest of the paragraph implies about calculations described in the paper and results displayed in the figures. It sounds as if different repertoires may have been subjected to different calculations.

“Average population frequencies” was replaced by “population frequencies”. The usage of population frequencies is optional. The users may need it when they face repertoires obtained using multiplex PCR with extremely broad primer efficacies. Analysis in the main text was performed without this option. Since eLife rules do not allow Supplementary notes, we decided to not discard this section from the paper to show to potential users the way of overcoming this possible pitfall.

Lines 157 – 158: The last sentence of this paragraph suggests that V- or J-genes with large or small relative frequencies may have been excluded. More details are needed. Were only large clones excluded? How were they identified for exclusion? By what method and what exclusion threshold was used? How many were excluded from each repertoire? If more than "a few", then before and after data need to be shown for some repertoires to show how this impacts the OAR distributions and averages.

In the current iROAR version user can define manually how many top non-functional clonotypes has to be excluded from OAR calculation. For most normal samples exclusion of one top clonotype for each V and each J-segment is sufficient to prevent inaccuracy in OAR calculation. To illustrate it we added suggested plots to Figure 4.

Lines 167 – 169: how is "adequate" defined? Within the 10% error range discussed in the context of Supplemental Figure 3?

Corrected:

“…adequate OAR values with the acceptable error rate of ~10%”.

Lines 257 – 262: These sentences suggest that the proposed approach is more reliable than biological spike-ins. This should be substantiated before being claimed. The authors state that this is because "the impact of CDR3 structure" is minimized. Why is this? Because the number of out-of-frame rearrangements is much much larger than the number of cell-line or synthetic spike-ins and so captures more CDR3 diversity? Or some other reason? More explanation of why this is being suggested should be given and the wording should make clear that this is an untested hypothesis.

We corrected the discussion sentence: “It makes this method potentially more reliable due to the ability to minimize the impact of CDR3 structure on PCR bias calculation since out-of-frame captures significantly higher CDR3 diversity”.

https://doi.org/10.7554/eLife.69157.sa2

Article and author information

Author details

Anastasia O Smirnova
1. Skolkovo Institute of Science and Technology, Moscow, Russian Federation
2. Department of Genomics of Adaptive Immunity, Shemyakin-Ovchinnikov Institute of Bioorganic Chemistry, Moscow, Russian Federation
Contribution
Data curation, Software, Formal analysis, Investigation

Competing interests
No competing interests declared
Anna M Miroshnichenkova
1. Department of Genomics of Adaptive Immunity, Shemyakin-Ovchinnikov Institute of Bioorganic Chemistry, Moscow, Russian Federation
2. Dmitry Rogachev National Medical Research Center of Pediatric Hematology, Oncology and Immunology, Moscow, Russian Federation
Contribution
Validation, Writing – review and editing

Competing interests
No competing interests declared
Yulia V Olshanskaya
1. Department of Genomics of Adaptive Immunity, Shemyakin-Ovchinnikov Institute of Bioorganic Chemistry, Moscow, Russian Federation
2. Dmitry Rogachev National Medical Research Center of Pediatric Hematology, Oncology and Immunology, Moscow, Russian Federation
Contribution
Resources, Supervision, Writing – review and editing

Competing interests
No competing interests declared
Michael A Maschan

Dmitry Rogachev National Medical Research Center of Pediatric Hematology, Oncology and Immunology, Moscow, Russian Federation

Contribution
Resources, Supervision, Writing – review and editing

Competing interests
No competing interests declared
Yuri B Lebedev
1. Department of Genomics of Adaptive Immunity, Shemyakin-Ovchinnikov Institute of Bioorganic Chemistry, Moscow, Russian Federation
2. Pirogov Russian National Research Medical University, Moscow, Russian Federation
Contribution
Resources, Supervision, Writing – review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-4554-4733
Dmitriy M Chudakov
1. Skolkovo Institute of Science and Technology, Moscow, Russian Federation
2. Department of Genomics of Adaptive Immunity, Shemyakin-Ovchinnikov Institute of Bioorganic Chemistry, Moscow, Russian Federation
3. Pirogov Russian National Research Medical University, Moscow, Russian Federation
4. Abu Dhabi Stem Cells Center, Abu Dhabi, United Arab Emirates
Contribution
Resources, Supervision, Writing – review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-0430-790X
Ilgar Z Mamedov
1. Department of Genomics of Adaptive Immunity, Shemyakin-Ovchinnikov Institute of Bioorganic Chemistry, Moscow, Russian Federation
2. Pirogov Russian National Research Medical University, Moscow, Russian Federation
Contribution
Resources, Supervision, Funding acquisition, Writing – original draft, Project administration, Writing – review and editing

Competing interests
No competing interests declared
Alexander Komkov
1. Department of Genomics of Adaptive Immunity, Shemyakin-Ovchinnikov Institute of Bioorganic Chemistry, Moscow, Russian Federation
2. Dmitry Rogachev National Medical Research Center of Pediatric Hematology, Oncology and Immunology, Moscow, Russian Federation
Contribution
Conceptualization, Formal analysis, Supervision, Funding acquisition, Investigation, Methodology, Writing – original draft, Project administration, Writing – review and editing

For correspondence
alexandrkomkov@yandex.ru

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-9113-698X

Funding

Russian Science Foundation (20-75-10091)

Alexander Komkov

Russian Foundation for Basic Research (20-015-00462)

Alexander Komkov

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank Grigory Armeev and Valery Novoseletsky for their help with data storage and processing. This work was supported by Russian Science Foundation (grant 20-75-10091 to A.K.), the analysis of TRA MPlex data was supported by Russian Foundation for Basic Research (grant 20-015-00462 to A.K.).

Senior Editor

Tadatsugu Taniguchi, University of Tokyo, Japan

Reviewing Editor

Gabrielle T Belz, University of Queensland, Australia

Reviewer

Lindsay Cowell

Version history

Preprint posted: March 25, 2021 (view preprint)
Received: April 6, 2021
Accepted: January 22, 2023
Accepted Manuscript published: January 24, 2023 (version 1)
Accepted Manuscript updated: January 26, 2023 (version 2)
Version of Record published: February 6, 2023 (version 3)
Version of Record updated: February 9, 2023 (version 4)

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.