Development and Validation of a Bronchoalveolar Lavage Genomic Classifier for Acute Cellular Rejection
Description
Background: Acute cellular rejection (ACR) is the main risk factor for chronic lung allograft dysfunction (CLAD), but diagnosis requires invasive transbronchial biopsy (TBB). We previously demonstrated the feasibility of bronchoalveolar lavage cell pellet (BAL-cp) gene expression for ACR diagnosis. We sought to develop and validate a genomic classifier for ACR in a multicenter cohort. Methods: We performed RNA-seq on 806 BAL-cp from 181 lung transplant recipients enrolled in CTOT-20. Differential expression was based on fold difference >2.0 and False Discovery adjusted p-value <0.05. Samples were randomly split 80:20 into training and testing sets. A Random Forest model was optimized for area under the curve (AUC), and the threshold for genomic ACR was selected for classification accuracy. We validated performance in an independent single-center cohort. Cox models evaluated risk for CLAD. Findings: From 37 cases and 151 controls, we identified 62 ACR genes, indicating upregulation of T-cell receptor signaling, and downregulation of CTLA4 signaling in cytotoxic lymphocytes, among other enriched pathways. A 31-gene Random Forest model’s AUC was 0.99 (SE 0.0053) in the training set, and 0.72 (SE 0.0874) in the test set. At a probability threshold of 0.396, accuracy for distinguishing clinically significant ACR cases from stable controls was 93.1% (specificity 95.4%, sensitivity 83.8%). In the independent validation cohort, accuracy was 82.1% (specificity 87.5%, sensitivity 73.3%). The model classified 138 (17.1%) CTOT-20 samples as genomic ACR. Late genomic ACR (≥90 days posttransplant) associated with increased CLAD risk (HR 2.52, 95% CI 1.47 - 4.34, p<0.001). Interpretation: A BAL-cp genomic classifier can identify ACR, predict CLAD risk, and may be a less invasive alternative to TBB after lung transplant.
Files
Steps to reproduce
Differentially expressed genes were detected using LIMMA and voom packages, We only included genes with >2 mean reads per million mapped reads. LIMMA was used in conjunction with voom, which weighs the mean-variance relationship of the log-counts needed for accurate generalized linear modeling. Empirical Bayes estimator was adopted, with adjustment for age, self-reported sex, days to sample collection (log2 transformed), transplant type (bilateral or single), primary graft dysfunction grade 3 at 72 hours, and use of induction. Correlation between repeated measures were considered and adjusted in differential expression analysis. Candidate genes were identified based on an absolute fold change >2.0 and Benjamini−Hochberg adjusted p-value less than 0.05. Models were trained to distinguish csACR from SC. Differentially expressed genes were included as features for classifier modeling. For modeling, we treated repeated measures as independent observations. Several machine learning methods were evaluated, including simple logistic regression, Elastic-Net with 10-fold Cross-Validation, and Random Forest (RF) with 10-fold Cross-Validation. We selected a RF model as they are generally robust to overfitting and capable of learning non-linear relationships. We tuned the number of variables available for splitting at each tree node (mtry range 1 to 15) and the number of trees (ntree range 100 to 5000) optimized for area under the curve (AUC). We used the Boruta variable selection method to reduce the list of genes and derive a final model. The final locked model threshold was selected based on optimal accuracy for distinguishing csACR from SC in the training set and validated in the test set. Classification performance was reported as sensitivity, specificity, and accuracy, as well as negative predictive value (NPV) and positive predictive value (PPV) where appropriate. We applied the classifier to the entire dataset and describe the proportion of non-csACR phenotypes classified as genomic ACR (gACR). We used Cox models to determine risk factors for CLAD. Time was calculated from day 90 posttransplant to either CLAD onset or censor date (last PFT). Time-independent covariates included baseline characteristics and “early” events occurring up to 90 days post-transplant, including gACR. “Late” events occurring after 90 days was considered time-dependent covariates. Class I and class II DSA were included in models only if occurring >90 days (not early DSA), as these events more reflect humoral responses that persist after transplant. Analyses were conducted using the Bioconductor suite of packages in the R statistical software environment version 4.2.1. Packages “randomForest” and “glmnet” were used for classifier modeling.
Institutions
- University of California Los AngelesCA, Los Angeles