A putative role for genome-wide epigenetic regulatory mechanisms in Huntington’s disease: A computational assessment

Concept profile analysis

In Concept Profile Analysis (CPA), the vector space model is used to associate two concepts mined from the literature with each other. Advantages of this model include efficient and transparent comparisons, and the possibility of attaching a weight to the association²⁰. The CPA algorithms have previously been used for a range of different gene expression data analysis purposes such as functional annotation²², comparison of studies²³, prediction of novel interactions²⁴, generation of gene sets^19,25, and association with chemical structures²⁶). The methodology has been described previously²⁷. In short: In our database, every concept is associated with PubMed records using the indexing engine Peregrine (https://trac.nbic.nl/data-mining/) which is equipped with an in-house thesaurus of biomedical and chemical concepts that have been prepared for text mining^28,29. For all concepts except genes and Gene Ontology (GO) terms the PubMed records are comprised of the texts in which the concept is mentioned. For genes, only a subset of PubMed records are used in order to limit the impact of ambiguous terms and distant homologs. GO terms are sometimes given as words or phrases that are infrequently found in the normal texts. To still provide broad coverage of GO terms, the PubMed records that were used as evidence for annotating genes with the GO term are added. For every concept in the thesaurus that is associated to at least five PubMed records, a vector containing all concepts related to the main concept (direct co-occurrence), weighted by the symmetric uncertainty coefficient is created. We call this a "Concept Profile". Concept profiles are matched to identify similarities via their shared concepts (indirect relations). Any distance measure can be used for this matching such as the mutual information, inner product, cosine angle, Euclidean distance or Pearson’s correlation. The CPA Web Services that we used for our analysis use an inner product measure³⁰. These web services can be found in the BioCatalogue web service registry https://www.biocatalogue.org/services/3559.

Data analysis and interpretation

For data analysis and interpretation we implemented a series of workflows using the Taverna workbench (Taverna workbench version 2.4.0)^31,32. Taverna is an open source software for the development and execution of workflows. Our workflows are deposited online on the Zenodo repository (https://doi.org/10.5281/zenodo.164201³³ and https://doi.org/10.5281/zenodo.164198³⁴) and the myExperiment platform (http://www.myexperiment.org/packs/553).

Workflows for data analysis

The first Taverna workflow load_data_identify_DE_genes_Array_A, was implemented to examine differential gene expression between control and HD samples. Required workflow inputs were two data files with gene expression values and the phenotype information that describes the samples from the microarray experiment. Differential expression was computed using moderated t statistics with the package limma³⁵ (version 3.14.1), which is provided by the bioconductor project³⁶, R version https://www.bioconductor.org/. We analyzed each brain region separately, because previous analysis³⁷ revealed regional patterns in gene expression. The workflow maps the expression data from probes to entrez gene ids using the Affymetrix Human Genome U133 set annotation data, (packages hgu133a and hgu133b, version 2.8.0; https://bioconductor.org/packages/hgu133a.db; https://bioconductor.org/packages/hgu133b.db/). When multiple probe names map to the same gene id, the ones exhibiting the most significant changes were used for further analysis. Final outcome of this workflow is a report. Each row is composed of a gene id, a fold change and its corresponding P-value indicating the significance of every change in gene expression, between HD and controls for each brain region. Adjusted P-values, generated by Benjamini and Hochberg’s method for multiple testing correction, are also included³⁸. This workflow can be adjusted to compute differential gene expression between other variables such as male/female or the grade of disease pathology by editing the nested workflow “compute_DE_limma” within the R workflow component. An additional workflow create_exprs_obj_download_files is included in the myExperiment pack that was used to download data from the ArrayExpress repository. The workflow saves the gene expression data and the corresponding phenotype file in the directory indicated in the workflow input.

We note that this particular microarray experiment was composed of two microarrays, Human Genome U133A and U133B. For convenience we added the workflow: Get differentially expressed genes for Array B one brain region, but in principle the first workflow could be reused by adjusting the “libraries” component.

The second workflow map_genes_on_chromosome uses the output from the first workflow in order to map genes to their corresponding genomic location. The workflow uses the Biomart³⁹ service within R, to obtain information regarding the position of each gene at the chromosome, HGNC gene symbols⁴⁰, transcription start and end site and the transcription strand. The database that was used was the Ensembl genes 68 from Sanger institute and the Homo Sapiens dataset GRCh37.p8. The mouse assembly that was used to map genes to their chromosomal location was Dec 2011, GRCm38mm10.

The last workflow get_promoter_region_calculate_overlaps, first computes a promoter region for each gene and then operates on genomic intervals to compute gene promoters that overlap with a genomic region. The promoter region is computed for each gene, according to prespecified values, indicating the number of base pairs (bp) upstream and downstream of the transcriptional start site (TSS); for the CpG island analysis 5000bp upstream and 2000bp downstream and for the chromatin states analysis 50bp upstream and 50bp downstream was used. The decision for the promoter size in each case was taken after discussing with the domain experts and from knowledge acquired from previous experiments using the ENCODE data from Ernst et al.⁴¹. Using those data we performed multiple runs with different input values for the “upstream” and “downstream” variables, and the overlap parameter (Supplementary File 1).

When genes had multiple transcription start sites, we computed a promoter region for every TSS. Next part of this workflow is to compute overlapping regions between the input datasets. It includes a two sample Kolmogorov – Smirnov test to compare the empirical cumulative distribution functions (ecdf) of the P-values between the gene promoters that overlap with a specific genomic region and the ones that do not. The null hypothesis tested here was that there is no difference between the two groups. The workflow returns two lists of genes, one for the genes that overlap with a particular genomic region and another that does not. Furthermore, the results of the statistical test are reported: the ks test statistic (maximum distance D between the ecdf of the two samples) and the P-value of the test. If P-value < 0.05 we reject the null hypothesis.

Workflows for data interpretation

The workflows that we implemented for gene interpretation and gene prioritization are based on the workflow pack at Zenodo, https://doi.org/10.5281/zenodo.164198 (see also, at myExperiment: http://www.myexperiment.org/packs/368), and are implemented using CPA web services³⁰.

The CPA workflow Annotate gene list with top ranking concepts annotates a gene list with top ranking concepts by matching concept profiles of genes with for example in our case the concept set of Biological Processes. The web services part of this workflow query the Anni database²⁷ that stores the concept profiles for each concept of interest. The first web service mapDatabaseIDListToConceptIDs maps a list of concepts, in our case Entrez gene identifiers, to their corresponding concept profile ids. Necessary inputs are a concept list (gene list) in a comma separated file, and the database identifier of the gene list necessary for the mapping (EG for Entrez Gene, see here: https://www.biocatalogue.org/soap_operations/41197 for more details on database identifiers). The next web service getSimilarConceptProfilesPredefined matches our gene list with the predefined concept set “Biological Process” (ID = “5”), and gives the top scoring biological processes that describe our gene list. For a complete list of predefined concept sets, the workflow List Concept Sets (provided in the current pack) can be run to choose the ID of the predefined concept set of interest. The web service getConceptName, gives the complete (human readable) names of the top matching biological processes. Lastly, the workflow Explain score between two concepts can be included to the analysis to provide evidence for the association between each gene and the annotations of the biological processes. The evidence reported is a list of concepts that link one concept with another and the contribution to the overall strength of the association. In addition, the corresponding concept ids are reported.

The workflow Prioritize gene list can prioritize a set of genes with respect to their association with particular concepts, in our case the HTT concept and epigenetics (concept profile: “epigene”). In order to obtain the concept profile identifiers the workflow getConceptSuggestionsFromTerm needs to run first.

Data obtained from public sources

Human brain data. The HD human brain data that was used in this analysis was originally produced and analyzed by A. Hodges and co-workers³⁷. This experiment contains 44 HD positive cases and 36 age and sex matched controls. The processed data are available from the public repository NCBI Gene Expression Omnibus, entry GSE3790. Three brain regions were included and analyzed; the caudate nucleus, frontal cortex and cerebellum, with an Affymetrix Microarray GeneChip (Human Genome U133A and U133B). Furthermore, the HD positive cases were further classified based on whether symptoms were present or absent and according to Vonsattel grade of disease pathology (scale = 0 – 4). In our analysis we used the processed data and performed our own differential gene expression analysis (Dataset 1⁴²) with the workflow that performs differential gene expression analysis that was described previously in Methods.

CpG island data. CpG island information in the human genome was obtained from UCSC genome browser⁴³, hg19 assembly FEB 2009 1 (Dataset 2⁴⁴). Here, CpG islands are marked as the DNA regions where the following conditions hold:

The ratio observed/expected (Obs/Exp) CpG was calculated as follows:

where N is the total amount of nucleotides in the sequence that is being analyzed. For CpG island information of the mouse genome the assembly dec 2011 (GRCm38/mm10) was used.

Chromatin states data. The chromatin marks were obtained from the encode project⁴⁵. The chromatin states were part of an integrative analysis of 111 reference human epigenomes profiled for histone modification patterns based on DNA accessibility, DNA methylation and RNA expression. We used the two cell types that were more suitable for our analysis; the anterior caudate and dorsolateral prefrontal cortex (Dataset 3⁴⁶ and Dataset 4⁴⁷, respectively). The chromatin states we used for the current analysis were: active transcription start site proximal promoter (TssA), bivalent regulatory region (TssBiv), heterochromatin (Het) and repressed Polycomb (ReprPC).

Mouse data. The mouse brain data was taken from a published study¹⁰. There, total RNA was extracted from cortex, striatum, and cerebellum from WT and R62 transgenic mice. This study examines the effect of the HDACi 4b inhibitor on the disease phenotype. However only the data from animals treated with vehicle was used (Dataset 1⁴²). This study used Illumina Mouse Mouseref-8 Expression Beadchips v1. The raw data were analyzed using the Bioconductor packages and contrast analysis of differential expression was performed by using the LIMMA package. The differential expression values are available in the supplemental material of that publication.

Chromatin state analysis and semantic interpretation

To test if and how differential gene expression in HD is associated with particular chromatin states we used publicly available datasets from the GEO and ENCODE public repositories. We selected HD gene expression data from three regions of the brain made available by Hodges and coworkers³⁷, and data from the ENCODE consortium carrying information about four chromatin states. These were: active TSS proximal promoter (TssA), bivalent regulatory region (TssBiv), heterochromatin (Het) and repressed Polycomb (ReprPC)⁴⁵. Briefly, the first state pertains to active genes, the second to repressed genes that are ready to be activated, and the latter two represent repressed genes. The chromatin states were part of an integrative analysis of 111 reference human epigenomes profiled for histone modification patterns, DNA accessibility, DNA methylation and RNA expression. The different chromatin states were defined by a computational model that is based on a multivariate Hidden Markov Model⁴⁸.

First we confirmed the results reported by the previous study³⁷: large and numerous changes in transcriptional activity in the caudate nucleus brain region and notably smaller changes in the frontal cortex and cerebellum. More specifically, the number of differentially expressed genes was 5219, 127 and 96 for each brain region respectively with a FDR of 0.05. These results confirm previous observations illustrating that specific HD-affected brain regions exhibit defined changes in gene expression, in line with observable physiological effects³⁷.

Secondly, we paired the data describing gene expression changes in the caudate nucleus and frontal cortex from the Hodges study to chromatin state data from the anterior caudate and dorsolateral prefrontal cortex from ENCODE. We considered these brain regions as the most comparable between the two studies. The cerebellum was excluded from this part of our analysis due to the absence of chromatin state data at ENCODE for the cerebellum. We determined the reference chromatin state of the genes from the gene expression study by determining the overlap of their promoter regions with the start and end positions of the chromatin state from the aforementioned brain regions used by ENCODE. The effect of a chromatin state on differential gene expression in HD was assessed by comparing the distribution of expression levels of genes overlapping with a particular chromatin state with the distribution of expression of non-overlapping genes. If chromatin state has substantial functional impact on the genes that are differentially expressed in HD, then these are expected to be significantly different. More details on the promoter region calculation and the overlap between the genomic regions can be found in the Methods section.

Specifically, we compared the distribution of the P-values for differential expression between these two groups of genes and assessed the difference by a KolmogorovSmirnov test, in the two brain regions and for each chromatin state (Figure 1). In the caudate region, we found an enrichment of genes overlapping with the active TSS and a depletion of genes overlapping with ReprPC (Figure 1). Heterochromatic or bivalent chromatin were not significantly associated with a genomic location of differentially expressed genes. A similar, but less pronounced pattern was observed for the frontal cortex for both the active TSS and ReprPC. In addition, in the frontal cortex the number of genes overlapping with the bivalent state was reduced. The magnitude of the association of a brain region with each chromatin state is in line with the HD neurodegeneration pattern, where the caudate nucleus exhibits the largest gene expression changes while the frontal cortex exhibits an intermediate to low pathology. In summary, we found enrichment of differentially expressed genes in the active chromatin state of reference tissue, but no strong evidence for significant enrichment in chromatin regions that are associated with repressed gene expression activity.

Figure 1. Chromatin states.

For each brain region, the P-value distribution for differential expression in HD patients compared to controls was compared between genes overlapping with a specific chromatin state and all other genes, in the anterior caudate and dorsolateral prefrontal cortex cells. The plot displays the maximum absolute distance D between the cumulative distribution function of the P-values for each group of genes (overlapping and non overlapping), as reflected by the KS test statistic. Positive values correspond to an enrichment of differentially expressed genes in a chromatin state, negative values with a depletion. Stars indicate a significant enrichment/depletion cf the KS test (P-value < 0.05). TssA= active TSS, TssBiv= bivalent Tss, Heterochr= Heterochromatin, ReprPC= repressed Polycomb state. For caudate: TssA D= 0.31, TssBiv D=0.04, Heterochr D=0.17, ReprPC D= 0.27, for frontal cortex TssA D=0.168, TssBiv D=0.097, Heterochr D=0.095, ReprPC D=0.147. The corresponding p-values were for the caudate : TssA pval < 2.2e-16, TssBiv pval= 0.59, Heterochr pval= 0.10, ReprPC pval=4.5167e-12, for frontal cortex TssA pval < 2.2e-16, TssBiv pval= 0.0002, Heterochr pval= 0.872, ReprPC pval= 1.4596e-05.

To further interpret the results from the chromatin state analysis, we used literature information (CPA^19,27) to assess the biological processes within the lists of genes per chromatin state that were associated with gene deregulation (Dataset 5⁴⁹).

Using CPA, we found that for the caudate nucleus, genes in the active TSS state are mainly associated with processes related to transcription, cell cycle, protein transport and modification, RNA splicing, chromatin modifications and RNA processing. Genes in the repressed polycomb chromatin state are mainly associated with (brain) developmental processes and responses related to zinc and cadmium stimulus.

For the frontal cortex, genes in the active TSS state are associated with protein transport and modification, cell cycle, RNA splicing and signal transduction pathways (MAPK, notch, smoothened). Genes in the bivalent TSS are associated with brain development, neurogenesis, synapse and responses related to zinc and cadmium stimulus. We found that the genes that are part of the polycomb repressed group to be associated with similar functions as the genes in the bivalent state: (brain) developmental processes, responses related to zinc, cadmium and copper stimuli.

We note that at this point of our analysis, we filtered out genes of which the promoters were labelled with more than one state, using the criterion of at least 50 bps overlap in a promoter region of 100 bps. Interpretation of this set of genes would be ambivalent.

Overlap of HD deregulated genes with CpG islands and semantic interpretation

Because CpG island methylation is a known epigenetic regulatory mechanism that is ubiquitous in the human genome and may be a target for genome-wide regulation^50,51, we also applied our approach to test if genes within CpG islands are overrepresented among differentially expressed genes in HD. We measured this by similar KS test statistics as in the previous section. We found that genes overlapping with CpG-islands in their promoter region, were significantly enriched in the group of HD-deregulated genes in all three brain regions (Figure 2A; P-value < 0.05), which supports taking this mechanism into account to formulate hypotheses when studying gene deregulation in HD.

Figure 2. CpG island effect in human and mouse data.

Maximum distance between the cumulative distribution function of the P-values between genes containing a CpG island in their promoter and genes that do not. The maximum distance (D) for each brain region was plotted for both datasets (human (A) and mouse (B) datasets). The gradient adjacent to each plot indicates the extent of neurodegeneration in each brain region. Black represents severe and white mild neurodegeneration. A: The analysis performed in the human data. Caudate nucleus is exhibiting the largest differences with a distance D between the two distributions of 0.2808, frontal cortex follows with D = 0.1482 and cerebellum with D = 0.147. B: The analysis performed in the mouse data. The results are analogous to the human data with striatum showing the largest differences with D= 0.118, cortex following with D=0.0739 and cerebellum with an insignificant difference (P-value = 0.596 >> 0.05) of D = 0.0383. * : depicts significance.

To further analyze and verify the association of gene deregulation in HD with the presence of CpG islands, we also analysed overrepresentation in data obtained from a HD transgenic mouse model¹⁰. Our results based on the mouse data corroborated those from the human data (Figure 2B). Comparable to the results for caudate nucleus in human data, the most highly affected mouse brain region, the striatum, showed the biggest overrepresentation of CpG islands among deregulated genes.

We next examined the biological processes that are associated with the differentially expressed genes that do or do not overlap with CpG islands by CPA. We inspected the top 50 annotations of each group. We identified many similarities between these two groups, but also annotations that were specific to each group (Dataset 6⁵²). For example, concepts related to chromatin alterations, such as chromatin remodelling and chromatin (dis)assembly, and histone post-translational modifications, such as (de)acetylation, and (de)methylation, were only found with a high rank in the list of genes containing CpG islands in their promoters.

Conversely, lymphocyte activation, angiogenesis, antigen presentation and neurogenesis were only high ranking associations for the non-CpG containing genes.

Some annotations were found in both groups but in a different ranking order. For example, the rankings of gene silencing, RNA splicing and phosphorylation were increased for genes within CpG islands while transcriptional activation and mitotic cell cycle rankings were higher for non-CpG containing genes (Dataset 6⁵²).

Semantic analysis to identify proteins associated with HTT and epigenetics

CPA can also be used to prioritize findings by a specific biological interest. In our case, we aimed to prioritise the list of genes overlapping with CpG islands that were also differentially expressed in HD caudate nucleus by their association with HTT and epigenetics. Here, we prioritised 100 proteins based on their association with HTT and epigenetics and absence of a direct relation with HTT by CPA (Dataset 7⁵³). Such relations are typically novel, or lost in tables or figures. If a novel relation is found (i.e. the relation is not found in our database of relations in MedLine abstracts), then CPA also provides intermediate concepts that link the two concepts.

In Table 1, we present the top 5 novel proteins that have the strongest association, and the intermediate concepts that link each prioritized protein to HTT and epigenetics. The intermediate concepts are grouped under the semantic categories “General”, “Biological Processes”, “Disease or Syndrome”, “Homo Sapiens proteins” and “Molecular Functions”.

Grouping the evidence provides a more complete insight into the processes involved in gene deregulation in HD mediated by CpG islands. For example, for one of our candidates, the amyloid beta (A4) precursor protein-binding, family A, member 1 (APBA1), we found nerve tissue (General), endocytosis (Biological Process), Alzheimer’s Disease (disease or syndrome), GRIN2B (H. sapiens Genes) and kinesin activity (Molecular Function) as intermediate links with HTT. This suggests that mechanisms involving APBA1 in HD share common components with the mechanisms involved in endocytosis and Alzheimer’s Disease and those involving GRIN2B and kinesin activity. Intermediate links with epigenetics were respectively: epigene (General), methylation (Biological Process), hypomyelination and congenital cataract (disease or syndrome), CDKN2A (H. sapiens Genes) and MGMT - O6-alkylguanine DNA alkyltransferase - (Molecular Function). Accordingly, these concepts provide suggestions about the epigenetic role of APBA1, which can be taken into account when further studying the role of APBA1 in HD.

Validation of concept profile gene prioritization

We next investigated whether the prioritized genes reflect valid biological knowledge. Figure 3 shows that CPA is able to prioritize true associations with huntingtin as measured by a gene expression experiment, but that combining experimental (differential expression) measurements and literature evidence enables to select even more specific HD signatures. We used our concept profile technology to match all genes in our database that have a concept profile (12,391 genes) to the “huntingtin” concept profile (black line). We then compared the distribution of those CPA scores to the CPA scores of genes that were found to be differentially expressed in the caudate nucleus (p value < 0.05). We included two gene lists in our analysis: the top 100 most differentially expressed genes (red line) and the top 1000 (green line). The shift in the distribution of CPA match scores between the differentially expressed genes (top 100 and top 1000) and the scores of all genes reflects the added value of CPA (CPA scores of top 100 and top 1000 can be found in Dataset 8⁵⁴). We found a significant shift in the scores when comparing all CPA scores from our database with the top 100 differentially expressed genes: p = 2.67e −08 and the top 1000 p < 2.2e− 16.

Figure 3. Combination of CPA with differential gene expression for effective gene prioritization.

Cumulative distribution of match scores of the concept profiles (CP) between differentially expressed genes in HD with the concept profile of htt: match scores of all genes with a concept profile (black), match scores of the top 1000 differentially expressed genes (green), and the match scores of the top 100 of differentially expressed genes (red).

Also the top 100 and top 1000 differ significantly (p = 0.03184), showing that it is useful to narrow down on the top ranks for follow-up research. In principle, more extreme p-values are associated with higher CPA scores. In addition, to show that our list of 100 prioritized gene-HTT CPA match scores would not have been found by chance, we assessed the percentile score of our list when compared to the frequency distribution of 100 match scores of randomly sampled gene-HTT concept pairs out of the 12,391 genes in our concept profile database (Dataset 9⁵⁵). All genes were in the top 95 percentile, except NTRK3 (55 percentile).