We collected the CDR3 of TCRβ from TCR-seq data to study. Considering that low-quality CDR3 sequences will affect the downstream analysis, the following types of sequences were removed as described in the previous study (Beshnova et al., 2020): 1) too short (<10) or too long (>24) sequences; 2) sequences containing special characters (X, +, ∗, etc.); 3) incomplete sequences, according to the ImMunoGeneTics (IMGT) nomenclature (Lefranc et al., 2015), not starting with the cysteine (C) or not ending with the phenylalanine (F); and 4) sequences where variable gene locus was not solved. There were some overlapping TCR sequences between healthy individuals and cancer patients, which were considered irrelevant to cancer and therefore also needed to be excluded. We used known training data to generate a reference dataset (denoted as D_R) containing sequences that appeared at a high frequency in both healthy individuals and cancer patients (the top 20,000 sequences with the highest cloning frequency in each TCR-seq sample). Any sequences in each sample appearing in the D_R were removed. After the above sequences were removed, the remaining TCR sequences were sorted in descending order of cloning frequency and the top k sequences were extracted for downstream analysis.

The raw TCR sequences were not directly input to the CNN because their antigen-binding ability was not well represented. AAs in TCR sequences can be represented by biochemical features, and a TCR sequence with length l is able to be encoded by a 20 × d feature matrix for 20 AAs into an l × d TCR matrix. The AA index database (Kawashima and Kanehisa, 2000) documented 566 AA indices containing rich biochemical information based on previous literature. Since the original 566 AA indices are very large, directly utilizing them to characterize AAs may lead to a large input data size and too many parameters of CNN, which may cause problems such as high computational complexity and overfitting. In addition, many of the original AA indices are highly correlated with each other. As a result, we considered dimensionality reduction of the original AA indices. Currently, many studies have extracted low-dimensional orthogonal features from high-dimensional AA indices, reducing the dimensionality of a large number of AA indices with minimal information loss. Kidera et al. derived a 20 × 10 feature matrix from 188 AA indices (Kidera et al., 1985), and Atchley et al. derived a 20 × 5 feature matrix from 494 AA indices (Atchley et al., 2005). Beshnova et al. employed principal component analysis (PCA) to generate a 20 × 15 feature matrix from 531 AA indices to characterize AAs (Beshnova et al., 2020). Considering that the Beshnova matrix was obtained from the largest number of AA indices and encompassed the most biochemical information (explaining more than 95% of the variation in the data), we used the Beshnova matrix in our experiments (d = 15) to encode the TCR sequences into matrixes.

The CNN is able to predict whether a TCR sequence is associated with cancer. TextCNN in natural language processing (NLP) firstly applied the CNN model to sequence analysis (Zhang and Wallace, 2015). Referring to the idea of TextCNN, the model consisted of a convolutional layer, a pooling layer, and a linear layer, as shown in Figure 1A. Different from TextCNN, we developed innovative convolution filters in the convolutional layer to handle the amino acid fragments with different lengths according to the X-ray crystal structure analysis (Ostmeyer et al., 2019). For a TCR matrix with dimension l × d, the model extracted its features by a set of convolution filters with various sizes respectively, performed 1-max pooling operations for the low-dimensional feature maps obtained by each convolution filter, and concatenated the pooled results to generate a TCR feature vector. The TCR score, the probability of being a caTCR, was obtained by a one-layer linear classifier finally.

The convolutional layer of the model was designed with multiple convolution filters to extract the key biochemical motifs with distinct lengths in TCRs. By analyzing the X-ray crystal structure, Ostmeyer et al. identified the CDR3 residues in contact with peptide, which were thought to make the greatest contribution to the antigen-binding specificity of the TCRs, and determined the length of the AA fragments according to the analysis (Ostmeyer et al., 2019). We counted the contiguous CDR3 residue regions (size ≥2) in the 55 CDR3 sequences used for analysis (Table 1) and observed that the regions ranging in size from 2 to 8 were present, with regions of sizes 2–4 occurring more frequently and regions of size eight occurring least frequently at 0.017. These contiguous regions were considered as potential cancer-specific motifs that contribute to the antigen-binding ability of TCRs. After excluding the regions with low frequencies (<0.05), we designed a set of various convolution filters and specified the number of corresponding convolution filters according to the occurrence frequencies of the regions (Table 1). As shown in Table 1 and Figure 1A, the convolution filters with six different sizes were designed to extract features from TCR matrixes, and the number of convolution filters with each size was positively correlated with the occurrence frequency of the corresponding region, for a total of 14 convolution filters. Every complete convolution operation is defined as: 𝓸 a i , F j = σ ( W F j ⋅ M i [ a : a + h − 1 ] + b F j ) ∈ 𝓸 i , F j = [ 𝓸 1 i , F j , … , 𝓸 l − h + 1 i , F j ] T , (1) where the output sequence 𝓸 i , F j ∈ ℝ l − h + 1 is composed of the results of each convolution 𝓸 a i , F j ( a = 1 , … , l − h + 1 ) , the activation function σ ( x ) = max ( 0 ,   x ) is the Rectified Linear Unit (ReLU), W F j ∈ ℝ h × d and b F j ∈ ℝ are respectively the weight matrix and bias of the jth convolution filter F j with size h × d , M i ∈ ℝ l × d is the ith TCR matrix, M i [ a : b ] denotes the submatrix from row a to row b in M i , and ∙ denotes the dot product (a sum of element-wise multiplications) between the submatrix and the filter. Equation 1 showed the convolution operation that one convolution filter performed on the TCR matrix, and then the corresponding feature map, o i , F j , was obtained. Taking the 11 × d TCR matrix shown in Figure 1B as an example, when the convolution filter with size 2 × d performed a complete convolution operation on the matrix from top to bottom, it could be regarded as extracting the biochemical features of the 2-mers in TCRs such as "CA", "AS", etc., and finally obtaining the 10 × 1 feature map, the feature set of all 2-mers, which probably included cancer-specific motifs. The filters with other sizes performed similar convolution operations. Because the frequencies of key motifs with different lengths varied, the numbers of convolution filters with various sizes were adjusted to give them appropriate weights in the model. ReLU was chosen as the activation function following the convolution operations due to its low computational complexity. Additionally, it can avoid the vanishing gradient or exploding gradient problems that Sigmoid and Tanh can cause. And the disadvantage of ReLU, the dead ReLU problem, was mitigated by using the Xavier initialization (Glorot and Bengio, 2010).

The 1-max pooling function was adopted to pool the feature maps generated after convolution operations of each convolution filter, reducing the mapping dimension to 1. The following pooling functions are commonly used: 1) the element-wise maximum function ( f Max ); 2) the element-wise average function ( f Avg ); 3) the log-sum-exp (LSE) function ( f LSE ) (Sahasrabudhe et al., 2021). These are defined as: f Max ( 𝓸 i , F j ) = max ( { 𝓸 1 i , F j , … , 𝓸 l − h + 1 i , F j } ) , (2) f Avg ( 𝓸 i , F j ) = 1 l − h + 1 ∑ a = 1 l − h + 1 𝓸 a i , F j ,   and (3) f LSE ( 𝓸 i , F j ) = log ( 1 l − h + 1 ∑ a = 1 l − h + 1 exp ( 𝓸 a i , F j ) ) , (4) where 𝓸 i , F j ∈ ℝ l − h + 1 is the output sequence of F j with size h × d after the convolution operation on M i . The feature maps obtained could be viewed as the feature sets of all z-mers ( z = 2 , … , 7 ) . Regarding that the goal of CNN was distilling the potential cancer-specific motifs in TCRs, we focused on the most contributing features in the feature sets, which are most likely from the key motifs, and ignored the features of other z-mers. With the 1-max pooling function, the most contributing features of each feature set were extracted and other non-key motifs’ features were discarded, whereas the features pooled by the average function and the LSE function were affected by other non-key motifs, which caused adverse effects on the classification ability of the model. Considering that one TCR may contain multiple short key motifs with the same length, the drawback that the 1-max pooling operation can only extract the feature of one key motif from one sequence can be compensated by the design of multiple convolution filters with the same size in the convolutional layer. As shown in Figure 1B, the most contributing features (marked with blue boxes) extracted by pooling were interconnected to generate a 14 × 1 TCR feature vector as: p i , F j = f Max ( 𝓸 i , F j ) ∈ p i = [ p i , F 1 ,   … , p i ,     F 14 ] T . (5) Equation 5 showed the 1-max pooling operation performed on the feature map, and then the TCR feature vector, p i , was obtained, which consisted of the most contributing features from the feature maps.

Ultimately, a one-layer linear classifier L was applied to aggregate the features extracted from each convolution kernel and predict the score for that TCR sequence. L that assigned scores to the TCRs is given by: y ˜ i = P ( y i = 1 | M i ) = σ ′ ( W L T p i + b L ) , (6) where P ( y i = 1 | M i ) denotes the probability that the TCR is associated with cancer, the activation function σ ′ ( x ) = 1 / ( 1 + exp ( − x ) ) is the sigmoid function to normalize the scores, and W L ∈ ℝ 14 and b L ∈ ℝ are respectively the weight matrix and bias of L. Equation 6 showed the operation process in L, and the probability that the TCR was associated with cancer was obtained. The TCR was predicted to be caTCR when y ˜ i > 0.5 , and was otherwise predicted to be noncancerous. A multi-layer linear classifier is capable of fitting the data better, but it also makes the structure of CNN more complicated and introduces the risk of overfitting. To reduce overfitting, a one-layer linear classifier is applied to the model to predict the TCR scores.

The CNN model jointly learns the various convolution filters and L so that it is end-to-end trainable and the preprocessed TCRs and the corresponding labels are needed for model training. Because the probability that a TCR is a caTCR obeys a Bernoulli distribution, the log-likelihood function (also known as the cross-entropy function) was used as the loss function to train the model, which is defined as: ℒ CNN = − [ y ˜ i ln y ˜ i + ( 1 − y ˜ i ) ln ( 1 − y ˜ i ) ] . (7)

During the training process, random dropouts at a rate of 40% were applied to L to mitigate overfitting.

Because zero-padding will alter the distribution of input data, DeepCAT (Beshnova et al., 2020) built five distinct CNN models for TCR sequences from 12 to 16 in length, which made itself cumbersome and unable to utilize the information of TCRs with other lengths, whereas TCRs with different lengths were processed by one model in this method. In the actual training process, we added zeros to the end of the shorter sequences to achieve the length of the longest sequences (l = 24) to ensure the consistency of the dimensionality of the input TCR matrixes. Because the features of all motifs containing zero viewed as non-key motifs were discarded after the 1-max pooling operations in the model, the classification ability of the model did not deteriorate.

Predicting whether a repertoire is cancer-associated from the TCRs in every repertoire can be described as MIL, where the TCRs are instances and the repertoires are bags. The standard MIL assumption assumes that each instance in a bag can be classified as either positive (1) or negative (0), and the label of a bag is 1 when including one or more positive instances (Dietterich et al., 1997; Foulds and Frank, 2010). Ostmeyer et al. applied this assumption in their study (Ostmeyer et al., 2019). However, it is inappropriate to predict the repertoire as having cancer by the presence of a non-normal TCR because a cancer patient usually contains many caTCRs, which are somehow related to each other. In addition, the needed labels of TCRs are unknown, which means that it is difficult to know whether a TCR is associated with cancer or not. Although Beshnova et al. obtained potential caTCRs used for model training from TCGA tumor RNA-seq samples in advance in their study by TRUST (Li et al., 2017) and sequence filtering based on a reference database (Beshnova et al., 2020), the evidence that all sequences were positive was lacking. The definition of cancer score for the repertoire by averaging the predictions for TCRs in a repertoire was also inaccurate because we could not prove that all TCRs enjoyed the same weight. Therefore, we designed a one-layer linear classifier L ′ as an aggregating function to combine the scores of k TCRs collected from repertoires to predict the repertoire. L ′ is defined as: Y ˜ = P ( Y = 1 | { M 1 , … , M k } ) = σ ′ ( W L ′ T [ y ˜ 1 , … , y ˜ k ] T + b L ′ ) , (8) where P ( Y = 1 | { M 1 , … , M k } ) denotes the probability that the repertoire has cancer, the activation function σ ′ ( x ) is the sigmoid function, and W L ′ ∈ ℝ k and b L ′ ∈ ℝ are respectively the weight matrix and bias of L ′ . Equation 8 showed the operation process in L ′ , which was similar to Eq. 6, and the probability that the repertoire was associated with cancer was obtained. The repertoire was predicted to be cancer-associated when Y ˜ > 0.5 , and was otherwise predicted to be noncancerous. A MIL model consisting of the CNN and L ′ is also end-to-end trainable, whose loss function is the log-likelihood function defined as: ℒ MIL = − [ Y ˜ ln Y ˜ + ( 1 − Y ˜ ) ln ( 1 − Y ˜ ) ] . (9)

Instead of simply taking the max value (Ostmeyer et al., 2019) or the average (Beshnova et al., 2020) value of all TCR scores as the cancer score of the repertoire, L ′ was capable of learning the interrelationships among TCRs and assigning the appropriate weights to each TCR after model training. Similar to the idea of L, the multi-layer linear classifier was replaced by L ′ and random dropouts at a rate of 40% were applied to L ′ during training in order to alleviate overfitting.

We used the publicly available TCR-seq data from Adaptive Biotechnologies immuneACCESS online database (IA), which contains eight groups of peripheral blood mononuclear cell (PBMC) samples with diverse cancer types and one group of non-cancer PBMC samples. To validate the performance of the models on Asian patients, the TCR-seq data from the clinical database of Geneplus Technology Ltd. in Shenzhen (Geneplus) were also used (Lan et al., 2020; Li et al., 2021), which included PBMC and tumor-infiltrating T lymphocyte (TIL) samples from patients with thyroid cancer (THCA) and lung cancer, as well as non-cancer PBMC samples. Table 2 shows the specifics of the various datasets that were used in the experiments.

In section 3.2, we gathered the training data of DeepCAT to train the CNN framework in DeepLION. The label-encoded training data was composed of cancer (n = 30,000) and non-cancer (n ∼ 60,000) TCR CDR3 sequences. The cancer sequences were derived from The Cancer Genome Atlas (TCGA) (Tomczak et al., 2015) tumor RNA-sequencing (RNA-seq) samples covering multiple cancers using TRUST algorithm (Li et al., 2017), whereas non-cancer sequences were derived from the healthy individuals H₂ (n = 120) (Emerson et al., 2017), which are independent of healthy individuals H₁ (n = 666). The test data consisted of two datasets (T₁ and T₂). T₁ consisted of eight groups of samples, each of which contained cancer and H₁ PBMC samples from IA. All TCR-seq data in T₂ were obtained from Geneplus, which were collected from Asian populations. T₂ was composed of two groups of samples: the first group contained THCA PBMC & TIL samples (n = 170) and the PBMC samples from the healthy individuals H₃ (n = 260), and the second group contained lung cancer PBMC and TIL samples (n = 184) and H₃ PBMC samples. Table 3 shows the training and test data in section 3.2.

In section 3.3, we used T₂ to train and test the entire DeepLION due to the larger sizes of cancer samples in T₂ (Table 3). The samples in each group were randomly divided into five equal parts, three of which were used as the training set, one as the validation set, and one as the test set.

DeepCAT (M^CAT) (Beshnova et al., 2020) is currently a preferred method for de novo caTCR prediction from the peripheral blood. To validate the performance of the CNN framework in DeepLION (M^CNN) when predicting the caTCRs, we conducted experiments to compare M^CNN with M^CAT. M^CAT and M^CNN shared the same training data (Table 3), and we processed the data as described in M^CAT.

Before the training process, all the training sequences were encoded into l × 15 TCR matrixes (d = 15) by the Beshnova matrix (Beshnova et al., 2020). To gain the final model, we trained M^CNN five times independently. In each training process, we randomly selected two-thirds of the training data (20,000 cancer and 40,000 non-cancer samples) as the training set and used them to train the model for 1,000 epochs at a learning rate of 0.001 with the assurance that the model reached convergence. The remaining training data were used as the validation set and AUCs were estimated to evaluate the trained models. Ultimately, the five trained models had similar AUCs, and the model with the highest AUC (0.85) was selected as the final model. The trained M^CAT was obtained from Github.

All the test data (Table 3) were processed in several steps. First, the top 10,000 most abundant sequences (k = 10,000) of each sample were extracted after low-quality sequences were removed. Second, we used iSMART (Zhang et al., 2020) with default parameters to cluster the sequences, and then the antigen-specific sequences were selected. Third, all the sequences were encoded into TCR matrixes for the downstream analysis.

The trained M^CAT, as well as the trained M^CNN, was applied to the processed test data, and the cancer score of each sample was defined by averaging all the input sequence scores in the sample. The sensitivities (at 0.9 specificities) and AUCs of both models were estimated for evaluating the models (Table 4). The results showed that M^CNN performed significantly better than M^CAT in terms of both the sensitivities and AUCs on each group of samples except for the melanoma, ovarian cancer, and colorectal cancer samples in T₁, where M^CNN’s performance was close to that of M^CAT.

To further compare the performance of these two models, M^CAT and M^CNN were also applied to the combined cancer samples of T₁ and T₂, which contained all various cancer samples as well as the control samples, and the ROC curves were generated based on the prediction results (Figure 2). As shown in Figure 2, the AUCs of M^CNN were higher than those of M^CAT on both T₁ and T₂ (T1: 0.83 > 0.78; T2: 0.73 > 0.71), indicating that it had the better feature extraction and prediction ability for TCRs.

In comparison with M^CAT and M^CNN, we expected that the entire DeepLION (M^LION) as a MIL method was capable of achieving more accurate caTCR prediction after correctly modeling the correlations among TCRs. Because the logistic regression model combined with MIL (denoted as M^LOG) (Ostmeyer et al., 2019) is a classical MIL method to distinguish tumor tissues from healthy tissues accurately by identifying the cancer-specific motifs in TCRs, we also compared M^LION with it.

We trained two models for THCA and lung cancer (denoted as M T LION and M L LION ) respectively (Table 3). To obtain M T LION , we first extracted the top 100 most abundant sequences (k = 100) from each sample and encoded them into l × 15 TCR matrixes (d = 15) by the Beshnova matrix (Beshnova et al., 2020). Then, we trained the model five times independently, in each of which the model was trained with the training set for 700 epochs at a learning rate of 0.001 with the assurance that the model reached convergence. The AUCs of the trained models on the validation set were estimated for model selection, and the model with the highest AUC (0.95) was selected as the final model. M L LION (AUC = 0.86) was obtained by the same training procedure. After using the same training sets as M^LION to train M^LOG with default training parameters, we also obtained M T LOG and M L LOG for THCA and lung cancer respectively.

Four trained models were applied to the corresponding test sets and the accuracies, sensitivities, specificities, and AUCs were calculated based on the sample predictions to evaluate their performances (Table 5). The results of M^CAT and M^CNN on the samples were also shown in Table 5 for comparison. For a fair comparison, the classification thresholds of M^CAT and M^CNN for the samples (M^CAT: 0.336 for THCA and 0.321 for lung cancer; M^CNN: 0.433 for THCA and 0.419 for lung cancer) were determined by the Youden index (Fluss et al., 2005) enabling the selection of an optimal threshold value for classification, whereas the classification threshold of both M^LOG and M^LION for any sample is a fixed value (0.5). The ROC curves were generated based on the predicted probabilities of all models on the two samples (Figure 3). The results demonstrated that despite the specificity of M^LION being slightly lower than M^CNN’s on THCA sample, its accuracies, sensitivities, and AUCs were significantly better than the other models on both THCA and lung cancer samples, which indicated that M^LION can accurately predict the samples.