Cell‐type <scp>eQTL</scp> deconvolution of bronchial epithelium through integration of single‐cell and bulk <scp>RNA</scp>‐seq
Authors: Cancan Qi, Marijn Berg, Xiaojing Chu, Wim Timens, Tessa Kole, Maarten van den Berge, Cheng‐Jian Xu, Gerard H. Koppelman, Martijn C. Nawijn, Yang Li
Journal: Allergy (2022)

Expression quantitative trait loci (eQTL) analyses may help to understand the function of disease-associated genetic risk factors. However, eQTL studies obtained from bulk RNA-seq data are limited by the lack of cell-type resolution. In this study, we hypothesized that eQTL deconvolution1 of bulk RNA-seq datasets may help to identify cell-type-associated regulation of gene expression by asthma-associated genetic variants in bronchial epithelial cell subsets. We present a tool for predicting bronchial cell proportions in bulk bronchial RNA-seq datasets, called “BronchiCellPred” (https://git.web.rug.nl/GRIAC/BronchiCellPred). The description of the methods can be found in the Appendix S1. Relative expression of the signature genes in each cell type of this tool is shown in Figure S1. We used two single-cell (sc-) RNA-seq datasets as validations to evaluate the performance of this tool, and found that the estimated cell proportions were significantly correlated with the measured frequencies for most of the cell types with spearman correlation R > .7 and p values <.05, except for the rare cell-type ionocyte (Figure 1B,C). We then applied eQTL deconvolution1 to bulk RNA-seq data of bronchial biopsies of 146 subjects from the INDURAIN cohort (Figure S2), to identify cell-type-associated cis-eQTLs of asthma-associated genetic variants,2, 3 in three main epithelial cell types (ciliated, basal and secretory cells). In total, we identified 14 decon-eQTLs (p value <.001, Table S1). For 11/14 eQTLs, decon-QTL associations were more significant than bulk-eQTL (Table S2). These decon-eQTLs were replicated using an independent bronchial epithelial single-cell-eQTL dataset (N = 28). Three eQTLs were replicated with p value <.05 (rs4749894 ~ CALML5 in basal cells, rs8103278 ~ FOXA3 in secretory cells, and rs10876864 ~ RPS26 in basal cells) (Figure 2, Table S2). All the study protocols were approved by the local medical ethics committee. All subjects gave their written informed consent. The eQTL analysis of asthma SNPs in bronchial tissue is still limited, especially at the cell-type resolution.2 To the best of our knowledge, this is report of the first cell-type-associated eQTL of asthma-related SNPs from bronchial biopsies by a deconvolution method, confirmed by single-cell eQTL. These approaches may help us to identify the potential targeted genes and cell types for drug discovery and improve our understanding of the regulation of genetic variants on the disease development. Two of the three replicated decon-eQTLs have previously been identified as bulk-eQTLs by GTEx (https://gtexportal.org/), including rs10876864 ~ RPS26 which were identified in multiple tissues including lung tissue, and rs8103278 ~ FOXA3 which were identified in minor salivary gland. Our study revealed the eQTL effect in cell types relevant to asthma. FOXA3 is a transcription factor expressed in goblet cells in the airways, which induces mucus production and goblet cell differentiation and may be involved in antiviral responses of epithelial cells.4 In our study, we showed that the asthma-associated risk allele was associated with increased FOXA3 expression in secretory epithelial cells, a cell-type annotation that encompasses both club and goblet cells,5 thereby providing a molecular mechanism for goblet cell hyperplasia and increased mucus production in asthma. CALML5, expressed in bronchial epithelial cells and keratinocytes, is important for epithelial differentiation and barrier function.6 We identified a basal cell-specific increased expression of CALML5 associated with the asthma-associated variant, which may play a role in epithelial repair and differentiation, and thus, in maintaining and/or restoring epithelial integrity in asthma. The accuracy of cell-type estimation is important for eQTL deconvolution.1 Most of the cell types can be well estimated using our methods, which was supported by the validation analyses. However, the method may not have enough power to predict rare cell type such as ionocytes. This should be kept in mind when filtering on the estimated cell-type proportion to be applied in downstream analysis. Many of the decon-eQTL genes showed low expression in scRNA-seq data (Figure S3), which makes it difficult to identify these cell-type eQTL by sc-eQTL. This indicates that the deconvolution methods may even be more powerful to identify these cell-type eQTL effects. As shown in our results, we were able to identify the cell-type eQTL effect with a limited sample size, some of which can be replicated by a sc-eQTL method for genes that had sufficiently high expression in the scRNA-seq datasets. However, there are still some limitations in this study: 1) the sample size of bulk data is small. We were only able to assess the three main epithelial cell types, and we applied a suggestive p value threshold considering the limited power. Here, we demonstrated the eQTL deconvolution analysis by focusing on the three abundant cell types. The method can be applied to include less abundant cell subpopulations when the sample size increase. 2) The identified eQTL genes were lowly expressed in scRNA-seq data in the replication dataset, leading to low power in our validation design. In conclusion, this study provides the first proof of concept of using eQTL deconvolution method to identify cell-type eQTL for asthma-associated genetic variants in bronchial biopsies. The method to deconvolute cell types from the bronchial airways is openly available for future use. The eQTL deconvolution analyses are likely to be a powerful tool to better understanding cell type-specific effects of asthma-associated genetic variants, thereby providing a broader overview of the molecular mechanisms of asthma development and can be applied on other large bulk RNA-seq datasets in respiratory disease. YL wasti supported by an ERC starting Grant (948207) and a Radboud University Medical Centre Hypatia Grant (2018). MCN is supported by the Lung Foundation Netherlands (grant no. 5.1.14.020) and the European Union's H2020 Research and Innovation Program under grant agreement no. 874656 (discovAIR). GHK is supported by a ZON-MW VICI grant. CQ is supported by a grant from the China Scholarship Council. CJX is supported by the Helmholtz Initiative and Networking Fund (1800167). The authors wish to acknowledge the patients and participants involved in the study. Open Access funding enabled and organized by Projekt DEAL. The authors declare no conflict of interest. Figure S1 Figure S2 Figure S3 Table S1-S2 Appendix S1 Please note: The publisher is not responsible for the content or functionality of any supporting information supplied by the authors. Any queries (other than missing content) should be directed to the corresponding author for the article.