Copy number variation is highly correlated with differential gene expression: a pan-cancer study

1Genetic structural variation in the human genome can be present in many forms, ranging from single nucleotide polymorphisms (SNPs) to large chromosome aberrance. In the past, SNPs are regarded as the predominant form of structural variation and account for much phenotypic variation. However, recent studies show the widespread existence of copy number variation (CNV) in individuals, and since that these observations have been extremely appreciated and expanded. In general, CNV is defined as an amplifying or decreasing number of DNA segments that is 1 kb or larger in the human genome, which accounts for an important part of genetic structural variation. Currently great efforts in science community have been directed to catalog and characterize somatic CNV in a comprehensive manner, which provides key knowledge on how they impact biological function, evolution and human diseases on genomic level.

2It is generally accepted that somatic CNV is highly associated with the development and progression of numerous cancers by impacting gene expression level. Samulin Erdem et al. found that Neurofascin (NFASC) gene is significantly amplified and overexpressed in non-small cell lung cancer (NSCLC) patients and the novel role of NFASC is identified in the regulation of cell motility and NSCLC migration. Dong et al. analyzed the copy number alterations and differentially transcribed genes in esophageal cancer and observed a noteworthy association between CNV and differential gene expression for FAM60A, TFDP1, CDC25B and MCM2. Subsequently, FAM60A was identified as a potential prognostic factor with a striking correlation to overall survival and clinical-pathological parameters. Lines of evidence support differential gene expression might be a vital intermediate mechanism for CNV to exert effect on the downstream phenotype.

3Despite a number of studies have explored CNV and differential gene expression of several classical oncogenes or tumor suppressor genes in different cancers, there has been no systematic study about the relationship between CNV and differential gene expression across a broader spectrum of cancer types and cell lines. It is unclear to what extent the expression level is affected by CNV for the whole genomics. Previous observations from single gene or single cancer type may not be representative for other genes or other types of cancer. Here we aimed to systematically investigate the specific relationship between somatic CNV and differential gene expression across cell lines and different cancer types for known genes. This study may help us better understand the correlation between CNV and differential gene expression and provide new insights into the mechanism of development and progression of cancer.

4The copy number and mRNA expression data of Broad-Novartis Cancer Cell Line Encyclopedia (CCLE), NCI-60 and the Cancer Genome Atlas (TCGA) were collected from the cBioportal for Cancer Genomics (http://www.cbioportal.org/). The cell lines datasets of CCLE and NCI-60 contained 966 and 59 cell lines respectively and the TCGA datasets contained 31 types of cancer of 9159 samples (see Additional file 2: Table S1). Putative copy number calls were determined by using GISTIC 2.0, while expression levels were quantified by RSEM from RNA-Seq data for TCGA. Another independent dataset were curated from COSMIC Cell Lines Project (CCLP, v81, http://cancer.sanger.ac.uk/cell_lines) containing 1020 cell lines. The copy number was obtained from Affymetrix SNP6.0 array data with PICNIC.

5For CCLE, NCI60 and TCGA datasets, gene-wise homozygous deletion or high level amplification were regarded as copy number amplified or deleted gene (Copy number values: − 2 = homozygous deletion; − 1 = hemizygous deletion; 0 = neutral / no change; 1 = gain; 2 = high level amplification). For CCLP, copy number amplification was obtained by the following criteria: (the average genome ploidy < =2.7 AND total DNA segment copy number > =5) OR (average genome ploidy > 2.7 AND total DNA segment copy number > =9). While the criteria for copy number deletion was: (the average genome ploidy < =2.7 AND total DNA segment copy number = 0) OR (average genome ploidy > 2.7 AND total DNA segment copy number < (average genome ploidy – 2.7)). Gene expression levels were quantified by RSEM from RNA-Seq data and mRNA Z scores were computed using the tumors samples that are diploid for the corresponding gene. For each gene, Z score = (x-u)/o, where u, o represent the average expression and standard deviation of this gene across samples, respectively; x represents the specific expression of this gene in a specific sample. Differentially expressed genes (DEGs) were further filtered out as Z scores more than 2 (upregulated genes) or less than − 2 (downregulated genes). Thus, the four variation tendencies were defined as follows: amplified AND upregulated; amplified AND downregulated; deleted AND upregulated; deleted AND downregulated.

6copy number amplification and expression level upregulation ratio against copy number deletion and expression level downregulation, copy number deletion and expression level downregulation ratio against copy number amplification and expression level upregulation more than 50 % on each gene across 9159 tumor samples were applied to identify AUGs and DDGs with a higher ρ (> 0.4) and a higher number (> 146.5) of copy number amplification and expression level upregulation than the median level of 30 most popular oncogenes (Additional file 2: Table S4) or a higher ρ (> 0.41) and a higher number (> 18.5) of copy number deletion and expression level downregulation than the median level of 10 tumor suppressor genes (Additional file 2: Table S5). AUGs and DDGs were sorted by the amount of copy number amplification and expression level upregulation, copy number deletion and expression level downregulation respectively and 31 representative AUGs and 29 representative DDGs matched with KEGG genes were identified from top 100 highly concordant genes.

7Search Tool for the Retrieval of Interacting Genes (STRING, http://www.string-db.org/) was used to construct the PPI network of FYTTD1. The line thickness indicates the strength of data support from the sources of text mining and experiments with a cutoff value of medium confidence (0.4). Then the functional enrichment results of KEGG pathways and Gene Ontology (GO) biological process were applied to the PPI network with false discovery rate less than 0.05.

8R (version 3.4.2) including R packages of data.table for data cleaning and management, survival for survival analysis, ggplot2 for data visualization and GraphPad Prism 7.0 were used for the statistical analysis. ρ is equal to the Pearson correlation between the rank values of those two variables to assesses how well the relationship between two variables can be described using a monotonic function and were calculated by the function of cor in R. Differences between two groups were determined using the Welch’s t-test (significant with p < 0.05).

9It is known that CNV is related to alteration on gene expression. However, the correlations between CNV and gene expression change on a global cellular scale remain to be elucidated. To comprehensively exploit this question, we collected available genomic datasets from the Broad-Novartis Cancer Cell Line Encyclopedia (CCLE), NCI-60, and the Cancer Genome Atlas (TCGA). Integrated analysis of CNV and differential gene expression was performed across 31 cancer types and 2 cancer cell lines resources (Additional file 2: Table S1). Our results showed an apparent effect of copy number on gene expression (Fig. 1). Interestingly, gene-wise copy number amplification appears to harbor a higher median Z score, indicating enhanced gene expression, while copy number deletion usually led to decreased gene expression among 966 cell lines for CCLE datasets (Fig. 1a). Even for NCI-60 project and mesothelioma (MESO) in TCGA, which have smaller sample size, this consistent tendency still is kept (Fig. 1a; results of other 30 cancer types in TCGA found in Additional file 1: Figure S1 and Figure S2). For all cancer types under study, median Z scores of gene expression with copy number amplified were strikingly higher than those with copy number deleted (Fig. 1B). In order to further detect the relation between copy number and transcript level, the linear regression model was taken for each cancer type of TCGA datasets. Fitting results indicated gene expression level was prominently linear-correlated with copy number (r = 0.93, p < 0.0001; Fig. 1c). In addition, the spearman’s correlation coefficient (ρ) test and linear regression analysis among cell lines and tumors were implemented to uncover the close association between differential gene expression and CNV with ρ ranging from 0.075 to 0.53 (mean r = 0.97; Fig. 1d), indicating the positive correlation of differentially transcriptional level and CNV in general.

10We next analyzed the relationship between copy number and expression level at the basis of individual gene. Remarkably, for substantial number of genes, their copy numbers exert a positive correlation with the corresponding gene expression level among 1025 cell lines dataset (88.11%; Additional file 1: Figure S3A; Figure S3B) and 9159 tumor samples (94.15%; Fig. 2a). This trend is particularly obvious for a vast number of oncogenes (e.g. MYC, AKT1, CDK9, KRAS, and MDM2) and tumor suppressor genes (e.g. CDKN2A, RB1, PTEN, and TP53). Linear regression was taken to explore whether copy number had a positive linear influence on gene expression. Plenty of genes presented a high degree of fitting with r intensively ranging from 0.8 to 1 (Fig. 2b; Additional file 1: Figure S3C), genomic structural variation may generate a direct effect on gene transcriptional level for the majority of genes. Concordantly, a recent study by GTEx consortium also associated genetic variants with gene expression levels across 44 human healthy tissues and gene expression levels were found to be affected by local genetic variation for most genes based on eQTL analysis. In addition, we selected 10 genes with a significant correlation between CNV and differential gene expression through TCGA and cell lines datasets and validated this relationship through literature mining. For 9 genes except the weak evidence for NFASC, a strong correlation between CNV and differential gene expression was confirmed combining ρ and Pearson’s r of fitting (Table 1), in accordance with previous experimental findings.

11Interestingly, our analysis shows a fraction of genes present low Pearson’s r of fitting and low ρ (Additional file 1: Figure S3D; Figure S3E), suggesting stable expression of these genes despite of CNV. To gain further insight into those genes (top 1000 genes with Pearson’s r of fitting from the lowest to the highest), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was carried out by utilizing the Database for Annotation, Visualization and Integrated Discovery (DAVID). Results demonstrated a strong enrichment of retinol metabolism, olfactory transduction, calcium signaling pathway, neuroactive ligand-receptor interaction, etc. (Additional file 2: Table S2; Table S3). Moreover, for a total number of 16,639 shared genes from cell lines and TCGA datasets, a high level of agreement on ρ was observed (Fig. 2c). Cell lines and tumors seem to exert similar trend. Wherein, 86.42% of shared genes showed a double-positive ρ among cell lines and TCGA datasets, while genes with a double-negative ρ only take up 1.63% in shared genes.

12Based on two variation trend of copy number (amplification and deletion) and expression level (upregulation and downregulation), we divided variant trend into four classes: copy number amplification and expression level upregulation, copy number amplification and expression level downregulation, copy number deletion and expression level downregulation, copy number deletion and expression level upregulation (see Methods). We analyzed the frequencies of four variant trends across 24 primary sites of cell lines in CCLE database and 31 cancer types in TCGA. Intriguingly, the amount of copy number amplification and expression level upregulation, copy number deletion and expression level downregulation, representing our so-called highly concordant genes on CNV and differential gene expression, occupied 15 and 20% of the total variant copy number count on average for cell lines and tumor samples, respectively. Nevertheless, the mean proportion of copy number amplification and expression level downregulation, copy number deletion and expression level upregulation count only took up 0.7 and 0.5% for both cell lines and tumor samples (Additional file 1: Figure S4A and S4B), which indicated that the copy number amplification barely causes gene expression downregulation and the copy number deletion hardly promotes gene expression upregulation. It is obvious that the frequency of copy number amplification and expression level upregulation totally exceeded copy number deletion and expression level downregulation for both cell lines and tumor samples (Fig. 3a and b), which may be a result of selection on deletions.

13Gene names in literature were transformed into official symbols. ρ: spearman’s correlation coefficient.

14In order to identify highly concordant genes on CNV and differential gene expression, the sum of copy number amplification and expression level upregulation, copy number deletion and expression level downregulation for all genes across 9159 tumor samples in TCGA were counted (Fig. 3c). The red points close to x-axis owning a high count of copy number amplification and expression level upregulation and a low count of copy number deletion and expression level downregulation are preferentially defined as copy number amplified and expression level upregulated genes (AUGs), while the blue points close to y-axis owning a high count of copy number deletion and expression level downregulation and a low count of copy number amplification and expression level upregulation are regarded as copy number deleted and expression level downregulated genes (DDGs). Moreover, we observed that a majority of genes (83%, ratio < 0.1) belongs to AUGs (Fig. 3d). It may be caused by unknown selective pressures on amplification. We attempt to revalidate the ten highly concordant genes through literature mining (9 AUGs, 1 DDG; shown in Table 2). Nine AUGs obtained preponderant count of copy number amplification and expression level upregulation compared to copy number deletion and expression level downregulation across 9159 tumor samples, while the only DDG (STK11) showed an overwhelming degree of copy number deletion and expression level downregulation (73%) against copy number amplification and expression level upregulation, in accordance with previous reports in literature (Table 2).

15Gene names in literature were transformed into official symbols. A&U represents the frequency of a gene with copy number amplified and expression level upregulated across 31 cancer types in TCGA, while D&D represents the frequency of a gene with copy number deleted and expression level downregulated. ↑ means the gene with copy number amplification and expression level upregulation in tumor, AUG; ↓ means the gene with copy number deletion and expression level downregulation in tumor, DDG.

16Thus, to identify AUGs and DDGs, the criteria applied in this work is the copy number amplification and expression level upregulation ratio against copy number deletion and expression level downregulation, copy number deletion and expression level downregulation ratio against copy number amplification and expression level upregulation with the cutoff value of 50 %. Further by filtering parameter included a higher ρ and a higher number of copy number amplification and expression level upregulation than the median level of 30 most popular oncogenes (Additional file 2: Table S4) or a higher number of copy number deletion and expression level downregulation than the median level of 10 tumor suppressor genes (Additional file 2:Table S5), which ultimately led to 560 AUGs such as DERL1, DVL3, FADD and 365 DDGs (e.g. MTAP, KLHL9, PTEN) (complete list in Additional file 2: Table S6). For representative AUGs and DDGs matched with KEGG pathway-related genes, their aberrant rate of copy number amplification and expression level upregulation, copy number deletion and expression level downregulation in corresponding tumor samples across 31 cancer types were analyzed and shown in Fig. 4a and Fig. 4b, respectively. Taken the AUG DERL1 as an example, this gene is amplified and overexpressed in a number of cancers, such as breast invasive carcinoma (BRCA, 16.84% patients), esophageal carcinoma (ESCA, 12.02% patients), liver hepatocellular carcinoma (LIHC, 10.17% patients), uterine carcinosarcoma (UCS, 10.71% patients), uveal melanoma (UVM, 16.25% patients), especially in ovarian serous cystadenocarcinoma (OV) with copy number amplification and expression level upregulation rate reaching 32.67% in 300 OV patients (Fig. 4a). Previous studies demonstrate DERL1 overexpression is significantly related to cancer cell proliferation, invasion and poor prognosis. Our study suggests this dysregulation may be driven by copy number amplification based on the extreme concordant correlation between copy number and gene expression. Figure 4b shows the distribution of copy number deletion and expression level downregulation rate across 31 cancers of 38 representative DDGs involving in many known tumor suppressor genes such as MTAP, KLHL9, PTEN, SMAD4, RB1, etc.

17Additionally, we analyzed the distribution of AUGs and DDGs across 22 chromosomes (Additional file 2: Table S7). The maximal proportion of highly concordant genes of AUGs and DDGs (21.62%) were located in chromosome 8, followed by chromosome 1 (16.54%), which is in coherence with the previous finding published on nature that chromosome 1 possesses the largest number of genes exerting a strong relevance to numerous diseases such as cancer, Alzheimer’s disease, Parkinson’s disease. As expected, the smallest chromosome 21 acquired the fewest percentage of highly concordant genes with 0.11%.

18To gain further insights into the function of AUGs and DDGs identified in this work, pathway enrichment analysis was performed on AUGs and DDGs by using DAVID. According to KEGG enrichment analysis results, a large proportion of AUGs were affiliated with metabolic pathways related to Oxidative phosphorylation and Glycosylphosphatidylinositol (GPI)-anchor biosynthesis (Fig. 5a), which manifest well-known gained function of metabolism-related proteins in tumors. In contrast, DDGs were significantly related with ubiquitin mediated proteolysis, wnt signaling pathway (Fig. 5a), whose dysfunction tend to lead to tumorigenesis, metastasis, resistance, etc. Interestingly, both AUGs and DDGs are remarkably enriched in the ubiquitin-proteasome system, which has been increasingly reported to be highly related to cancer recently.

19In addition, we found 9 highly concordant genes with a strong relation to patient overall survival including 5 AUGs and 1 DDG (Table 3). For examples, the copy number amplification and expression level upregulation of FYTTD1 and CTTN would remarkably cause a poor prognosis in tumor patients of ESCA and head and neck squamous cell carcinoma (HNSC) respectively (Fig. 5b), while the copy number deletion and expression level downregulation of MTAP was significantly associated with a worse prognosis (Fig. 5b). Wherein, it has been widely reported in literatures that five of these highly concordant genes are highly related with the development and progression of numerous cancers (See Additional file 2: Table S8 references) but FYTTD1. Thus, we used STRING platform to analyze the protein-protein interaction (PPI) relationships of FYTTD1 (Fig. 5c). By integrated the functional enrichment results of KEGG pathways and Gene Ontology (GO) biological process of the PPI network (Fig. 5d), we found FYTTD1 plays an important role in RNA splicing in close relation to tumorigenesis, which may be a potential prognostic marker in ESCA.

20To validate our previous findings, the copy number and gene expression profiles were integrated across 1020 cell lines from CCLP, an independent dataset. Although CCLP utilized a different algorithm to calculate copy number variation, it still showed a positive correlation between copy number and expression level (Fig. 6a) by using a linear regression fitting on median Z score versus copy number (r = 0.89, p = 9.089e-6; Additional file 1: Fig. S5). In addition, gene expression levels with gene-wise copy number amplification exhibited a significant higher Z score than those with copy number deletion (Fig. 6b), in accordance with our results. Similarly, the amount of copy number amplification and expression level upregulation, copy number deletion and expression level downregulation occupied 28% of the total variant copy number count, while the proportion of copy number amplification and expression level downregulation, copy number deletion and expression level upregulation count only took up 1%. Thus, we counted the sum of copy number amplification and expression level upregulation, copy number deletion and expression level downregulation for the copy number aberrant genes among 1020 cell lines. Concordantly, most genes showed an overwhelming level of either copy number amplification and expression level upregulation or copy number deletion and expression level downregulation (93%, ratio > 0.9). Very fewer genes were with both high level of copy number amplification and expression level upregulation, copy number deletion and expression level downregulation (Fig. 6c). In addition, we overlapped the AUGs and DDGs identified by the ratio of copy number amplification and expression level upregulation versus copy number deletion and expression level downregulation between CCLP and TCGA respectively. Quite a number of highly concordant genes of AUGs (e.g. CTTN, GRB7, NSMCE2, KIAA0196, etc.) and DDGs (e.g. KLHL9, PTEN, PLAA, SMAD4, etc.) were shared within these two independent datasets. Collectively, these data suggested the close correlation between CNV and differential gene expression for most genes.

21In this study, we provided strong evidences to support the high correlation between CNV and differential gene expression. This finding reveals the qualitative relationship between genetic variation and its downstream effect, especially for oncogenes and tumor suppressor genes, which is of a critical importance for prevention, diagnosis and treatment of cancer. First, by integrated analysis of CNV and differential gene expression of CCLE, NCI-60 and TCGA, it revealed a positive association between copy number and expression level with high Pearson’s r of fitting, positive ρ and significant p values. Besides, not only in cell lines but in patients copy number amplification strikingly harbored a higher expressed level compared to copy number deletion. Secondly, we investigated every gene over the relationship of CNV and differential gene expression across 9139 tumor samples and 1025 cell lines respectively. Our results showed the majority of genes the copy number displayed a positive linear influence on gene expression, indicating that genetic variation generated a direct effect on gene transcriptional level. In addition, we validated 10 genes with a significant correlation between CNV and differential gene expression through literature (Table 1). A strong correlation was confirmed combining ρ or Pearson’s r of fitting for 9 genes except the weak evidence for NFASC, possibly due to the difference of analytical method.

22A recent study by GTEx consortium associated genetic variants with gene expression levels across 44 human healthy tissues and gene expression levels are found to be affected by local genetic variation for most genes based on eQTL analysis. Meanwhile, it was reported that copy number and expression levels had a strong positive correlation for 99% abundantly expressed human genes by integrating predicted copy number and corrected expression level from 77,840 expression profiles. Moreover, it has been widely reported that copy number is remarkably correlated with expression of protein in literature such as FGFR1, HER2, MET, FADD, EGFR. Message RNA, as intermediates between genes and functional proteins, plays a vital role in proteins production. Thus, we speculated gene expression might be correlated with protein expression as well for the high correlation and concordance between CNV and differential gene expression across cell lines and TCGA datasets and between CNV and differential protein expression of these five genes in literature (Additional file 2: Table S9, Table S10). Notably, FGFR1, known as fibroblast growth factor receptor 1, has been discovered that its copy number amplification is strikingly correlated with FGFR1 gene upregulation and FGFR1 protein upregulation in tumor samples. It indicated that the dysregulation of protein might attribute to original copy number aberrance through the concordant differential gene expression.

23However, a fraction of genes’ expression level did have nothing to do with copy number keeping in a stable expression level over various copy number. We think these genes might be involved in the maintenance of the basal cellular function such as metabolism and signal transduction by the results of significant KEGG pathway enrichment including retinol metabolism, olfactory transduction, calcium signaling pathway, neuroactive ligand-receptor interaction, etc. (Additional file 2: Table S2 and Table S3). Otherwise, it has been well documented that 24% of the 575 housekeeping (HK) genes accounted for the metabolic proteins and 19% for RNA-interacting proteins. Thus, we focus on the whole small nucleolar RNAs and found most genes were indeed expressed very stably versus CNV (Additional file 2: Table S11 and S12). Third, our results revealed the little existence of highly inconsonant genes of copy number amplification and expression level downregulation, copy number deletion and expression level upregulation (Additional file 1: Figure S4), which indicated that the copy number amplification barely causes gene expression downregulation and the copy number deletion hardly promotes gene expression upregulation. Otherwise, among the highly concordant genes with copy number amplification and expression level upregulation, copy number deletion and expression level downregulation, the frequency of copy number amplification and expression level upregulation evidently exceeded copy number deletion and expression level downregulation (Fig. 3b) possibly as a result of selection on deletions for it is unknown of the selective pressures on amplification. We attempt to revalidate the ten highly concordant genes in literature (9 AUGs, 1 DDG; Table 2), whose results was highly consistent with the variation trend in literature.

24Note that although the sample sizes of CCLE, NCI-60 and 31 cancers in TCGA were discrepant (Additional file 2: Table S1), they still showed a similar tendency of the association between CNV and differential gene expression (Fig. 1a; Additional file 1: Figure S1 and Figure. S2). Moreover, we observed a high level of agreement between cell lines and TCGA datasets which showed a consistent distribution of genes in Fig. 2a and Additional file 1: Fig. 3a including the ρ for 16,639 shared genes from (Fig. 2c) and a comparable Pearson’s r of fitting (Fig. 2b; Additional file 1: Figure S3C and Figure S3D). Our results suggested that this phenomenon was well conserved within cell lines and tissues.

25In total, we identified 925 highly concordant genes including 560 AUGs and 365 DDGs. For examples, numerous studies reported that DERL1 overexpression was significantly related to cancer cell proliferation, invasion and poor prognosis, which might be driven by copy number amplification for DERL1 obtained the majority of copy number amplification and expression level upregulation in many cancers (Fig. 4a). Obviously, CNV-driven differentially expressed genes (DEGs) might broaden our insights into the mechanism of tumorigenesis, migration, resistance, poor prognosis, etc. for the increasing studies on CNV-driven DEGs. In our study, a large proportion of AUGs were affiliated with metabolic pathways especially in terms of Oxidative phosphorylation and GPI-anchor biosynthesis (Fig. 5a), which suggested the gained function of metabolism-related proteins in tumors to provide more energy for cancer cells. In contrast, DDGs were significantly related with ubiquitin mediated proteolysis and wnt signaling pathway (Fig. 5b), whose dysfunction tend to lead to tumorigenesis, metastasis, resistance, etc. With respect to wnt signaling, lost function of DDGs such as inhibitory SMAD4 and APC would definitely enhance the function of wnt signaling leading to tumorigenesis, while attenuated function of ubiquitin mediated proteolysis facilitate proliferation. Wherein, we found 10 highly concordant genes with a strong relation to patient overall survival including 5 AUGs and 1 DDG (Table 3), while FYTTD1 has been hardly reported to be associated with cancer. By further integrated analysis of CNV and differential gene expression of FYTTD1 in ESCA patients, we observed that 24.73% patients showed a high level of copy number amplification with a median Z score of 4.15 which means FYTTD1 was strikingly overexpressed (Additional file 1: Figure S6). Wherein, AU patients occupied an overwhelming part among these 48 copy number amplified patients (84.44%) indicating global effect of CNV on FYTTD1 gene expression, which may be a potential driver gene or prognostic marker in ESCA. Therefore, highly concordant genes of AUGs and DDGs may provide new insights into the development and progression of cancer.

26Additionally, we utilized another independent dataset (CCLP) to revalidate the relationship between CNV and differential gene expression. Although CCLP applied a different algorithm to calculate copy number variation, it also showed a positive correlation between copy number and expression level (Fig. 6a). Our results demonstrated that gene expression levels of copy number amplification substantially surpassed gene expression levels of copy number deletion (Fig. 6b). Besides, copy number amplification and expression level downregulation, copy number deletion and expression level upregulation versus copy number aberrant counts took up the smallest part (1%). Concordantly, most genes showed an overwhelming level of either copy number amplification and expression level upregulation or copy number deletion and expression level downregulation (93%, ratio > 0.9), and it was hardly existed of genes with both high level of copy number amplification and expression level upregulation, copy number deletion and expression level downregulation (Fig. 6c).

27In conclusion, this study demonstrated the close correlation between CNV and differential gene expression. Moreover, this trend is consistent across cell lines and patient samples. For the majority of genes, copy number shows a positive linear influence on gene expression (copy number amplification and expression level upregulation, copy number deletion and expression level downregulation), while copy number amplification barely causes gene expression downregulation and the copy number deletion hardly promotes gene expression upregulation. Furthermore, both AUGs and DDGs are remarkably enriched in the ubiquitin-proteasome system. In addition, we identified amplification and overexpression of FYTTD1 is highly related with poor prognosis in ESCA, which may be a potential prognostic marker in ESCA. Whereas, more in-depth studies are needed to further reveal molecular mechanisms between CNV and differential gene expression. Overall, it is of a critical importance to elucidate the relationship between copy number variation and gene expression for prevention, diagnosis and treatment of cancer.