Gene expression of Mouse Organoids was retrieved from the NCBI GEO database (https://www.ncbi.nlm.nih.gov/geo/). Six RNA-seq datasets from two studies, GSE150425 (15) and GSE180360 (11), were downloaded from the GEO website. Differential expression analysis was performed using the “edgeR” R package (version 3.40.0). Genes with adjusted P value < 0.05 and |log2FC| > 0.5 were identified as differentially expressed genes (DEGs).

TIMER (http://timer.cistrome.org/) was used to analyze the differential expression of IRAK1BP1 between normal tissues and tumor types ( Figure 1G ). [IRAK1BP1 expression levels in various tumor types and the corresponding normal tissues were analyzed by TIMER (http://timer.cistrome.org/) ( Figure 1G )]. The clinical characteristics, somatic mutations and RNA-seq data in TPM (transcripts per million reads) of lung cancer patients, as well as the corresponding normal lung tissue data in the cancer genome atlas (TCGA) (https://portal.gdc.cancer.gov/) were downloaded using TCGAbiolinks (16). The correlation between IRAK1BP1 expression and overall survival (OS) was analyzed by the “survival” R package (version 3.4-0). Patient samples were divided into four groups based on 25%, 50% and 75% quartiles ranked by IRAK1BP1 expression, and the Kaplan–Meier survival analysis for OS was performed in the top and bottom groups. Statistical significance of survival curves was assessed by a log-rank test.

We calculated the percentage of mutation for each gene across all the LUAD samples and selected the top 20 genes with the highest mutation rates. We then ranked the samples based on the expression of IRAK1BP1 from high to low, and split all cases into three groups (high, < 25%; medium, ≥ 25% and ≤ 75%; low, > 75%). A waterfall plot was displayed using the “maftools” R package (version 2.14.0).

Gene expression data of LUAD cell lines were retrieved from the cancer cell line encyclopedia (CCLE) (https://sites.broadinstitute.org/ccle/). Differentially expressed genes in LUAD samples from TCGA and LUAD cell lines from CCLE were identified with adjusted P value < 0.05 and |log2FC| > 1. Gene Ontology enrichment analyses were performed by the “clusterProfiler” R package (version 4.5.2). The “fgsea” R package (version 1.22.0) was used to perform gene set enrichment analysis (GSEA).

To search for novel biomarkers for lung cancer that predict disease outcome, we analyzed RNA-seq datasets from two studies using murine tumor organoids derived from murine alveolar epithelial progenitor (AT2) cells and/or bronchioalveolar stem cells (BASCs), which were postulated to give rise to LUAD (10). Figure 1A demonstrates the workflow of RNA-seq data collection from mouse organoids with or without the in vitro induced activation of KRAS^G12D (10) or EGFR^{T790M / L858R} (11), which are two of the most predominant oncogenic mutations in NSCLC (17); this was then followed by DEG analysis to find genes aberrantly expressed in tumor organoids by comparing the transcriptomic profiles of tumor versus normal (control) organoids. Data collected for a total of six organoid types were analyzed for our study ( Figure 1A ). Lists of DEGs for each dataset were then compared to find 9 DEGs that were common to all six organoid models of LUAD ( Figures 1A-C ). Genes that were consistently up-regulated across all LUAD organoid types include Wnt7a, Unc13c, Rap1gap, and Ank1 ( Figure 1B , upper heatmap), while genes consistently down-regulated were Gstm1, Gpx2, Acox2, Irak1bp1, and Erich2 ( Figure 1B , bottom heatmap). The log fold change (log2FC) for each DEG was displayed in Figure 1D , showing that these 9 genes were indeed consistently up- or down-regulated in all 6 organoid types that were analyzed.

To find candidate biomarkers that may be valuable to the prediction of LUAD progression, we retrieved and analyzed LUAD patient transcriptomic and prognostic data from the cancer genome atlas (TCGA) as well as 3 published tumor microarray datasets (GSE19188, GSE30219, GSE43458) to obtain our list of shared DEGs. To see how well these candidate biomarkers predict patient outcome, we measured the hazard ratio (hr) of the 9 DEGs discovered, and found that one gene, in particular, called IRAK1BP1 had a lower hr value with statistical significance, suggesting that it had greater predictive power than the rest of the DEGs ( Figure 1E ). As a result, we decided to further investigate the expression of IRAK1BP1 in patient samples from TCGA as well as 3 independent cohorts; our findings suggested that IRAK1BP1 expression in LUAD patients resembled what was observed in the organoid models, showing reduced expression in tumor tissue relative to normal tissue ( Figure 1F ).

Interestingly, the differential expression of IRAK1BP1 is also observed in other cancer types: we assessed IRAK1BP1 expression for all cancer types with available data ( Figure 1G ). The differential expression of IRAK1BP1 is again evident in LUAD and LUSC, the two main NSCLC subtypes (red box); IRAK1BP1 expression is lower in tumor relative to normal tissue in many other cancer types such as breast cancer (BRCA), kidney cancers (KICH, KIRC), and prostate adenocarcinoma (PRAD); while in some cancer types, such as cholaniocarcinoma (CHOL) and colon adenocarcinoma (COAD), the opposite is true ( Figure 1G ). The differential expression of IRAK1BP1 in various cancer types warrants further survival analysis in cancers other than lung cancer.

We further investigated whether differential IRAK1BP1 expression correlated with differential LUAD prognosis by plotting the survival curves for LUAD patients divided into two groups (top and bottom quartiles), based on the level of IRAK1BP1 expression ( Figure 2A ). The overall survival (OS) and early-stage survival for LUAD patients in the TCGA cohorts (TCGA-LUAD), as well as the survival curves calculated from patient microarray data, suggested that LUAD patients with high IRAK1BP1 expression (red) had higher survival rates hence better prognosis compared with patients with low IRAK1BP1 (blue) ( Figure 2A , respectively). To see how well IRAK1BP1 alone performs as a prognostic predictor relative to other known biomarkers of NSCLC prognosis, we measured the coefficients of correlation between IRAK1BP1 expression and well-studied indicators of lung cancer prognosis such as patient age, gender, and tumor stage at diagnosis, and found that IRAK1BP1 expression had no correlation with any of these factors; our analysis showed that IRAK1BP1 alone was able to robustly predict LUAD patient prognosis ( Figure 2B ). In addition, pair-wise correlation between each of a collective set of lung cancer prognostic marker genes and IRAK1BP1 expression was measured and shown in Figure 2D , where the color represents the coefficient of correlation and the size of the dot represents statistical significance; again, our analysis suggests that IRAK1BP1 had very weak correlation with most known biomarkers and was able to act alone as a prognostic predictor ( Figure 2D ).

To compare the performance of IRAK1BP1 as a prognostic predictor with known biomarkers of LUAD, hazard ratios (hr) were calculated for IRAK1BP1 along with previously described lung cancer biomarkers listed in Figure 2C : hr values lower than 1 indicates low hazards (i.e. better prognosis); the lower the hr value, the better the prognosis it predicts. Our analysis showed that IRAK1BP1 was able to robustly predict LUAD patient prognosis, outperforming each of the known biomarkers in the list ( Figure 2C ).

Tumor mutation burden (TMB) and mutation rates are well-known biomarkers of lung cancer patient treatment response and prognosis. To assess whether the correlation between IRAK1BP1 expression and patient survival is a function of TMB and mutation rates, we evaluated the correlation between IRAK1BP1 expression and a number of factors including TMB, the mutation status of the top 20 genes with the highest mutation rates, age, gender, OS, smoking status, and pathological stage of LUAD ( Figure 3 ). Patient samples were divided into three groups based on the level of IRAK1BP1 expression and were aligned horizontally from high to low expression ( Figure 3A ). Here, we did not observe any significant correlation between IRAK1BP1 expression and most of the factors evaluated; mutation rates were not dramatically altered for the 20 most mutated genes selected, though mutation rates in the high-expression group were slightly lower than those in low- and/or medium-expression groups ( Figure 3B ). The frequencies of different types of mutations found in LUAD patients were also fairly even between the High-IRAK1BP1 and Low-IRAK1BP1 groups ( Figure 3C ). Our results suggest that IRAK1BP1 expression correlates with LUAD prognosis and is able to act independently as a candidate predictor of disease outcome.

To find out how IRAK1BP1 may assist in tumor formation and progression, we performed gene-set enrichment analysis (GSEA) as well as gene ontology analysis for LUAD samples in the high-IRAK1BP1 group ( Figure 4 ). By comparing transcriptomes between LUAD patients with high- and low-IRAK1BP1 expression, we found genes that were up-regulated or down-regulated in patients with high IRAK1BP1 ( Figure 4A ). GSEA analysis in the high-IRAK1BP1 transcriptomes from the TCGA-LUAD patient cohorts suggested that some of the most well-known oncogenic pathways, such as the Notch, canonical Wnt, as well as the TGF-β pathways, were down-regulated in LUAD tumors from the high-IRAK1BP1 group ( Figure 4B ). Further validating the predictive power of IRAK1BP1, our GSEA analysis of high-IRAK1BP1 expressers using RNA-seq data from lung cancer cell lines in the cancer cell line encyclopedia (CCLE) database showed remarkable similarity to the GSEA results in patient data ( Figures 4C, D ). Interestingly, in addition to down-regulated oncogenic pathways, signatures of inactive KRAS pathway were also present in the high-IRAK1BP1 group ( Figure 4D ). GO term analysis of up-regulated genes ( Figure 4E ) and down-regulated genes ( Figure 4F ) in TCGA-LUAD samples with high IRAK1BP1 suggested that high IRAK1BP1 expression was associated with suppression of pathways associated with immune/inflammatory response as well as developmental/wound healing, possibly reflecting reduced TMB in this group ( Figure 4F ); in contrast, pathways associated with cytoskeleton and cell adhesion were up-regulated in the high-IRAK1BP1 group, suggesting that these cells may have better structural stability ( Figure 4E ). These findings prompt further studies to better characterize the association of IRAK1BP1 with tumorigenesis and cancer progression.