Smoking-dependent expression alterations in nasal epithelium reveal immune impairment linked to germline variation and lung cancer risk

Study subjects

We recruited 487 subjects among which were 114 healthy volunteers from the Cambridge Bioresource (https://www.cambridgebioresource.group.cam.ac.uk/) and 373 patients referred to the outpatient clinic at Royal Papworth Hospital (Cambridge, UK) or Peterborough City Hospital (Peterborough, UK) with symptoms or imaging suspicious for lung cancer (clinic group). Healthy volunteers are here defined as individuals without any prior history or current suspicion for airway malignancy. Within the clinic group, 301 patients were diagnosed with cancer and 72 patients had a final diagnosis of a benign condition, the majority of which were inflammatory or infectious diseases (Fig. 1, Sup. Tab. 1). From these donors we collected a total of 649 samples: 413 nasal epithelial samples by mini-curette from 114 healthy donors and 299 clinic patients, and 236 bronchial brushings from clinic patients (Fig. 1, Methods). For 162 clinic patients both nasal and bronchial samples were collected (Sup. Tab. 2). Samples from healthy volunteers and clinic patients were collected and processed by the same staff using identical experimental protocols.

Smoking history was obtained for all subjects, confirmed by cotinine test, and recorded as never smokers (NV, n=45), current smokers (CS, n=153) and former smokers (FS, n=289). Former smokers were stratified into 3 categories based on their time from smoking cessation: former smokers who had quit less than one month (n=10), 1 to 12 months (n=45), or more than 1 year (n=234, median=168 months) prior to sample collection (Fig. 1, Methods). Cumulative smoke exposure was measured in pack-years, and stratified into 4 categories: none, 0-10, 11-30, >31 pack-years. In addition to smoking status, sex, age, lung cancer subtype and stage and presence of chronic obstructive pulmonary disease (COPD) were recorded according to the GOLD criteria (Vogelmeier et al. 2017) (Sup. Tab. 2). While most clinic patients with cancer were diagnosed with non-small cell lung cancer (NSCLC; n=245), 56 subjects presented with metastatic disease from an extra-thoracic primary (n=8), small-cell lung cancer (SCLC, n=31), or a rare pulmonary cancer e.g. carcinoid (n=17). Given the different underlying biology between NSCLC and other types of tumours, these subjects (with cancer status marked as Ineligible in Sup. Tab. 2) were included in all analyses investigating smoke injury response, but were excluded for lung cancer risk prediction. Clinic patients with a final diagnosis of a benign condition were followed up for a minimum of 1 year to confirm absence of cancer.

Airway samples underwent RNA sequencing using standard protocols. Blood samples were taken from 467 subjects for germline genotyping with Illumina Infinium Oncoarray pla{orm at 450K tagging germline variants. Total gene expression was quantified as variance-stabilised counts and corrected for batch effects in all downstream analyses (Methods).

Healthy volunteers and clinic patients show widespread differences in gene expression

To investigate overall gene expression patterns, we first tested for gene expression differences between all clinic patients (benign and cancer diagnoses) and healthy volunteers using nasal epithelium samples from both current and former smokers correcting for smoking status, pack-years, sex and age. We found extensive differences in gene expression between the healthy volunteer and clinic groups, with 5359 genes differentially expressed (FDR < .05, Methods). Genes showing increased expression in clinic patients were enriched for cilium assembly and organization, while genes showing reduced expression were enriched for oxidative phosphorylation and several immune-related pathways, such as neutrophil activation, antigen processing and presentation and response to interferon gamma (Sup. Tab. 3). When performing the same comparison in current smokers only, similar enrichment was found in the genes with increased and reduced expression. In former smokers who had quit for more than 1 year, there was no increased expression compared to healthy volunteers for genes related to ciliary function, but there was reduced expression of genes related to immune pathways such as inflammatory response, neutrophil activation and response to interferon gamma. These analyses demonstrate widespread expression differences betweenhealthy volunteers and clinic patients not solely aPributable to differences in smoke exposure, and suggest that an immunosuppressed state can be detected in the nasal epithelium of subjects from the clinic group during active smoking and for years after smoking cessation.

In contrast, comparing gene expression between patients with and without cancer in the clinic group and accounting for the same confounding (analysing current and former smokers together) yielded only 28 significantly altered genes (Padj < .05, Methods) in the bronchus, and no significantly differentially expressed genes in the nose. Among the 28 differentially expressed genes in the bronchus, 3 were up-regulated in patients with cancer: MMP13, a metalloproteinase known to increase lung cancer invasion and metastasis (19), EDA2R, a member of the tumour necrosis factor (TNF) receptor superfamily, members of which modulate immune response in the tumour microenvironment (20), and CTSL, a lysosomal cysteine protease involved in epithelial-mesenchymal transition (21). The 25 genes down-regulated in cancer patients were enriched in immune related GO terms, in particular neutrophil-mediated immunity (Sup. Tab. 4), consistent with our finding in the comparison between clinic patients and healthy volunteers in nasal tissue.

In summary, we observe major gene expression differences in nasal epithelium between healthy volunteers and clinic patients. However, we find no significant signal when comparing patients with cancer with those with a benign diagnosis within the clinic group. This result is in contrast to that obtained in the AEGIS study (16), which reported a notable difference in nasal gene expression between clinic-referred cancer and benign patients. However, we found a significant overlap between the set of differentially expressed genes between cancer and no-cancer in AEGIS and the set of differentially expressed genes between our clinic and healthy groups (P = 1.44×10⁻⁵). These results may be explained by subtle differences in the nature of the benign (non-cancer) diagnoses between the two studies. In our study, the majority of patients in the clinic group had clinical symptoms/imaging highly suspicious for lung cancer. Patients with a final benign diagnosis were predominantly due to significant typical bacterial infection/inflammation (pneumonia). However in the AEGIS cohorts many of the benign diagnoses, where known, were due to sarcoidosis, fibrosis, benign tumours or atypical infections (fungal and mycobacterial). Therefore, in our cohort the pre-test probablility for malignancy in the benign group was higher than in the AEGIS benign group.

Gene expression response to smoke injury differs between healthy volunteers and clinic patients

Intrigued by these overall expression differences between volunteers and clinic patients we investigated the post-cessation dynamics of individual genes using a population-based approach. We first employed a Bayesian linear regression model to predict nasal gene expression in healthy volunteers as a function of smoking status, accounting for sex and age (Methods). This model classified genes as either unaffected by smoking (US), rapidly reversible (RR; no difference between former and never smokers), slowly reversible (SR; intermediate expression levels in former smokers compared to never and current), or irreversible (IR; no difference between former and current smokers). Additionally, genes were classified as cessation-associated (CA) if no difference was present between current and never smokers, but elevated or reduced expression was observed in former smokers (see Sup. Fig. 1 for a schematic).

In healthy volunteers 5755 genes were found to be affected by smoking status, out of which 513 genes show a strong effect (effect size > 0.4 for rapidly reversible, slowly reversible, irreversible genes, > 0.25 for cessation activated genes, Methods, Sup. Tab. 5). Most genes (485/513) were found to be rapidly reversible, in line with previous findings in bronchial tissue (9). GO pathway analysis of these genes revealed up-regulation of cellular detoxification, response to oxidative stress (e.g. CYP1A1, CYP1B1, AHRR, NQO1, GPX2, ALDH3A1) and keratinization (e.g. KRT6A, KRT13, KRT17, SPRR1A, SPRR1B, CSTA) pathways, and down-regulation of cilium organization (e.g FOXJ1, DNAH6, IFT81, CEP290, UBXN10), extracellular matrix organization (e.g. FN1, COL3A1, COL5A1, COL9A2) and interferon signaling (e.g IFI6, IFIT1, IFI44, RSAD2) in current compared to never smokers. Genes involved in inflammatory response were found both among the up-regulated (IL36A, IL36G, S100A8, S100A9, CLU) and down-regulated (SAA1, SAA2, IL33) genes. Principal components analysis using the rapidly reversible genes showed a clear separation of current smokers from all other subjects. In contrast, slowly reversible and irreversible genes placed patients on a trajectory from never smokers to current smokers, as expected (Sup. Fig. 2a).

We next repeated the above analysis on the clinic subgroup. In the absence of clinic never smokers, and since no technical or biological covariates could explain the observed overall expression differences between the groups (see Methods), we considered the healthy volunteer never smokers as a bona fide reference group for this analysis. We found 4112 genes with smoking-dependent expression changes, 584 of which showed a strong effect (same effect size thresholds as above, see Methods and Sup. Tab. 5). We evaluated this classification with a principal components analysis on the clinic subjects, similar to what was done for healthy volunteers, and found that patients clustered according to their smoking status, as expected (Sup. Fig. 2b). Of the 584 genes identified as dysregulated by smoke in the clinic patients, 233 were also found in the healthy volunteer analysis (P < .001, chi-squared test). However, while 227/233 were rapidly reversible in the healthy volunteers, only 112 of these 227 genes were also classified as rapidly reversible in the clinic group (Fig. 2a). Of the remaining 115 genes, one gene (BPIFA2) was now classified as irreversible and 22 genes as slowly reversible, including CYP1B1, a well-known detoxification gene, and BMP7, a gene previously shown to have a role in immunoregulation (22) (Fig. 2b). Ninety-two of 115 genes that were classified as rapidly reversible in the healthy volunteer group and as cessation-associated in the clinic group (e.g. UBXN10, Fig. 2b) showed a strong enrichment for cilia structure and function (Sup. Tab. 6). While cilia-associated genes were down-regulated in current smokers in both groups (consistent with cigarette smoke damaging airway cilia), the same genes showed increased expression in current and former smokers in the clinic group compared to the healthy volunteers. This observation in the clinic group might be linked to the decreased expression of interferon gamma-related genes in the clinic group, as it has been shown that interferon gamma suppresses ciliogenesis and ciliary movement (26).

• Download figure
• Open in new tab

Figure 2: Smoke injury dynamics.

(a) Plot showing a change of reversibility dynamics for the same genes in the healthy volunteer (left) and clinic (right) donor groups (genes classified as unaffected by smoking in both donor groups were removed); (b) Normalized gene expression over smoking status for 4 exemplar response genes with different post-cessation dynamics in the clinic and healthy groups, with linetype and shape representing donor status and colors representing the genes’ assigned reversibility classes.

Lastly, the 351 genes that showed smoking-dependent expression changes in the clinic group but not in the healthy volunteers (Fig. 2a) were strongly enriched in extracellular matrix organization and immune-related genes (including response to interferon gamma, neutrophil activation, chemotaxis and inflammation). For example, GBP6 showed down-regulation and slow reversibility in the clinic group (Fig. 2b) and is known to be associated with reduced overall survival in squamous cell carcinoma of the head and neck (27).

Overall, we observe striking differences in smoke-dependent gene expression in the clinic patients compared to volunteers that could not be explained by comorbidities or other covariates, with generally slower reversibility post-cessation in the clinic group. We hypothesise that some of the 749 genes with differences in smoke-dependent expression might reflect individual responses to the smoke injury and thus refer to them as response genes.

Response gene expression levels predict disease status and may improve risk stratification for population screening

We postulated that the smoke-injury response genes we identified might provide evidence for a personalised smoke injury response and be candidate genes for a molecular biomarker of lung cancer risk. In the clinic group, where patients already show evidence of lung disease, such a biomarker would help identify patients with the highest need for further investigation. In the general smoker and former smoker population it could be added to existing methods of risk stratification to improve the identification of individuals who would most benefit from lung cancer screening thereby sparing those at lowest risk who would have least to benefit from screening.

Therefore, we trained two independent classifiers: a ‘clinic classifier’ that predicts the cancer status of each sample (cancer vs clinic benign and healthy volunteers: potentially of use in the clinic), and a ‘population classifier’ that predicts the donor group that the samples were taken from (clinic benign or clinic cancer vs healthy volunteers: potentially of use in risk stratification for population screening). For both classifiers, we used gene expression data from the 749 response genes together with clinical information (sex, age, smoking status and pack-years; see Methods) in a lasso-penalized multivariate logistic regression, and derived a log-odds score from each classifier. In line with the observed strong expression differences between healthy volunteers and clinic patients, the ‘population’ score clearly separates healthy volunteers from clinic subjects (Fig. 3a). Interestingly, the ‘clinic’ score (Fig. 3b) additionally distinguishes the benign and cancer patients within the clinic group, placing benign subjects between healthy volunteers and cancer subjects. As expected, the two scores are highly correlated (Pearson correlation = 0.8, P < .001, Sup. Fig. 3a). Both scores yielded high area under the curve (AUC) values for both precision-recall (clinic score: mean AUC-PR=0.83 after 10-fold cross validation; population score: mean AUC-PR=0.85, 10-fold cross validation, Fig. 3c-d) and receiver-op-erator characteristics (clinic score: mean AUC-ROC=0.84, 10-fold CV; population score: mean AUCROC=0.92, 10 fold CV, see also Methods) and performed significantly better than a model using the same number of randomly selected genes (Sup. Fig. 4). In practice, to reach a sensitivity of 95% for the population score, one would use a score threshold of 2.69, that would result in an average false positive rate of 42.8%, while to reach a similar sensitivity using clinical data alone would result in a false positive rate of 74.5%. For the clinic score, a score threshold of -1.46 gives a 95% sensitivity and false positive rate of 62.1%, while similar sensitivity with clinical data alone would result in a false positive rate of 67.8% (Fig. 3c-d). These results indicate that models incorporating gene expression data of the response genes defined above performed significantly better than models built on clinical covariates alone (see also inset of Fig. 3c-d for a comparison of the performance of models based on gene expression data alone, clinical covariates alone or a combination of gene expression data and clinical covariates). In addition, both scores retained their ability to separate the patient groups after regressing out all potential confounders, confirming that gene expression data improves classification compared to using clinical covariates alone (Sup .Fig. 3b-c).

• Download figure
• Open in new tab

Figure 3: Disease status prediction based on response genes.

(a & b) Risk score distribution for the population test (a) and the clinic test (b) predicted from the clinical variables and the expression of the response genes using a penalized regression (See Methods). The risk distributions are presented separately for healthy volunteers (green), clinic patients without cancer (orange) and clinic patients with cancer (purple) (c & d) ROC curves for the population (c) and clinic (d) scores. For each case, we present the ROC curve for the model trained on clinical data (triangles) or on gene expression and clinical data (squares). Each curve is an average obtained across 100 cross validation (CV) experiments and the grey area surrounding the curve gives the standard error. The color of the curve represents the test threshold corresponding to the represented Sensitivity / False Positive Rate compromise. (Inset) Area under the ROC curve, in 100 CV rounds, for a clinical-only model (red) the model constructed on the response genes (blue) and a model constructed on a combination of clinical information and response genes (green) for the population (c) and clinic (d) classifiers. (e) The population and clinic classifiers applied to nasal samples from the AEGIS cohort.

We also assessed the performance of the trained population and clinic risk score models separately on current and former smokers. We found that the population risk score is equally applicable to current and former smokers: a significant difference in the risk score of the healthy volunteers and clinic subjects can be observed, even after regressing out clinical covariates and confounding (Sup. Fig. 5). While the clinic risk score performs well on both groups, the added value from gene expression data appears less important in the clinic score, in particular in former smokers (Sup. Fig. 5). Finally, we find that our classifiers are efficient at separating subjects regardless of their cancer stage, cancer type (squamous carcinoma or adenocarcinoma), and COPD status (Sup. Fig. 6), and that our classifiers capture differences in risk that persist for more than 10 years after smoking cessation (Sup. Fig. 7).

Finally, we validated our classifiers by applying them to an independent cohort. No publicly available cohort matches the composition of our cohort, in particular because of the absence of a healthy group of current and former smokers distinct from the clinic-referred patient group. However, the AEGIS cohort (28) includes nasal samples from clinic-referred patients with pulmonary nodules and a diagnosis of lung cancer or benign disease. We applied our two classifiers to this cohort, and found a good separation between subjects with and without cancer, despite the different gene expression quantification technologies and populations of origin of the patients (Fig. 3e, Sup. Fig. 8). We found a stronger separation between patients with and without cancer using the AEGIS nasal classifier from Perez-Rogers et al (2017) (16) on the AEGIS data (Sup. Fig. 8a). However, we note that the AEGIS classifier (16), when applied to our data, mostly differentiates healthy volunteers and clinic patients while the difference between the scores of the cancer and no-cancer patients is only modest (Sup. Fig. 8b). These results confirm the ability of our classifier to stratify patients, even when applied to patients from different clinical contexts.

Overall, our results demonstrate that classifiers based on nasal gene expression have the potential to improve risk stratification of current and ex-smokers in both a population screening context and a clinic context.

Alterations in immune pathways underlie the lung cancer risk classification

To gain insights into the mechanisms of risk, we asked which genes robustly contributed most to the classifiers by identifying genes selected in more than 80% of the cross validation (CV) rounds (Sup Fig. 9). Among the 46 genes selected most often in either of the risk prediction models, we found genes that were previously identified as important players in lung cancer development, e.g. SAA2 [18], HAS2 (29), (30, 31) or TGM3 (32–35), in line with the current literature.

However, the genes used as predictors of risk in our model reflect a wide variety of smoking-associated alterations. In order to gain some mechanistic insight, we investigated risk contribution at the pathway level. First, we performed GO enrichment analysis on the list of smoke injury genes (both the ones identified in the healthy volunteers and in the clinic group) to identify the main pathways affected by smoke. We found that the smoke injury genes are mainly involved in xenobiotic metabolism and response to oxidative stress, extracellular matrix organization, keratinization, ciliary structure and mobility, and immune response (Sup. Tab. 7). We then chose 8 GO terms as represen-tatives of these alterations: Keratinization, Extracellular matrix organization, XenobioIc metabolism, Cilium organization, Inflammatory response, Neutrophil mediated immunity, Response to interferon gamma, and Antigen processing and presentation. We calculated geneset metascores for each of these GO terms (Fig. 4a and Sup. Fig. 10). For some of these pathways, such as Keratinization, we observed a similar, rapidly reversible dynamic in healthy volunteers and clinic patients (Sup. Fig. 10a). For most pathways, however, the dynamics were different in the two donor groups. Cilium organization appeared to be rapidly reversible in healthy volunteers, while in clinic patients it showed increased expression in former smokers, with no difference between current and never smokers. XenobioIc metabolism showed a slower reversibility in clinic patients than healthy volunteers (Sup. Fig. 10a). For all immune-related pathways, we observed reduced expression in current smokers, and a slow reversibility dynamic, uniquely in clinic patients (Fig. 4a); we also observed that their activity does not revert to healthy never-smoker level even long after smoking cessation (Sup. Fig. 10b).

• Download figure
• Open in new tab

Figure 4: Pathway analysis and contribution to risk.

(a) Comparison of geneset metascore over smoking status for 4 immune-related GO terms in healthy and clinic subjects. (b) Correlation between the population or clinic risk score and geneset metascore for the 8 gene sets representing biological functions altered by smoking; correlation is shown separately for current and former smokers (> 12 months); Spearman correlation values are reported, as well as the associated p-values (*: P <= 0.05, **: P <= 0.01, ***: P <= 0.001).

To identify which of these pathways contributed most to increased risk, we then calculated the correlation between geneset metascore in each subject and subjects’ risk scores from the population and clinic classifiers. We calculated these correlations for current and former smokers (> 12 months) separately, to be able to identify differences in geneset contribution to risk in the two groups that might reflect differences between acute smoke injury response and the long-term consequences of past smoke exposure (Fig. 4b). In current smokers, while Keratinization and Extracellular matrix organization did not significantly correlate with either risk score, the remaining four genesets tested showed moderate but significant correlation with both risk scores, pointing to alterations of the xenobiotic detoxification pathways, ciliary function and immune response as major contributors to patient-specific differences in risk. In former smokers, the population risk score correlated with the same 4 GO terms indicating that detoxification pathways, ciliary function and immune response are the main contributors to overall risk of lung disease. In contrast, only pathways related to immune alterations (Response to interferon gamma, Neutrophil-mediated immunity, Antigen processing and presentation) correlated with the clinic risk score in former smokers, while no correlation was observed with XenobioIc metabolism, and only a very weak correlation with Cilium organization (Fig. 4b). These results indicate that immune alterations are significant contributors to the risk of cancer in both current and former smokers in the clinic group.

Patient-specific genetic background modulates the smoke injury response

Germline genetic variation may influence individual differences in response to airway smoke injury, and hence, risk of smoking-related lung cancer. To investigate this, we first conducted an eQTL analysis on nasal and bronchial epithelium separately and jointly to identify variants that affect the expression of neighbouring genes (Methods). We obtained 990 (bronchial), 1316 (nasal) and 1695 (combined) eQTL effect genes (e-genes) at 1% FDR. We found a significant overlap between the nasal and bronchial e-genes (Sup. Fig. 11a), with 574 genes in common (corresponding to 58% and 44% of the bronchial and nasal eQTL respectively, Fisher’s exact test P < .001). Similarly, we found a correlation of 0.56 between the adjusted p-values of the lead variants between both sets (Sup. Fig. 11b), confirming shared cis-regulation between the nasal and bronchial epithelium.

To further study the interaction between subject-specific genetic background and environmental factors we next leveraged this eQTL catalogue to search for genetic variants within the 749 response genes that might modulate gene expression differently depending on subjects’ smoke exposure. We identified 78/749 genes with at least one lead eQTL variant with genome wide significance at 10% FDR, (Sup. Tab. 8). We then tested for an interaction effect between smoking status and genotype for all 78 lead eQTL variants on gene expression. We identified 11 genes (CH25H, LHX6, WNT5A, DRAM1, SULF1, LGALS7B, HAPLN4, FXYD5, EFCAB2, TOX and SPRR1A, see Sup. Fig. 12) whose expression changes in response to smoke are modulated by the presence of genetic variants (nominal P < .1, Sup. Tab. 8), suggesting that those genetic variants might modulate the response to smoke injury and to lung cancer risk. For example, up regulation of FXYD5 has been shown to correlate with tumor size (36) and poor survival (37) in NSCLC and to be implicated in many cancer types as FXYD5 enhances NFκ-B transcriptional activity, promotes angiogenesis and increases tumor cell’s migration and invasion abilities (38). Finally, this protein also promotes inflammation in epithelial cells, notably in lung tissues (39). Analysing the expression of this gene in our cohort, we find that subjects with a homozygous reference genotype at the 19:35660670:G:A locus have similar levels of expression both in never, ex, and current smokers (Fig. 5a). On the contrary, subjects that have a heterozygous or homozygous alternative genotype present higher levels of expression of this gene in response to smoke (Fig. 5a), which might increase their lung cancer risk. We observe similar trends for the 10 other response genes stated above (Sup. Fig. 12, Sup. Tab. 8). This finding demonstrates how subjects’ specific genetic background can influence their reaction to cigarette smoke and in turn might affect their risk of developing lung cancer.

• Download figure
• Open in new tab

Figure 5: Genotype background influences lung cancer risk

(a) Combined environmental and genetic effect on the expression of the FXYD5 gene in nasal tissues. For each nasal sample, we present the expression level of the gene FXYD5 separately for never (pink), former (green) and current (blue) smokers. Samples are further stratified depending on the genotype of the subject at the 19:35660670:G:A locus (Ref/Ref: homozygous reference; Ref/Alt: heterozygous; Alt/Alt homozygous Alternative). The p-value gives the significance level of an interaction effect of the smoking status and the genotype at the 19:35660670:G:A on the expression of the FXYD5 gene (see Methods). GWAS enrichment analysis: (b) Network representation of the 4 bronchial regulons enriched in GWAS genes. The 4 TFs are shown as squares and their targets in the bronchial network as circles. Color of the nodes indicates whether the gene/TF is a smoke injury risk gene (blue), a gene that co-localizes with a GWAS hit (i.e. no threshold on eQTL significance) (red) or both (green). The level of overrepresentation for genes in the network of those TFs can be found in Table 2 (for the GWAS) and Table 3 (for the response genes); (c) Activity level of each of the 4 TFs in nasal tissue, depending on the disease status of the patient.

Common germline variants regulate interferon gamma genes and are linked to known lung cancer risk loci

We next identified GWAS hits that were in strong linkage disequilibrium in the UK population to SNPs that we found to be regulating the expression of nearby genes in our eQTL analyses (Methods). Among the 1261 GWAS lung cancer risk loci, our analysis identified 63 GWAS risk loci from 13 different studies with variants that significantly affect the expression of a nearby gene at a 5% FWER threshold (Sup. Tab. 10). These 63 eQTL/GWAS variants were linked to the expression of 41 genes, notably including 10 genes implicated in the interferon gamma signalling pathway. Pathway enrichment confirmed a strong enrichment for genes involved in response to interferon gamma (hypergeometric test, Padj = 7 × 10⁻¹³), as well as for other immune-related functions (e.g. innate immune response, anIgen processing and presentation of exogenous pepIde anIgen, regulation of immune response, T cell receptor signalling pathway; see Sup. Tab. 11 for the full list of enriched GO terms).

To better understand the mechanisms by which GWAS variants might increase lung cancer risk, we looked for a link between 41 genes linked to a GWAS risk locus and transcriptional regulatory network in bronchial tissue. To do so, we inferred a TF-target interaction network from bronchial expression data (see Methods) and searched for TFs whose targets were enriched for the 41 genes. We found 4 TFs showing a strong enrichment (hypergeometric test, Padj < .05, see Methods), ETV7, SPIB, IRF1 and CIITA (Fig. 5b), all of which are known players in the interferon gamma mediated signalling pathway (40–43). We further confirmed the enrichment of GWAS genes in these 4 TFs by using a wider list of GWAS genes with a relaxed eQTL cut-off (nominal P < .05), and still found a 2- to 3-fold enrichment in all 4 TFs (Table 1). Analyzing the activity of those 4 TFs in the nasal samples, we found significant differences between healthy volunteers and clinic patients, in particular a lower activity in clinic patients, confirming the importance of these 4 TFs in the progression toward a disease status (Fig. 5c, and see Sup. Fig. 13 for the activity of the same 4 TFs in the bronchial samples of clinic patients with and without cancer). In contrast, we found that the levels of activity of those 4 TFs were similar in clinic patients with and without cancer (Fig. 5c and Sup. Fig. 13). We further tested whether our set of response genes was enriched within the targets of those TFs and indeed found that all 4 TFs are enriched for response genes (2 to 3-fold enrichment, nominal P < .05, Table 2).

Altogether, these findings suggest that the effects of inherited variation on lung cancer risk may be exerted in part through a different immune response following smoke injury, creating an immunosuppressed environment that favours the final steps to the emergence of a cancer.

Availability of the code and data

Scripts that were used to conduct the analysis presented in this paper are available in the bitbucket repository accessible with the following link:

https://bitbucket.org/schwarzlab/paper-debiase-massip-2021

Raw and processed data will be available upon acceptance of the paper on the EGA website.

Cohort and sample collection

487 donors were recruited into the CRUKPAP cohort at Royal Papworth Hospital, Cambridge (UK), including 114 healthy volunteers (HV) and 337 patients being investigated for suspicion of lung can-cer. The eligibility criteria for healthy volunteers were: age 18 or above; current or former smokers must have smoked at least 100 cigarettes in their lifetime; exclusion of individuals with previous history or current suspicion of airway or lung cancer; exclusion of individuals with ‘bleeding disorders’.

All participants were stratified into smoking cessation categories as follows: 45 never smokers (NV), 289 former smokers (FS) and 153 current smokers (CS). Former smokers were further divided into categories: > 1 year after cessation (FS1, n=234), 1-12 months after cessation (FS2, n=45) and < 1 month after cessation (FS3, n=10). Smoking status for all subjects was confirmed via blood cotinine test. Cumulative smoke exposure measured in pack years was recorded and stratified into four categories: ‘none’ (PY1), < 10 years (PY2), 10-30 years (PY3) and > 30 years (PY4). For suspected lung cancer patients, both COPD status and final cancer diagnosis (lung cancer / no lung cancer) were recorded.

From these donors 413 nasal epithelial curettages were collected using Arlington Scientific ASI Rhinopro nasal curettes. Briefly, the nostril is opened with a nasal speculum to identify the inferior turbinate. Under direct vision the tip of the nasal curette is gently scraped over the turbinate to obtain a ‘peel or curl’ of epithelial tissue. The curl of tissue is then removed by flicking the curette while the tip is submerged in RNAlater™ collection medium and presence of the curl floating in the medium is confirmed by visual inspection. This procedure is repeated twice for each nostril per donor. RNA integrity (RIN) was checked for all samples and we found >80 % of samples to have RIN 6 or better.

Bronchial brushings were collected using 2.0mm brush diameter cytology brushes (Olympus Medical, UK) from 236 patients undergoing flexible bronchoscopy as part of investigations for suspected lung cancer.

For 162 donors, both nasal and bronchial samples were available. Sample collection and diagnosis took place contemporaneously. All samples underwent short-read RNA sequencing using Illumina TruSeq library generation for the Illumina HiSeq 2500 pla{orm. Blood samples were taken from 467 donors and germline genotyped using the Illumina Infinium Oncoarray pla{orm at 450K tagging germline variants. Total gene expression levels (TPM and variance stabilised) were determined for 18,072 protein coding genes for all samples using DeSeq2.

Research ethics approvals for sample collection from participants in this study were given by East of England Cambridge Central REC 13/EE/0012 and the National Research Ethics Service Committee South East Coast – Surrey 13/LO/0889.

RNA extraction and sequencing

Tissue samples from bronchial brushings and nasal curettes were stored in 500μl RNALater overnight at 4oC, and then at -80oC for longer-term storage. RNA was extracted using Qiagen MiRNeasy columns according to the manufacturer’s protocols. Briefly, bronchial brushes were rinsed in PBS, brushes transferred into 700μl Qiazol and cells lysed by vortexing twice for 30 seconds. For nasal samples the RNALater containing nasal tissue (500μl) was diluted with 2ml of PBS and spun at 10,000 rpm for 10 min. The cell pellet was lysed by resuspension in 700μl Qiazol. For both types of samples, the Qiazol lysate was applied to a QiaShredder tube (#217004) and spun at 13,000 rpm for 2 mins. The homogenate was kept at room temperature for 5 mins, followed by chloroform extraction using PhaseLock tubes. Nucleic acids in the aqueous phase were precipitated using 1.5 volumes of 100% ethanol and DNA was digested using DNAse I. Finally, RNA was isolated from the mixture using RNAeasy mini spin columns. RNA was quantified using a Qbit measurement and quality assessed us-ing an Agilent Bioanalyzer. For samples with a RIN greater than 7, a total of 500ng of RNA was used for Illumina TruSeq Library generation. Sequencing was carried out on HiSeq 2500 Illumina sequencers. Sequencing was carried out in two separate multiplexed experiments.

RNA sequencing data processing

Quality control using FastQC showed good sequence quality and no adapter contamination for all samples. Alignment was carried out with TopHat2, using as reference the human genome version GRCh37. Read counts were computed for all protein-coding genes with subread featureCounts v1.6.0. The data was produced in 2 experimental batches, producing a strong batch effect that can be observed in the raw data. Moreover, a group of samples from one of the batches has lower total counts compared to the other samples (Sup. Fig 14a).

Raw counts were normalized using DESeq2’s variance-stabilizing transformation, which had the advantage of partly correcting the previously mentioned batch effects (Sup. Fig 14a). Genes with across-samples log variance smaller than -4 were discarded from further analysis. Total gene expression levels (variance stabilised) were determined for 18,072 protein coding genes for all samples. To ascertain that experimental batch did not covary with any clinical covariate, we computed the strength and significance of association between batch and the other covariates using Cramer’s V and chi-square test. We did not observe a significant association between batch and age, sex, COPD, smoking status, pack-years and donor population of origin (healthy volunteer/clinic patient). We only observed a weak but significant association with cancer status (Sup. Fig 14b).

To assess the overall contribution of clinical and environmental variables to gene expression in the nasal epithelium, we also extracted variance components using a linear model, regressing donor population of origin (healthy volunteer/clinic patient), cancer status, smoking status, pack-years, sex, age, COPD and experimental batch against total gene expression across all genes (Sup. Fig 14c). We found that donor population of origin and smoking status contribute most to gene expression vari-ability (28.8 and 25.4% of total explained variance. Notably, donor population of origin still con-tributes significantly to the explained variance after accounting for all other clinical and technical covariates.

Differen3al expression analysis

All differential expression analyses were performed with DESeq2 v1.26.0. Age, experimental batch, sex and pack-years were included as confounding variables. Genes with multiple-testing-adjusted (Benjamini-Hochberg) p-values < 0.05 were considered differentially expressed. For differential ex-pression between clinic cancer and clinic benign in bronchial samples, 8 genes had artificially high (>20) absolute fold-change, due to their very low average expression across samples. These genes were removed from the list of differentially expressed genes.

Modelling time-dependent dynamics of smoke injury in nasal 3ssue

Derivation of population and clinic risk scores

L1-penalized multivariate logit regression was performed with R package glmnet 3.0-2 using only the nasal gene expression data. Patients status was encoded with a binary variable (cancer: 1; no cancer 0 for the clinic classifier; Clinic patient: 1, Healthy Volunteer: 0 for the population classifier), and patients with Ineligible status were excluded from the analysis. In the gene expression classifiers, the status of each patient was predicted based on the expression of the 749 response genes and 4 clinical covariates, namely sex, age, smoking status and packyears, all of which were encoded as numerical variables (smoking status encoding: Never smokers: 0, Ex >1year: 1; Ex 1-12months: 2; Ex <1m: 3, current smokers: 4). For the clinical classifier we also used a lasso regression, using only sex, age, smoking status and packyears as predictors. The lasso shrinkage parameter (λ) was chosen to minimize the mean cross-validated error (“lambda-min” option in the cv.glmnet function). Area under the receiver operating characteristic curve and precision recall curves were computed using the PRROC package, after 10 rounds of 10-fold cross validation experiments. To compare performances of the response genes to performances on random genes, we randomly drew 20 sets of 749 genes among the 18,072 protein coding genes retained for all analyses, and cross validations experiments were conducted on the same test and training set as the one used with the response genes.

Gene ontology analysis and pathway analysis

All Gene Ontology (GO) enrichment analyses were performed using clusterProfiler v3.14.3. GO terms with adjusted (Benjamini-Hochberg) p-values < 0.05 were considered enriched.

Pathway metascores were calculated by averaging vst-normalized gene expression of genes belonging to the selected genesets, after regressing out experimental batch effect.

Genotyping data pre-processing

eQTL Analysis

We computed the eQTL tests for the set of 18,072 protein coding genes for which we have sufficient coverage (see filter criteria for RNAseq data above). For each gene, we tested all SNPs in a 500kb cis window (500kb upstream from the TSS, 500kb downstream from the transcription termination site).

For each test, we model the effect of known clinical and technical covariates (Sex, Age, Batch, Smok-ing Status and PacksYear) using a fix effect. All clinical covariates were encoded as numerical values (0-4 for smoking status, 0-3 for age and packYears, and binary 0-1 for sex and batch), and genotypes are encoded as a numeric variable (0: Ref/Ref; 1: Alt/Ref; 2: Alt/Alt). P-values were computed using the R package Matrix eQTL (49). We used a two step multiple testing correction procedure, as described in (50). First, for each gene, we correct for the number of tests using Bonferroni correction. Second, we performed a global correction across the lead variants, that is, the most significant SNPs, per eQTL, using a Benjamini-Hochberg procedure.

Gene environment interaction test

To test for a combined effect of genotype and environment on the gene expression level of the smoke injury gene, we conducted an interaction test between the genotype background and the smoking status of the patient, encoded in a 0/1/2 form (Never/Ex/Current). For each of the 749 smoke-injury gene, we retrieved the lead eQTL variant identified in the genome-wide eQTL analysis, and tested for an interaction effect between the genotype encoded in a 0/1/2 numeric and the smoking status, correcting for the effect of age, sex, smoking status, packyears, and genotype.

Identification of GWAS-linked genes

To study the mechanisms by which germline genotype background influences the lung cancer risk, we adopted the approach developed by (51). We downloaded a curated set of 1261 GWAS lung cancer risk loci from the GWAS catalog (52) (see Sup. Tab. 9) and mapped genotyped and imputed SNPs of all patients to the nearest GWAS risk locus as follows. For each GWAS risk locus, we retrieved a list of variants in our cohort within a 500kb cis-window using a linkage disequilibrium (LD) cutoff of R^2>0.8 in the UK population using the Linkage Desiquilibrium Calculator of the ensembl website (53), yielding 9,739 candidate variants and 135,513 gene-SNP pairs. 3,455 of those 9,739 variants had a significant effect to their corresponding e-gene at a 5% FDR threshold. Many of those 3,230 hits were in LD with the same GWAS variant, such that all eQTL variants mapped to 67 unique GWAS risk loci (Sup. Tab. 8) from 10 different studies and were linked to the expression of 44 genes.

Transcription factor network and activity

A context-specific protein-protein interaction network for nasal and bronchial epithelium was built using ARACNe-AP (54) on the vst-normalized expression data and a list of 1988 human transcriptional regulators, compiled using information available on public databases, from (55). ARACNe-AP was able to infer context-specific interactions across 1548 nasal and 1535 bronchial regulators. Activity of each of these regulators in each nasal and bronchial sample was inferred using VIPER v1.20.0 (56). Network representations of TF-TF and TF-targets interactions were produced with Cytoscape v3.8.1. To find TFs that had an overrepresentation of GWAS genes in their target network, we used a context-specific TF-TF interaction network built using ARACNe-AP on bronchial vst-normalized gene expression data and a list of 1988 human transcriptional regulators (see above). For each TF i, we first counted the number (N_G(i)) of genes in its target network that were identified as a GWAS gene. We then compared the proportion of GWAS genes in each TF target network to the expected number that would be found for a similar number of randomly selected genes with a one-tailed hypergeometric test using the phyper function in R with the following parameters:

m: total number of genes in the network of TF i; n: 18,062 m; k= the number identified of GWAS genes and q = N_G(i), the number of GWAS genes in the target network of the TF i. Obtained p-values where adjusted for multiple testing using a Benjamini-Hochberg correction. We applied the same procedure to test for the enrichment of response genes in the 4 identified GWAS TFs, although we did not correct the p-values for multiple testing this time since we conducted only 4 tests.