
==== Front
Brief Bioinform
Brief Bioinform
bib
Briefings in Bioinformatics
1467-5463
1477-4054
Oxford University Press

10.1093/bib/bbae431
bbae431
Problem Solving Protocol
AcademicSubjects/SCI01060
Robust detection of infectious disease, autoimmunity, and cancer from the paratope networks of adaptive immune receptors
Xu Zichang Department of Systems Immunology, Immunology Frontier Research Institute (IFReC), Osaka University, 3-1 Yamadaoka, Suita 565-0871, Japan

Ismanto Hendra S Department of Systems Immunology, Immunology Frontier Research Institute (IFReC), Osaka University, 3-1 Yamadaoka, Suita 565-0871, Japan
Department of Genome Informatics, Research Institute for Microbial Diseases (RIMD), Osaka University, 3-1 Yamadaoka, Suita 565-0871, Japan

Saputri Dianita S Department of Systems Immunology, Immunology Frontier Research Institute (IFReC), Osaka University, 3-1 Yamadaoka, Suita 565-0871, Japan
Department of Genome Informatics, Research Institute for Microbial Diseases (RIMD), Osaka University, 3-1 Yamadaoka, Suita 565-0871, Japan

Haruna Soichiro Department of Genome Informatics, Research Institute for Microbial Diseases (RIMD), Osaka University, 3-1 Yamadaoka, Suita 565-0871, Japan

https://orcid.org/0009-0008-4704-7072
Sun Guanqun School of information Science, Japan Advanced Institute of Science and Technology, 1-1 Asahidai, Nomi, Ishikawa 923-1292, Japan

Wilamowski Jan Department of Genome Informatics, Research Institute for Microbial Diseases (RIMD), Osaka University, 3-1 Yamadaoka, Suita 565-0871, Japan

Teraguchi Shunsuke Department of Genome Informatics, Research Institute for Microbial Diseases (RIMD), Osaka University, 3-1 Yamadaoka, Suita 565-0871, Japan
Faculty of Data Science, Shiga University 1-1-1 Banba, Hikone, Shiga 522-8522, Japan

Sengupta Ayan Cogent Labs, 3-2-1 Roppongi, Minato-ku, Tokyo 106-6122, Japan

Li Songling Department of Systems Immunology, Immunology Frontier Research Institute (IFReC), Osaka University, 3-1 Yamadaoka, Suita 565-0871, Japan
Department of Genome Informatics, Research Institute for Microbial Diseases (RIMD), Osaka University, 3-1 Yamadaoka, Suita 565-0871, Japan

https://orcid.org/0000-0003-4078-0817
Standley Daron M Department of Systems Immunology, Immunology Frontier Research Institute (IFReC), Osaka University, 3-1 Yamadaoka, Suita 565-0871, Japan
Department of Genome Informatics, Research Institute for Microbial Diseases (RIMD), Osaka University, 3-1 Yamadaoka, Suita 565-0871, Japan

Corresponding author. Department of Systems Immunology, Immunology Frontier Research Institute (IFReC), Osaka University, 3-1 Yamadaoka, Suita 565-0871, Japan. E-mail: standley@biken.osaka-u.ac.jp
9 2024
03 9 2024
03 9 2024
25 5 bbae43122 5 2024
19 7 2024
21 8 2024
© The Author(s) 2024. Published by Oxford University Press.
2024
https://creativecommons.org/licenses/by-nc/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (https://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com

Abstract

Liquid biopsies based on peripheral blood offer a minimally invasive alternative to solid tissue biopsies for the detection of diseases, primarily cancers. However, such tests currently consider only the serum component of blood, overlooking a potentially rich source of biomarkers: adaptive immune receptors (AIRs) expressed on circulating B and T cells. Machine learning–based classifiers trained on AIRs have been reported to accurately identify not only cancers but also autoimmune and infectious diseases as well. However, when using the conventional “clonotype cluster” representation of AIRs, individuals within a disease or healthy cohort exhibit vastly different features, limiting the generalizability of these classifiers. This study aimed to address the challenge of classifying specific diseases from circulating B or T cells by developing a novel representation of AIRs based on similarity networks constructed from their antigen-binding regions (paratopes). Features based on this novel representation, paratope cluster occupancies (PCOs), significantly improved disease classification performance for infectious disease, autoimmune disease, and cancer. Under identical methodological conditions, classifiers trained on PCOs achieved a mean AUC of 0.893 when applied to new individuals, outperforming clonotype cluster–based classifiers (AUC 0.714) and the best-performing published classifier (AUC 0.777). Surprisingly, for cancer patients, we observed that “healthy-biased” AIRs were predicted to target known cancer-associated antigens at dramatically higher rates than healthy AIRs as a whole (Z scores >75), suggesting an overlooked reservoir of cancer-targeting immune cells that could be identified by PCOs.

adaptive immune receptor
diagnosis
machine learning
liquid biopsy
Japan Agency for Medical Research and Development 10.13039/100009619 JP20am0101108
==== Body
pmcIntroduction

Liquid biopsies, which detect the presence of disease by analyzing biomarkers in bodily fluids, offer a promising approach to early disease detection [1]. Their low invasiveness and the ability to leverage trace amounts of biomarkers through machine learning represent a paradigm shift in medical diagnosis [2]. In the case of blood-based biopsies, both cellular and noncellular components carry information about disease status. Among blood cells, B and T lymphocytes are of particular interest, as they express adaptive immune receptors (AIRs), which recognize antigens associated with many diseases. Because of the critical importance of lymphocytes in both disease detection and prevention, a great deal of effort has been devoted to the development of bioinformatics methods to associate AIR coding sequences with their target antigens [3–6].

Furthermore, the rapid growth in publicly available AIR sequence data from diverse populations, along with the emergence of powerful machine learning methodologies for extracting biological knowledge from such data, has expanded the focus of AIR analysis beyond the molecular level to the donor level in order to predict disease states in individuals [3, 7–14]. Such AIR-based disease detection methods represent a new class of liquid biopsy [15]. Until now, most blood-based liquid biopsies have utilized only the noncellular components of blood—in particular extracellular RNAs or DNAs—and have been limited primarily to the detection of cancers [16–18]. By extending liquid biopsies to include AIRs, the breadth of detectable diseases can, in principle, be expanded beyond cancer to include infectious, autoimmune, and neurodegenerative diseases. This aligns with broader trends in utilizing immune-related markers for disease assessment, as seen in recent population studies on systemic immune-inflammation indicators [19].

The computational challenge in utilizing AIRs as a disease biomarker arises from their diversity: because each individual draws only a small fraction of gene rearrangements from a vast pool of possibilities (Fig. S1A), no two individuals are likely to share more than a few % of their B-cell receptors (BCRs) or T-cell receptors (TCRs) [20, 21]. Because of this “donor sharing problem,” it is necessary to group similar but nonidentical AIR sequences using some form of clustering in order to improve the overlap between AIRs from different individuals [5, 9, 11, 22]. As a consequence, some degree of information loss inevitably occurs upon clustering nonidentical AIRs from different donors.

In the current study, we explore two approaches to clustering AIRs in order to overcome the donor sharing problem: the conventional and more conservative “clonotype” approach and an alternative and more permissive “paratope” approach. A clonotype is defined by the V gene name, J gene name, and amino acid sequence of the third complementarity-determining region (CDR3), where the V and J coding sequences join. Paratopes, by contrast, are defined in terms of the amino acid sequences of all three CDRs: CDR1 + CDR2 + CDR3 (Fig. S1B). Because most physical contacts between AIRs and antigens occur within or near the CDRs (Fig. S1C), paratope sequences have been used to functionally cluster both BCRs [23] and TCRs [3]. Unlike clonotypes, which are comprised of categorical (V gene, J gene) and quasi-continuous (CDR3 sequence) data types, paratopes can be represented as a single sequence, which can easily be compared using conventional sequence alignment methods.

We next take advantage of the paratope representation to extend AIR clusters. To this end, we compute the network of paratope sequence similarities. This, in turn, allows an individual AIR to “occupy” more than one cluster if has edges to such cluster members in its similarity network. The rationale of paratope occupancies is that they will, to some extent, mitigate the information loss due to clustering.

To test this hypothesis, we constructed a benchmark of six diseases spanning three broad categories: infection, autoimmunity, and cancer. We found that the use of paratope networks significantly improved the accuracy of the resulting classifiers over clonotype cluster–based features, or previously published classifier methods, under identical conditions. In addition, we found that, in the case of cancer, paratope network–based features allowed the identification of a previously unnoticed reservoir of potentially therapeutic T cells.

Methods

Adaptive immune receptor datasets

AIR data were collected from published studies covering six diseases from three categories (Table S1): BCRs of 44 coronavirus disease 2019 (COVID-19) patients with 58 healthy donors [1]; BCRs of 94 human immunodeficiency virus (HIV) patients with 128 healthy donors [2]; TCRs of 59 autoimmune hepatitis (AIH) patients and 59 healthy donors [3]; TCRs of 34 type 1 diabetes (T1D) patients and 11 healthy donors [4]; TCRs of 204 nonsmall cell lung cancer (NSCLC) patients and 294 healthy donors from seven studies [3, 5–10]; and TCRs of 20 colorectal cancer (CRC) patients and 88 healthy donors [10].

Adaptive immune receptor representation

We represented AIRs as either clonotypes or paratopes. Clonotypes were constructed based on V and J gene names and the CDR3 amino acid sequence, as previously described [7, 11, 22, 24, 25]. Paratopes were constructed as pseudo sequences by concatenating the three CDR segments.

Clustering

Clustering was used to merge non-identical but highly similar AIRs. We used two different methods (clonotype and paratope) for clustering. To cluster clonotypes, AIRs with identical V genes, J genes, and CDR3 lengths were first grouped together; subsequently, the AIRs within these groups with identical length were clustered according to their CDR3 amino acid sequences. Here, we used 80% sequence identity and 100% alignment coverage, as implemented in CD-HIT [26]. Paratope pseudo sequence clustering was carried out using the MMseqs2 (version 13.45111) cluster method [27] with a sequence identity threshold of 80% and a coverage threshold of 90%, as described previously [28].

Adaptive immune receptor adjacency matrices

AIR adjacency matrices were computed from paratope pseudo sequences using the MMseqs2 search method [27], with a sequence identity threshold of 80% and a coverage threshold of 90%. Let us assume that there is a total of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $M$\end{document} distinct AIR sequences for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $D$\end{document} donors. For any given pair of AIRs denoted by\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${a}_i$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${a}_j$\end{document}, let us denote the sequence identity by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${p}_{ij}$\end{document} and the coverage by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${q}_{ij}$\end{document}. Then, the elements of the adjacency matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${A}^{M\times M}$\end{document} can be defined by equation (1):

(1) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{equation*} {A}_{ij}=\left\{\begin{array}{@{}l}1\ if\ {p}_{ij}>0.8\ and\ {q}_{ij}>0.9\\{}0\ \kern4pc otherwise\end{array}\right.\textrm{where}\ i\in \left[1,M\right]\ and\ j\in \left[1,M\right] \end{equation*}\end{document}

and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${A}_{ij}$\end{document} denotes the element in row \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $i$\end{document} and column \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $j$\end{document}. By definition, the diagonal elements of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $A$\end{document} are equal to unity.

Cluster frequencies

In order to build a feature vector that can be used to compare donors, we compute cluster frequencies. To this end, we first clustered all of the AIR data for a given disease dataset and described each donor by the frequency of AIRs in the resulting clusters. Hence, the feature vector for each donor is essentially a numerical representation that captures the distribution of their AIRs across the different clusters in the dataset.

Let us denote an AIR dataset by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${X}_{AIR}$\end{document}, and let us assume that there are \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $N$\end{document} clusters defined as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\left\{{S}_1\cdots{S}_N\right\}$\end{document} after performing the clustering as defined above. Then, each of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $M$\end{document} total AIRs belong to one of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $N$\end{document} clusters. Formally, we can define a set of distinct AIRs as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${X}_{AIR}=\left\{{a}_1\cdots{a}_M\right\}$\end{document}, where,

For all \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${a}_i\in{X}_{AIR}$\end{document}, there exists a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $k\in \left[1,N\right]$\end{document} such that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${a}_i\in{S}_k$\end{document}.

For all \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${a}_i,{a}_j\in{X}_{AIR}$\end{document}, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${a}_i={a}_j$\end{document} if and only if \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $i=j.$\end{document}

Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${C}^{M\times N}$\end{document} be an \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $M\times N$\end{document} matrix that contains the one-hot encoding of the cluster assignment for each AIR, such that each element \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${C}_{ik}\in{C}^{M\times N}$\end{document} is defined as follows:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{equation*} {C}_{ik}=\left\{\begin{array}{@{}cc}1& if\ {a}_i\in{S}_{k,}\ \\{}0& otherwise\end{array}\right. \quad{\textrm{where}}\ i\in \left[1,M\right]\ and\ k\in \left[1,N\right]. \end{equation*}\end{document}

That is, if AIR \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${a}_i$\end{document} belongs to cluster \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${S}_k$\end{document}, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${C}_{ik}$\end{document} is equal to one and zero for all the other values of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $k$\end{document}. If the set of rows corresponding to donor \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $D$\end{document} is defined as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${X}_{AIR}^D$\end{document}, the cluster frequency \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $CF$\end{document} of a given cluster \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${S}_k$\end{document} for donor \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $D$\end{document} can be defined as follows:

(2) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{equation*} {CF}_k^D={\sum}_{i\in{X}_{AIR}^D}{C}_{ik}\quad \textrm{ where}\ k\in \left[1,N\right]\end{equation*}\end{document}

Now, the feature vector for a donor \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $D$\end{document} is defined as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\{{CF}_1^D,\cdots, {CF}_k^D\}$\end{document}. The procedure above is illustrated in Fig. 2A.

Figure 1 Paratope networks enable donor sharing. BCR heavy-chain networks were constructed for 10 COVID-19 patients; in these networks, nodes are colored by the donor, and edges represent BCR sharing of the same clonotype or similar paratope. Only networks larger than 10 are shown for simplicity. (A) Clonotype networks are typically small and do not connect different donors. (B) Paratope networks are much larger and generally connect many donors. (C) network size and frequency were distinct even on a log scale. (D) A close-up view of one of the larger paratope networks [circled in (B)], which is composed of many different clonotypes and includes all 10 donors.

Figure 2 Cluster frequency and occupancy strategy in the classification algorithm and method workflow. (A) Cluster frequencies are illustrated with two example donors (A and B) whose AIRs form six clusters, which are one-hot encoded in matrix C. The cluster frequency feature of a donor is obtained by summing the rows of this matrix that belong to that donor. (B) Cluster occupancies are derived similarly to cluster frequencies except that we introduce an adjacency matrix describing the pairwise similarities between each AIR. If there are similarities between AIRs belonging to different clusters, the occupancy matrix O will differ from C. Summing the rows for each donor yields a feature vector that is generally less sparse than that derived from cluster frequencies. (C) Flowchart of the classification pipeline. Starting from a set of training and test donors, the first step is feature generation using one of four methods, which differ in the type of clustering (clonotype or paratope) and whether the clusters are transformed into features based on frequencies or occupancies. Thereafter, the remaining steps are identical and consist of training/hyperparameter selection and testing the classifier on the test features.

As mentioned above, we can perform the clustering in two different ways; hence, we can produce two different representations of feature vectors for each donor. We refer to features based on clonotype clusters as clonotype cluster frequency (CCF) and those based on paratope clusters as paratope cluster frequency (PCF).

Cluster occupancies

We utilized the adjacency matrix to identify the “occupancy” of a given AIR in a given cluster. The occupancy matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${O}^{M\times N}$\end{document} is defined as the matrix multiplication of the adjacency matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${A}^{M\times M}$\end{document} and cluster matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${C}^{M\times N}$\end{document}containing one-hot cluster encoding. Formally, it can be written as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${O}^{M\times N}={A}^{M\times M}$\end{document}\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${C}^{M\times N}$\end{document}, where each element \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${O}_{ik}\in{O}^{M\times N}$\end{document} can be defined as:

(3) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{equation*} {O}_{ik}=\sum_j^M{A}_{ij}\ {C}_{jk} \end{equation*}\end{document}

and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${A}_{ij}$\end{document}is defined by equation (1). If the adjacency matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $A$\end{document} contains nonzero elements that connect AIRs from different clusters, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${O}^{M\times N}$\end{document}will be less sparse than \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${C}^{M\times N}$\end{document}. Moreover, similar to equation (2), we can compute the overall cluster occupancy \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${CO}_k^D$\end{document} of cluster \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${S}_k$\end{document}for a given donor \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $D$\end{document} as follows:

(4) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} \begin{equation*} {CO}_k^D=\sum_{i\in{X}_{AIR}^D}{O}_{ik} \end{equation*}\end{document}

The cluster occupancy–based feature vector for a donor \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $D$\end{document} is defined as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\{{CO}_1^D,\cdots, {CO}_k^D\}$\end{document}.

We refer to cluster occupancy–based features derived from clonotype clusters and paratope clusters as clonotype cluster occupancy (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $CCO$\end{document}) and paratope cluster occupancy (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $PCO$\end{document}), respectively. The occupancy calculation is illustrated in Fig. 2B.

Cluster downsizing

Using all the cluster information for creating the feature vectors will often result in a high-dimensional feature space. This is mainly because the cluster sizes follow a long-tailed distribution, containing many clusters with very few members (Fig. S1D). Additionally, the time complexity of calculating the adjacency matrix is quadratic in nature; this results in very long computational times for large numbers of AIRs, making it increasingly difficult to calculate the full matrix for increasing dataset sizes. To address both problems, we propose an algorithm for decreasing the number of clusters and thus reduce the number of AIRs and features for each donor, i.e. we seek to downsize the number of clusters while retaining potentially meaningful AIRs for feature generation. To this end, we retain only those clusters that meet a donor diversity criterion, as described below.

We first define a set of binary classes \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $Y=\left\{{y}_1,{y}_2\right\}$\end{document}, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${y}_1$\end{document} denotes the class of healthy donors and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${y}_2$\end{document} denotes the class of donors with the disease. Next, we let the number of donors that belong to class \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${y}_p$\end{document} be defined by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${N}_p$\end{document}, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $p\in \left[1,2\right]$\end{document}. The number of donors that belong to class \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${y}_p$\end{document}and who have AIRs in cluster \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${S}_k$\end{document} can be defined by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${n}_p\left({S}_k\right).$\end{document} Then, we define the cluster ratio tuple of a cluster \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${S}_k$\end{document} as (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${r}_k^{y_1},{r}_k^{y_2}$\end{document}), where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${r}_k^{y_p}=\frac{n_p\left({S}_k\right)}{N_p}$\end{document}. In other words, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${r}_k^{y_p}$\end{document} is the ratio of the number of donors in class \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${y}_p$\end{document} in cluster \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${S}_k$\end{document} to the total number of donors in class \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${y}_p.$\end{document} We next introduce a threshold parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${d}_{min}$\end{document}, which defines a minimum ratio, below which all clusters are dropped and the respective AIRs ignored. For a cluster to be selected, at least one value of the cluster ratio tuple should be more than \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${d}_{min}$\end{document}. Specifically, for a cluster \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${S}_k$\end{document}, if \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} $\mathit{\max}\left({r}_k^{y_1},{r}_k^{y_2}\right)<{d}_{min}$\end{document}, then it will be dropped. This downsizing is applied to all the datasets in the context of training, with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${d}_{min}$\end{document} used as a hyperparameter. In addition, for the very large HIV dataset, we used this downsizing as a preprocessing step since the adjacency calculation was prohibitively time-consuming for the entire dataset. Importantly, such downsizing was performed after splitting the data into training and test sets, i.e. only the labels of the training AIRs were used to drop clusters.

Classification

For the binary classification task, we used the XGBoost algorithm via the XGBClassifier (version 1.6.2) Python API. The training was performed on a randomly selected 70% of the donors, and aside from weighting positive classes by the ratio of true to false training donors, all parameters were set to default values. Classification performance was evaluated on the remaining 30% of AIR data (the test set) using the area under the receiver operating characteristic curve (ROC AUC) or the area under the precision–recall curve (PR AUC) as implemented in the Scikit-Learn Python package (version 1.0.2). Four workflows were constructed, differing only in terms of the clustering method (clonotype or paratope) and feature construction (cluster frequency or occupancy), as shown in Fig. S1C.

Comparison with other repertoire classification methods

To systematically evaluate the performance of repertoire classification and disease diagnosis methods, two previously published algorithms, DeepRC [14] and immuneML [29], were evaluated using identical training and test data.

Healthy-biased and disease-biased clusters

After choosing the best \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{upgreek} \usepackage{mathrsfs} \setlength{\oddsidemargin}{-69pt} \begin{document} ${d}_{min}$\end{document} value in the training classification, we investigated the clusters to identify bias to either the patient or healthy donor classes. We divided clusters into those with significantly more TCRs in one class (healthy or disease) than in the other using a t test at a P-value cutoff of .001. We then ranked the healthy- or disease-biased clusters based on their P values to identify the 100 most significant healthy- or disease-biased clusters.

Cancer-associated T-cell receptor datasets

We collected TCR sequences that were annotated to target cancer-associated antigens in VDJdb (https://vdjdb.cdr3.net/) or McPAS (http://friedmanlab.weizmann.ac.il/McPAS-TCR/). After preprocessing, 3654 and 951 cancer-associated TCR beta chains were independently prepared from VDJdb and McPAS, respectively. These 4605 TCR beta chains were used to search healthy- or disease-biased clusters with the following criteria: sequence identity of CDR1 ≥ 90%; CDR2 ≥ 90%; CDR3 ≥ 80%; and alignment coverage ≥90%.

Comparison of cancer-associated T-cell receptors and background T-cell receptors

To fairly compare the hit rate of cancer-associated TCRs from a given set of clusters with background TCRs from the same class of donors (healthy or disease), we constructed a set of queries by randomly selecting an identical number of sequences from the corresponding donors. The selection was repeated 100 times. We calculated the Z score based on the mean and standard deviation of the hit ratio of these 100 background database queries and the hit rate of TCRs from the clusters of interest.

Results

Paratope adjacencies connect different clonotypes and donors

To address the limitations of conventional AIR-based disease classification methods, we developed a novel paratope-based classifier that represents AIRs using their paratope sequences and connects them by constructing paratope networks. To illustrate this approach, we constructed clonotype (Fig. 1A) and paratope (Fig. 1B) networks using the BCR data of 10 COVID-19 patients from a previous study [30]. Notably, none of the clonotype networks, but many of the paratope networks, connected different donors, consistent with the previously reported low sharing of BCR clones [21, 31]. The frequency of cluster sizes was also systematically lower for clonotypes than for paratopes (Fig. 1C). Indeed, the largest paratope network connected all 10 donors (Fig. 1D). These observations on clonotype and paratope networks form the basis of the features introduced here. The key observation is that AIRs can be “paratope adjacent” even if they come from different clonotypes or different donors. Therefore, features based on such paratope adjacency are likely to have elements that are shared across donors. This is essential to grouping different individuals belonging to a common disease class.

Table 1 Summary of performance metrics. Performance of previously published methods (DeepRC, immuneML) using various settings, along with the classifiers developed in this study, each trained on one of the four features (CCF, PCF, CCO, and PCO) and applied to six diseases (COVID-19, HIV, AIH, T1D, CRC, and NSCLC). The numbers of donors are listed under each disease. ROC AUC and PR AUC curves are given. For the PR AUC values, the expected performance of a random predictor (Rand) is given by the ratio of positive to negative donors. The abbreviations in the second column are as follows: convolutional neural network (CNN), long short-term memory (LSTM), logistic regression (LR), support vector machine (SVM), K-nearest neighbor (KNN), random forest (RF), and probabilistic binary classifier (PBC).

		COVID-19	HIV	AIH	T1D	CRC	NSCLC	Average
AUC	
Donors	Healthy	58	128	59	11	88	294	
Disease	44	94	59	34	20	204	
ROC AUC	DeepRC (CNN)	0.407	0.888	0.66	0.5	0.541	0.715	0.619	
DeepRC (LSTM)	0.813	0.861	0.5	0.5	0.458	0.741	0.646	
immuneML (LR)	0.717	0.946	0.688	0.65	0.780	0.802	0.764	
immuneML (SVM)	0.781	0.744	0.5	0.5	0.750	0.836	0.685	
immuneML (KNN)	0.565	0.851	0.656	0.525	0.760	0.795	0.692	
immuneML (RF)	0.838	0.876	0.706	0.5	0.850	0.895	0.777	
immuneML (PBC)	0.533	0.5	0.713	0.5	0.660	0.317	0.537	
CCF	0.648	0.941	0.653	0.500	0.779	0.762	0.714	
PCF	0.867	0.968	0.787	0.500	0.800	0.967	0.815	
CCO	0.863	0.975	0.828	0.625	0.700	0.409	0.733	
PCO	0.896	0.985	0.947	0.725	0.821	0.985	0.893	
PR AUC	DeepRC (CNN)	0.706	0.796	0.636	0.714	0.451	0.651	0.659	
DeepRC (LSTM)	0.928	0.856	0.444	0.714	0.416	0.696	0.676	
immuneML (LR)	0.796	0.914	0.729	0.871	0.718	0.802	0.805	
immuneML (SVM)	0.894	0.774	0.722	0.857	0.719	0.905	0.812	
immuneML (KNN)	0.71	0.785	0.7	0.835	0.782	0.827	0.773	
immuneML (RF)	0.878	0.831	0.757	0.857	0.810	0.851	0.831	
immuneML (PBC)	0.767	0.723	0.722	0.857	0.560	0.31	0.657	
CCF	0.818	0.936	0.544	0.857	0.332	0.605	0.682	
PCF	0.893	0.930	0.740	0.857	0.565	0.965	0.825	
CCO	0.890	0.947	0.808	0.829	0.390	0.371	0.706	
PCO	0.932	0.966	0.953	0.902	0.468	0.982	0.867	
Rand	0.516	0.236	0.444	0.714	0.152	0.421	0.413	
Bold values indicate the highest ROC AUC or PR AUC for each disease classification, representing the best-performing method for that particular disease.

Assessment across infectious disease, autoimmune disease, and cancer

As described in Methods, the clonotype cluster frequency (CCF) feature vector represents the conventional and more conservative definition of AIR data, while paratope cluster occupancy (PCO) represents a more permissive definition wherein a given AIR can occupy multiple clusters. As controls, two intermediate features, paratope cluster frequency (PCF) and clonotype cluster occupancy (CCO), were assessed. In addition, two classifiers based on DeepRC and five classifiers based on immuneML, were assessed. All features were applied to the same training and test datasets under the same conditions. We found that the PCO-based classifier clearly exhibited the best overall performance on independent test donors (Table 1). The PCO-based classifier achieved the highest ROC AUC among all the classifiers for each of the six diseases except for CRC, the disease with the fewest patients (20) and for which the PCO ROC AUC (0.821) was rather close to the highest value (0.850). The mean ROC AUC of the PCO-based classifier (0.893) was substantially higher than that of the CCF-based classifier (0.714) or the previously published methods (0.537–0.777). Three of the PCO ROC AUCs—those for HIV (0.985), AIH (0.947), and NSCLC (0.985)—were close to perfect. The remaining three diseases—COVID-19 (0.896), T1D (0.725), and CRC (0.821)—were represented by fewer disease donors, suggesting that a threshold number of training donors is needed for robust classifier performance. In contrast to PCO, the ROC AUCs of the CCF-based classifier exceeded 0.9 for only one disease (HIV), for which the sequencing depth was much greater than that of the remaining diseases and was no better than chance for another disease (T1D). Notably, removal of the paratope features reduced the CRC performance (0.714) to a value similar to that of immuneML (RF) (0.777), suggesting that the improvement was likely due to the use of paratope networks. As additional points of comparison, we investigated three repertoire encodings: the Kmer frequency of CDR3 sequences, V/J gene usage frequency, and evenness profiles using the same cross-validation framework as our original PCO-based classifier. The PCO-based classifiers outperformed these three feature-based classifiers on balance in terms of both ROC and PR AUC (Table S2). The above results demonstrate that PCO-based classifiers can successfully distinguish disease patients from healthy donors in infectious disease, autoimmunity, and cancer. We discuss the performance of the classifiers in the six diseases in detail in the Supplementary Results section.

Nonsmall cell lung cancer classifiers trained on paratope cluster occupancies were robust against batch effects

Because our large NSCLC dataset was composed of data from multiple studies, it provided the opportunity to systematically explore the robustness of the different features against batch effects, which arise from differences between samples that are not rooted in the experimental design [32]. To this end, we sampled all possible splits of the datasets where the test set consisted of the data from one healthy study and one NSCLC study, while the data from the remaining studies were used for training. This process resulted in 12 splits with various numbers of AIRs in the training and test sets. These results indicated that, in terms of the ROC AUC (Fig. 3A) and PR AUC (Fig. 3B), the classifiers based on PCF and PCO features performed well overall, but those based on the CCF and CCO features performed inconsistently, in agreement with the data shown in Figs S7C and D. Although 2 of the 12 splits showed that it was possible to combine the NSCLC studies in such a way that the PCO model failed to generalize to the test data, in these two splits, the size of the training dataset was much smaller than that of the test dataset. This result was consistent with the findings for the other under-sampled diseases such as T1D and CRC, again emphasizing the need for sufficient training data. Overall, this result demonstrated that the classifiers trained on PCOs were robust against batch effects, especially when compared to the purely cluster-based classifiers.

Figure 3 Holdout testing demonstrates the robustness of PCO against batch effects in NSCLC diagnosis. We collected a diverse dataset from seven studies: four from healthy individuals and three from NSCLC individuals. We exclusively split the holdout test set consisting of the data from one healthy study and one NSCLC study, while the data from the remaining studies were used for training. Twelve combinations of training and test sets were performed the identical training process. (A) Summary of 12 ROC AUC values in the test sets of NSCLC using different combinations of studies for training and testing. (B) PR AUC values in the same 12 test sets.

Surprising relationship between healthy-biased clusters and cancer antigen specificity

To interpret the underlying biological meaning of paratope networks, we investigated the clusters corresponding to features with high importance values, as determined by XGBoost. Because each feature corresponds to an AIR cluster, we hypothesized that AIRs within clusters corresponding to important features might be more likely to target antigens that are specific to the disease in question. To test this hypothesis, we again examined the NSCLC dataset, for which we had the largest amount of data. First, we ranked the PCO clusters by their feature importance (Fig. 4A) and examined the proportions of healthy- and NSCLC-derived TCRs in each. We observed both clusters with significantly more TCRs from healthy donors (“healthy-biased” clusters) and with significantly more TCRs from disease donors (“NSCLC-biased” clusters) (Fig. 4B). The features of these class-imbalanced clusters are shown as a heatmap in Fig. 4C. Next, using these class-imbalanced clusters, we searched two antigen-annotated TCR databases: McPAS [33] and VDJdb [34]. Interestingly, we identified strong hits to cancer-associated antigens using the TCRs from either the healthy- or NSCLC-biased clusters. Indeed, there were ‘more’ hits overall from the healthy-biased clusters than from the NSCLC-biased clusters (Fig. 5A). Figure 5B illustrates some of these hits, which included antigens such as MLANA [35], EphA2 [36], BST2 [37], TKT [38], and IGF2BP2 [39], all of which have been reported as prognostic markers for NSCLC. These findings demonstrate the potential for PCO-based features to identify cancer-fighting T cells.

Figure 4 Assessment of cluster features in NSCLC patient classification. (A) Feature importance. The Y-axis shows the importance of the XGBoost features in descending order. (B) The six most important features. Each dot represents the ratio value of healthy (HD) and disease (NSCLC) donors exhibiting those features. A t test was performed to compare the proportions between the HD and NSCLC cohorts (***P ≤ .001). (C) Heatmap of the top 100 significant disease and HD cluster features. P values were calculated with the t test; by sorting the P values in ascending order, the ranks of the corresponding features are determined. The value of each cell is the normalized ratio from 0 to 1. The left Y-axis shows each donor from the healthy and NSCLC cohorts grouped by the corresponding study.

Figure 5 Functional annotation of healthy- and NSCLC-biased clusters. (A) Histogram of hits to TCRs targeting the indicated antigens from healthy- and NSCLC-biased clusters after querying the McPAS and VDJdb cancer databases. Each count represents one database hit. (B) Sequences from important clusters matching known cancer-targeting antigens, as shown by their alignments and sequence logos based on all members of the same cluster. (C) Hit rates of TCRs from healthy-biased clusters to cancer-targeting TCRs (vertical line) and the corresponding hit rates of 100 randomly selected sets of queries from the same donors (histogram with vertical lines representing the means). The horizontal arrow indicates the Z score. (D) Hit rates of TCRs from NSCLC-biased clusters to cancer-targeting TCRs (vertical line) and the corresponding hit rates of 100 randomly selected sets of queries from the same donors (histogram with vertical lines representing the means). The horizontal arrow indicates the Z score.

We speculated that T cells that target cancer-associated antigens might be more common in healthy donors because they migrate from the peripheral blood toward tumor-presenting tissues (i.e. the lungs) in cancer patients. To test this idea, we repeated the database queries using randomly selected TCRs from healthy donors. Compared to TCRs from healthy-biased clusters, randomly selected TCRs from healthy donors resulted in dramatically fewer hits (Z score = 75), as shown in Fig. 5C. In contrast, TCRs from NSCLC-biased clusters exhibited a similar level of cancer-associated hits to randomly selected TCRs from donors with NSCLC (Z score = 1.41), as shown in Fig. 5D. These observations were not qualitatively sensitive to the similarity threshold used to define a hit (Fig. S8). Taken together, they show that the healthy-biased TCRs are indeed distinct from typical healthy donor–derived TCRs, supporting the notion that they would have the potential to migrate out of the periphery toward tumors in cancer patients. To investigate whether the presence of healthy-biased clusters with cancer-associated antigens is a general phenomenon, we performed an analogous analysis on CRC data. These results qualitatively agreed with the NSCLC observations above (Figs S9–S11).

Discussion

Liquid biopsies that can detect specific diseases from blood have the potential to reshape the future of medical diagnosis. Adaptive immune cells, which circulate in peripheral blood, are highly sensitive to a broad range of diseases. However, harnessing this sensitivity in a reproducible and general manner has presented a challenge due to the high diversity and low interdonor sharing of AIRs. Moreover, the natural inclination to gather more extensive datasets for machine learning purposes is at odds with achieving widespread mutual sharing among all donors. This problem is exacerbated by the use of the clonotype nomenclature, which separates AIRs into discrete clusters by their V and J gene names. Our results demonstrate that an AIR representation that incorporates the extended networks of paratopes allows donors to be connected and improves classifier performance compared to clonotype-based classifiers.

We found that there was obvious improvement with additional training data, which is encouraging given the rapid growth of publicly available AIR data (Fig. S12). Beyond a threshold value of 200 donors, the PCO-based classifiers were very accurate, achieving ROC AUC values of 0.985 (HIV and NSCLC). This suggests that initial clinical validation of specific diseases could be performed with smaller numbers of donors (e.g. 20 patients) and then scaled up to larger numbers (e.g. 100 patients) based on initial results.

In addition to these general observations, we found that in the two cancer datasets, clusters that corresponded to important features included both healthy- and disease-biased clusters. Typically, in repertoire analysis, there is a tendency to focus on clusters that are overrepresented in patients [40, 41]. However, we observed that healthy-biased clusters were important for disease classification as well. We hypnotized that these clusters are relevant to cancer but have migrated out of the peripheral blood (e.g. to the site of the tumor) in patients. A recent report showed that several cancer-specific TCRs could target multiple tumor types via HLA A*02:01-restricted epitopes EAAGIGILTV, LLLGIGILVL, and NLSALGIFST from Melan A (MLANA), BST2, and IMP2 (IGF2BP2) and that PBMCs from healthy donors expressed such TCRs [42]. Consistently, many TCRs from healthy-biased clusters were predicted to target these antigens in both NSCLC (Fig. S13) and CRC (Fig. S14) datasets. Our observations on healthy-biased T-cell clusters may relate to recent research on T-cell exhaustion in cancer immunotherapy [43], suggesting complex interactions between circulating T cells and tumor microenvironments. These findings highlight the interpretability of our Machine learning (ML) classifiers along with their ability to identify potential therapeutic immune cells against tumors.

AIR-based liquid biopsies may have diverse applications in combination with other advanced medical technologies including: integration with advanced natural language processing techniques for medical records [44]; combining with advanced sampling techniques, such as robotic bronchoscope systems [45]; and integrating with emerging tissue analysis techniques, like virtual staining methods [46]. The above applications could further improve disease detection accuracy, potentially enhance AIR data interpretation, and provide a more comprehensive view of immune responses in various diseases.

Our study represents the first use of paratope networks as a feature for disease detection. The results indicate that use of this novel feature improves the robustness of disease classifiers across a range of disease types. Although paratope clustering is computationally more complex than clonotype clustering due to the lack of partitioning into V-J-length groups, the complexity is offset by the use of mmseqs to efficiently perform the adjacency and clustering steps.

Since most liquid biopsies utilizing peripheral blood utilize only the serum component, our technology may add breadth and sensitivity to existing tests. As the field of AIR repertoire analysis continues to evolve, we anticipate that further refinement will enhance the performance and generalizability of this approach, ultimately leading to improved liquid biopsies for a wide range of diseases.

Key Points

Novel Representation with Improved Classification: This study introduces a novel representation of AIRs based on the occupancy of paratope similarity networks (PCO), which significantly enhances disease classification performance for infectious diseases, autoimmune diseases, and cancer compared to traditional clonotype-based methods.

Robust Performance Across Diseases: Classifiers trained on PCOs achieved high mean AUC values, outperforming both clonotype cluster–based classifiers and previously published methods, demonstrating robust performance in distinguishing disease patients from healthy donors across multiple diseases, including COVID-19, HIV, AIH, T1D, NSCLC, and CRC.

Identification of Therapeutic Immune Cells: The study reveals the existence of a reservoir of cancer-targeting immune cells in healthy individuals, which are identifiable through paratope networks and may have therapeutic potential.

Addressing Donor-Sharing Problem: The paratope-based approach mitigates the information loss due to the vast diversity and low interdonor sharing of AIRs, enhancing the ability to group different individuals belonging to a common disease class and improving classifier generalizability.

Future Implications for Liquid Biopsies: This study highlights the potential of using paratope networks to improve the breadth and sensitivity of liquid biopsies, which traditionally focus on noncellular blood components.

Supplementary Material

Suppl_Figures-r2_bbae431

TableS1_bbae431

TableS2_bbae431

Supplementary_Material-0517_bbae431

Supplementary_Material-0830_bbae431

Acknowledgements

The authors thank all members of the Genome Informatics and System Immunology labs for valuable discussions in the development of this project.

Funding

This work was funded by the Japan Agency for Medical Research and Development (AMED), Platform Project for Supporting Drug Discovery and Life Science Research (Basis for Supporting Innovative Drug Discovery and Life Science Research), under JP20am0101108.

Conflict of interest: D.M.S, S.L. and Z.X. have a patent pending for Generalizable features for the diagnosisof disease from adaptive immune receptor repertoires. D.M.S and S.L are shareholders in MyImmune Corporation.

Data availability

All AIR sequence data of this study are already published and available. We have uploaded the input and output data files as well as the training model scripts and documentation on GitLab at: https://gitlab.com/sysimm/2024-bib-data.

Author contributions

Z.X., S.L., and D.M.S. contributed to the conceptualization of the study and its design. They also conducted the investigation by performing the experiments and contributed to writing the original draft of the manuscript. Z.X., S.L., K.W., S.T., A.S., and D.M.S. were involved in developing the methodology, specifically the method pipeline. G.S. provided resources in the form of GPU computational support. H.S.I., D.S.S., G.S., and S.H. participated in critical project discussions and contributed to the writing process through review and editing. D.M.S. took on the role of supervision for the study.
==== Refs
References

1. Lone SN , NisarS, MasoodiT. et al. Liquid biopsy: a step closer to transform diagnosis, prognosis and future of cancer treatments. Mol Cancer 2022;21 :79.35303879
2. Ko J , BaldassanoSN, LohPL. et al. Machine learning to detect signatures of disease in liquid biopsies - a user's guide. Lab Chip 2018;18 :395–405. 10.1039/c7lc00955k.29192299
3. Dash P , Fiore-GartlandAJ, HertzT. et al. Quantifiable predictive features define epitope-specific T cell receptor repertoires. Nature 2017;547 :89–93. 10.1038/nature22383.28636592
4. Glanville J , HuangH, NauA. et al. Identifying specificity groups in the T cell receptor repertoire. Nature 2017;547 :94–8. 10.1038/nature22976.28636589
5. Sidhom JW , BarasAS. Deep learning identifies antigenic determinants of severe SARS-CoV-2 infection within T-cell repertoires. Sci Rep 2021;11 :14275.34253751
6. Xu Z , LiS, RozewickiJ. et al. Functional clustering of B cell receptors using sequence and structural features. Mol Syst Des Eng 2019;4 :769–78.
7. Chen Y , YeZ, ZhangY. et al. A deep learning model for accurate diagnosis of infection using antibody repertoires. J Immunol 2022;208 :2675–85.35606050
8. Foers AD , ShoukatMS, WelshOE. et al. Classification of intestinal T-cell receptor repertoires using machine learning methods can identify patients with coeliac disease regardless of dietary gluten status. J Pathol 2021;253 :279–91. 10.1002/path.5592.33225446
9. Ostrovsky-Berman M , FrankelB, PolakP. et al. Immune2vec: embedding B/T cell receptor sequences in R (N) using natural language processing. Front Immunol 2021;12 :680687.34367141
10. Park JJ , LeeKAV, LamSZ. et al. Machine learning identifies T cell receptor repertoire signatures associated with COVID-19 severity. Commun Biol 2023;6 :76.36670287
11. Shemesh O , PolakP, LundinKEA. et al. Machine learning analysis of naive B-cell receptor repertoires stratifies celiac disease patients and controls. Front Immunol 2021;12 :627813.33790900
12. Cinelli M , SunY, BestK. et al. Feature selection using a one dimensional naive Bayes' classifier increases the accuracy of support vector machine classification of CDR3 repertoires. Bioinformatics 2017;33 :951–5. 10.1093/bioinformatics/btw771.28073756
13. Eliyahu S , SharabiO, ElmedviS. et al. Antibody repertoire analysis of hepatitis C virus infections identifies immune signatures associated with spontaneous clearance. Front Immunol 2018;9 :3004.30622532
14. Widrich M , SchäflB, PavlovićM. et al. Modern hopfield networks and attention for immune repertoire classification. Advances in neural information processing systems 2020;33 :18832–45.
15. Snir T , EfroniS. T cell repertoire sequencing as a cancer's liquid biopsy—can we decode what the immune system is coding? Curr Opin Syst Biol 2020;24 :135–41.
16. Cescon DW , BratmanSV, ChanSM. et al. Circulating tumor DNA and liquid biopsy in oncology. Nat Cancer 2020;1 :276–90. 10.1038/s43018-020-0043-5.35122035
17. Ignatiadis M , SledgeGW, JeffreySS. Liquid biopsy enters the clinic - implementation issues and future challenges. Nat Rev Clin Oncol 2021;18 :297–312. 10.1038/s41571-020-00457-x.33473219
18. Tomasik B , SkrzypskiM, BienkowskiM. et al. Current and future applications of liquid biopsy in non-small-cell lung cancer-a narrative review. Transl Lung Cancer Res 2023;12 :594–614. 10.21037/tlcr-22-742.37057121
19. Zhang Y , MengY, ChenM. et al. Correlation between the systemic immune-inflammation indicator (SII) and serum ferritin in US adults: a cross-sectional study based on NHANES 2015-2018. Ann Med 2023;55 :2275148.37883981
20. Robins HS , SrivastavaSK, CampregherPV. et al. Overlap and effective size of the human CD8+ T cell receptor repertoire. Sci Transl Med 2010;2 :47ra64.
21. Soto C , BombardiRG, BranchizioA. et al. High frequency of shared clonotypes in human B cell receptor repertoires. Nature 2019;566 :398–402. 10.1038/s41586-019-0934-8.30760926
22. Roskin KM , JacksonKJL, LeeJY. et al. Aberrant B cell repertoire selection associated with HIV neutralizing antibody breadth. Nat Immunol 2020;21 :199–209. 10.1038/s41590-019-0581-0.31959979
23. Richardson E , GalsonJD, KellamP. et al. A computational method for immune repertoire mining that identifies novel binders from different clonotypes, demonstrated by identifying anti-pertussis toxoid antibodies. MAbs 2021;13 :1869406.33427589
24. Miho E , RoskarR, GreiffV. et al. Large-scale network analysis reveals the sequence space architecture of antibody repertoires. Nat Commun 2019;10 :1321.30899025
25. Ruiz Ortega M , SpisakN, MoraT. et al. Modeling and predicting the overlap of B- and T-cell receptor repertoires in healthy and SARS-CoV-2 infected individuals. PLoS Genet 2023;19 :e1010652.36827454
26. Fu L , NiuB, ZhuZ. et al. CD-HIT: accelerated for clustering the next-generation sequencing data. Bioinformatics 2012;28 :3150–2. 10.1093/bioinformatics/bts565.23060610
27. Steinegger M , SodingJ. MMseqs2 enables sensitive protein sequence searching for the analysis of massive data sets. Nat Biotechnol 2017;35 :1026–8. 10.1038/nbt.3988.29035372
28. Saputri DS , IsmantoHS, NugrahaDK. et al. Deciphering the antigen specificities of antibodies by clustering their complementarity determining region sequences. mSystems 2023;8 :e0072223.37975681
29. Pavlovic M , SchefferL, MotwaniK. et al. The immuneML ecosystem for machine learning analysis of adaptive immune receptor repertoires. Nat Mach Intell 2021;3 :936–44. 10.1038/s42256-021-00413-z.37396030
30. Edahiro R , ShiraiY, TakeshimaY. et al. Single-cell analyses and host genetics highlight the role of innate immune cells in COVID-19 severity. Nat Genet 2023;55 :753–67. 10.1038/s41588-023-01375-1.37095364
31. Briney B , InderbitzinA, JoyceC. et al. Commonality despite exceptional diversity in the baseline human antibody repertoire. Nature 2019;566 :393–7. 10.1038/s41586-019-0879-y.30664748
32. Sprang M , Andrade-NavarroMA, FontaineJF. Batch effect detection and correction in RNA-seq data using machine-learning-based automated assessment of quality. BMC Bioinformatics 2022;23 :279.35836114
33. Tickotsky N , SagivT, PriluskyJ. et al. McPAS-TCR: a manually curated catalogue of pathology-associated T cell receptor sequences. Bioinformatics 2017;33 :2924–9. 10.1093/bioinformatics/btx286.28481982
34. Shugay M , BagaevDV, ZvyaginIV. et al. VDJdb: a curated database of T-cell receptor sequences with known antigen specificity. Nucleic Acids Res 2018;46 :D419–27. 10.1093/nar/gkx760.28977646
35. Der SD , SykesJ, PintilieM. et al. Validation of a histology-independent prognostic gene signature for early-stage, non-small-cell lung cancer including stage IA patients. J Thorac Oncol 2014;9 :59–64. 10.1097/JTO.0000000000000042.24305008
36. Brannan JM , SenB, SaigalB. et al. EphA2 in the early pathogenesis and progression of non-small cell lung cancer. Cancer Prev Res (Phila) 2009;2 :1039–49. 10.1158/1940-6207.CAPR-09-0212.19934338
37. Suzuki K , KachalaSS, KadotaK. et al. Prognostic immune markers in non-small cell lung cancer. Clin Cancer Res 2011;17 :5247–56. 10.1158/1078-0432.CCR-10-2805.21659461
38. Niu C , QiuW, LiX. et al. Transketolase serves as a biomarker for poor prognosis in human lung adenocarcinoma. J Cancer 2022;13 :2584–93. 10.7150/jca.69583.35711845
39. Han L , LeiG, ChenZ. et al. IGF2BP2 regulates MALAT1 by serving as an N6-Methyladenosine reader to promote NSCLC proliferation. Front Mol Biosci 2021;8 :780089.35111811
40. Huang C , LiX, WuJ. et al. The landscape and diagnostic potential of T and B cell repertoire in immunoglobulin a nephropathy. J Autoimmun 2019;97 :100–7. 10.1016/j.jaut.2018.10.018.30385082
41. Liu X , ZhangW, ZhaoM. et al. T cell receptor beta repertoires as novel diagnostic markers for systemic lupus erythematosus and rheumatoid arthritis. Ann Rheum Dis 2019;78 :1070–8. 10.1136/annrheumdis-2019-215442.31101603
42. Dolton G , RiusC, WallA. et al. Targeting of multiple tumor-associated antigens by individual T cell receptors during successful cancer immunotherapy. Cell 2023;186 :3333, e3327–49.37490916
43. Wang X , YangT, ShiS. et al. Heterogeneity-induced NGF-NGFR communication inefficiency promotes mitotic spindle disorganization in exhausted T cells through PREX1 suppression to impair the anti-tumor immunotherapy with PD-1 mAb in hepatocellular carcinoma. Cancer Med 2024;13 :e6736.38204220
44. Li Q , YouT, ChenJ. et al. LI-EMRSQL: linking information enhanced Text2SQL parsing on complex electronic medical records. IEEE Trans Reliab 2024;73 :1280–90.
45. Duan X , XieD, ZhangR. et al. A novel robotic bronchoscope system for navigation and biopsy of pulmonary lesions. Cyborg Bionic Syst 2023;4 :0013.36951809
46. Liu Z , ChenL, ChengH. et al. Virtual formalin-fixed and paraffin-embedded staining of fresh brain tissue via stimulated Raman CycleGAN model. Sci Adv 2024;10 :eadn3426.38536925
