Volume 19 Supplement 5

## Selected articles from the Biological Ontologies and Knowledge bases workshop 2017

- Research
- Open Access

# Higher-order partial least squares for predicting gene expression levels from chromatin states

- Shiquan Sun†
^{1, 3}Email author, - Xifang Sun†
^{2}and - Yan Zheng†
^{1}

**19 (Suppl 5)**:113

https://doi.org/10.1186/s12859-018-2100-y

© The Author(s) 2018

**Published:**11 April 2018

## Abstract

### Background

Extensive studies have shown that gene expression levels are strongly affected by chromatin mark combinations via at least two mechanisms, i.e., activation or repression. But their combinatorial patterns are still unclear. To further understand the relationship between histone modifications and gene expression levels, here in this paper, we introduce a purely geometric higher-order representation, *tensor* (also called multidimensional array), which might borrow more unknown interactions in chromatin states to predicting gene expression levels.

### Results

The prediction models were learned from regions around upstream 10k base pairs and downstream 10k base pairs of the transcriptional start sites (TSSs) on three species (i.e., Human, Rhesus Macaque, and Chimpanzee) with five histone modifications (i.e., H3K4me1, H3K4me3, H3K27ac, H3K27me3, and Pol II). Experimental results demonstrate that the proposed method is more powerful to predicting gene expression levels than several other popular methods. Specifically, our method enable to get more powerful performance on both commonly used criteria, R and RMSE, as high as 1.7% and 11%, respectively.

### Conclusions

The overall aim of this work is to show that the higher-order representation is able to include more unknown interaction information between histone modifications across different species.

## Keywords

- Higher-order partial least squares
- Chromatin states
- Tensor decomposition
- Gene expression levels
- Histone modification

## Background

In epigenetics, histone modifications like methylation, acetylation, and phosphorylation play critical roles in transcriptional regulation process [1]. Specifically, during gene expression process, each unit of chromatin like beads wrapping around DNA subsequences (about 147 base pairs) is highly impact the process of gene expression by chemical modification of chromatin condensation and DNA accessibility when genetic information are converted into gene products [2]. These modifications are shown to regulate gene transcription with active or repressive manners [3]. For example, tri-methylation on K4 of histone H3 (i.e., H3K4me3) is primarily associated with transcriptional activation [4, 5], while tri-methylation on K27 of histone H3 (i.e., H3K27me3) are primarily associated with transcriptional repression [6, 7].

One of challenges in this study is to discover or characterize what chromatin mark combinatorial patterns can affect the process of gene expression, further revealing complex gene expression mechanisms in downstream analysis [8–13]. This topic have attracted extensive attentions [14–16], however, up to now it is still limited knowledge to understand the degree of complexity of “histone code”. Recent studies have shown that machine learning-based methods can statistically offer higher prediction power to predict gene expression levels, and it can be considered as a promising way to reveal some interesting results in many cases ([17, 18], Devadas L, Yen A, Kellis M. Various localized epigenetic marks.predict expression across 54 samples and reveal underlying chromatin state enrichments. 2015; bioRxiv 030478, unpublished). For example, Chen et al., utilized support vector machine to train a prediction model for each bin. The results demonstrated that all bins are useful to predict gene expression levels, but they are not equally informative. In order to investigate the higher-order interactive relationship between chromatin features, they modeled an interaction model \( y \sim {\sum \nolimits }_{i} {{x_{i}}} + {\sum \nolimits }_{i < j} {{x_{i}}{x_{j}}} \) to predict gene expression levels, where the expression level *y* as a linear combinations of the interactions between individual histone modification features *x*_{
i
} and their products *x*_{
i
}*x*_{
j
} [19]; Dong et al., established a two-step model using linear regression model and random forest method to reveal the relationship between chromatin features and gene expression levels across various cellular contexts. The best bin was selected to represent the remaining signals for each histone modification. The predictor matrix was formed from the best bin for each histone and the whole gene expression levels [20]; Zhou et al., developed a linear mixed model to evaluate the association of each and joint contribution of the five marks with gene expression levels. The marginal effects of each mark are the summation of all window size. The higher-order interactions between markers were also studied by considering them as the covariates in linear mixed model [21].

To naturally characterize higher-order interactions between different markers, tensor representation (also called multilinear or *N*-way) are frequently introduced to model higher-order interactions in different research fields [22, 23]. More recent studies ([24], Khan SA, Ammad-ud-din M. tensorBF: an R package for Bayesian tensor factorization. 2016; bioRxiv 097048, unpublished) leveraged tensor representation to integrate different omics, environmental, and phenotypic data sets to uncover unclear biological problems; Also, our previous work [25] used tensor representation to identify transcription factor binding sites. All results from these applications are demonstrated that tensor representation enable to achieve a powerful performance.

## Methods

### Data sets and pre-processing

In this study, we used the real data sets which are from lymphoblastoid cell lines (LCLs) over three species, namely Human (GSE47991 and GSE19480), Rhesus Macaque (GSE60269), and Chimpanzee (GSE60269), and these data set are all available in Gene Expression Omnibus (GEO). For each species, eight individuals were considered, and for each individual, 26,115 genes were considered in our experiments.

- (i)
We used COVERAGEBED tool [28] to convert the reads into the given window for each mark and each individual, and then we normalized the peak read counts for each individual of each mark by subtracting the number of mapped reads divided by total number of mapped reads and input reads divided by total number of input reads (the detailed procedure see the reference [21]);

- (ii)
We used logarithmic transformation log2(

*x*+*θ*) to normalize the data. In order to obtain the optimal parameter*θ*in prediction phrase, we divided the whole data set into two parts. One-third of data set was used to optimize the parameters*θ*, then the optimized*θ*^{∗}was added to the same modification of the remaining two-thirds of data set to train the prediction model and test their performance.

### High-order representation

In this section, we give more detailed description how the original data were represented by a higher-order representation. The first step is to divide the each gene body, which flanking TSSs both sides with 10k base pairs, into different bins for each individual and each mark (the first two steps of Fig. 1). For each gene, the data for each bin and each histone marker was reformulated into a three-dimension data structure (genes × marks × bins), therefore, each gene was represented by a matrix instead of a vector. The gene expression levels were used the averaged values across 8 individuals. The detailed process to form tensor data is showed in Fig. 1.

As shown in Fig. 1, the first step is to show the different signal patterns near the TSSs over different marks. The histone marks H3K4me3 and Pol II show more informative, while H3K27me3 and H3K3me1 show weaker informative. Each mark was represented by multiple bins (e.g., 41 bins). Therefore, we first combined the five marks into a matrix for each gene. In third step (Fig. 1), we used the contour of distribution of each gene to represent its signals. Assume we have 26115 genes. Finally, we collected all genes to form a tensor data \(\mathcal {X}\).

### High-order model and algorithm

Higher-order partial least squares (or *N*-way partial least squares, NPLS) was proposed by Bro et al. [29]. It is adapted to high-order data
. Here, *N* is the number of order for high-order data \({\mathcal {X}}\) (in our case, *N*=3), and the variable *I*_{
i
} represents the dimensionality of the mode *i*. The response variable
is the averaged gene expression levels across eight individuals for all marks.

**p**

_{1}and

**p**

_{2}such that:

**t**

_{1}and the response variable

**Y**. The and is the loading vector for mode 2 and 3, respectively and \({\mathcal {E}}\) is the residual of data \(\mathcal {X}\) after the first extraction. For the given number of factors

*f*, the predicted variable can be estimated by the equation,

**T**=(

**t**

_{1},

**t**

_{2},⋯

**t**

_{ f }),

**b**is the regression coefficient with respect to

**T**. The MATLAB code of

*N*-way partial least squares is freely available at: http://www.models.life.ku.dk/source/nwaytoolbox.

*x*+

*θ*

_{ k }) with respect to the parameter

*θ*

_{ k }(the parameter for

*k*th bin). In order to determine the optimal \(\theta ^{*}_{k}\) for each bin, we divided the whole data into two parts: one-third of dataset was used for finding the optimal parameter \(\theta ^{*}_{k}\) and then the same \(\theta ^{*}_{k}\) was added to the corresponding bins in the remaining data set. The gene expression levels

**Y**was also logarithm-scaled using the equation log2(

**Y**). A high-order multivariate regression model was developed using the logarithm-scaled training data set and the 10-fold cross validations was used to avoid the over-fitting in training model phrase. Finally, the performance of methods were measured by the Pearson’s correlation coefficient (R) and root mean square error (RMSE).

## Results

In our experiments, to avoid the risk of over-fitting to training prediction model and obtain the reliable results, we used 10-fold cross-validation (with 10 random splitting replicates) in which nine parts are used for training the prediction model, while the remaining part for testing the performance of learned model (similar to previous study [30]).

**Y**is the real gene expression levels, while \({ \hat {\mathbf {Y}}}\) is the predicted the expression levels;

*I*

_{1}is the number of genes;

*μ*

_{ Y }is the mean of gene expression levels

**Y**.

### Experiments on simulation data sets

Our first experiments were conducted on a series of simulation data sets. Interestingly, we found that the distribution of each mark is similar to the normal distribution but does not exactly the same as it. Therefore, in this simulation experiments, we simulated the five different histone marks according to their expression levels, which are described the height *h* of their distributions, *h*=3,2,1,0.5,0.4. 10,000 genes were simulated and the distributions of expression levels for each gene were divided into 100 bins to represent the expression levels of corresponding simulated histone. Similar to the PLS model, we also introduced a latent variable *Z* to simulate the data in simulation experiments.

*Z*∼

*N*(0,

*β*), then

*ε*

^{ y }∼

*N*(0,0.2) and \(\varepsilon _{i}^{\sigma } \sim N\left ({0,{\gamma _{i}}} \right), i=1,2,3,4,5\).

The simulation data \(\mathcal {X}\) with 100 bins can be obtained by partitioning the density function \(N \left ({0,\sigma _{i}^{2}} \right)\) into 100 intervals and calculating their area of corresponding intervals for each mark.

*β*. As shown in Fig. 3, our method (npls) steadily outperforms others on both criteria, i.e., R (left) and RMSE (right). The second-best performance is achieved by random forest regression. Compared with other three methods, linear regression method linearly increases on RMSE.

### Experiments on real data sets

Our second experiments were conducted on real data sets. In this section, we investigated the relationship between gene expression level and chromatin features based on three species (Humans, Chimpanzees, and Rhesus Macaques). The results from the previous work [21] have shown that five marks are significantly enrich near TSSs regardless of species, and the enrichments pattern is robust with respect to the choice of the size of the TSS regions.

The performance of different models on three species data sets

Linear model | Random forest | Support vector machine | NPLS(41bins) | NPLS(21bins) | |
---|---|---|---|---|---|

Hum | 0.769(2.43) | 0.775(2.46) | 0.774(2.46) | 0.784(2.37) | 0.787(2.35) |

Chi | 0.756(2.52) | 0.767(2.47) | 0.765(2.53) | 0.780(2.41) | 0.784(2.39) |

Rhe | 0.760(2.52) | 0.761(2.51) | 0.765(2.54) | 0.774(2.46) | 0.778(2.43) |

## Discussion and conclusion

In this paper, we proposed a higher-order representation method for predicting gene expression levels from chromatin state enrichments. The effectiveness of proposed method was validated by a series of simulation and real data sets. Our method can outperforms others, most likely because higher-order representation method can integrate more unknown interaction information than standard representation method. These results again demonstrate that the gene expression levels are strongly correlated with the combination of chromatin markers.

Each chromatin state shows specific functional and annotation. Our study provides a way to study genomic annotation via chromatin mark combinations, which can extend the epigenetic functional interpretation of the human genome. Therefore, our further work is to incorporate epigenetic factors into the downstream analysis, such as gene expression analysis [9, 31], GO ontologies [32, 33], and disease-related ncRNAs [34].

## Declarations

### Acknowledgements

We would like to acknowledge Dr. Xiang Zhou at University of Michigan, Ann Arbor for very helpful guide to pre-process the real data sets and simulation data sets.

### Funding

Publication costs were funded by the Fundamental Research Funds for the Central Universities (Grant No. 3102017OQD098).

### Availability of data and materials

All data set used in this work are available in Gene Expression Omnibus (GEO): Human (GSE47991;https://doi.org/https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE47991 and GSE19480;https://doi.org/https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE19480), Rhesus and Chimpanzee Macaque (GSE60269;https://doi.org/https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE60269). The MATLAB code of higher-order partial least squares is freely available at: https://doi.org/http://www.models.life.ku.dk/source/nwaytoolbox.

### About this supplement

This article has been published as part of *BMC Bioinformatics* Volume 19 Supplement 5, 2018: Selected articles from the Biological Ontologies and Knowledge bases workshop 2017. The full contents of the supplement are available online at https://doi.org/https://bmcbioinformatics.biomedcentral.com/articles/supplements/volume-19-supplement-5.

### Authors’ contributions

SS and XF conceived and wrote the manuscript. SS implemented the software and analyzed the data, and YZ analyzed the data. All authors read and approved the final manuscript.

### Ethics approval and consent to participate

Not applicable.

### Consent for publication

Not applicable.

### Competing interests

The authors declare that they have no competing interests.

### Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

## Authors’ Affiliations

## References

- Bannister AJ, Kouzarides T. Regulation of chromatin by histone modifications. Cell Res. 2011; 21:381–95.View ArticlePubMedPubMed CentralGoogle Scholar
- Kouzarides T. Chromatin modifications and their function. Cell. 2007; 128:693–705.View ArticlePubMedGoogle Scholar
- Karlic R, Chung HR, Lasserre J, Vlahovicek K, Vingron M. Histone modification levels are predictive for gene expression. Proc Natl Acad Sci USA. 2010; 107:2926–31.View ArticlePubMedPubMed CentralGoogle Scholar
- Santos-Rosa H, Schneider R, Bannister AJ, Sherriff J, Bernstein BE, Emre NCT, et al. Active genes are tri-methylated at K4 of histone H3. Nature. 2002; 419:407–11.View ArticlePubMedGoogle Scholar
- Ruthenburg AJ, Allis CD, Wysocka J. Methylation of lysine 4 on histone H3: Intricacy of writing and reading a single epigenetic mark. Mol Cell. 2007; 25:15–30.View ArticlePubMedGoogle Scholar
- Mikkelsen TS, Ku MC, Jaffe DB, Issac B, Lieberman E, Giannoukos G, et al. Genome-wide maps of chromatin state in pluripotent and lineage-committed cells. Nature. 2007; 448:553–60.View ArticlePubMedPubMed CentralGoogle Scholar
- Barski A, Cuddapah S, Cui KR, Roh TY, Schones DE, Wang ZB, et al. High-resolution profiling of histone methylations in the human genome. Cell. 2007; 129:823–37.View ArticlePubMedGoogle Scholar
- Sun SQ, Peng QK, Shakoor A, Vol. 9. A Kernel-Based Multivariate Feature Selection Method for Microarray Data Classification; 2014, p. e102541.Google Scholar
- Sun SQ, Hood M, Scott L, Peng QK, Mukherjee S, Tung J, Zhou X. Differential expression analysis for RNAseq using Poisson mixed models. Nucleic Acids Res. 2017; 45:e106.View ArticlePubMedPubMed CentralGoogle Scholar
- Peng J, Xue H, Shao Y, et al.A novel method to measure the semantic similarity of HPO terms. Int J Data Min Bioinforma. 2017; 17:173.View ArticleGoogle Scholar
- Chen L, Jiang Y, et al.DisSim: an online system for exploring significant similar diseases and exhibiting potential therapeutic drugs. Sci Rep. 2016; 5:30024.View ArticleGoogle Scholar
- Peng J, Lu J, Shang X, et al. Identifying consistent disease subnetworks using DNet. Methods. 2017; 131:104–10.View ArticlePubMedGoogle Scholar
- Hu Y, Zhou M, et al. DisSetSim: an online system for calculating similarity between disease sets. J Biomed Semant. 2017; 28:71.Google Scholar
- Strahl BD, Allis CD. The language of covalent histone modifications. Nature. 2000; 403:41–5.View ArticlePubMedGoogle Scholar
- Ernst Jason, Kellis Manolis. Discovery and Characterization of Chromatin States for Systematic Annotation of the Human Genome. Nat Biotechnol. 2010; 28:817–25.View ArticlePubMedPubMed CentralGoogle Scholar
- Wu SH, et al. Independent regulation of gene expression level and noise by histone modifications. Plos Comput Biol. 2017; 13:e1005585.View ArticlePubMedPubMed CentralGoogle Scholar
- Daemen A, Griffith OL, Heiser LM, Wang NJ, Enache OM, Sanborn Z, et al. Modeling precision treatment of breast cancer. Genome Biol. 2013; 14:R110.View ArticlePubMedPubMed CentralGoogle Scholar
- Kim K, Bolotin E, Theusch E, Huang HY, Medina MW, Krauss RM. Prediction of LDL cholesterol response to statin using transcriptomic and genetic variation. Genome Biol. 2014; 15:460.View ArticlePubMedPubMed CentralGoogle Scholar
- Cheng C, Yan KK, Yip KY, Rozowsky J, Alexander R, Shou C, et al. A statistical framework for modeling gene expression using chromatin features and application to modENCODE datasets. Genome Biol. 2011; 12:R15.View ArticlePubMedPubMed CentralGoogle Scholar
- Dong XJ, Greven MC, Kundaje A, Djebali S, Brown JB, Cheng C, et al. Modeling gene expression using chromatin features in various cellular contexts. Genome Biol. 2012; 13:R53.View ArticlePubMedPubMed CentralGoogle Scholar
- Zhou X, Cain CE, Myrthil M, Lewellen N, Michelini K, Davenport ER, et al. Epigenetic modifications are associated with inter-species gene expression variation in primates. Genome Biol. 2014; 15:547.View ArticlePubMedPubMed CentralGoogle Scholar
- Freitas MP, da Cunha EFF, Ramalho TC, Goodarzi M. Multimode Methods Applied on MIA Descriptors in QSAR. Curr Comput Aided Drug Des. 2008; 4:273–82.View ArticleGoogle Scholar
- Guzman E, Baeten V, Pierna JAF, Garcia-Mesa JA. Evaluation of the overall quality of olive oil using fluorescence spectroscopy. Food Chem. 2015; 173:927–34.View ArticlePubMedGoogle Scholar
- Hore V, Vinuela A, Buil A, Knight J, McCarthy MI, Small K, et al. Tensor decomposition for multiple-tissue gene expression experiments. Nat Genet. 2016; 48:1094–100.View ArticlePubMedPubMed CentralGoogle Scholar
- Sun SQ, Zhang XP, Peng QK. A high-order representation and classification method for transcription factor binding sites recognition in Escherichia coli. Artif Intell Med. 2017; 75:16–23.View ArticlePubMedGoogle Scholar
- Liaw A, Wiener M. Classification and regression by randomForest. R News. 2002; 2:18–22.Google Scholar
- Lu ZQJ, Vol. 173. The Elements of Statistical Learning: Data Mining, Inference, and Prediction, 2nd edition; 2010, pp. 693–4.Google Scholar
- Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010; 26:841–2.View ArticlePubMedPubMed CentralGoogle Scholar
- Bro R. Multiway calibration. Multilinear PLS. J Chemometr. 1996; 10:47–61.View ArticleGoogle Scholar
- Sun SQ, Peng QK, Zhang XK. Global feature selection from microarray data using Lagrange multipliers. Knowl-Based Syst. 2016; 110:267–74.View ArticleGoogle Scholar
- Sun SQ, Peng QK. A hybrid PSO-GSA strategy for high-dimensional optimization and microarray data clustering. In: IEEE International Conference on Information and Automation, vol. 105. Hailar: IEEE: 2014. p. 41–6.Google Scholar
- Cheng L, Sun J, Xu W, et al. OAHG: an integrated resource for annotating human genes with multi-level ontologies. Sci Rep. 2016; 6:34820.View ArticlePubMedPubMed CentralGoogle Scholar
- Peng J, Wang H, Lu J, et al. Identifying term relations cross different gene ontology categories. BMC Bioinformatics. 2017; 18:573.View ArticlePubMedPubMed CentralGoogle Scholar
- Hu Y, Zhou M, et al. Measuring disease similarity and predicting disease-related ncRNAs by a novel method. BMC Med Genomics. 2017; 10:71.View ArticlePubMedPubMed CentralGoogle Scholar