Consensus label propagation with graph convolutional networks for single-cell RNA sequencing cell type annotation

Abstract Motivation Single-cell RNA sequencing (scRNA-seq) data, annotated by cell type, is useful in a variety of downstream biological applications, such as profiling gene expression at the single-cell level. However, manually assigning these annotations with known marker genes is both time-consuming and subjective. Results We present a Graph Convolutional Network (GCN)-based approach to automate the annotation process. Our process builds upon existing labeling approaches, using state-of-the-art tools to find cells with highly confident label assignments through consensus and spreading these confident labels with a semi-supervised GCN. Using simulated data and two scRNA-seq datasets from different tissues, we show that our method improves accuracy over a simple consensus algorithm and the average of the underlying tools. We also compare our method to a nonparametric neighbor majority approach, showing comparable results. We then demonstrate that our GCN method allows for feature interpretation, identifying important genes for cell type classification. We present our completed pipeline, written in PyTorch, as an end-to-end tool for automating and interpreting the classification of scRNA-seq data. Availability and implementation Our code for conducting the experiments in this paper and using our model is available at https://github.com/lewinsohndp/scSHARP.


Introduction
Single-cell RNA sequencing (scRNA-seq) measures the amount of RNA from each gene present in an individual cell, serving as a proxy for gene expression. High-quality labels of cell type (based on the transcriptional profile produced by scRNA-seq) have proven valuable for characterizing gene expression of cells, and for discovering cell types and genetic drivers of disease. Traditionally, these labels are produced by unsupervised clustering followed by labeling clusters with known marker genes. However, unsupervised clustering is limited by issues such as the size of scRNA-seq datasets as well as subjectivity in reclustering and biological interpretation of clusters (Kiselev et al. 2019).
The limitations of traditional cell type annotation methods have necessitated the development of automated methods for cell labeling. Three main categories of tools have emerged: marker gene-based, correlation-based, and supervised classification-based (Pasquini et al. 2021). Marker-based approaches employ known marker genes for labeling, while correlation and supervised learning-based approaches require manually labeled scRNA-seq datasets with the cell types of interest. Within these broad categories, the performance of individual tools varies widely across datasets (Abdelaal et al. 2019). As a result, using the consensus of multiple classification tools could yield higher accuracy. However, there currently exist no tools for researchers to easily apply multiple classification algorithms to their scRNA-seq data. Even once different predictions have been collected, deciding how to handle cells with ambiguous predictions between tools remains an open question.
We address these issues in two ways ( Fig. 1). First, we provide a pipeline for annotation of scRNA-seq data with multiple state-of-the-art annotation algorithms. Second, we implement and test a semi-supervised Graph Convolutional Network (GCN) as a mechanism to propagate labels from confidently labeled cells to unconfidently labeled cells. We compare our method to five state-of-the-art cell classification tools, along with a simple max consensus and a nonparametric neighbor majority approach. We show our method improves overall classification accuracy (and specifically, classification accuracy on unconfidently labeled cells) compared to component tools and the simple consensus of the tools. We also demonstrate our approach shows comparable performance with the nonparametric neighbor majority approach.
Following accurate classification, an important aspect of model usage within the biological sciences is interpretability. There has been significant development around neural network interpretation, resulting in a variety of approaches for assigning attribution scores to input features. Gradient-based attribution methods for neural networks have proven valuable given their efficiency and ease of implementation (Ancona et al. 2019). Importantly, trained neural network models can be queried for important features for classification that are not otherwise obvious. One effective algorithm for gradient-based interpretation is DeepLIFT. This method decomposes the output prediction of a neural network on a specific feature by backpropagating through the network. DeepLIFT tries to explain to the difference between the output and a neutral reference output, generated by a neutral reference input (Shrikumar et al. 2017).
We demonstrate the use of DeepLIFT (Shrikumar et al. 2017) as an effective interpretation tool for our GCN model, allowing researchers insight into classification decisions and important cell type gene markers. The ability to apply DeepLIFT to our model is a valuable feature compared to other baseline methods used in this study. Specifically, we show interpretation of our trained model identifies differentially expressed genes with biological significance not explicitly used as markers. Furthermore, we show adding these highly attributed genes as markers for marker-based tools can improve accuracy.

Materials and methods
2.1 Our model 2.1.1 Selected component methods Our method relies on picking confident labels for a subset of cells in a dataset. To do this, our pipeline currently includes five different state-of-the-art annotation methods: SCINA (Zhang et al. 2019), ScType (Ianevski et al. 2022), scSorter (Guo and Li 2021), SingleR (Aran et al. 2019), and scPred (Alquicira-Hernandez et al. 2019). We briefly describe each of these tools and refer the reader to these tools' original papers for more details. ScType employs a clustering-based approach that inputs the scRNA-seq gene expression matrix along with cell type gene markers, and outputs predictions (Ianevski et al. 2022). scSorter is another clustering-based approach that uses the scRNA-seq gene expression matrix along with cell type gene markers. This method recognizes that over-expression of certain marker genes is not present in populations of many cell types and attempts to address this problem (Guo and Li 2021). SCINA is another state-of-the-art approach that uses the same information as ScType and scSorter, with the added benefit of being much faster. SCINA uses an expectation-maximization algorithm to assign labels (Zhang et al. 2019). SingleR requires both the scRNA-seq gene expression matrix and a labeled reference gene expression matrix. This method uses correlation between the reference and query sets to extend labels (Aran et al. 2019). scPred requires the same information as SingleR. This method uses feature reduction to pull out important cell type features and then a probabilistic machine learning algorithm for label prediction (Alquicira-Hernandez et al. 2019). We note that our pipeline could be applied to a broader set of tools.

Picking confident labels via consensus
We briefly describe our pipeline for generating "confident" i.e. consensus cell labels. Gene counts are initially filtered to first remove cells expressed in fewer than 200 genes and then filtered to remove genes expressed in fewer than 3 cells. Tools make predictions on further preprocessed gene counts as specified by the given tool, returning an n by t matrix R where t is the number of annotation tools and R i;j is the prediction from tool j for cell i. We then one-hot encode the predictions from each column of R resulting in t matrices of size n by the number of cell types. We then add these matrices to construct a new matrix E where the rows correspond to cells and columns correspond to the amount of votes for each cell type. Finally, we create confident labels c 2 R n where c n corresponds to the predicted cell type for cell n if a majority of tools agreed and -1 if not. Selecting confident labels entails setting a threshold of confidence which we set at 0.51 for a simple majority. For cell x i we find m i ¼ maxðE i =sumðE i ÞÞ and assign label -1 if m i < :51 and the label of the corresponding column of m i if m i >¼ :51. This results in c 2 R n containing the confidently labeled and unconfidently labeled cells.  Figure 1 Starting from an scRNA-seq dataset, a user can apply any number of tools which classify cell types (e.g. classification from a prelabeled reference dataset or from reference genetic markers). If a majority of these underlying tools assign the same label to a cell, we say the ensemble is confident in this label. Our GCN learns to propagate labels from confidently labeled cells to the rest of the cells ("unconfident cells") via message passing in a K-nearest neighbor graph.

Semi-supervised GCN
Our data are presented to our model as a preprocessed scRNA-seq matrix for n cells, denoted by X ¼ fx 1 ; . . . ; x n g R F where F is the number of features for each cell. We use Principal Components Analysis (PCA) to project the cell-bygene matrix X into a 500D vector space. To train our GCN, we draw random minibatches of size b at a time. We construct G as the k-nearest neighbor (k-NN) graph of the cells included in each minibatch with k neighbors, including selfloops. We experimented with different values of k, finding lower neighbors and batches performed better. We input G into the following GCN architecture. EdgeConv (Wang et al. 2019) is a graph convolution operator developed for prediction from point clouds that captures local structure in graphs while remaining permutation invariant. We construct a GCN with l EdgeConv layers with SiLU (Elfwing et al. 2018) activation function and summation aggregation in PyTorch (Paszke et al. 2019). The SiLU activation function is a common neural network activation function that we found to work well for our GCN architecture. A final linear layer projects node embeddings into label space (whose dimension is the number of cell types in our dataset). For architecture details see Supplementary Appendix SA. We train our GCN for 150 epochs, with the Adam optimizer (Kingma and Ba 2014) at an initial learning rate of 0.0001. Our training loss is Cross-Entropy loss on the set of confidently labeled cells. Nearest neighbors are generated separately for each batch of each epoch. We consider multiple choices for each of the hyperparameters batch size b, neighbors k, number of message passing steps l, and final embedding layer size e, selecting the hyperparameters that result in best performance on held-out validation data. See Experiment Settings for details.

Interpretation with DeepLIFT
We employ DeepLIFT (Shrikumar et al. 2017) with the Rescale rule as implemented by Captum (Kokhlikyan et al. 2020), a model interpretability and understanding library for PyTorch. The Rescale rule is one option for assigning the contribution scores of the immediate inputs to each neuron. We use the same hyperparameters batch size b, neighbors k, number of message passing steps l, and final embedding layer size e as used during training. Please refer to the Experiment Settings section for how these hyperparameters are chosen. DeepLIFT requires a reference input to determine the reference activation of a neuron. We use a reference count matrix of all zeros and the same size as our scRNA-seq count matrix X. We then use our tool's predictions for each cell type as its target neuron for calculating the gradient. This results in attribution matrix A with the same size as X. We use A to calculate the mean feature attribution scores for each gene, for all cells of each type. We then take the absolute value of these attributions to compute a general importance score per feature/gene. Note that our transformed DeepLIFT scores do not indicate negative or positive importance for each feature, just overall importance. Finally, we scale each cell type to have unit variance to allow comparison of gene importance scores across cell type. Calculating these attribution scores is only possible for a differentiable model, like our GCN-not possible for any of the underlying tools. We note here that the pipeline we analyzed with DeepLIFT included the PCA step, resulting in attribution scores per gene for each cell.

Nonparametric neighbor majority approach
In addition to the max consensus baseline described below in Experiment Settings, we implemented a nonparametric neighbor majority approach. This is similar to the message passing step in our GCN model, with the difference that this method does not use a neural network to encode/decode messages. This method operates on the 500D vectors produced as the principal components of the preprocessed gene expression matrices for each dataset. We first construct a k nearest neighbors graph in Euclidean distance, using a batch size of 1000. During each round of message passing, the unconfidently labeled nodes' labels are updated as the majority label of their k nearest neighbors in Euclidean distance (only considering those neighbors who have been labeled thus far). We test three strategies during hyperparameter optimization for convergence of this method: • one round of label propagation; • iterating until <5% of labels change between epochs; and • iterating until all cells are labeled or 50 epochs have gone by.
We show accuracy scores across a variety of k values and convergence strategies in Supplementary Appendix SD.

Preprocessing data
The initial input to our pipeline is a scRNA-seq count matrix X where X ng corresponds to observed gene counts for gene g in cell n. Cells expressing < 200 genes and genes expressed in < 3 cells are removed, and X is row-normalized according to x i;: ¼ logð1 þ ððx i;: Ã 10000Þ= P j x ij ÞÞ. Both of these steps are common scRNA-seq preprocessing steps (He et al. 2020, Collin et al. 2021). Finally, we use Principal Components Analysis (PCA) to project X down to 500 features per cell.

Prototyping data
Separate simulated data that were generated with a different seed than the test simulated data was used for developing our model. We also used a subset of the PBMC data that was not included in the final test set to develop our model. We approached this task as a transductive semi-supervised graph learning task where each scRNA-seq dataset is its own graph. However, no ground truth labels are provided for the learning task. We did not run our model on our four test datasets until architecture details were finalized. In addition, only confidently labeled nodes are used for hyperparameter optimization in these test sets, meaning no ground truth labels are involved. It is also important to note no trained model weights carried over from the prototyping datasets to the testing datasets.

Simulated data
We generated our simulated datasets with Splatter (Zappia et al. 2017) and parameters estimated from 4000 Pan T Cells from a healthy donor (Available Here). Each simulated dataset contains 1000 cells, evenly split between four cell types with different transcriptomic profiles. To demonstrate that our pipeline can also incorporate reference-based tools (like SingleR and scPred), we also generated 1000 reference cells for each dataset with the same gene profiles (separated from the original data with a batch.facScale of 0.5 to simulate batch effects). Simulated datasets vary by the de.facScale Consensus label propagation with GCNs for scRNA-seq cell type annotation 3 parameter which determines the magnitude of differential expression in differentially expressed genes between cell type groups. Five markers were selected randomly from the top ten differentially expressed genes from each cell type. The 0.7 de.facScale and 0.8 de.facScale simulated datasets had 156 and 85 unconfidently labeled cells respectively.

Real data
It is not usually feasible to acquire ground truth labels for scRNA-seq data. An alternative gold standard is Fluorescenceactivated Cell Sorting (FACS), which pre-sorts cells by markers prior to conduction of scRNA-seq (Tung et al. 2007). We test our model on two FACS-labeled datasets. First, we use a scRNA-seq dataset generated from mouse testis cells (Jung et al. 2019). This data contains three cell types: 292 Spermatogonia, 244 Spermatocytes, and 156 Spermatids after filtering. Cell type markers were selected from relevant literature (Green et al. 2018) and a full list of markers is included in Supplementary Appendix SI. In addition, no reference dataset was used for this data. There were 46 unconfidently labeled cells after prediction tool voting on the dataset.
Second, we use an scRNA-seq dataset generated from Human peripheral blood mononuclear cells (PBMCs) (Zheng et al. 2017). This data contains ten cell types, however, we removed cell types not purely sorted by FACS, combined CD4þ T cells, and combined CD8þ T Cells. This resulted in five cell types: 9106 B Cells, 2341 Monocytes, 7572 Natural killer (NK) Cells, 38 006 CD4þ T Cells, and 19 856 CD8þ T Cells. We used the same markers as scSorter (Guo and Li 2021) with the addition of CD8A, CD8B, and CD4 to distinguish CD4 and CD8 T Cells. A full list of markers is included in Supplementary Appendix SI. We used the 10X PBMC 3k dataset as a reference as in https://satijalab.org/seurat/articles/ pbmc3k\_tutorial.html. This dataset contained 2310 unconfidently labeled cells.

Accuracy on test sets
3.1.1 Experiment settings For the simulated and testis datasets, 20% of confidently labeled cells were masked and held out as a validation set. We performed a hyperparameter optimization search (see Supplementary Appendix SA for details) for options of batch size b, neighbors k, layers l, and embedding layer size e, selecting the GCN architecture with the highest validation accuracy. For the PBMC dataset, 5% of the dataset was used for this same hyperparameter optimization process as using the entire dataset was computationally intractable. Each GCN used EdgeConv feature propagation between each node and its k closest neighbors, with distance determined dynamically between node features (including the PCA features at the first layer). Five random initializations of the optimal model were then trained for 150 epochs as described above, and mean accuracy was recorded. We use max consensus, the nonparametric neighbor majority approach, and other tool accuracies as baselines for our model. Max consensus simply chooses the cell type with the most votes. In the event of a tie, this method returns "unknown." The k hyperparameter and convergence method for the nonparametric neighbor majority approach was chosen with the same validation set used for GCN hyperparameter optimization. We report total accuracy and unconfident cell accuracy for each dataset. "Unconfident cell accuracy" refers to the accuracy on only the cells where the underlying tools did not find consensus.

Simulated datasets
For both simulated datasets, the optimal model has batch size 20, 2 nearest neighbors, and 2 EdgeConv layers. For the simulation with 0.7 de.facScale, 25-dimensional embedding space was optimal, whereas for the simulation with 0.8 de.facScale, the optimal value was 40. In addition, 15.6% and 8.5% of cells were unconfidently labeled in the 0.7 de.facScale and 0.8 de.facScale datasets respectively. Table 1 shows our model outperforms all other methods for total accuracy and slightly under performs SingleR for unconfident cell accuracy on the 0.7 de.facScale data.
The nonparametric neighbor majority approach proves extremely ineffective in propagating confident labels to unconfidently labeled cells in both simulated datasets. This could be due to the transcriptional similarity of all four cell types in the simulated datasets. We compare the average distance between each cell type cluster in Supplementary Appendix SH, finding all cell types in each simulated dataset are very transcriptionally similar. In contrast, the testis and PBMC datasets show varying levels of transcriptional dissimilarity depending on which cell types are compared. This suggests that in datasets with many similar cell types, the nonparametric neighbor majority approach would fail.

Testis dataset
Only marker-based prediction tools were used for this dataset as no labeled reference was easily available, and we wanted to test our method without reference-based tools. The optimal model for this dataset used batch size 20, 2 nearest neighbors, 2 EdgeConv layers, and embedding layer size of 25. Importantly, 6.6% of cells were unconfidently labeled by the component tools. Table 1 shows accuracy results, demonstrating our GCN model outperforms all other methods for both total and unconfident cell accuracy, except for the nonparametric approach. The accuracy of the nonparametric approach is within the standard deviation of our model for both the total and unconfident accuracies.
These results suggest both our GCN model and the nonparametric neighbor majority effectively propagate consensus labels to unconfidently labeled cells in this testis scRNA-seq data. This is clear given the improved accuracy over all component methods and significant improvement of our model over the tool average accuracy by 14.1%. Our model therefore demonstrates improved accuracy over component methods and comparable accuracy to the nonparametric neighbor majority approach, with the major advantage of being interpretable.

PBMC dataset
The optimal model for this dataset used batch size 50, 2 nearest neighbors, 2 EdgeConv layers, and embedding layer size of 25. 3.0% of cells were unconfidently labeled by the component tools. Table 1 shows accuracy results. For accuracy on unconfident cells, the GCN model places third behind ScType and the nonparametric approach. For overall accuracy, our model ties with the nonparametric approach to the third significant digit. In addition, the accuracy of the nonparametric approach is well within the standard deviation of our model accuracy on the unconfident set.
These accuracy results highlight the ability of our GCN model and the nonparametric neighbor majority approach to propagate confident labels in the PBMC data. We again show a clear improvement over all component methods, demonstrating the value of this approach in labeling scRNA-seq cells. Figure 2A shows the five most highly attributed genes by DeepLIFT for each cell type in the testis dataset. Interestingly, all of these top genes have uniquely high attribution in their important cell type. The highly attributed genes for a cell type also have relatively high gene expression in that cell type as indicated by Fig. 2B. We also observe high expression of Spermatocyte genes in Spermatid cells. It is important to note DeepLIFT interpretation indicates importance of both input marker genes and genes not explicitly given to marker-based tools. For example, DeepLIFT indicates genes Ncl, Ldhc, and Tnp1 that are differentially expressed in those cell types, but not explicitly included as marker genes. Tsga8 and Acrv1 are both Spermatid marker genes that appear in the top ten highest attributed DeepLIFT genes. In addition, Spermatogonia marker gene, Dazl, is both a used marker gene and highly attributed DeepLIFT gene. A full list of marker genes used and a list of the top ten highest attributed DeepLIFT genes for each cell type are included in Supplementary Appendices SI and SG, respectively. Figure 2C shows the DeepLIFT attribution scores for the top five highly attributed genes by cell type for the PBMC data. Figure 2D shows actual gene expression of these highly attributed genes in the dataset. For B Cells, Monocytes, and NK Cells we see a clear connection between the genes picked out as important by DeepLIFT and the genes expressed by those cell types. However, for CD4 and CD8 T Cells, the expression is not clearly higher for all genes. Importantly, we do observe CD8B as the most important gene for CD8 T Cell classification, a key marker for the cell type. We also observe CD3E (another important marker for all T Cells) as an important gene for both subtypes of T Cells. However, we do not observe CD4 as highly attributed in CD4 T Cells. This result is not totally unexpected given that CD4 is only measured by scRNA-seq in a relatively low percent of CD4 T Cells (Wang et al. 2022). In addition, the DeepLIFT scores for CD4 and CD8 T Cells may be less informative because the GCN often misclassifies CD8 T Cells as CD4 T Cells. This is likely due to the large heterogeneity and general transcriptional similarity of T Cell subtypes. Finally, DeepLIFT again reveals both used marker genes and novel genes. Specifically, DeepLIFT highly attributes CD74, FTH1, GNLY, and more genes that are differentially expressed, but not originally included as marker genes. DeepLIFT also highly attributes original marker genes such as CD79A, CD8B, CD3E, and CD3D. Again, a full list of PBMC marker genes is included in Supplementary Appendix SI.

Feature interpretation
We also compared highly attributed DeepLIFT genes in both datasets to genes identified in differential expression analysis. This experiment and the top ten highest ranked genes for each method can be found in Supplementary Appendix SG. In both datasets, both DeepLIFT and differential expression analysis revealed many common genes. However, there are also many genes not jointly identified by both methods. We suggest the use of DeepLIFT as a complimentary analysis to differential expression analysis, providing another method for identifying candidate cell type marker genes.
Overall, DeepLIFT interpretation proves effective in both the testis and PBMC datasets, revealing novel genes not originally included as markers. Importantly, the GCN is the only one of our tested methods that can be interpreted using DeepLIFT. This adds a valuable feature for researchers to both further understand classification decisions and get insight into novel cell type markers.

Highly attributed genes improve marker-based methods
We added the top three highly attributed genes from DeepLIFT for each cell type not already included as a marker gene to the list of marker genes. We then re-ran the markerbased tools ScType, scSorter, and SCINA with these markers on the testis and PBMC datasets to compare tool accuracy. This experiment shows the value of highly attributed DeepLIFT genes in assisting future classification. Table 2 shows the change in accuracy before and after adding highly attributed genes in the testis dataset. Both scSorter and SCINA show significant improvement in accuracy with the addition of highly attributed DeepLIFT genes. scSorter shows large in improvement in Spermatocytes while SCINA shows extreme improvement in Spermatids, demonstrating the effectiveness of using DeepLIFT attributed genes across different cell types and tools. ScType accuracy remains the same across all cell types. Although ScType does not improve, it is important to note accuracy does not go down with the addition of these new marker genes.   FTH1  TYROBP  S100A9  S100A8  FTL  GNLY  NKG7  TYROBP  RPS21  GZMB  CD8B  CD3E  CTSW  TPT1  EEF1A1  LTB  FTH1  CD3E MT.CO2 RPL34  maintained improvement in B Cells, Monocytes, and NK Cells, indicating the effective use of DeepLIFT attributed genes as markers in these cell types. However, both T Cell subtypes show varying performance, with decreases in accuracy in both scSorter and SCINA. We believe this is due to the inaccuracy of our model in distinguishing these T Cell subtypes as previously discussed. This suggests using DeepLIFT attributed genes falls short in cell types with poor model accuracy. We suggest further validation of all genes with direct expression of genes in the scRNA-seq data to confirm new markers.

Discussion
scRNA-seq cell type annotation remains an open and essential problem in bioinformatics. Accurate cell type labels are critical for downstream analysis and robust hypothesis testing with single-cell data. In this work, we propose a novel framework for scRNA-seq cell type annotation that improves upon existing approaches. Building upon existing annotation tools, we implement an EdgeConv-based GCN model to propagate consensus-based confident labels to the remaining unlabeled cells. We show an improvement in accuracy over a baseline max consensus algorithm and the average tool accuracy. We also show similar, if not improved, performance over a nonparametric neighbor majority propagation approach. We then demonstrate the ability to identify important genes for classification via model interpretation with DeepLIFT. The model interpretation is especially valuable for researchers as it has the potential to uncover novel gene markers and provide insight into the model's decisions.
We provide evidence for the accuracy of our GCN approach through comparison to state-of-the-art methods ScType, scSorter, SCINA, SingleR, and scPred on a variety of datasets, both simulated and real, including both Mouse and Human tissues. Our method improves accuracy in all datasets over all of these methods, providing strong evidence for the value of our method. A key takeaway from our research is the value of using a consensus of multiple annotation tools to improve general classification accuracy in scRNA-seq data.
We further provide evidence for the interpretability of our trained model with DeepLIFT. Highly attributed DeepLIFT genes include both given markers and novel genes that are differentially expressed in the actual data. This interpretation method is only possible for our method out of all compared methods, and provides a valuable tool for researchers when analyzing their scRNA-seq data.
One interesting aspect of our research is the similarity in classification accuracy of the nonparametric neighbor majority approach and our GCN approach for both real datasets. This suggests for these two datasets, the majority of valuable information for classification is simply in the existing labels of each cell's neighbors. However, the GCN approach does improve performance when multiple cell types are transcriptionally similar (as in the simulated datasets). This suggests the value of using our GCN approach emerges when transcriptional profiles are similar across multiple cell types. These findings are relevant to both this project and the general function of message passing in graphs. As a whole, we find the most value in our GCN approach given its consistency in classification accuracy across all tested datasets and interpretability.
Overall, our method will serve as a valuable tool for researchers. They will be able to easily apply five state-of-theart labeling tools and obtain accurate labels with our GCN propagation. Importantly, our pipeline accommodates any set of predictions, allowing our tool to be applied to any cell type with or without a reference dataset. a Accuracy difference is recorded for the following cell types: B Cells, Monocytes, Natural Killer Cells (NK Cells), CD4 T Cells (CD4), and CD8 T Cells (CD8).