Cis-eQTL analysis of Korean asthma patients
A total of 436 Korean asthma patients were used for the eQTL analysis, and their basic characteristics are shown in Table 1. The mean age of these asthma patients, at diagnosis, was 45.1 ± 16.4 years, and the proportion of males was 40.6%. Among them, 42.9% were current or ever smokers and 45.4% were asthma patients with atopy. The FEV1 and FEV1/forced vital capacity (FVC) were 75.2 ± 18.0 and 72.7 ± 13.4, respectively.
Table 1
Characteristics of 436 Korean asthma patients
Characteristics | Korean asthma patients mean ± SD / N (%) |
Age, years | 49.5 ± 15.5 |
Sex (%) | |
Male | 177 (40.6) |
Age at diagnosis, yearsa | 45.1 ± 16.4 |
Age group at diagnosisa | |
0 ~ 19 | 25 |
≥ 20 | 341 |
Smoking status (%) | |
Current | 53 (12.2) |
Ever | 134 (30.7) |
Never | 249 (57.1) |
Atopy (%) | 198 (45.4) |
FEV1 pred (%) | 75.2 ± 18.0 |
FEV1/FVC (%) | 72.7 ± 13.4 |
All data are presented as the mean ± standard deviation (SD) or the number of participants (%). aparticipants available with this data were analyzed. FEV1 pred, predicted forced expiratory volume in 1 second; FVC, forced vital capacity. |
We performed the cis-eQTL analysis with RNA-seq of PBMCs from the asthma patients using FastQTL [39]. Following QC for genotype and gene expressions, 6,671,573 SNPs and 24,059 genes were chosen for the analysis. After applying FDR < 0.05 as an eQTL P-value threshold, 268,719 SNPs (eSNPs) had at least one target gene (eGene), and 2,931 genes (eGenes) had at least one expression-associated SNP (eSNP); consequently 342,054 eGene-eSNP pairs (cis-eQTLs) were discovered (Additional File 1: Table S1). Of the 268,719 eSNPs, 224,087 (83.4%) had one target gene, 30,143 (11.2%) had two target genes, 8,490 (3.2%) had three target genes, and 5,999 (2.2%) had more than four target genes (Fig. 1A). Among the 342,054 eGene-eSNP pairs, the distance from the TSS of the target gene to the eSNP was within 500 kb in 95.5% (326,493 pairs), within 250 kb in 85.3% (291,795 pairs), within 100 kb in 59.7% (204,359 pairs), and within 50 kb in 38.7% (132,445 pairs) (Fig. 1B). The 2,931 eGenes were composed of 2,604 protein-coding genes (88.8%), 268 pseudogenes (9.1%), 54 non-coding RNAs (1.8%), and 5 other types of genes (0.2%) (Fig. 1C).
We compared our blood eQTL results from asthmatic patients with healthy subjects. For this comparison, we used two blood eQTL datasets from a GTEx v7 study [41] and a Japanese study [42]. The total numbers of eGenes from the GTEx and Japanese datasets were 8,663 and 3,386, respectively. Among the 2,931 eGenes from the Korean asthmatic eQTL, 63.2% (1,853 of 2,931 eGenes) overlapped with the GTEx eGenes and 38.5% (1,128 of 2,931 eGenes) overlapped with the Japanese eGenes. For reference, 69.6% of the eGenes from the Japanese data (2,355/3,386) overlapped with the GTEx eQTL data. We also compared the direction of eSNP allelic effects on the target gene expressions between the asthmatic and normal eQTL data. There were 1,052,542 eGene-eSNP pairs in the GTEx and 335,813 eGene-eSNP pairs in the Japanese population. Among the 342,054 eGene-eSNP pairs from the Korean asthmatic eQTL data, 98,666 eGene-eSNP pairs were shared with GTEx and 75,745 eGene-eSNP pairs were shared with the Japanese pairs. Among the shared eGene-eSNP pairs, approximately 93.0% (91,721/98,666 pairs) for GTEx and 95.7% (72,449/75,745 pairs) for Japanese showed the same direction of allelic effects on gene expression as the Korean asthmatic eQTL (Fig. 2). As a reference, 142,923 eGene-eSNP pairs were shared between the GTEx and Japanese eQTL data, and 99.2% (141,842 pairs) showed allelic effects in the same direction (Fig. 2). As shown in Additional File 2: Table S2, we also found eQTL results showing allelic effects in the opposite direction, between the asthmatic and normal eQTL data as follows: 7.0% (6,945 pairs from 158 eGenes) of 98,666 pairs for GTEx, and 4.3% (3,296 pairs from 24 eGenes) of 75,745 pairs for the Japanese eQTL. Among the 158 eGenes from the GTEx data and 24 eGenes from the Japanese data, showing the opposite direction of the eQTL effect, 17 genes, including CEACAM21, NDE1, LGALS8, and AC124944.3, overlapped between them (Additional File 2: Table S2).
Identification of the asthma-susceptibility genes from the Korean asthmatic blood eQTL data using COLOC and SMR analysis
We performed a colocalization analysis to identify the target genes corresponding to the genetic variants from the asthma GWAS using the asthmatic eQTL data. We applied the COLOC method [15] to the three blood cis-eQTL datasets (Korean asthma, GTEx, and Japanese) and two asthma-GWAS summary statistics using European participants [4] and East Asian ancestries [51]. The European asthma-GWAS analyzed 64,538 asthma cases and 239,321 controls from the UK Biobank and identified 148 lead SNPs [4], whereas the Japanese asthma-GWAS analyzed 8,216 cases and 201,592 controls from a Japanese cohort and identified 11 lead SNPs [51] (Additional File 2: Table S3 and S4). We conducted the colocalization analysis by applying the criteria described in the Methods section and defined the GWAS colocalized genes with a posterior probability of colocalization (PP4) greater than 0.5.
We identified 46 genes in the Korean asthmatic eQTL data with evidence of colocalization between the eQTL and European asthma-GWAS signals (43 genes) or Japanese asthma-GWAS signals (3 genes) (Additional File 2: Table S5). In addition, 59 and 39 genes were identified as having significant colocalization, from the normal eQTLs of GTEx and the Japanese study, respectively, using both the European and Japanese asthma-GWAS summaries (Additional File 2: Table S5). Of the 46 colocalized genes identified in the Korean asthmatic eQTL data, 23 genes were also significant in the colocalization results of the GTEx or Japanese eQTL datasets, 13 genes were not available in the GTEx or Japanese eQTL datasets, and the remaining 10 genes were found to be significant only in the Korean asthmatic eQTL dataset (Table 2, Fig. 3, and Additional File 3: Fig. S1-S9). For example, an association signal around rs2070901 located on chromosome 1 showed a clear colocalization pattern with the FCRG3A gene only in the Korean asthmatic eQTL data (Fig. 3). In the normal eQTL data of GTEx and the Japanese study, strong eQTL lead variants were identified in the LD block, different from that of the asthmatic eQTL data, with modest GWAS P-values. As a reference, the colocalization results for GTEx and the Japanese study showed 17 overlapping genes (Additional File 2: Table S6). Among these 17 overlapping genes, 11 were also colocalized in the Korean asthmatic eQTL data. Among the remaining six genes, four were not available, and two did not show significant colocalization in the Korean asthmatic eQTL data.
Table 2
Korean asthma-specific colocalization genes
Study ID | Ensemble ID | Gene symbol | Var num | PP4 Korean asthma | PP4 GTEx | PP4 Japanese |
GCST010042 | ENSG00000124920 | MYRF | 796 | 0.921 | 3.25 × 10− 10 | 0.304 |
GCST010042 | ENSG00000185800 | DMWD | 1181 | 0.823 | 0.113 | 0.109 |
GCST010042 | ENSG00000186088 | GSAP | 1129 | 0.953 | 0.002 | 0.096 |
GCST010042 | ENSG00000173917 | HOXB2 | 961 | 0.91 | 0.116 | 0.074 |
GCST010042 | ENSG00000104147 | OIP5 | 991 | 0.745 | 0.211 | 0.069 |
GCST010042 | ENSG00000163608 | NEPRO | 1249 | 0.634 | 3.40 × 10− 6 | 0.063 |
GCST010042 | ENSG00000203747 | FCGR3A | 1281 | 0.648 | 9.39 × 10− 6 | 0.053 |
GCST010042 | ENSG00000197070 | ARRDC1 | 711 | 0.981 | 2.00 × 10− 4 | 0.016 |
GCST010042 | ENSG00000260261 | AC124944.3 | 1162 | 0.994 | 0.143 | 6.58 × 10− 4 |
GCST010042 | ENSG00000116688 | MFN2 | 1032 | 0.982 | 1.38 × 10− 8 | 1.13 × 10− 5 |
Study ID, ID of the paper which provided asthma-GWAS summary; Var num, number of genetic variants used in the colocalization analysis; PP4, posterior probability of colocalization; NA, not available because gene was out of analysis criteria or does not exist in the dataset; |
As another way of identifying asthma genes using the asthmatic eQTL data, we performed an SMR analysis using the asthmatic eQTL data and two asthma-GWAS summaries from Europe and Japan. SMR is generally used to discover genes satisfying the hypothesis that genetic variants affecting the risk of disease may act via changes in the expression levels of certain genes [16]. We applied the Bonferroni correction and HEIDI method for the SMR analysis, to reduce the likelihood of false-positive results.
After applying the Bonferroni correction (Psmr <2.08×10− 6) and HEIDI threshold (PHEIDI >0.01), a total of five SMR genes were identified: four SMR genes with the European asthma-GWAS summary and one SMR gene with the Japanese asthma-GWAS summary (Additional File 2: Table S7). We also performed an SMR analysis using the GTEx and Japanese eQTL datasets by applying the Bonferroni correction (Psmr <2.57×10− 6 for GTEx and < 2.27×10− 6 for Japanese) and HEIDI threshold (PHEIDI >0.01), and found 18 and one SMR genes, respectively (Additional File 2: Table S7).
Among the five SMR genes in the asthmatic eQTL data, three genes, DDX5, MED24, and NDFIP1, were not identified in the SMR analysis of the two normal eQTL data. In addition, these three genes were significant in both the colocalization and SMR analyses (Additional File 2: Table S5 and S7). In a previous SMR study [18], NDFIP1 was reported to be associated with asthma in the SMR analysis using skeletal muscle tissue (GTEx) and whole blood (eQTLgen) (Additional File 2: Table S8). However, the direction of the NDFIP1 effect on asthma was different: it was positive in the Korean asthmatic eQTL and GTEx, and negative in the eQTLgen. The other two genes, DDX5 and MED24, have not been reported to be significantly associated with asthma.