Cis- and Trans-Acting Expression Quantitative Trait Loci of Long Non-Coding RNA in 2,549 Cancers With Potential Clinical and Therapeutic Implications

Many cancer risk loci act as expression quantitative trait loci (eQTLs) of transcripts including non-coding RNA. Long non-coding RNAs (lncRNAs) are implicated in various human cancers. However, the pathological and clinical impacts of the genetic determinants of lncRNAs in cancers remain largely unknown. In this study, we performed eQTL mapping of lncRNA expression (elncRNA) in 11 TCGA cancer types and characterized the biological processes of elncRNAs in the setting of genomic location, cancer treatment responses, and immune microenvironment. As a result, 10.86% of the cis-eQTLs and 1.67% of the trans-eQTLs of lncRNA were related to known genome-wide association studies (GWAS) cancer risk loci. The elncRNAs are significantly enriched for those which are previously annotated as predictive of drug sensitivities in cancer cell lines. We further revealed the downstream transcriptomic effectors of eQTL-elncRNA pairs. Our data specifically suggested that the genes affected by eQTL-elncRNA associations are enriched in the immune system processes and eQTL-elncRNA associations influence the constitution of tumor infiltrating lymphocytes. In ovarian cancer, the “rs34631313-AC092580.4” pair was associated with increased fraction of CD8+ T cells and M1 Macrophage; whereas in KIRC, the “rs9546285-LINC00426” pair was associated with increased fraction of CD8+ T cells and a decreased fraction of M2 macrophages. Our findings provide a systematic view of the transcriptomic impacts of the eQTL landscape of lncRNA in human cancers and suggest its strong potential relevance to cancer immunity and treatment.


INTRODUCTION
Transcript abundance is an inheritable quantitative trait which serves as a major intermediate variable to explain the functional background of intergenic trait associated loci (TAL) (1). The germline variants which are associated with the transcript levels are known as "expression quantitative trait loci" (eQTL) (2). The eQTL may act either in-cis or in-trans, depending on the relative positioning between the loci and the target transcripts (1). The biological mechanisms underlying the functions of cis-and trans-acting eQTLs for mRNA have been addressed by a plethora of studies in various cell lines and tissues (3,4). Mapping of eQTLs in tumor tissues provides fundamental clues the functional implications of non-coding, risk-associated variants (5)(6)(7)(8). Quite a few cancer risk loci have been proven to be genetic determinants of gene expression. For example, the breast cancer risk loci 1p13.2 (rs11552449) acts in-cis to influence the activity of DCLRE1B, an evolutionarily conserved gene involved in DNA stability and the repair mechanism of inter-strand cross-link (9)(10)(11). Recent studies have reported associations between cancer risk loci and lncRNA expression and empirically verified the functional phenotypical effects with clinical significance (12).
Non-coding RNA (ncRNA) comprises diverse RNA transcripts including microRNA (miRNA), long non-coding RNA (lncRNA), circular RNA (circRNA), and piwi-interacting RNA (piRNA), many of which are known to function as regulators of transcription (13). The altered abundance of ncRNAs affects the cancer transcriptome and the consequent processes of carcinogenesis and tumor development (14). Among the ncRNA species, lncRNAs participate in diverse regulatory activities in the cell, including chromatin structuring and reprogramming cis-regulation of enhancers, competing endogenous RNA (ceRNA) networks, immune response, and post-transcriptional regulation of mRNA processing (15)(16)(17)(18)(19). The regulatory activities of lncRNAs play a critical role in the immune-related disease and the biological processes which determine the pathological and clinical phenotypes of cancers (20,21). Growing evidence indicated that dysregulation of lncRNAs involved in the regulation of immune system (22). Yongsheng Li et al. demonstrated the immune-associated lncRNAs (ImmLnc) show expression perturbation in cancers and are significantly correlated with immune cell infiltration (23), which suggested that lncRNA play an active role in cancer immunity.
Studying the genetic determinants of lncRNAs is crucial for understanding the etiology of carcinogenesis. For example, recent evidence has suggested that risk-associated loci in prostate cancer act in-cis on specific lncRNAs, thereby predisposing cells towards malignancy (24). However, eQTL mapping in cancer is more challenging than that in established cell lines due to tissue heterogeneity. In addition, other somatic changes in the genome, such as somatic copy number alterations and CpG methylation might confound the effect of germline determinants. To control for the confounders in the expression data, several studies used multivariate methods to adjust the eQTL mapping results for the effects of heterogeneity, sample purity, and somatic alterations on transcription activity (25). These methods have yielded valuable associations between germline variants and gene expression in cancer. Recent studies have reported eQTLs for various non-coding transcripts (lncRNA and miRNA) in human cancers (11,26). Moreover, study of Xia et al. demonstrated that lncRNAs are active participants in immune regulation in 33 cancer types (23). However, the extensive biological impacts of cancer eQTLs of lncRNA on the whole cancer transcriptome and the tumor microenvironment and the consequential therapeutic implications have not been thoroughly investigated. In addition, differences in the landscape of genetic determinants of lncRNAs and mRNAs in human cancers have not yet been fully addressed. Hence, we investigated the eQTL landscape of lncRNAs in eleven cancer populations from TCGA and revealed the potential pathological and clinical significance of elncRNAs in cancer immunity and treatment responses using Instrumental Variable Analysis (IV) or called Mendelian Randomization (MR) (27).
Our findings suggest that the downstream targets of eQTL-elncRNA pairs are enriched for immune system pathways and are consistently associated with varied fractions of immune cell types and patient clinical outcomes. In addition, we show that the elncRNAs are significantly enriched for known predictors of the treatment responses. Altogether, our data confirm that elncRNAs are active intermediates between genetic variants and the transcriptional activities in cancers, which further influence the tumor immune microenvironment, treatment responses and clinical outcomes.

Genotype Data Collection, Imputation, and Processing
To identify the eQTLs across eleven cancer types, we obtained genotype data from the TCGA portal (https://portal.gdc.cancer. gov/), which contains Affymetrix genome-wide human SNP 6.0 array-based genotype data for 898,620 single nucleotide polymorphisms (SNPs). To increase the power for eQTL discovery, we imputed the variants on autosomal chromosomes for all samples in each cancer type using IMPUTE2 algorithm with 1000 Genomes Phase 3 (28) as the reference panel. Then, a two-step procedure of IMPUTE2 were performed, which consists of pre-phasing and imputation of the phased data. After that, the following criteria were used to select SNPs: (i) imputation confidence score, INFO ≥0.5, (ii) minor allele frequency (MAF) ≥5%, (iii) SNP missing rate <5% for bestguess genotypes at posterior probability ≥0.7, (iv) Hardy-Weinberg Equilibrium P-value >1×10 −6 estimated using the Hardy-Weinberg R package (29).

Transcript Expression Data Collection and Processing
The lncRNA expression profiles of eleven TCGA cancer types were directly obtained from 'The Atlas of Noncoding RNAs in Cancer' (TANRIC) database (2018/06/11 -Version 1.3) (30). The expression levels of 12,727 lncRNAs were quantified as reads per kilobase per million mapped reads (RPKM). To ensure the detection power and reduce the noise, we applied a three-step filtering according to a previous study (31): (1) The lncRNAs with the 10 th -percentile RPKM value is equal to 0 were excluded, and those with the 90 th -percentile RPKM value is greater than 0.1 were retained for further analysis. (2) Quantile-normalization of the lncRNA expression matrix ensured that the underlying gene expression distribution was similar for all samples. (3) Finally, we applied an inverse normal transformation to the gene expression matrix, such that the expression levels of each lncRNA across samples was matched to the quantiles of standard normal distribution. The expression value of each lncRNA was also transformed based on log 2 (RPKM + 1).

Correction for Technical Confounders
The tumor microenvironment is the non-cancerous cells present in and around a tumor. These non-cancerous components of the tumor may play an important role in cancer biology. Therefore, we removed samples whose tumor purity less than 0.6 based on the criteria of prior study (32). In order to account for the possible latent confounding effects on transcript levels, we used the probabilistic estimation of expression residuals (PEER) method to estimate a set of latent covariates for the gene expression levels in each cancer type (25).

Covariates of Expression
Previous studies have shown that various factors affecting global gene expression confound eQTL-mapping (33). To exclude the effect of population structure on gene expression, we performed principal component analysis (PCA) for each cancer type and selected the CEU (Northern Europeans from Utah) population for further analysis. To regulate the other potential confounders, age, gender, gene-level somatic copy number alterations (SCNA), and CpG methylation levels were included as additional covariates. The copy number changes of a given lncRNA were determined by the segmented copy number scores of the tumor sample and paired-normal tissues obtained from the UCSC Xena database (https://xena.ucsc.edu/). Then, the CpG methylation status of the promoter area of each lncRNA (TSS 200/ TSS 1500) was determined by the respective discretized values obtained using Agilent HumanMethylation450K, with cutoff values 0.2 and 0.6.

Identification of eQTLs
For each cancer type, the genotype data, expression data, and covariates were processed to three N (genotype, expression, or covariates) × S (samples) matrices with matched sample identifiers. The lncRNA annotation file (hg19) was obtained from the GENCODE database (https://www.gencodegenes.org/). The SNP location (hg19) was downloaded from dbSNP148 (https://www.ncbi.nlm.nih.gov/projects/SNP/) (v148). eQTL analysis was performed according to the Genotype-Tissue Expression (GTEx) eQTL workflow (34) using the R package of "MatrixEQTL" in the linear regression model (35). The cis/trans-eQTL analysis was based on MatrixEQTL with PEER factors and other covariates: where ϵ i ∼ N(0, s 2 ) is a Gaussian error term; G i is the genotype of the i th sample; A i is the age of the i th sample; S i is the sex of the i th sample. In the cis-eQTL analysis, we evaluated the associations between lncRNA expression levels and the genotype at a given SNP locus located within 1 Mb upstream and 1 Mb downstream of lncRNA. trans-eQTLs are defined as associations beyond the 2-Mbp interval. The PEER factors were computed using the R package PEER to exclude the effects of confounding variables (such as the "batch" effects) from the gene expression matrix.
To reveal the complete landscape of eQTLs in cancer, we also performed trans-analysis for the lncRNAs. A false discovery rate (FDR) cutoff of 0.1 was applied for the significant cis-eQTL-elncRNA pairs, and, due to the vast number of tests performed in the trans-analysis, we set a threshold for the test P-values as 1 × 10 −7 for trans-eSNP-elncRNA pairs. In addition, Plink's clump function was used to shorten the list of cis/trans-eQTLs. Then, only the most significant cis/trans-eSNPs per haplotype block (r 2 =0.2, clump distance =500 kb) was selected to perform the second round of validation based on multivariate regression, wherein DNA methylation and SCNA were adjusted.
where ϵ i ∼ N(0, s 2 ) is a Gaussian error term; cis -G i is the genotype of the most significant cis-eQTL per LD block of the i th sample; trans -G i is the genotype of the most significant trans-eQTL per LD block of the i th sample; A i is the age of the i th sample; S i is the sex of the i th sample; SCNA i is the somatic copy number alteration of the i th sample; M i is the median of TSS200/ 1500 DNA methylation status of the nearest gene in the i th sample. A complete eQTL analysis results are available at https:// xcqxcq.github.io/lncRNA/LncRNA.html.

Functional Annotation of eQTLs
The eQTLs implicated in the analysis were functionally annotated by FUMA, a web-based bioinformatics tool that uses a combination of positional, eQTL and chromatin interaction mapping to prioritize likely causal variants and genes. The functional annotations included ANNOVAR categories and RegulomeDB scores. Moreover, the ANNOVAR categories identified the genomic positions of the SNPs (for example, intron, exon, intergenic) and the associated function. The RegulomeDB is a categorical score based on the information from normal tissue eQTLs and chromatin markers; 1a-7 with low scores indicate a decreased regulatory function. The scores are stated in (Supplementary Table S1).

Meta-Analysis of Epigenetic Markers
We used chromatin-immunoprecipitation sequencing (ChIPseq) peak data from the ENCODE database. For cis-eQTL overlapped with epigenetics modifications, Bedtools was applied (https://bedtools.readthedocs.io/en/latest/content/tools/ fisher.html) to calculate the enrichment of the cis-eQTLs for epigenetic markers. Then, meta-analysis was performed using the R package (meta) to integrate the fold-enrichment of each epigenetic marker.

Instrumental Variable Analysis
Instrumental variable (IV) analysis was employed to identify the cis-regulated lncRNAs that affect the transcription of proteincoding genes in-cis/trans according to the study performed by McDowell et al. (37). The genotype of an SNP as the genetic instrument across n individuals (eQTLs) directly affects the independent variable elncRNA in-cis, which in turn, affects the dependent variable mRNA in-cis/trans. Thus, the IV analysis was performed using the most significant eQTL-lncRNA pairs from the cis-analyses on the "AER" package in R-3.5.0, wherein the eQTLs are the genetic instrument, cis-lncRNA expression is the "mediator," and the cis/trans-mRNA expression is the "outcome" (38). We obtained the mRNA expression data of the corresponding TCGA samples from UCSC Xena (https://xena.ucsc.edu). To reduce the burden of multiple statistical tests, we filtered the mRNA in each cancer type by the following criteria: (1) Coefficient of variation (CV) >0.5, the 30 th percentile log 2 (count + 1) value >0 and the 90 th percentile log 2 (count + 1) >5.
(2). The absolute value of Spearman's correlation between lncRNA and mRNA >0.3. IV analysis of tumor infiltrating immune cells was conducted for the significant eQTL-elncRNA associations identified above. The immune cell fraction measures in TCGA cancer samples are based on published data from CIBERSORT (39). We picked the significant eQTL-elncRNA pairs (FDR<0.1, adjusted R 2 >0.1). Again, eQTL is considered as a genetic instrument and elncRNA as an independent variable. Since the dependent variable immune cell fraction is not normally distributed, it was log 2 (immune cell fraction + 1) transformed.

Enrichment Test of elncRNAs for lncRNA Predictors of Anticancer Drug Sensitivity
We obtained drug sensitivity data for lncRNAs from Nath et al.'s (21) study from https://osf.io/m2qja/. We analyzed the following datasets: Effects_CTRP.csv, Pval_CTRP.csv, Effects_GDSC.csv, and Pval_GDSC.csv. We performed Fisher-exact test to assess the enrichment of cis-elncRNAs/trans-elncRNAs for annotated lncRNA predictors of anti-cancer drug sensitivity. The magnitude of enrichment using odds ratios.

RESULTS eQTL Mapping for 11 TCGA Cancer Types
We performed pan-cancer eQTL mapping for lncRNAs (30) Table S2). The cis/trans-eQTL analyses were performed using a two-step regression model. In the first step, we used the MatrixEQTL software (35) to screen for significant expression SNPs (eSNP) (FDR<0.1). Next, to account for linkage disequilibrium (LD), we collapsed the significant eSNPs by the linkage disequilibrium into independent eQTL (r 2 =0.2, clump distance =500 kb). Each eQTL was represented by a tag-eSNP with maximal significance of association. Finally, we verified the significant cis/trans-eQTLs using multivariate model adjusting for somatic copy number alteration (SCNA) and DNA methylation in poise cis-regulatory regions (TSS-200/TSS-1500, Figure 1A).
A majority of the cis-eQTLs (65.86%) were cancer type specific. However, a set of cis-eQTLs were noted in multiple cancer types. For example, 2.4% cis-eQTLs were present in more than eight cancer types, which are termed as "pan-cancer" cis-eQTLs hereafter ( Figure 1D). As for the elncRNAs in-cis, 1,150 (24.25%) were cancer-specific, and 237 (5.0%) occurred in more than eight cancer types ( Figure 1D).
In the trans-eQTL analysis, we report 10,721 significant trans-eQTL-elncRNA association pairs with genome-wide P-values < 1 × 10 −7 ; these pairs were mapped to 9,996 unique trans-acting eQTLs and 3,284 unique elncRNAs ( Figure 1C and Table 1). Among the trans-eQTLs, 91.26% were cancer-specific, only 0.5% of the trans-eQTLs were observed in seven cancer types, suggesting that the trans-eQTLs of lncRNA are more cancer-specific ( Figure 1D). As for the elncRNAs in-trans, 1,730 (52.67%) of elncRNAs in-trans were cancer-specific ( Figure 1D), and only 1 (0.03%) elncRNA in-trans (RP11-667M19.2) were observed in eight cancer types. Together, these results suggested that that the overall eQTL landscape of elncRNAs are highly specific to the cancer types. However, trans-eQTL and their elncRNAs were more cancer-specific than their counterparts in-cis.
The eQTL-elncRNA association in cancers features a manyto-many relationship, where one lncRNA can associate with different eQTLs or one eQTL can associate with different elncRNAs. For example, a total of 3,107 elncRNAs were associated with both cis-and trans-eQTLs in the same cancer type (Supplementary Table S4). Moreover, a total of 1,039 eQTLs act both in-cis and in-trans on different lncRNAs in the same type of cancer (Supplementary Table S5). These results suggested that a complex network is involved in the genetic determinants of lncRNA transcription.
Consistent with previous eQTL studies for mRNA (40), the number of significant cis-eQTLs for lncRNA identified in the current study was strongly associated with the size of the cohort (Spearman's correlation coefficient, rho=0.87, P=0.00095, Figure  1E). A similar tendency was observed in trans-eQTLs for lncRNA (rho=0.79, P=0.0061, Spearman's correlation, Figure 1E).   Effect Size of the Determinants of lncRNA Transcript Levels The fractions of variation attributed to the major factors of lncRNA expression, including age, sex, PEER factors, gene-level somatic copy number, DNA methylation and cis-/trans-eQTLs are calculated (Figure 2A). A total of five PEER factors were selected to remove the latent confounding effects and maximize the cis-elncRNA discovery (Supplementary Figure S2). In the 11 cancer types analyzed, cis-eQTLs accounted for 7.52% (OV) to 20.15% (STAD) of the total variance in the lncRNA transcript levels, whereas trans-eQTLs accounted for 2.16% (LUAD) to 27.76% (STAD) (Supplementary Figure S3). The other factors, such as SCNA, explain 1.09% to 11.82% of the total variance, which is much larger than the effect size of CpG methylation (0.01-0.42%) (Supplementary Figure S3). However, 9.05% to 26.20% of the variation in lncRNA expression was attributed to non-specific, random effects represented by the PEER factors (Supplementary Figure S3). Finally, the effects of age and sex were much smaller than the other factors ( Figure 2A). Overall, the fractions of variation in lncRNA expression explained by the major factors in cancers are highly similar to that of mRNA, as shown by prior studies (40).
As for the transcript-wise effect sizes, on average, 10.93% of the variance of the elncRNAs expression was attributed to cis-eQTLs, compared to 4.32% attributed to trans-eQTLs ( Figure  2B). Notably, the average effect size of cis-eQTLs of lncRNA is twice as much as that of the mRNA (5.3%, based on prior eQTL data) (40). For the majority of the elncRNAs, 85.50% cis-and 98.16% trans-eQTLs accounted for less than 20% of the total variance. However, for a small fraction of the elncRNAs, the effect of genetic determinants mounts to 69.72% (cis) and 75.34% (trans) ( Figure 2B). Furthermore, the average allelic effect size of cis-eQTLs of lncRNA (fold-change of elncRNA with every unit increase of the effect allele) was 0.54 ( Figure 2C), which is larger than that of mRNA based on prior study (0.37, Supplementary  Table S6) (40). For all elncRNA in cancers, the allelic effect sizes of the cis-eQTLs ranged from 0.002 to 2.71, compared to 3.0×10 −6 to 3.04 for the trans-eQTLs ( Figure 2C). Taken together, these results suggest that the effect size of genetic determinants on the expression of lncRNA is larger than that of mRNA.

Functional Characteristics of eQTLs of lncRNAs in Cancers
The majority of genetic determinants of transcripts are localized in intergenic or intronic regions. Thus, we sought to reveal functional characteristics of the eQTLs of lncRNA from three different aspects: the genomic location, the epigenetic landscape and the association with cancer risk-associated loci.
The genomic locations of the cis-eQTLs and trans-eQTLs were annotated by Functional Mapping and Annotation of Genome-Wide Association Studies (FUMA) (41). The results demonstrated that the majority of cis-eQTLs of lncRNA are located in the intergenic regions (45.25%) and intronic regions (25.91%) ( Figure 3A). Only a small fraction of cis-eQTLs was localized in the exonic (0.74%) or ncRNA exonic (3.66%) regions. The genomic distribution of cis-eQTLs of mRNA was similar to that of the cis-eQTLs of lncRNA (Supplementary Figure S4A). On the other hand, 62.42% of the trans-eQTLs of lncRNA were localized in the intergenic regions and 23.27% in the intronic regions ( Figure 3A); while some were located in the exonic (0.44%) or ncRNA exonic (1.01%) regions. Overall, the genomic distribution of the eQTLs is highly comparable between lncRNA and mRNA. A number of non-coding QTLs interfered with the regulatory elements of the genome (42). The present study showed that 38.99% of the non-coding cis-eQTLs of lncRNAs were colocalized in the transcription factor (TF) binding sites and/or DNase I hypersensitive sites (RegulomeDB score <6) Figures 3B, C and Supplementary Table S1. Whereas 48.63% of the noncoding cis-eQTLs of mRNA were co-localized in the TF binding sites and/or DNase I hypersensitive sites, and 4.95% had a putative regulatory function (Supplementary Figure S4B).
Since chromatin structures are tightly associated with the cisregulatory activities, we further evaluated the histone modification landscape of the eQTLs of lncRNA in matched cancer cell lines ( Figure 3D) and normal tissues (Supplementary Figure S5A). As a result, the eQTLs for lncRNA in the matched cancer cell lines were significantly enriched in the active promoter or enhancer markers, such as, H3K4me3 (OR=2.78, 95% CI: 1.60 to 4.83) and H3K27ac (OR=1.91, 95% CI: 0.93-3.92). The eQTLs were also negatively enriched in transcription repressor markers, such as H3K27me3 (OR=0.31, 95% CI: 0.13-0.77) ( Figure 3D). In the matched normal tissues, the eQTLs of lncRNA showed a consistent tendency of enrichment for the histone markers (Supplementary Figure S5A). About 10.86% independent cis-eQTLs of lncRNAs was intercepted with known cancer risk loci (r 2 >0.5) ( Figure  3E), which is much higher than that of cis-eQTLs of mRNAs (3.07%) reported recently (40), suggesting lncRNA is strongly related to GWAS cancer risk loci (Supplementary Table S6 and Figure S4C). Herein, the cis-eQTLs of lncRNA were significantly enriched for GWAS cancer risk loci than for non-eQTLs (fold of enrichment =3.51, P<2.2 × 10 −16 ) (For example, in ER-positive BRCA, we identified 12 cis-eQTLs which were in LD (r 2 >0.5) with breast cancer risk loci of rs62073257 (17q21.31), rs7104902 (11p15.5), and rs9393716 (6p22.2). In the case of trans-eQTLs, about 1.67% independent trans-eQTLs were intercepted with known cancer risk loci ( Figure 3E) which corresponds to 7.74-fold of enrichment (P<2.2 × 10 −16 ) compared with non-eQTLs. In the case of TAL, we observed that an average of 40.82% of cis-eQTLs of lncRNA were overlapped with TAL, while 27.89% of that of mRNA were intercepted with TAL (Supplementary Table S6).
Recent studies suggest that lncRNA are associated with the treatment responses of cancers (21). Therefore, we evaluated the predictive power of the elncRNA in public databases contains drug IC50 for various cancer cell lines (21). As a result, both   (Figures 4C, E). In particular, we reported top drugs from each database with the highest significance of association to elncRNAs, namely, PX.12, decitabine, cytabarine hydrochloride in CTRP, and PIK-93 and I-BET-762 in GDSC (Figures 4D, F). There are also drugs which are specifically associated with elncRNAs in-cis, such as CH542802 (ALK inhibitor) in GDSC, and we found that except for CH542802 (ALK inhibitor), both the elncRNA incis/trans were significantly enriched in the predictive lncRNA sets of PX.12, decitabine, cytabarine hydrochloride in CTRP, and PIK-93 and I-BET-762 in GDSC (Figures 4D, F).
As prior studies show lncRNA strongly influence cancer immunity, we move on to evaluated the enrichment of elncRNAs in the annotated immunoreactive lncRNA database (23). As a result, both elncRNAs in-cis/trans were significantly overrepresented in the ImmLnc database (FDR<0.01). Moreover, the magnitude of enrichment of the cis-elncRNAs range from 1.20 fold (STAD) to 1.44 fold (THCA), while the enrichment for the trans-elncRNAs range from 1.18 fold (KIRC) to 1.47 fold (COAD) (Supplementary Figure S5B).

The Downstream Transcripts of eQTL-elncRNA Influence Tumor Immune Microenvironment
We further investigated the impact of eQTLs-lncRNA associations on the downstream transcripts and the consequent pathological clinical phenotypes, such as the tumor immune microenvironment and the disease outcomes.
We conducted IV analysis to identify the eQTL-elncRNA-mRNA regulatory axes in each cancer type. The IV regression takes the eQTL as an instrument to the independent variable (elncRNA in-cis), which in turn, affects the dependent variable (mRNA). Thus, we identified a total of 191 unique regulatory axes in 11 cancer types. (FDR<0.1; adjusted R 2 >0.1, Supplementary Table S7 and Figure 5A). We then evaluated the biological processes enriched in the downstream transcripts of the significant regulatory axes. As a result, we noticed that a set of the downstream target genes such as IFNG, CALCA, CXCL1, NLRP6, BLK, CD79, FASLG, CCR4, GZMM, were significantly enriched in the immune system process and immune response pathways (FDR<0.01 Figure 5B). This result prompted us to further study the impacts of eQTL-elncRNA on cancer immunity. We evaluated the effects of the 191 eQTL-elncRNA pairs on the fractions of tumor infiltrating immune cells in each TCGA cancer type estimated by prior studies (39). As a result, we identified 23 eQTL-elncRNA pairs in different cancer types (10 elncRNAs in 5 cancer types) which were significantly associated with the fraction of at least one immune cell subtype (P<0.05 and R 2 >0, Supplementary Figure S7). The eQTL-elncRNA pairs showed two opposite ways of impact (Supplementary Figure  S7). In one way, AC092580.4 in OV and LINC00426 in KIRC were associated with lower fraction of Macrophage M2 and higher fraction of Macrophage M1 and CD8+ T cells (Supplementary Figure S7). M1 macrophages and CD8+ T cells are well known for their anti-tumoral effect, while M2 macrophages are reported as anti-inflammatory and associated with pro-tumor phenotypes (43). In the other way, RP4-756H11.3 in LIHC and RP11-677M14.7 in THCA were associated with increased fraction of M2 macrophages and decreased percentage of CD8+ T cells hence a pro-cancer effect (Supplementary Figure S7).
In another case, rs9546285 (13q12.3) was a robust genetic instrument for LINC00426 and thereby influences an immune gene set of TNIP3 (b=1. 58 Figure S6). These genes are involved in the antiviral and antitumor effects (48,49). In KIRC, the same pair is associated with higher fraction of CD8+ T cells (b= 0.72, P = 6.34 × 10 −4 ) and lower fraction of M2 macrophages (b = −0.39, P = 0.04) ( Figures 5G, H). These results suggest that eQTLs of lncRNA influence the immune processes by regulating relevant genes and thus influence the constitutions of the tumor immune microenvironment.

DISCUSSION
The present study aimed to provide an in-depth view of the eQTL landscape of lncRNA in cancers by systematic assessment the downstream effects of the eQTL-elncRNA associations (40). Our data suggested that the genetic determinants of lncRNA expression are more active in cancer than those of mRNA, as the latter is directly associated with cancer phenotypes hence subject to higher selection pressure. Nevertheless, our results confirmed that the eQTL-elncRNA associations show more diverse but indirect effects on the phenotypes via the transcriptome such as the treatment responses (50) and immune microenvironments (51,52). While mapping the eQTLs in cancer, we accounted for different confounding effects (such as PEERs, SCNA, DNA methylations). Thus, our model is more stringent than those reported previously (11,26), and the eQTLs directly interfered with the regulome of cancers. In addition, we demonstrated that the genomic landscape of cis-eQTLs of lncRNA was similar to that of the cis-eQTLs of mRNA. However, the cis-eQTL of lncRNA was greater than that of mRNAs in both cis-acting gene numbers per eQTL and effect sizes, which suggested that lncRNA expression is more susceptible to the genetic regulators than mRNA expression in cancers. Intriguingly, a number of cancer-related studies have demonstrated that germline risk variants effectuate via specific lncRNAs to impact tumor phenotypes (24,53,54). Such differences might be attributed to the different selection pressure in mammals (55) or different regulatory mechanisms underlying these two types of transcripts (56). Furthermore, the current data showed a significant higher enrichment of cancer risk loci in the eQTLs of lncRNAs than in those of mRNA in cancer, which is consistent with a previous study (26).
LncRNAs are known for a variety of regulatory activities in gene expression, many of which might exert functional and phenotypical impact in cancers. Therefore, we hypothesized that these elncRNAs can further interact with various downstream mRNA targets to influence the whole cancer transcriptome. Therefore, IV analysis was performed to identify regulatory axes in different cancer types. Especially, we report that the downstream transcripts of the elncRNAs in cancers were significantly enriched for the immunoregulatory processes. These results suggested that the eQTL act on immune genes in trans via the elncRNA. A previous study of cis-eQTL analysis of mRNA in cancers has also demonstrated that the cis-eQTL regulated mRNA are enriched in immunity pathways (40). We further conducted a second round of IV analysis, and 23 eQTLs-elncRNAs pairs were found to be significantly associated with at least one immune cell type. We closely analyzed the top representative eQTL-elncRNA pairs to determine the target immune genes and immune cell types of these pairs. These observations suggest germline variants have a strong impact on cancer immunity and may serve as predictive markers for immune therapy.
Notably, the top immune related eQTL-elncRNA pairs reported in the present study are also reported for association with drug sensitivity in prior independent studies (21). For example, LINC00426 is a strong predictor for Dexamethasone sensitivity in CCLE cell lines (21). Dexamethasone is a glucocorticoid (GC) steroid that is used as a supportive care co-medication for cancer patients undergoing standard pemetrexed/platinum doublet chemotherapy. It is used to reduce inflammation and suppress the body's immune response by inhibition of IL-2, IL-12, and IFN-g of signaling activities in NK, Th1, and CD8+ T cells (57). LINC000426 is associated with IFNG transcription and the activities of tumor infiltrating CD8+ T cells, hence LINC000426 expressing tumors are more sensitive to Dexamethasone. Since rs9546285 is a robust genetic regulator of LINC00426, this eQTL-elncRNA pair can inform the efficacy of Dexamethasone in cancer patients and also for cancer immunotherapy.
In addition, we identified several regulatory axes involving eQTL, elncRNA and mRNA by taking eQTLs as instrument variables. These regulatory axes cannot be easily identified by conventional trans association analysis due to the huge number of tests needed and the stringent control of false discovery rate. Many studies have shown that instrument variable regression, which is also known as Mendelian Randomization, is more efficient in identification of potential causal relationships between germline variants and different traits.
In certain cancer types, different eQTL-lncRNAs pairs influence the same set of mRNA targets; for example, in ovarian cancer, two eQTL (2p25.2 and 7q34) eventually regulate the same set of mRNAs (FASLG, GZMM, PYHIN1 and TRAT1) but through different elncRNAs (AC092580.4 and TRBV11-2). These results suggested that lncRNA act as a flexible intermediate between germline variants and gene expression intrans and hence may play more active roles in cancers.
Nevertheless, the present study has several limitations. First, the inferred interaction patterns were based on multi-omics data, and further cellular and molecular experiments are needed to validate the findings. Then, the results were derived from 11 cancer types, and, thus, our findings may not apply to other cancer types. Third, the mapping of eQTLs of lncRNA was limited by the sample size, missing data, and confounding factors. The current results for the trans-eQTL hold limited statistical power due to a relatively small sample size (2549 samples). As a consequence, many eQTLs with small effect sizes cannot be identified. Also, we used the IV analyses to identify several regulatory axes with functional implications. However, these associations can occur in different cell types in the tumor tissue, or even via more complex intercellular interactions (58). Further functional analysis is needed to reveal the underlying biological mechanisms.
Taken together, the current findings provided insight into the genetic regulation of lncRNAs in 11 cancer types and explored the biological role as well as clinical phenotypes of eQTL-elncRNA in cancer immunology. Our findings may provide valuable genetic or lncRNA biomarkers for drug sensitivity and cancer immune therapy.

CONCLUSIONS
This study investigated the eQTL landscape of lncRNAs in cancers and the potential biological function of elncRNAs in cancer microenvironment and drug sensitivity prediction. We performed instrumental analysis (Mendelian Randomization) to identify genes, pathways and immune cell types influenced by eQTL-elncRNA associations.
Our findings suggest that the downstream targets of eQTL-elncRNA pairs are enriched for immune system pathways and are consistently associated with varied fractions of immune cell types and patient clinical outcomes. Our data confirm that elncRNAs are active intermediates of non-coding genetic variants in cancer immunology and clinical outcome of cancers and provide valuable genetic and lncRNA biomarkers for drug sensitivity and cancer immune therapy.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: The Cancer Genome Atlas (https://portal.gdc. cancer.gov/). The original contributions presented in the study are publicly available. This data can be found here: https:// xcqxcq.github.io/lncRNA/LncRNA.html.