TCellR2Vec: efficient feature selection for TCR sequences for cancer classification

Features of TCR protein sequences

The amino acid sequence of a TCR is highly variable, reflecting the diversity of the TCR repertoire that is necessary for effective immune function. Given a TCR (protein) sequence as input, various features can provide insights into its structure, function, and diversity. These features include (i) amino acid sequence composition, (ii) the length of the CDR3 region, (iii) hydrophobicity, (iv) charge distribution, (v) average motif similarity score, (vi) Shannon entropy, and (vii) Simpson index. We will now discuss those features one by one.

Amino acid sequence composition

Amino acid sequence composition in the shape of a numerical vector refers to the frequencies or proportions of the different amino acids in a protein sequence. The numerical vector has a length of 20 because 20 common amino acids can be found in protein sequences. In the context of TCR sequencing data, amino acid sequence composition can be used to characterize the diversity and functional properties of the TCR repertoire (Izraelson et al., 2018). Moreover, the amino acid sequence composition can provide information about the functional properties of the TCR repertoire. For example, certain amino acids may be enriched in TCR sequences that are specific for tumor-associated antigens or other disease-associated antigens (Pan & Li, 2022).

Length of the CDR3 region

The complementarity-determining region 3 (CDR3) is highly variable in the antigen receptor genes of B-cells and T-cells, which is critical for the specificity of antigen recognition. The length of the CDR3 region is typically measured in terms of the number of amino acids, forming a scalar integer.

The length of the CDR3 region can be associated with the affinity and specificity of TCR recognition (Glanville et al., 2017).

Hydrophobicity

In TCR biology, hydrophobicity refers to the tendency of certain amino acid residues in the TCR to interact with hydrophobic residues in the major histocompatibility complex (MHC) molecule and/or the peptide antigen that is presented by the MHC(Chowell et al., 2015). It forms a scalar value represented by a float. Hydrophobic interactions play an essential role in the binding of the TCR to the MHC-peptide complex, and the strength of these interactions can affect the affinity and specificity of the TCR for the antigen (Maruyama et al., 1993).

Charge distribution

Charge distribution in TCR refers to the distribution of charged amino acid residues(i.e., positively charged lysine and arginine and negatively charged aspartic acid and glutamic acid) in the complementarity-determining regions (CDRs) of the TCR (Robbins et al., 2008). The charge distribution in the TCR can affect the binding affinity and specificity of the receptor for the antigen (Davis et al., 1998). Therefore, the charge distribution in the TCR CDRs can provide information about the antigen-binding properties of the TCR repertoire.

Average motif similarity score

The average motif similarity score in TCR shows as a scalar float and refers to the average similarity score between the TCR CDR amino acid sequences and a set of predefined motifs or patterns that are known to be associated with certain antigens or disease (Gupta et al., 2007). The average motif similarity score in TCR can provide information about the antigen specificity of the TCR repertoire. TCRs that have high similarity scores with known antigen motifs may be more likely to recognize and respond to the corresponding antigen (Wadie et al., 2022).

Shannon entropy

Shannon entropy is a measure of diversity or uncertainty in a set of data (Shannon, 1948). In the context of TCR sequencing data, Shannon entropy can be used to determine the frequency of each amino acid at each position in the sequence and then use that information to compute the entropy. A high Shannon entropy implies that multiple different amino acids are observed at a particular position in the TCR protein sequence. This can be an indication of functional flexibility or adaptability, as different amino acids at that position may confer different properties or functions to the TCR protein (Krishna et al., 2020).

Simpson index

The Simpson index is a scalar float value between 0 and 1 and is a measure of diversity commonly used in ecology to quantify the evenness of species abundance in a community (Simpson, 1949). In the context of a single TCR (T-cell receptor) protein sequence, the Simpson index can be used to assess the clonality or dominance of specific amino acids in the sequence. A Simpson index closer to 0 indicates higher diversity or richness (Leinster & Cobbold, 2012), meaning that the amino acids in the sequence are evenly distributed or relatively equal in abundance. Conversely, a Simpson index closer to 1 suggests lower diversity or richness, indicating that a few amino acids are dominant or highly abundant in the sequence (Choudhury et al., 2016).

Shannon entropy and Simpson index are measures of diversity or uncertainty within a sequence or a set of sequences. While they are commonly used to assess diversity across a repertoire or a population of sequences, they can also be computed for individual sequences. It’s important to note that when considering a repertoire or a set of sequences, Shannon entropy and Simpson index can provide a broader picture of diversity within the entire population. However, computing these measures for individual sequences allows for a more detailed analysis at the sequence level.

TCellR2Vec generation

The seven features (discussed above) are used to generate our TCellR2Vec embedding. The dimensionality of all features in our approach is equal to one, except for amino acid sequence compositions, which is 20. Therefore, for each TCR sequence, we generate a numerical feature vector with a length of 26. It’s important to note that each of these features is calculated individually for each TCR sequence, rather than being derived from the entire TCR protein sequence. Algorithm 1 presents the pseudocode for computing TCellR2Vec features and its final embedding. The input for this algorithm is a set of sequences S = {s₁, s₂, …, s_n} where S is the set of sequences, n is the number of sequences, and s₁ to s_n represent the first sequence to the n^th sequence in the sequence set. One sequence is shown as an example in Fig. 2A. The algorithm operates by iterating over each sequence in S and calculating features such as CDR3, aa_comp, hydrophobicity, charge, similarity, Shanon, and Simpson as shown in lines 4 to 10, also shown in Fig. 2B. The first feature, CDR3, is computed by counting the number of unique amino acids in each sequence. The second feature, aa_comp, is generated by analyzing the sequence of amino acids present in each protein sequence and counting the frequency of occurrence of each amino acid, line 5 of Algorithm 1 . For example, if a protein sequence contains 15 amino acids and the amino acid alanine appears four times in the sequence, then the frequency of alanine in that sequence would be 0.26. Similarly, the frequency of occurrence of all the other amino acids is also calculated to generate the aa_comp feature vector. We followed the same procedure to calculate the frequency distribution of each amino acid in the sequence “CASSATGNEQFF”, as depicted in the Fig. 2 and the result is the following: [0, 0.166, 0, 0, 0.083, 0, 0.166, 0, 0.083, 0.083, 0, 0, 0.083, 0.166, 0, 0.083, 0.083, 0, 0, 0].

The similarity score feature is calculated by using the BLOcks SUbstitution Matrix 62 (Blosum62) which is a substitution matrix used in bioinformatics for sequence alignment of protein (Henikoff & Henikoff, 1992). It assigns a score to each possible substitution of one amino acid with another, based on the observed frequency of that substitution in the set of aligned sequences. The value calculated for the specific sequence in Fig. 2B is 1. Finally, to calculate hydrophobicity, charge, Shanon, and Simpson features we used specific libraries. For example, to calculate the hydrophobicity, we used the Kyte-Doolittle hydrophobicity scale, which assigns a score to each amino acid based on its hydrophobicity (Kyte & Doolittle, 1982). The hydrophobicity value calculated for the sequence “CASSATGNEQFF” in Fig. 2 is −0.125. The charge distribution was calculated by counting the number of positively and negatively charged amino acids in each sequence. The charge distribution value computed for the sequence “CASSATGNEQFF” in Fig. 2 is −1.465.

For the Shannon feature, we used the Shannon entropy algorithm, which measures the diversity of amino acids in each sequence. This algorithm takes into account both the frequency of each amino acid and the number of different amino acids present in each sequence. The Simpson diversity feature was calculated using the Simpson index, which measures the probability that two amino acids selected at random from a sequence will be different. This index ranges from 0 to 1, with 1 indicating the highest diversity and 0 indicating no diversity. Figure 2 displays the calculated values for Shannon entropy (3.084) and Simpson diversity (0.875). After calculating these features for each sequence in the input set, the algorithm concatenates all the calculated features to obtain the final embedding vector which is shown in Algorithm 1 , line 11 and Fig. 2C. The resulting embeddings can be merged with baseline embedding generator methods and utilized as an input for classification methods in supervised analysis.

In Algorithm 1 , the loop runs for each sequence in |S|, which has n sequences. Thus, the complexity is O(n). Within each sequence of average length m, finding the CDR3 length and calculating amino acid composition, hydrophobicity/charge scores, charge distribution, and the CDR3 motif all take O(m) time. Since these operations are performed for all n sequences, their overall complexity becomes O(m * n). Loading the Blosum62 matrix has a constant time complexity of O(1) as it’s independent of sequence size. Calculating the similarity score for the motif involves iterating over its length (k), which is much smaller than m, resulting in O(k) complexity. Finally, calculating Shanon Entropy and Simpson Index based on sequence length have a complexity of O(m). Overall, the code’s complexity is dominated by O(m * n) operations due to its dependence on both sequence length and the number of sequences.

 
_______________________ 
Algorithm 1 TCellR2Vec____________________________________________________________________ 
     Input:  TCR Sequences S 
    Output:  TCellR2Vec Embedding ϕ 
 1:  ϕ ← []                                                ⊳ Initialize Embedding 2D array 
  2:  for i ← 0 To |S| do 
 3:       seq ← S[i] 
  4:       CDR3 ← CDR3LENGHT(seq) 
  5:       aa_comp ← AMINOACIDCOMPOSITIONS(seq) 
  6:       hydrophobicity ← HYDROPHOBICITY(seq) 
  7:       charge ← CHARGEDISTRIBUTION(seq) 
  8:       similarity ← MOTIFSIMILARITY(seq) 
  9:       Shannon ← SHANNONINDEX(seq) 
10:       Simpson ← SIMPSONINDEX(seq) 
11:       Vec ← CONCATENATE(CDR3, aa_comp, hydrophobicity, charge, similar- 
     ity, Shannon, Simpson) 
12:       ϕ.append(Vec) 
13:  end for 
14:  return ϕ             ⊳ Returning list of embeddings for all TCR sequences    

Concatenating TCellR2Vec with baseline methods

After generating TCellR2Vec embedding using the features extracted from the protein sequences, we aggregate (concat) those embeddings with recent embedding methods from the literature to enhance their performance. To this end, the embedding methods used from the literature are (i) one-hot encoding, (ii) Spike2Vec, (iii) PWM2Vec, (iv) spaced K-mers, (v) autoencoder, and (vi) WDGRL. We will now discuss these methods in more detail and the process of their concatenation with TCellR2Vec using a TCR protein sequence associated with lung cancer as an example.

One-hot encoding (OHE) and TCellR2Vec

OHE is a method for generating numerical embeddings of sequences by creating binary vectors for each character. Each binary vector has a size equal to the number of possible characters and assigns a value of 1 to the location corresponding to the character and 0 to all others. The resulting binary vectors are then concatenated to form the final numerical embedding of the sequence (Kuzmin et al., 2020). Figure 3 illustrates a flowchart outlining the process of generating an OHE feature embedding merged with embeddings generated through our feature selection method, TCellR2Vec. First, we generate OHE for each amino acid in the TCR sequence, Figs. 3A–3B, and then combine these vectors together to get the final numerical representation of the TCR sequence, in Fig. 3C. Eventually, we merged them with the numerical vectors of TCellR2Vec, Fig. 3D, and used them as our input of classifiers to classify cancer types.

Spike2Vec and TCellR2Vec

Spike2Vec is a method for converting bio-sequences, such as DNA or protein sequences, into numerical embeddings that can be used in machine learning-based classification tasks. It does so by counting the occurrences of K-mers, which are consecutive substrings of length K that retain the ordering information of the sequence (Ali & Patterson, 2021). We set k equal to 3 in our experiments. The flowchart in Fig. 4 provides an overview of the process of a combination of embeddings generated by Spike2Vec and TCellR2Vec methods. First, in the Spike2Vec part of the approach, the TCR protein sequence is obtained and a set of K-mer, representing subsequences of length k, is generated from the TCR sequence while preserving the order of the sequence, in Figs. 4A to 4B. Next, each K-mer is transformed into a low-dimensional vector using the Spike2Vec embedding method, which learns representations of biological sequences based on cooccurrence patterns in large-scale sequence data, Fig. 4C. Finally, we merged the TCellR2Vec features with the Spike2Vec embeddings to get a single feature vector representation of the TCR sequence, which is used as an input for machine learning models to classify the TCR as associated with cancer target label, Fig. 4D.

PWM2Vec and TCellR2Vec

PWM2Vec is a technique for obtaining numerical embeddings of biological sequences using the K-mers concept. However, unlike other methods that use K-mer frequencies, PWM2Vec assigns weights to each amino acid in a K-mer and uses these weights to generate the embeddings (Ali et al., 2022). Figure 5 presents a flowchart outlining the steps involved in generating a feature embedding using PWM2Vec combined with embeddings generated through our feature selection method, TCellR2Vec. Initially, the TCR protein sequence is used and a set of K-mers are generated and the k is 3 in our approach, Figs. 5A to 5B. Next, a position frequency matrix (PFM) which shows occurrences of amino acids in each K-mer and subsequently a position probability matrix (PPM) is calculated by dividing the frequency of a character in the column by the sum of that column, Figs. 5C to 5D. Then, to avoid 0 values 0.1 is added to all values of PPM, and by using amino acids frequency tables a position weight matrix (PWM) is deliberated by taking the log-likelihood of each character in each cell of the matrix, Fig. 5E divided by its value in amino acids frequency tables, in Figs. 5E to 5G. Lastly from the PWM2Vec part, we generated the numerical representation of the TCR sequence which shows in Fig. 5H. Finally, we merged both, PWM2Vec and TCellR2Vec embeddings and used them as an input for our classifiers, in Fig. 5I.

Spaced K-mers and TCellR2Vec

In bioinformatics, generating embeddings using K-mers can be challenging due to the sparsity and high dimensionality of the resulting feature space. To address these issues, spaced K-mers have been introduced, which are non-contiguous substrings of length K called g-mers (Singh et al., 2017). The presented flowchart in Fig. 6 outlines the different stages involved in creating a feature embedding for a TCR protein sequence while using a combination of spaced K-mers and TCellR2Vec embeddings. As a first step, the TCR protein sequence is taken as input for the feature embedding generation process, Fig. 6A. the second step is generating g-mers (we considered g = 9 here) for the TCR sequence due to generating compact feature vectors with reduced sparsity and size, Fig. 6B. From those generated g-mers then we computed K-mers (K = 3) and produced the frequency of those K-mers which represent the numerical vector of spaced K-mers method, in Figs. 6C to 6D. Ultimately, we merged spaced k-mers and TCellR2Vec embeddings and used them as our classification models input, Fig. 6E.

Autoencoder and TCellR2Vec

The autoencoder approach uses a neural network to obtain numerical features from bio-sequence data. It follows an autoencoder architecture in which the encoder module is optimized to generate the embeddings. The encoder performs a non-linear transformation of the data from space X to a low dimensional numerical feature space Z (Xie, Girshick & Farhadi, 2016). In our experiments, we use a two-layered autoencoder network with an ADAM optimizer and MSE loss function to generate the embeddings. Figure 7 illustrates the process of generating embeddings of a TCR sequence, where we combine the outputs of the autoencoder method and TCellR2Vec. Initially, the TCR sequence is taken as input, and for each amino acid in the sequence, the OHE vector is computed, as shown in Figs. 7A and 7B. Subsequently, the OHE vectors are combined to generate the input for a two-layer autoencoder network that includes an encoder and a decoder to generate a low-dimensional numerical representation of the TCR sequence, as depicted in Fig. 7D. Finally, we combine the autoencoder output with TCellR2Vec embeddings to create the final embedding, which serves as input for the classification models, as shown in Fig. 7E.

Wasserstein distance guided representation learning (WDGRL) and TCellR2Vec

The WDGRL approach utilizes a neural network to extract numerical features by optimizing the Wasserstein distance (WD) between source and target distributions, making it an unsupervised domain adaptation technique (Shen et al., 2018). Figure 8 illustrates the process of embedding generation for a TCR protein sequence using WDGRL and the TCellR2Vec approach. The first step involved converting each amino acid in the sequence into separate one-hot encoded vectors, as shown in Figs. 8A and 8B. In the next step, we combined all the numerical vectors to form a single vector, which was then used as an input for the WDGRL network. This network was used to extract embeddings by optimizing the WD between the source and target features, as depicted in Fig. 8D. Finally, we combined the results of the WDGRL and TCellR2Vec embeddings to produce a final embedding, which served as the input for the classifier models, as shown in Fig. 8E.

Dataset statistics

We obtained our TCR beta chain sequence data from TCRdb, a comprehensive database for T-cell receptor sequences that offers a powerful search function (Chen et al., 2021) with more than 277 million sequences collected from over 8265 TCR-Seq samples derived from hundreds of tissues, clinical conditions, and cell types. This study focused on identifying and extracting data on the five most prevalent types of cancer, as determined by their incidence rates. To extract a subset from the original data while preserving the distribution of target labels(cancer types), we used the Stratified ShuffleSplit method and randomly extracted 50,731 TCR sequences for five different types of cancer (see Table 1). The Sequence Length Statistics in Table 1 indicate that the average sequence length is 15 for most cancer types, highlighting the challenge of working with sequences of very short lengths. Total unique TCR sequences compared with the total number of sequences for each cancer type shows the diversity of this dataset which is due to the importance of the uniqueness of the TCR sequences for the immune system’s ability to recognize and respond to a wide variety of pathogens.

Table 2 provides examples of TCR sequences, cancer names, and gene mutations for Five different cancer types. The “Gene Mutation” column presents the tumor suppressor gene mutations that are linked to an elevated risk of developing these cancers. For instance, BRAF and EGFR mutations are linked to a higher likelihood of developing lung cancer.

Data visualization

To explore the data visually, we employed t-distributed stochastic neighbor embedding (t-SNE) to map the input sequences into a 2D representation (Van der Maaten & Hinton, 2008). Figures 9A to 9F shows the t-SNE results for different embedding methods, including one-hot-encoding (OHE), Spike2Vec, PWM2Vec, spaced K-mer, autoencoder, and TCellR2Vec (our embeddings). Our observations indicate that OHE, Spike2Vec, PWM2Vec, and autoencoder exhibit a smaller number of groups for different cancer types, whereas spaced K-mer shows a more scattered pattern. Our method in Fig. 9F provides a cohesive representation, resulting in a cleaner overall structure. Additionally, we provided t-SNE plots for merging embeddings of different baseline methods with TCellR2Vec embeddings, as shown in Fig. 9G to 9K. The same pattern is repeated here, except for autoencoder, which shows a more scattered pattern. Overall, this improved representation enhances the interpretability and effectiveness of the results, facilitating a better understanding of the underlying patterns and relationships in the data.

Statistical significance

To ensure the credibility and reliability of the classification results, we performed p-value calculations for each method. These calculations were based on the average and standard deviation (the std. values are not included in the paper due to page limitation) of all evaluation metrics, obtained from five experiment runs. Notably, the p-values for all comparisons between the proposed model and the baselines were found to be less than the significance level of 0.05. However, it is important to highlight that the training runtime metric exhibited different behavior, revealing higher variability in the training runtime values, which led to certain p-values exceeding the 0.05 threshold. This discrepancy can be attributed to several factors, such as processor performance and the number of concurrent processes at any given time, which can impact the training time.

Statistical analysis

One approach to assess the efficacy of the feature embeddings is through an analysis of their compactness. To accomplish this, we conduct statistical analyses, specifically employing Pearson and Spearman correlation measures. These measures allow us to calculate the correlation values between the various features of the embeddings, which correspond to the target labels (cancer types). We then determine the proportion of attributes within each feature embedding that exhibit a strong correlation with the class labels. The Pearson correlation values for different thresholds are presented in Fig. 10, while Fig. 10 displays the Spearman correlation values for different thresholds spanning from −1 to 1, across multiple embeddings. We can observe that although WDGRL shows a higher correlation within the range < − 0.1 and >0.1, the autoEncoder and TCellR2Vec (only sequence features) show a comparable behavior. This shows that with only features (and not sequences) in the embeddings, TCellR2Vec is able to capture a correlation of features with the cancer-type labels that is comparable to the correlations of WDGRL and autoEncoder from the baseline. Compared to embedding methods like Spike2Vec and PWM2Vec, the TCellR2Vec’s embedding features are highly correlated with the class label, demonstrating that this embedding is more compact, hence preserving more valuable information.