powerTCR: A model-based approach to comparative analysis of the clone size distribution of the T cell receptor repertoire
[['Recent advances in high-throughput sequencing of the T cell receptor ( TCR ) repertoire provide a new , detailed characterization of the immune system .', 'T cells , each displaying a unique TCR , are capable of responding to presented antigens and initiating an adaptive immune response .', 'An immune response is described by rapid proliferation of T cell clonotypes whose TCRs are specific to the antigen .', 'In humans , it is estimated that the body is capable of producing more than 1018 different TCRs [1 , 2] , where high diversity of the TCR repertoire implies a greater range of pathogens that can be fought off .', 'A variety of studies have been published demonstrating the value in characterizing this immune response for purposes such as describing tumor cell origin [3] and predicting response to cancer therapy and infection [4] .', 'The applications of TCR sequencing are many , but this type of data presents new needs for analysis techniques not met by existing tools for other kinds of genomic experiments .', 'Several groups have identified that the distribution of larger clone sizes in a sample can be approximated by a power law [5–8] , which means that the number of clones of a given size decays approximately as a power of the clone size .', 'This heavy-tailed distribution comes as a consequence of extensively proliferated clones actively participating in an ongoing immune response .', 'More recent work has aimed to quantify statistically the diversity of the TCR repertoire , initially through the use of various estimators borrowed from ecology , such as species richness , Shannon entropy [9] , and clonality .', 'These estimators are known to be highly sensitive to sample size and missing observations .', 'Given that the TCR repertoire is mostly populated by rare clonotypes , many of the clonotypes in the system are absent from any one sample .', 'This presents a challenge to many of the ecological estimators .', 'Model-based approaches to approximating the clone size distribution have also been proposed , with the goal of providing added stability and consequently more statistical power .', 'Some examples are the Poisson-lognormal model [10] , Poisson mixture models [11 , 12] , and a heuristic ensemble method [13]; however , these models lack a biologically meaningful interpretation , and further do not sufficiently account for the power law-like nature of the data .', 'That is , power law distributions are heavier-tailed than the Poisson or even the lognormal distribution , leading to systematic bias in the model fit .', 'Previous research has also identified the imperfectness of the power law behavior for the clone size distribution below some clone size threshold [7 , 8] .', 'To handle the imperfectness , [7 , 8] proposed to model large clones above the threshold using a type-I Pareto distribution , which is a member of the power-law distribution family , and omitted the clones with frequency below that threshold .', 'The threshold is either user-specified or determined from the data based on a goodness-of-fit measure .', 'Indeed , this model has certain biological basis .', 'Through a stochastic differential equations setup that models the birth , death , selection , and antigen-recognition of cells active in the immune system , Desponds et al . [8] showed that the upper tail of the clone size distribution at equilibrium approximately follows a type-I Pareto distribution ( Fig 1A ) .', 'Unlike the Poisson and lognormal models , parameters in this model are related to relevant actors in the immune response , and can reveal certain biological insights into immune response , such as average T cell lifetime [8] .', 'Yet , the resulting model excludes all clones below a certain frequency threshold .', 'However , even small clones may provide information; for example , Desponds et al . [8] indicated that the generation of new T cells affects the landscape of smaller clone sizes .', 'Other studies have shown that low-frequency clones may support a diverse immune system and present a potential to mobilize against antigens , as in some cases having a clone size distribution highly dominated by a few clones has been correlated with unfavorable clinical outcome [14 , 15] .', 'With this in mind , we sought a means to exhaust all available data and consider modeling the complete clone size distribution .', 'To address this question , we propose a novel statistical tool , called powerTCR , to characterize the full distribution of the TCR repertoire .', 'Our method models large clones that are above the threshold , where the power law begins , using the generalized Pareto distribution ( GPD ) , which contains the type-I Pareto distribution as a special case , but provides a more flexible fit .', 'It also models the small clones below the threshold using a truncated Gamma distribution .', 'It determines the threshold in a data-driven manner simultaneously with the characterization of clone size distribution .', 'Our final model contains parameters that are analogous to those found in the type-I Pareto model of Desponds et al . [8] , relating our model to the biological interpretation of the dynamics of the immune system .', 'Altogether , this allows our model to more accurately describe the shape of the clone size distribution for both large and small clones .', 'Such a model is well suited for providing a global view of the state of the immune repertoire .', 'It can also be employed to perform comparative analysis of healthy and compromised individuals to identify descriptors of strengths and deficiencies in the immune system .'], ['Our goal is to model the clone size distribution of a sample immune repertoire .', 'Fig 1A shows a typical distribution plotted using the repertoire of a Sarcoidosis patient in [16] .', 'If the data are truly Pareto distributed , this plot would appear linear [17 , 18] .', 'However , noting the linear behavior is only true for the far upper tail of the data , this suggests that these data are a departure from the Pareto distribution .', 'This imperfect power law implicates the use of a heavy-tailed distribution above some threshold and a lighter-tailed distribution below that threshold .', 'Here , we model the tail part with a GPD .', 'The GPD , introduced by [19] , is a classical distribution typically used to model the values in the upper tail of a dataset .', 'This formulation results in a distribution with density', 'f ( x ) = 1 σ ( 1 + ξ x - u σ ) - ( 1 / ξ + 1 ) , ( 1 )', 'where u ∈ ( −∞ , ∞ ) is a threshold that typically needs to be prespecified , σ ∈ ( 0 , ∞ ) is a scale parameter , and ξ ∈ ( −∞ , ∞ ) is a shape parameter .', 'The GPD has support x ≥ u when ξ ≥ 0 and u ≤ x ≤ u − σ/ξ when ξ < 0 . We model the bulk part with a Gamma distribution with the upper tail truncated at the threshold .', 'The Gamma distribution has a flexible shape and can fit many different clone size distributions .', 'The threshold and the parameters in the two distributions are estimated from the data simultaneously .', 'This setup , where data above and below an unknown threshold are drawn from the “bulk” and “tail” distributions respectively , falls into a class of models called spliced threshold models .', 'The typical motivation for the model is the belief that the data above and below the threshold are driven by different underlying processes .', 'We refer the interested reader to [20] for a thorough review of the general spliced threshold model , and its applications in fields such as insurance , hydrology , and finance .', 'Denote the proportion of data above the threshold u as ϕ .', 'Let the bulk model distribution function be Hc ( x|θb ) and the tail model distribution function be Gc ( x|θt ) , where subscripts b and t denote the bulk and tail model parameter vectors , respectively .', 'Then the distribution function of the model is given by', 'F c ( x ) = { ( 1 - ϕ ) H c ( x | θ b ) H c ( u | θ b ) for x ≤ u 1 - ϕ + ϕ G c ( x | θ t , u ) for x > u ( 2 )', 'with corresponding density', 'f c ( x ) = { ( 1 - ϕ ) h c ( x | θ b ) H c ( u | θ b ) for x ≤ u ϕ g c ( x | θ t , u ) for x > u .', '( 3 ) Because the clone size distribution is count data that typically exhibit numerous ties in the less frequently observed clonotypes , it is appropriate to treat this as a discrete problem .', 'We modify the model in order to account for any quantized or censored data .', 'Let ψ and Ψ be the density and distribution function of a continuous distribution , and let d be the interval length at which the data are censored .', 'We obtain a quantized analog of ψ by letting', 'Pr ( X = x ) = Ψ ( x + d ) - Ψ ( x ) , x ∈ k · d , k ∈ Z .', 'This results in a discrete model with distribution function', 'F ( x ) = { ( 1 - ϕ ) H ( x | θ b ) H ( u - d | θ b ) for x ≤ u - d 1 - ϕ + ϕ G ( x | θ t , u ) for x ≥ u ( 4 )', 'and corresponding probability mass function', 'f ( x ) = { ( 1 - ϕ ) h ( x | θ b ) H ( u - d | θ b ) for x ≤ u - d ϕ g ( x | θ t , u ) for x ≥ u .', '( 5 )', 'where h ( x|θb ) ∼ discrete Gamma ( α , β ) , g ( x|θt ) ∼ discrete GPD ( u , σ , ξ ) , and d = 1 , which specifies that we model integer data ( see Methods for the functional form of the discrete Gamma distribution and the discrete GPD ) .', 'This discretization step turns out to be important for accurate estimation in our scenario .', 'See S2 Text for a comparison between the performance of the discrete and continuous models in settings resembling true clone size distributions .', 'The relationship between the discrete Gamma-GPD spliced threshold model and the type-I Pareto model in Desponds et al . [8] , hereafter referred to as the Desponds et al . model , allows us to draw connections between some of our parameters and the dynamics of immune response underpinning their approach .', 'First , results from [8] show that the threshold at which the power law begins is indicative of the point over which a clone’s large size can be attributed to active immune response , as opposed to noise in the body that arises from processes such as self-recognition .', 'The threshold fitted from the data provides an objective way to narrow down which clonotypes from a sample repertoire should be interrogated further .', 'This notion is convenient for studying factors such as CDR3 ( complementarity-determining region 3 ) amino acid motifs or specific V , D , and J genes important for combating certain antigens , which are typically determined based on a heuristic abundance cutoff .', 'For example , [21] studies the 1 , 000 most abundant CDR3 amino acid motifs across all sampled peripheral blood mononuclear cell ( PBMC ) libraries , while [16] determines CDR3 amino acid motifs from clones that are present with 10 or more reads in a sampled repertoire .', 'The threshold u estimated with our model , however , introduces a means to select motifs that does not rely on heuristics and automatically scales with sequencing depth .', 'Moreover , the shape parameter ξ of the GPD is inversely related to the shape parameter αd used in the Desponds et al . model ( see Methods ) .', 'As explained by Desponds et al . , a small αd , i . e . a large ξ , implies increased average T cell lifetime and antigenic noise strength .', 'They further show that antigenic noise strength grows as a consequence of a higher initial concentration of antigens and a higher rate at which new antigens are introduced .', 'Interestingly , ξ also positively correlates with the familiar clonality estimator ( 1-Pielou’s evenness [22] ) .', 'Indeed , as ξ increases , the clone size distribution becomes heavier-tailed—that is , more skewed towards dominating clones .', 'This trend is in line with that of the clonality estimator , which favors a more uniform clone size distribution as clonality approaches 0 and a distribution dominated by expanded clones as clonality approaches 1 . To numerically validate this relationship , we simulated the data from our model and computed the clonality ( see Methods ) .', 'We observed a high correlation between clonality and ξ ( Spearman’s ρ ≈ 0 . 9 ) , confirming that ξ reflects the skewness towards dominating clones ( see S3 Text ) .', 'It is worth noting that our model acquires a theoretical gain via the threshold stability property of the GPD [23] .', 'That is , for any generalized Pareto distributed data , the shape parameter ξ remains constant regardless of changes in u .', 'In our context , this means that at decreasing sequencing depths , though the threshold u would decrease due to fewer cells being sampled , the shape parameter ξ in principle would be stable against the variation in sequencing depth .', 'We will demonstrate this gain in stability on a murine tumor dataset .', 'See Methods for our extension of the threshold stability property to the case of the discrete GPD .', 'In the following sections , we inspect four different datasets using our model .', 'We compare our results to results from the Desponds et al . model to demonstrate the practical and theoretical benefits of our approach .', 'We also make comparisons to results from the widely used richness , Shannon entropy , and clonality estimators .', 'See Methods for information on computation of competing methods .', 'The expression of major histocompatability complex II ( MHC-II ) proteins in tumors correlates with boosted anti-tumor immunity .', 'As part of a study of how MHC-II expression impacts tumor progression and functional plasticity of T cells [24] , the CDR3 of TCRβ-chains of tumor infiltrating lymphocytes ( TILs ) were sequenced from breast cancer tumor tissue from six BALB/c mice [25] .', 'Three of the mice were grafted with MHC-II expressing tumor cells and three control animals received parental MHC-II-negative cells .', 'Samples were collected at 21 days after the date of treatment .', 'Table A in S4 Text summarizes the number of unique clonotypes observed and the total number of reads in each sample .', 'Sarcoidosis is an inflammatory disease that typically is accompanied by an accumulation of activated CD4+ T cells in the lungs .', 'A particularly acute form of Sarcoidosis , called Löfgren’s syndrom ( LS ) , occurs with additional , more severe symptoms .', 'A known signature of LS is the bombardment of the lungs with CD4+ T cells , which is expected to significantly alter the entire landscape of the TCR repertoire .', 'We applied our method to TCR repertoire data of LS and non-LS Sarcoidosis patients [29] , originally described in [16] .', 'In this study , bronchoscopy with the bronchoalveolar lavage was performed on a cohort of 9 LS and 4 non-LS individuals and prepared for TCR α− and β–chain sequencing .', 'We compared the TCR distribution between LS and non-LS Sarcoidosis patients using our method and the competing methods .', 'In order to visualize closeness of samples , we generated a distance matrix using JSD between fitted distributions using our method and the Desponds et al . model .', 'The estimated parameters are in Tables E and F in S4 Text respectively .', 'We then applied non-metric multidimensional scaling ( MDS ) to the distance matrix and plotted the first two coordinates .', 'For the ecological estimators , we simply plotted centered and scaled estimates .', 'As shown in Fig 4A , results from our model cluster LS patients into a tight group distinct from non-LS patients , bolstering the claim that LS patients exhibit a signature immune response .', 'On the other hand , competing methods fail to uncover any pattern ( Fig 4B and 4C ) .', 'We applied our method to data collected during a clinical trial of 13 glioblastoma patients receiving autologous tumor lysate-pulsed dendritic cell ( DC ) vaccine therapy [30] , first detailed in [31] .', 'Three intradermal injections were administered to patients at biweekly intervals .', 'TCRβ-chains from PBMC samples were sequenced for the patients prior to vaccinations and two weeks following the final injection .', 'Patients were followed up with and their time to progression ( TTP ) and overall survival ( OS ) were recorded .', 'TTP was defined as the time from the first DC vaccination until MRI-confirmed tumor progression .', 'OS was calculated as the time from the first DC vaccination until the patient’s death from any cause .', 'We investigated whether current tools using TCRs sequenced only from blood samples indicate anything about patients’ survival time and time to progression .', 'We first fit our model to the pre- and post-treatment samples .', 'In both cases , we classified the patients into two groups using the hierarchical clustering based on our model , the Desponds et al . model , and the richness , Shannon entropy , and clonality estimators .', 'No clear grouping with respect to either TTP or OS could be observed from any clustering on the pre-treatment samples , whether by the model-based methods or the selected estimators ( see S7 Text ) .', 'However , among post-treatment samples , our method tends to cluster together patients with better clinical outcome ( Fig 5A ) .', 'This may indicate that the DC therapy alters the landscape of the TCR repertoire into a form that promotes favorable clinical outcome .', 'We do , however , cluster one patient ( ID: 33296 ) with low TTP and OS in the group with overall higher TTP and OS .', 'Interestingly , this misplaced patient had the lowest estimated TIL count and tumor/PBMC overlap of the entire cohort ( S4 Text , Table G ) .', 'Tumor/PBMC overlap was defined as the total number of reads of shared CDR3s normalized by total reads in the tumor and PBMC samples .', 'Similarly , patient 17232 displayed among the best clinical outcome but clustered with lower-performing patients .', 'Patient 17232 had the highest TIL count and level of tumor/PBMC overlap in the whole cohort .', 'This information taken as a whole suggests that , while the clone size distribution found in blood may indicate something about a patient’s response to treatment , it still does not guarantee that T cells will infiltrate the tumor , an important factor for clinical benefit [32] .', 'S8 Text highlights the clone size distributions of these two patients against all others .', 'Notably , inferred thresholds ( minimum u = 4 , maximum u = 6 ) on this dataset are much lower than on other datasets .', 'This is likely because this dataset contains less deeply-sequenced samples than the others , which consequently reduces the threshold .', 'Noting that clones with size at or above the estimated threshold are considered active participators in the immune response , we sought to investigate whether any relationship existed between clinical outcome and the proportion of more highly stimulated cells .', 'We defined the proportion of highly stimulated cells to be the total number of reads at or above the threshold , normalized by the total number of reads in the entire repertoire ( S3 Text , Table I ) .', 'We found correlations between this measure and both TTP ( Spearman’s ρ = 0 . 54 ) and OS ( Spearman’s ρ = 0 . 80 ) .', 'Rank scatterplots for these correlations are in Fig 5B .', 'The positive correlation we uncovered suggests that this statistic could be a useful tool to quantify the antigen-specificity of the sample .', 'Risk factors for type 1 diabetes ( T1D ) are known to be heritable , yet genes alone are not sufficient explanation for drivers of the disease .', 'Studies of monozygotic twins have revealed that , given one twin has T1D , the other will only have it at most half of the time [33] .', 'The CD4+ T cell is viewed as the initiator of T1D as dysregulation of CD4+ antigen-recognition drives the autoimmune disease .', 'Seeking out apparently non-heritable determinants of T1D , [34] conducted a deeper investigation of the CD4+ T cell .', 'Briefly , the authors obtained PBMCs from 14 volunteer healthy donors ( HDs ) and 14 recently diagnosed patients with T1D .', 'The cells were sorted using flow cytometry into distinct T cell subsets ( true naïve; TN , central memory , CM; regulatory , Treg; and stem cell-like memory , Tscm ) and TCRβ-chains were sequenced .', 'The authors conducted a thorough analysis , finding shorter CDR3 sequence lengths and lower overall repertoire diversity among patients with T1D .', 'However , on a per-individual basis , the authors were unable to uncover a relationship between repertoire diversity and disease status .', 'Since the the spliced threshold model provides a new means to probe this complex data , we applied our approach to complement the original analyses .'], ['We have developed a model , the discrete Gamma-GPD spliced threshold model , and demonstrated its utility on several datasets .', 'As shown in our analyses , several biologically relevant descriptive features can be obtained from our model .', 'One is the tail shape parameter ξ , a measure of the weight of the upper tail of the clone size distribution , where a heavier tail of the fitted model implies a more dominated distribution of expanded clones .', 'Another is the proportion of total reads at or above the estimated threshold , a possible measure of intensity of the immune response .', 'The third is the estimated threshold , which is a useful guide to objectively identify CDR3 motifs for downstream analysis .', 'This could involve denoting motifs as only those CDR3s found in TCRs with frequencies at or above the estimated threshold for a given sample , or it could mean studying TCR gene usage among that same group of clonotypes .', 'Though the dynamics driving our model form a compelling argument for this interpretation of the threshold , we acknowledge that further biological validation on more datasets is still needed to confirm this .', 'Similar to other estimators , our model requires that a repertoire be adequately sampled .', 'Without adequate sampling , the differentiating features between TCR repertoires will be masked [8] , and the estimated model parameters will not be reliable .', 'Given the immense diversity of the TCR repertoire , one should in general be cautious about using any method to make inference about a sample TCR repertoire when few cells are sequenced .', 'With sufficient samples , though , the spliced threshold model provides the user a meaningful high-level view of the TCR repertoire .', 'The diversity of the TCR repertoire and its responsiveness to stimuli provide a high-dimensional biomarker for monitoring the immune system and its adaptivity .', 'Robust assessment of the clone size distribution through TCR sequencing is important for understanding this diversity .', 'The discrete Gamma-GPD spliced threshold model is a flexible model that effectively captures the shape of the clone size distribution .', 'It is especially appropriate since the heavy-tailed GPD is a good fit to model the highly expanded clones that dominate many TCR repertoire samples .', 'The method also provides a means to comparatively analyze a collection of TCR repertoire samples while maintaining convenient theoretical properties and interpretations .', 'Compared with existing approaches , our method is more flexible , utilizes the full clone size distribution , is less sensitive to sequencing depth , and identifies the threshold in a data-driven manner .', 'The parameters estimated from our method are biologically relevant and instructive to the dynamics of immune response .', 'Our results on multiple datasets also show that the spliced threshold model is powerful in a range of scenarios for comparing TCR repertoires across samples , revealing potential trends in the landscapes of clone size distributions of affected immune systems .'], ['We use maximum likelihood to estimate the parameters of our model .', 'First , we more explicitly specify the form of our distribution .', 'Letting x ∼ Gamma ( α , β ) , we write the probability mass function of a discrete Gamma distribution as', 'h ( x ) = 1 Γ ( α ) [ γ ( α , β ( x + 1 ) ) - γ ( α , β x ) ]', 'for α > 0 , β > 0 , x ∈ Z , and where γ ( α , βx ) is the lower incomplete gamma function', 'γ ( α , β x ) = ∫ 0 β x t α - 1 e - t d t .', 'If x ∼ GPD ( u , σ , ξ ) , we write the probability mass function of a discrete GPD as', 'g ( x ) = ( 1 + ξ x - u σ ) - 1 / ξ - ( 1 + ξ x + 1 - u σ ) - 1 / ξ', 'for u ∈ ( −∞ , ∞ ) , σ ∈ ( 0 , ∞ ) , and ξ ∈ ( −∞ , ∞ ) .', 'The discrete GPD has support x ≥ u when ξ ≥ 0 and u ≤ x ≤ u − σ/ξ when ξ < 0 , where x ∈ Z . In all analyses presented here , we make no assumptions on the sign of ξ , although empirically we tend to observe ξ > 0 .', 'To proceed , we employ a profile likelihood approach .', 'Let u be the threshold , θb be the bulk parameter vector {α , β} , θt be the tail parameter vector {σ , ξ} , and θ be the parameter vector {θb , θt} .', 'Let also h and H be the density and distribution function of a discrete Gamma distribution , respectively , and let g be the density of a discrete GPD .', 'Then the complete data likelihood is given by', 'L ( { θ , u } ∣ x ) = ∏ i = 1 n [ ( 1 - ϕ ) h ( x i ∣ θ b ) H ( u - 1 ∣ θ b ) 1 ( x i ≤ u - 1 ) + ϕ g ( x i ∣ θ t , u ) 1 ( x i ≥ u ) ]', 'and the profile likelihood of the model at u is denoted as', 'L p ( u ) = max θ L ( θ ∣ x , u ) .', 'A grid search over a suitable range of thresholds u⋆ = ( u1 , … , uk ) may be implemented to maximize the profile likelihood .', 'In this study , we adopted an approach similar to those of [19] and [40] , searching for thresholds at or above the 75% quantile of the sample .', 'The estimated parameters are then', 'u ^ = arg max u ⋆ ∈ u ⋆ L p ( u ⋆ ) , θ ^ = argmax θ L ( θ ∣ u = u ^ ) , and ϕ ^ = n u n ,', 'where n is the total number of clones and nu denotes the number of clones with size greater than or equal to the threshold .', 'The Desponds et al . model was fit as previously described [8] .', 'Briefly , the model has density', 'f ( x ) = α d u α d x α d + 1 ( 6 )', 'and distribution function', 'F ( x ) = 1 - ( u x ) α d ( 7 )', 'where u > 0 is the threshold and αd > 0 is a shape parameter .', 'For each sample TCR repertoire , a grid of potential thresholds u⋆ = ( u1 , … , uk ) was constructed by considering every unique clone size in the repertoire .', 'Then , for each ui , the shape parameter is estimated as', 'α ^ d = n i [ ∑ j = 1 n i ln x j u i ]- 1 ( 8 )', 'where ni is the number of clones with size larger than the threshold ui .', 'Once this value is computed for every threshold in u⋆ , the threshold and corresponding α ^ were chosen to minimize the Kolmogorov-Smirnov statistic .', 'The ecological estimators [9 , 22] were computed as follows .', 'For a sample X , let S ( X ) be the sample richness , defined as the number of unique clonotypes in X , and let pi be the number of cells of clonotype i normalized by the total number of cells in the sample .', 'Then , the Shannon entropy of X is', 'H ( X ) = - ∑ i = 1 S ( X ) p i ln p i ( 9 )', 'and the clonality of X is', 'C ( X ) = 1 - H ( X ) ln S ( X ) .', '( 10 )', 'The Desponds et al . model , which is a type-I Pareto distribution , and the “tail” part of our model , which is a GPD , are closely related .', 'In fact , the GPD contains the type-I Pareto distribution as a special case .', 'We can write the distribution function of y , where y ∼ GPD ( u , σ , ξ ) , as', 'F ( y ) = 1 - ( 1 + ξ y - u σ ) - 1 / ξ .', '( 11 )', 'Now , let x ∼ GPD ( u , u α d , 1 α d ) .', 'Then', 'F ( x ; u , u α d , 1 α d ) = 1 - [ 1 + 1 α d ( x - u u / α d ) ] - α d = 1 - [ 1 + ( x - u u ) ] - α d = 1 - ( u x ) α d', 'which is exactly the distribution function of a type-I Pareto distribution with threshold u and shape αd ( Eq 7 ) .', 'Of course , this exact relationship only holds when σ = u α d .', 'Nevertheless , αd and ξ perform the same function in their respective distributions , adjusting the weight of the tail .', 'This relationship always holds—a larger ξ ( smaller αd ) implies a heavier-tailed distribution , while a smaller ξ ( larger αd ) implies a lighter-tailed distribution .', 'We conjecture that ξ , the shape parameter of the GPD , positively correlates with clonality .', 'We numerically validated this claim using a simulated cohort of 48 clone size distributions .', 'That is , we generated samples of n = 20 , 000 clonotypes , where our 48 parameter settings were derived from every combination of α ∈ {3 , 5 , 10} , ξ ∈ { . 25 , . 5 , . 75 , 1 . 1} , and ϕ ∈ {0 . 1 , 0 . 15 , 0 . 2 , 0 . 25} .', 'We chose β = 0 . 15 , σ = α β , and u = ⌊Qα , β ( 1 − ϕ ) ⌋ in each simulation , where Qα , β is the quantile function of the Gamma distribution with mean α β .', 'To adjust for the effect of sample size on clonality , we downsampled the simulated data so that each sample contained the same number of reads ( 415 , 989 total reads per sample ) .', 'We computed the clonality of each simulated TCR repertoire on these adjusted datasets .', 'The relationship between a pair of TCR repertoires can be elucidated by evaluating the distance between their fitted spliced threshold models .', 'Several methods to compare densities are available .', 'We propose measuring the distance between each pair of distributions using Jensen-Shannon distance ( JSD ) [41] .', 'This metric is a symmetric and smoothed adaptation of the well-known Kullback-Leibler divergence that does not require the distributions under comparison to share the same support .', 'Given discrete distributions P and Q , the JSD between P and Q is', 'J S D ( P , Q ) = 1 2 [ ∑ i ( P i ln P i M i ) + ∑ i ( Q i ln Q i M i ) ] , ( 12 )', 'where M i = 1 2 ( P i + Q i ) .', 'The resulting distances allow analysis and visualization via MDS or hierarchical clustering of the samples .', 'Throughout our study , we use Ward’s method for hierarchical clustering .', 'The threshold stability property of the GPD is well-established [23] .', 'Here , we show that the property also holds for the discrete GPD .', 'Let X ∼ discrete GPD ( u , σ , ξ ) and denote its distribution function as F with Fc as its continuous analog .', 'Then we can write', 'P ( X - u ≤ x + 1 | X ≥ u ) = P ( u ≤ X ≤ x + u + 1 ) P ( X ≥ u ) = F ( x + u + 1 ; u , σ , ξ ) - F ( u ; u , σ , ξ ) 1 - F ( u ; u , σ , ξ ) = F c ( x + u + 2 ; u , σ , ξ ) - F c ( u + 1 ; u , σ , ξ ) 1 - F c ( u + 1 ; u , σ , ξ ) = ( 1 + ξ σ ) - 1 / ξ - ( 1 + ξ x + 2 σ ) - 1 / ξ ( 1 + ξ σ ) - 1 / ξ = 1 - ( 1 + ξ x + 1 σ + ξ ) - 1 / ξ = F c ( x + 1 ; 0 , σ + ξ , ξ ) = F ( x ; 0 , σ + ξ , ξ ) .', 'This states that if X ∼ discrete GPD ( u , σ , ξ ) , then X − u ∼ discrete GPD ( 0 , σ + ξ , ξ ) .', 'Or , for our application , consider a clone size distribution , where clones larger than some threshold u are distributed according to the discrete GPD .', 'At decreasing sequencing depths , this estimated u decreases , implying naturally that the size a clone in the sample must achieve to be considered “expanded” decreases .', 'Still , while u shrinks , the threshold stability property states that ξ remains constant .']]
['Introduction', 'Results', 'Discussion', 'Methods']
['Sequencing of the T cell receptor ( TCR ) repertoire is a powerful tool for deeper study of immune response , but the unique structure of this type of data makes its meaningful quantification challenging .', 'We introduce a new method , the Gamma-GPD spliced threshold model , to address this difficulty .', 'This biologically interpretable model captures the distribution of the TCR repertoire , demonstrates stability across varying sequencing depths , and permits comparative analysis across any number of sampled individuals .', 'We apply our method to several datasets and obtain insights regarding the differentiating features in the T cell receptor repertoire among sampled individuals across conditions .', 'We have implemented our method in the open-source R package powerTCR .']
['A more detailed understanding of the immune response can unlock critical information concerning diagnosis and treatment of disease .', 'Here , in particular , we study T cells through T cell receptor sequencing , as T cells play a vital role in immune response .', 'One important feature of T cell receptor sequencing data is the frequencies of each receptor in a given sample .', 'These frequencies harbor global information about the landscape of the immune response .', 'We introduce a flexible method that extracts this information by modeling the distribution of these frequencies , and show that it can be used to quantify differences in samples from individuals of different biological conditions .']
['blood cells', 'inflammatory diseases', 'rheumatology', 'medicine and health sciences', 'cluster analysis', 'immune cells', 'hierarchical clustering', 'immunology', 'blastomas', 'cancers and neoplasms', 'cloning', 'oncology', 'sarcoidosis', 'pareto distribution', 'neurological tumors', 'probability distribution', 'mathematics', 'molecular biology techniques', 'research and analysis methods', 'immune system proteins', 'white blood cells', 'animal cells', 'proteins', 'mathematical and statistical techniques', 't cells', 'molecular biology', 'immune response', 'probability theory', 'biochemistry', 'signal transduction', 't cell receptors', 'cell biology', 'neurology', 'biology and life sciences', 'cellular types', 'immune receptors', 'physical sciences', 'glioblastoma multiforme']