Identification of genomic functional hotspots with copy number alteration in liver cancer

Copy number alterations (CNAs) can be observed in most of cancer patients. Several oncogenes and tumor suppressor genes with CNAs have been identified in different kinds of tumor. However, the systematic survey of CNA-affected functions is still lack. By employing systems biology approaches, instead of examining individual genes, we directly identified the functional hotspots on human genome. A total of 838 hotspots on human genome with 540 enriched Gene Ontology functions were identified. Seventy-six aCGH array data of hepatocellular carcinoma (HCC) tumors were employed in this study. A total of 150 regions which putatively affected by CNAs and the encoded functions were identified. Our results indicate that two immune related hotspots had copy number alterations in most of patients. In addition, our data implied that these immune-related regions might be involved in HCC oncogenesis. Also, we identified 39 hotspots of which copy number status were associated with patient survival. Our data implied that copy number alterations of the regions may contribute in the dysregulation of the encoded functions. These results further demonstrated that our method enables researchers to survey biological functions of CNAs and to construct regulation hypothesis at pathway and functional levels.


Introduction
Chromosomal instability is one of the characteristics in cancer [1] and results in the numerical and structural alterations of DNA copy number variations (CNAs). Recently, some literatures have reported the association of CNAs and patient survival in different tumors [2][3][4]. Several important oncogenes or tumor suppressors were also showed with high frequency of gain or loss status in different cancers. For example, the copy number amplification of gene Her2, which is the addicted oncogene in the HER2+ subtype of breast cancer, was highly correlated with the gene overexpression [5]. However, in addition to focal amplification, most tumors display multiple and broad ranges of copy number change, where large number of genes are involved in and potentially to be induced or suppressed due to copy number amplifications or deletions. Some in vitro studies were performed to survey the affected functions of CNAs [6][7][8]. For example, Nicole et al. utilized the shRNA library to identify the GO and STOP genes which positively and negatively regulate proliferation to evaluate the effect of gene deletions [7], respectively. They also proposed a model called 'Cancer Gene Island', which encompasses high density of genes with the same function within a genomic region [7]. However, the in vitro studies were labor intensive if not cost prohibitive. Moreover, it is hard to perform a systematic analysis based on these approaches, thus, leaving the gene island model and their functions unexplored.
In conventional gene expression data analysis, several bioinformatics methods based on the concept of 'gene set enrichment analysis' (GSEA) have been successfully utilized to explore the underlying molecular pathways and Gene Ontology functions [9][10][11][12]. The GSEA method assesses the number of overlap genes between two gene sets: the differentially expressed genes of a certain functional annotation and genes from the entire genome with the same annotation, to estimate the probability of the overlapping through the statistical test. The procedure provides a high throughput and systematic analysis to explore the putative activated pathways or functions. Hepatocellular carcinoma (HCC) is one of the malignant cancers and the third leading cause of cancer death worldwide [13]. Major etiologies associated with HCC are hepatitis B virus (HBV) and hepatitis C virus (HCV) infection [14]. Previous studies have been reported in which comparative genomic hybridization by microarray (aCGH) was utilized to examine CNAs in HCC. Several regions with frequent copy number gain and loss were identified. The CNA-associated oncogenes and tumor suppressors, such as MYC, JAG1, TP53, and RB1, were also found [15][16][17][18]. The association between survival and CNAs has been investigated, and ten associated genes were reported [19]. However, the biological functions altered by CNAs remain unknown and thus need to be dissected.
According to the concept of Cancer Gene Island, here, we propose an algorithm to identify the spatial functional hotspots (SFHs) in human genome based on the enrichment analysis. The human genome is divided firstly into segments along the genomic sequence coordinate. Then, the tests of enrichment between the segments and whole genome functional categories are performed. Finally, a method which identifies the optimal regions of enriched functions between the segments was applied to examine putative SFHs. To demonstrate the ability of our method, we applied the method to an aCGH data set of HCC. The result showed several immune-related SFHs which showed gain and loss in HCC samples. Also, survival-associated SFHs were identified. The result also indicated that our system could serve as a useful method to understand the CNAs-affected functions.

Methods
To identify the SFHs in human genome, we proposed a novel enrichment analysis that compares the genes contained within a genomic segment with all genes belonging to the same function categories associated to the genes within the segment under consideration based on the concept of gene set enrichment. As shown in Figure 1A, two matrixes, B and P, were constructed first. The indicator matrix B contains information whether or not a gene belongs to a genomic region (spatial segment) determined by a sliding window along the genomic position of all chromosomes or B = (b k,i ) KxM , where M is the number of genes and K is the number of genome segments pre-determined and where b k,i = 1 when ith gene is in the kth segment, otherwise 0. The matrix P = {p i,l } MxL is also an indicator matrix of functional gene sets, where L is the number of functional gene sets and p i,l = 1 when ith gene is in the lth GO (Gene Ontology) function, otherwise 0. The enrichment is defined as scoring function C of the two matrixes B and P.
Here, we use Fisher's exact test as the score function C Based on the p values, we can determine if the function l was enriched at the genome segment k. Then, we merge and extend the enriched segments to a merged window to include all genes involved in the function l if the segments were located nearby and have position overlapping. As shown in Figure 1C, assuming qth to (q + R)th segments have enrichment for the function l, the genes involved in the merged windows can be expressed as vector d: where s i ¼ X qþR ð Þ t¼q b t;i . Assuming there are G genes (from eth to (e + G)th) located in the qth to (q + R)th enriched window, we defined the subsets of the G genes which exclude out genes gradually from left or right side according to the genome coordinate. Two parameters, pL and pR, which perform enrichment analysis (Fisher's exact test) between the subsets with the gene set of function l were introduced ( Figure 1C). pL and pR are defined as: where g = 1,…,G and C(.) is the enriched score function. Then, the optimal enriched region o of function l can be defined as: If the p value of the region o passed the selection threshold, o was defined as the SFHs of function l.

Gene sets of genome segments and biological functions
To define the genome segments, the detection window size was set as one million base pairs (Mbp) after the testing of three different conditions (Additional file 1: Figure S1). The sliding distance was set at 0.25 Mbp. The genomic position of each gene was obtained from Ensembl (version Homo Sapiens 65) [20], or equivalent to NCBI human genome GRCh37. Therefore, a total of 12,098 segments were defined. To construct the functional gene sets, we downloaded all records of Gene Ontology from the BioMart website of Ensembl 65 (http://useast.ensembl.org/info/ data/biomart.html) [20]. A total of 7,654 GO terms were downloaded. After excluding the gene sets containing fewer than 15 genes, 1,091, 404, and 275 gene sets associated to biological process (BP), molecular function (MF), and cellular component (CC) terms were utilized in this study, respectively.

aCGH arrays of hepatocellular carcinoma
To identify the functional effect of CNAs in HCC, the aCGH array data set, GSE14322, was downloaded from GEO/NCBI website. The data set contains 76 HCC samples. The determination of CNAs was through the NEXUS software (BioDiscovery, San Diego, CA, USA). The CBS segmentation algorithm was performed to identify the segments of CNAs [21] using the thresholds of log2 values of fold change larger or smaller than ±0.2.  (Figure 2A). The averaged SFH density is 0.43 SFHs per million base pairs (Mbp). Chromosome 6 has the highest SFH density (0.48 SFHs/Mbp) ( Figure 2B). For the 838 SFHs, the average length of SFHs was 0.56 Mbp ( Figure 2C) and the averaged 11.5 genes are in a SFH ( Figure 2D). The SFH of 'sugar binding' enrichment, which is located in the 7.88 to 10.6 Mbp region at chromosome 12, has the longest region length. The SFH of 'immune response' enrichment (31.2 to 33 Mbp at chromosome 6), which contains 93 genes, has the largest number of genes. The region located in 29.7 to 31.5 Mbp at chromosome 6 contained the most number of enriched gene sets (16) ( Figure 2E and Additional file 1: Table S1). The region includes lots of SFHs which have enrichment of immune-related gene sets, such like MHC class I protein complex, type I interferon-mediated signaling pathway, and immune response. Our finding indicated that the two regions are important for cell immunity.

The affected function of copy number variation in liver cancer
To evaluate the effect function of CNAs in liver cancer, the dataset GSE14322, which contains 76 aCGH arrays of HCC samples, was downloaded and analyzed. The percentage of CNA status of each SFH was calculated. There are 61 and 89 SFHs that contained copy number gain and loss in more than 30% patients (25). The result was showed in Table S2 in Additional file 1, and the top ten SFHs were listed in Table 1. One immune-related gene set had the gain status in most of the samples (innate immune response), and one had loss status (response to virus), since the major etiologies of HCC are the infection of HBV and HCV. We hypothesize that those immune-related SFHs that harbor CNAs may play a role in the HCC carcinogenesis.
We also analyzed the association between disease-free survival and the CNAs of the SFHs through log rank test. Using p < 0.01 as the threshold, a total of 20 and 19 SFHs of which gain and loss status were identified with survival association, respectively (see Table 2). The copy number gain status in the SFH which located at 41.1 to 41.9 Mbp at chromosome 19 had the smallest p value of the survival testing. The SFH had the enrichment of 'oxygen binding'. As shown in Figure 4A, the patients with copy number gain in the SFH had reduced survival comparing with neutral status. Interestingly, all the SFHs with survival-associated gain status were all located at were enriched in the region. The finding indicated that the immune functional island located at the region is sensitive to patient survival. The SFH located at 11.8 to 12.2 Mbp at chromosome 8, which has enrichment of 'defense response to bacterium' , has the smallest p value of copy number loss status ( Figure 4B). For SFHs with survivalassociated loss status, 12 of them were located at 55 to 76.9 Mbp region at chromosome 4, and 7 of them were located at 11.1 to 38.3 Mbp at chromosome 8.

Discussion
We introduced a system biology method, motivated by Cancer Gene Island, to identify the spatial functional hotspots in human genome. A statistical assay was presented to estimate the enrichment within genome regions to functional gene sets. By applying the terms of Gene Ontology into our method, the result provided the details of the function encoded in human genome. We set the two parameters of the algorithm, the length of window size and shift distance, as 1 and 0.25 Mbp, respectively. Although the setting of the parameters will affect the p value of enrichment testing for each segment, our algorithm performed an optimal procedure which merge the continual enriched segments and find the region with maximum p values by removing the gene one by one from both sides. Different settings of window size will not affect the results of final optimal regions. However, the detection of continual functional enriched segments could be missed under the condition of small window size because the windows contained no and less genes. To find out the workable parameters, we tested three conditions of window sizes, 0.5, 1, and 1.5 Mbp and found out that the condition of 0.5 Mbp contains large numbers of segments of which the gene number is less than three. The parameters of 1 and 1.5 Mbp contain fewer segments with low numbers of genes.
Through the testing, we set the window size as 1 Mbp to analyze the human genome. We applied the method in HCC data set to estimate the effect of hotspots in the genome. Using the data set GSE14322 as an example, a total of 150 SFHs have been identified with copy number alterations in most of the HCC patients, and the novelty of our analysis is to identify the functional hotspots in human genome. The region we identified is located with high density of genes that share the same biological function, and as we demonstrated in the HCC dataset, these functions may also be sensitive to CNAs. Two immune-related functional regions were identified with gain or loss in most of patients in the HCC dataset. The major carcinogenesis of HCC is the chronic and acute inflammation under HBV or HCV infection; thus, we hypothesize that these two regions we identified may also play a role in HCC oncogenesis.
We also identified 39 SFHs of which the copy number status was associated with patient survival. The result indicates that the copy number alterations in these regions may affect the function of tumor progression and then reflect on patient survival. For example, the patients who have copy number loss in the SFH which was enriched in inflammatory response have shorter survival. The chronic and acute inflammations induced by HBV and HCV infection have been proved to play an important role in HCC tumorgenesis [22,23]. Our data implied that copy number alterations may contribute in the inflammatory response. Also, other enriched functions in survival-related SFHs have been reported, such as regulation of immune response, cell growth, apoptosis, and caspase activator activity. The SFHs and enriched function we identified provided the clues of the association between CNAs and the regulations of the enriched functions. We expected that the SFHs we identified will provide further insight of affected functions of CNAs to uncover the mechanism of cancer.

Conclusions
In this paper, we systematically surveyed human genome and identified 838 functional hotspots based on Gene Ontology classification. To substantiate our findings, 76 HCC tumors and their DNA copy number gain/loss statuses were examined closely. Among the 838 hotspots, a total of 150 regions affected by CNAs, and the encoded enriched functions were identified. Our results indicate that two immune-related hotspots had copy number alterations in most of the patients and might be involved in HCC oncogenesis. In addition, 39 survival-related hotspots were identified. Taken together, our results demonstrated that the method presented in the paper is a powerful tool to survey biological functions of CNAs and to construct regulation hypothesis at pathway and functional levels.