Integrative Analysis of Epigenomics Data

CHAPTER

INTEGRATIVE ANALYSIS OF EPIGENOMICS DATA

7

Cenny Taslim1, Stephen L. Lessnick1, 2, Simon Lin3 1

Center for Childhood Cancer and Blood Diseases, Nationwide Children’s Hospital Research Institute, Columbus, OH, United States; 2Division of Pediatric Hematology/Oncology/BMT, The Ohio State University College of Medicine, Columbus, OH, United States; 3Research Information Solutions and Innovation, Nationwide Children’s Hospital, Columbus, OH, United States

INTRODUCTION Understanding the functional consequences of genome-wide variation and how it affects complex human diseases remains a critical challenge for the scientific community. With the advent of nextgeneration sequencing, a wealth of experimental information from diverse “omics” data has been generated enabling us to comprehensively study alterations in the sequences of DNA and how epigenetic modifications alter chromatin structure, modulating gene activation or repression. Gene regulation is a complex process that involves genetic and epigenetic mechanisms. A site-specific transcription factor (TF) is thought to interact with regulatory elements to govern the activity of target genes [1]. The binding of these regulatory factors is restricted by complex chromatin structure. Epigenetic modifications such as histone modifications alter chromatin structure and DNA accessibility, thereby regulating gene expression. For example, promoters are typically associated with histone H3 lysine 4 trimethylation (H3K4me3), while enhancers (distal regulatory elements) often show histone H3 lysine 4 monomethylation (H3K4me1) [2] (Fig. 7.1). State-of-the-art next-generation sequencing (NGS) technologies such as ChIP-seq (chromatin immunoprecipitation sequencing) can provide genome-wide profiles of the localization of transcription factors, global histone modification patterns, and chromatin accessibility depending on the antibody used to pull down the protein. In ChIP-seq experiment, a region in the genome enriched with many sequence reads (also called a peak) represents the location where the target protein binds on the DNA. ChIP-seq peaks can be identified using peak calling algorithms such as MACS2 [3]. For a review of ChIP-seq analysis workflow and other peak calling algorithms, see Ref. [4]. Another NGS technology, RNA-seq (RNA sequencing), provides genome-wide transcriptional activity (Fig. 7.1). Peaks from RNA-seq experiments represent gene expression level and can be analyzed using tools such as DESeq2 [5]. Conesa et al. have laid out a guideline for RNA-seq data analysis [6]. A rigorous mechanistic understanding of transcriptional regulation is needed to integrate the combinatorial effect of transcription factors and their target genes,

Computational Epigenetics and Diseases. https://doi.org/10.1016/B978-0-12-814513-5.00007-6 Copyright © 2019 Elsevier Inc. All rights reserved.

107

108

CHAPTER 7 INTEGRATIVE ANALYSIS OF EPIGENOMICS DATA

FIGURE 7.1 Chromatin structure and histone modification at regulatory elements. Each rectangle represents a genomic window (segment). Peaks in each row are regions showing reads enrichment at the genomic locations where the proteins are bound in ChIP-seq experiment or where the genes are expressed in RNA-seq experiment. Peak height represents the binding intensity or expression level in ChIP-seq and RNA-seq experiment, respectively. (A) Closed chromatin restricting binding of TF resulting in repression of target gene. (B) Accessible promoter regions can be marked by H3K4me3 histone modification, allowing binding of TF which activates target gene. (C) Active enhancer marked by H3K4me1 is bound by TF that activates the transcription of distal gene.

as well as the chromatin organization of their regulatory elements [1]. Several large research projects, including The Encyclopedia of DNA Elements (ENCODE) [7], The Cancer Genome Atlas (TCGA) [8], and the NIH Roadmap Epigenomics Mapping Consortium [9], provide expansive resources of large-scale multiomics data sets that are publicly available for researchers. Many computational and statistical methods have been implemented to perform integrative analysis of multiomics data [2,10e15]. Some studies focus on analyzing multiple histone modifications, others focus on integration of multiple transcription factor binding sites, and there are those that also integrate gene expression analysis. However, there is no universal method that is used to integrate all these multiomics data. This chapter provides a general guideline and overview of methodologies that have been implemented to jointly analyze multiple omics data generated using ChIP-seq and RNA-seq technologies. We will focus on methods for integrative analysis of transcription factors, histone modifications, and gene expression. This chapter does not intend to provide a detailed description of the methods or a comparison of the methods. Instead it attempts to provide general overview on how these machine learning/statistical methods were applied to jointly analyze multiple NGS data sets and answer different biological questions. For more detailed implementation of the analysis tools presented here, please refer to the following references [10e17]. It is important to note that there are many other methods developed to integrate multiple omics data that are not presented here. Table 7.1 lists all the computational tools described in this chapter.

QUALITY CONTROL AND DATA PREPROCESSING

109

Table 7.1 A Curated List of Computational Tools and Libraries Methods Used

Description

Multinomial logistics regression and hierarchical clustering Logistics regression

Semisupervised classification of differentially regulated genes into four classes Supervised classification of differentially regulated genes into two classes with distinct regulation mechanism Identification of novel and alternative promoter regions based on shapes of transcription factors and histone modification binding patterns Identification of differentially expressed genes by modeling both gene expression and histone modification binding using Bayesian mixture model Unsupervised method to learn and characterize chromatin states using hidden Markov model Discriminant random forest method used to identify genome-wide enhancers Method based on artificial neural network to perform high-dimensional clustering of multiple transcription factor binding sites studying patterns of colocalization Prediction of transcription factors and histone marks binding in multiple cell types and the effect of single nucleotide variants on binding using deep learning based algorithm

Finite mixture model

Bayesian mixture model

ChromHMM

RFECS Self-organizing map (SOM)

DeepSEA

Software Availability/URL

References

glmnet (R package)

[16]

glm (R package)

[17]

MATLAB (global optimization toolbox)

[14]

epigenomix (R package)

[12]

Java-based ChromHMM

[15]

In-house MATLAB scripts Kohonen (R package)

[13]

http://deepsea. princeton.edu/job/ analysis/create/

[10]

[11]

QUALITY CONTROL AND DATA PREPROCESSING Effective integration and analysis of multiple omics data are expected to provide insights into the complex molecular mechanisms of diseases. Unfortunately, technical heterogeneity (or batch effects) such as experiment times, reagents, sequencing depth, sequencing quality, etc. can confound the biological signal and lead to incorrect conclusions [18]. Therefore, all data need to be of high quality and processed uniformly. Uniformity in terms of read length and run type (single- or paired-end) is also important. Common quality control metrics such as number of reads, GC content, and base-calling quality scores are typically measured for all types of raw sequencing reads using FastQC (http://www. bioinformatics.babraham.ac.uk/projects/fastqc/). Postmapping quality metrics specific for ChIP-seq and RNA-seq experiments will be discussed below.

110

CHAPTER 7 INTEGRATIVE ANALYSIS OF EPIGENOMICS DATA

Chromatin immunoprecipitation followed by sequencing (ChIP-seq) has become the go-to method for identifying genome-wide proteineDNA interactions (e.g., transcription factor or histone protein binding). Eliminating or replacing poor quality data is important for successful data integration. Several postmapping quality metrics have been developed and standards are emerging [19,20]. Currently, it is recommended for each experiment to have two or more biological replicates with matching input control experiments to ensure reliability and reproducibility. Concordance across replicates can be assessed using irreproducible discovery rate (IDR) [21]. Several other criteria have been developed to assess the quality of experiments: (1) Library complexity which signifies low amount of DNA or problems with library construction can be assessed using nonredundant fraction (NRF) and/or PCR bottlenecking coefficients (PBC) and (2) Signal-to-noise ratio can be used to assess the degree of noise in the input control, while signal in the ChIP sample can be evaluated by crosscorrelation analysis and fraction of reads in peaks (FRiP) [4,19,20]. Different standards are applied depending on whether the target protein binds at specific locations that are punctate (i.e., transcription factor) or at longer chromatin domains (i.e., histone modification) [19]. Additionally, ChIP quality assessment in cells under which the protein is perturbed (i.e., performed in an siRNA, CRISPR, or genetic knockout background) or in cells where repressive marks are present should be interpreted differently [22]. The ENCODE Consortium has developed standards for these quality metrics including the recommended sequencing read depth and number of aligned reads that they keep up-todate to ensure high quality data (https://www.encodeproject.org/data-standards/). RNA-seq is a well-established approach to quantify gene expression genome-wide. Quality assessment of RNA-seq experiment is critical to ensure downstream analysis of biologically relevant data. As with ChIP-seq experiments, two or more biological replicates are recommended. Typically, concordance between replicates is measured at the gene level with a Spearman correlation [6]. Several quality metrics have been proposed to assess RNA-seq data: (1) Capture efficiency metric, which measures the percentage of reads that map to the intended target region, can be used to identify low sample quality, library contamination, or inefficient removal of ribosomal RNA; (2) GC content and 30 e50 bias can help identify nonuniform coverage of transcripts [23,24]. Other quality control metrics such as detection of aberrant splicing have also been implemented in RSeQC [25]. In addition to the quality metrics described above, when consolidating multiple ChIP-seq data for integrated analysis, it is also important to avoid artificial differences due to sequencing depth and mappability [26]. To avoid these artificial differences, all histone and transcription factor ChIP-seq data should have comparable sequencing depth, or uniformly subsampled. Read counts should be standardized. Additionally, reads map to multiple locations should be discarded and regions associated with repetitive regions should be excluded [2]. Finally, dimension reduction approaches such as principal component analysis (PCA), which decomposes high-dimensional data into components that explain most of the variation in the data, can be used to identify batch effects or outliers [27].

RELATIONSHIP BETWEEN HISTONE MODIFICATION PATTERN, TRANSCRIPTION FACTOR BINDING, AND MRNA EXPRESSION LEVEL A critical challenge in understanding disease mechanisms is elucidating the underlying regulation of gene expression. The widespread application of the next-generation sequencing method has generated large-scale genome-wide molecular profiles, enabling researchers to study the combinatorial effects of

RELATIONSHIP BETWEEN HISTONE MODIFICATION PATTERN

111

chromatin organization, regulatory elements, and transcription factor binding that govern gene expression and underlie cell state changes [7,28]. Transcription factors (TFs) bind to specific DNA motifs in the promoter or enhancer region and regulate transcription either through proximity to the transcription start site (TSS) or higher-order chromatin looping [29,30]. Epigenetic histone modifications are also known to regulate gene expression by either reorganizing the local chromatin structure to control accessibility of TF binding sites [31] or by recruiting other chromatin remodelers [32]. Thus, nucleosomes in the vicinity of transcriptional regulatory elements (e.g., enhancer or promoter regions) typically contain histones with specific posttranslational modifications (PTM), such as H3 lysine 4 mono- and trimethylation (H3K4me1/H3K4me3) (Fig. 7.1). Dysregulation of gene expression caused by aberrant TF or histone modification has been directly associated with many diseases including Ewing sarcoma [16], prostate cancer [17], and various other disorders [33]. Here, we focus on methods for integrative analysis of ChIP-seq (for measuring TF binding and histone modification) and RNA-seq (for measuring mRNA expression). In this section, we review several existing statistical methods used to investigate the relationship between histone modification, TF binding, and gene expression.

REGRESSION ANALYSIS To understand the effect of TF binding and chromatin regulation on gene expression, Tomazou et al. first identified the differentially regulated genes, comparing cell lines with two different conditions [16]. The expression of these genes was then correlated to the histone binding intensities at their promoter regions. For each histone mark, the maximum peak score within a 1e5-kb window (depending on the type of histone mark) around the TSS was assigned, resulting in a matrix of genes  histone modification. A weighted discretized version of this matrix was then used as input to perform unsupervised clustering. In a semisupervised approach, 40e50 genes from four distinct clusters were selected based on unsupervised clustering and expert classification. These genes were then used to build a multinomial logistic regression to model the four clusters. This model was used to predict the cluster assignment of remaining genes. Using this approach, the authors identified four clusters of genes: Cluster 1, genes associated with four active promoter marks but not repressive mark; Cluster 2, genes associated with two active promoter marks; Cluster 3, genes associated with one active promoter mark; Cluster 4, genes associated with repressive mark in the presence or absence of active promoter marks. Genes in the same group that respond differently to a transcription factor were associated with distinct molecular functions. For example, a gene identified in Cluster 1 whose expression is anticorrelated with a transcription factor is illustrated in Fig. 7.2. In their implementation, genes were clustered based on both supervised (i.e., expert classification) and unsupervised (i.e., hierarchical clustering) approaches. Taslim et al. [17] used a logistic regression model to perform a supervised classification to study the effect of distance and combinatorial binding of transcription factors, and active histone modification on regulation of gene expression. Specifically, a logistic model was built using variables such as distance from regulated genes to transcription factor binding sites and the presence of overlap between several transcription factors and active histone mark binding sites, with expression of regulated genes as response. Initially, models with up to six variables were considered, followed by selection of a model with the optimal specificity and sensitivity in a fivefold cross-validation. Bayes factor like criterion was used to perform classification on the regulated genes. Bayes factor is defined as a ratio of

112

CHAPTER 7 INTEGRATIVE ANALYSIS OF EPIGENOMICS DATA

FIGURE 7.2 Illustration of a gene characterized by four active promoter marks which has negative correlation with TF. Gene model is depicted on top with exons represented by blue boxes and introns as blue lines.

probability of null and alternative hypothesis. It provides an assessment that represents how well a hypothesis predicts the empirical data relative to the alternative [34]. Using this approach, the authors identified a group of treatment-responsive genes which has a relatively short long-range regulation mechanism facilitated by a combination of transcriptional factors (including a TF that regulates chromatin structure) and histone modifications.

MIXTURE MODEL In order to identify novel and alternative promoter regions, Taslim et al. investigated the shapes of RNA polymerase II and active histone modification binding patterns [14]. Combining binding patterns of histone modification that has been shown to be enriched in active promoter regions with transcriptional factor such as RNA polymerase II (Pol-II) and RNA-seq expression provide more specific identification of patterns within promoter regions. Furthermore, phenomenon such as RNA Pol-II stalling, in which the Pol-II binds at the promoter but the gene is not transcribed, may display binding patterns specific to their functions. Therefore, modeling the shapes of these binding patterns may provide additional insights. Specifically, normalized ChIP-seq read counts of Pol-II and active histone mark in the 10-kb regions around TSSs were segmented into 100-bp bins. These binding patterns were

RELATIONSHIP BETWEEN HISTONE MODIFICATION PATTERN

113

then grouped using K-means clustering to identify the four classes of promoter binding patterns. A mixture model with two double-exponential and uniform components was fitted to the averaged binding patterns within a group. The double-exponential component is designed to capture the shape of Pol-II and active histone mark that are unimodal and bimodal, respectively. The uniform component is used to capture the tails of Pol-II binding profile. The mixture model was fitted using a generalized pattern search algorithm [35] with KullbackeLeibler distance metrics [36]. This model of combined patterns representing binding signatures for promoter regions was then used to search for unknown genomic regions that may represent promoters of unknown genes or known genes with novel alternative promoters. The authors also showed that genes whose promoters display the combined binding patterns were associated with higher expression compared to those who do not. Fig. 7.3, box 2 shows an example of a gene which displays the combined binding patterns of active histone and TF. Joint analysis on histone ChIP-seq and RNA-seq expression data was performed by Klein et al. using Bayesian mixture model to identify genes with significant differences in both expression and histone modification bindings comparing two different conditions [12]. For each differential gene, expressions and histone modifications were integrated into a correlation variable (Z) which was

FIGURE 7.3 Illustration of the scanning process to find regions with combined TF and active histone patterns. The models of active histone and TF patterns are shown as blue and red curves, respectively. These fitted patterns were used to scan the genome. The actual binding data are shown as gray areas, top row for active histone and bottom row for TF. Box 1 shows a gene associated with TF but not active histone pattern. This gene is shown to have less expression than the gene in box 2 which displays both TF and active histone patterns.

114

CHAPTER 7 INTEGRATIVE ANALYSIS OF EPIGENOMICS DATA

calculated by multiplying the standardized difference of expressions and binding intensities in two distinct conditions. The distribution of Z is then fitted using Bayesian mixture model with normal and exponential distributions. Active genes associated with increased binding of active histone marks are considered to be in the same direction and called positive component. Repressive genes associated with increased binding of repressive histone marks are thought to be in the opposite direction and called negative component. Both of these components are modeled by exponential distribution. The normal component represents genes that do change in expressions and histone bindings. The proportions of the mixture model are modeled with a Dirichlet prior distribution. Klein et al. used a Bayesian approach in which informative prior distributions need to be selected [12]. The model is fitted using Gibbs sampling which is based on Markov Chain Monte Carlo (MCMC) techniques [37,38], and model estimates were obtained by averaging the posterior distribution. Using this approach, for example, activated (repressed) genes associated with increased binding of active (repressive) histone marks can be identified.

IDENTIFICATION OF FUNCTIONAL REGULATORY REGIONS Posttranslational histone modifications have been shown to change nucleosome dynamics allowing or blocking DNAeprotein binding and thereby affecting the transcriptional activity of target genes [39,40]. There are a wide variety of posttranslational histone modifications, and it was suggested that multiple histone modifications act in a combinatorial fashion to specify distinct functional regulatory regions [41]. Here, we describe two methods that have been developed to integrate multiple histone modification ChIP-seq data, to identify and define functional regulatory regions such as active promoters, enhancers, repressed promoters, etc. ChromHMM [15], which is based on hidden Markov model (HMM), was developed to model the combinatorial patterns of multiple chromatin elements. The goal is to capture the similarities of chromatin patterns across different cell types giving an unbiased classification of the genome into functional regions (i.e., chromatin states are modeled as latent variables in HMM). In the implementation of HMM, the genome was segmented into 200-bp regions (default value) that roughly approximate nucleosome sizes and are used as the “time” axis. Thus as input, multiple histone modification and general TF (and chromatin accessibility assays such as DNase-seq and FAIRE-seq) data were normalized and converted into a binary presence or absence call for each mark at 200-bp resolution [42]. ChromHMM was trained with two-stage nested parameter initialization using Euclidean distance for pruning the states. Models with up to 30 hidden states from multiple observed combinations of chromatin marks and 200 training iterations were completed to learn the parameters of the final HMM. Biologically these hidden states may represent different genomic features such as strong/weak enhancers, repressed/active promoters, open chromatin, etc. The number of states chosen was a compromise between capturing all the combinatorial complexity of chromatin marks and models that are useful for interpretation of biologically meaningful functional regions. The trained model was then used to compute the posterior probability of each state for each segment. Each region was then labeled using the state with the maximum posterior probability. The states discovered by the model were then annotated by the authors using large-scale systematic data-mining effort. An example of a state identified by ChromHMM is “poised promoter” regions that are associated with both activating and repressing histone modifications. In this state, genes are silenced but poised for activation.

MULTIPLE TRANSCRIPTION FACTORS USING SELF-ORGANIZING MAP

115

RFECS (Random Forest based Enhancer identification from Chromatin States) was developed to integrate multiple histone modification profiles to identify enhancers. Enhancers are distal regulatory elements where transcription factors bind to and regulate gene expression [43]. Transcriptional regulations via enhancer elements are thought to be cell type-specific and dynamically regulated during development [44,45]. Identification of enhancers has been challenging due to a lack of common features and their locations that are often distal from their target genes. Rajagopal et al. [13] implemented a discriminant random forest based algorithm to identify genome-wide enhancers. RFECS was also used to identify the most informative set of histone modifications that can accurately predict enhancers. For training, only histone modifications that provide largely nonredundant information were selected and used as features. RFECS was constructed using binary classification trees. Specifically, each histone profile in bins of 100-bp surrounding binding sites of a protein called p300 known to be recruited to enhancers was used to train enhancer class. Similarly, histone profile around non-p300 binding sites was used to train nonenhancer class. At each node in the tree, a random subset of features was selected and parametric multivariate linear discriminant analysis (LDA) was applied to split the tree nodes. The optimal number of trees was determined by examining the area under the curve (AUC) generated using fivefold cross-validation. The final class predictions (enhancer vs. nonenhancer) were decided by majority voting. Variable importance (i.e., importance of individual histone modifications) was estimated using out-of-bag measure implemented in MATLABÒ by Rajagopal et al. [13]. The most informative sets of histone modification for enhancer identification were chosen based on the ordering of variable importance across multiple replicates and cell types.

ASSOCIATION BETWEEN MULTIPLE TRANSCRIPTION FACTORS USING SELF-ORGANIZING MAP (SOM) Different TFs need to co-operate at specific loci to perform complex regulation of their target genes. Investigating the colocalization of a combination of transcription factors to achieve complex and accurate regulation of target genes is challenging due to the high dimensionality of the data. For 128 TFs, there are more than 1038 possible combinations. In order to take into account the full combinatorial co-binding of many TFs, Xie et al. [11] used a type of artificial neural network called self-organizing map (SOM) to perform high-dimensional clustering based on the combinations of up to 128 TF binding profiles. An artificial neural network is inspired by the neural network in the brain, and as such, consists of layers of fully connected neurons that nonlinearly transformed input [46]. SOM provides a way to reduce the dimensionality of data using a compression technique known as vector quantization that also retains the topological properties [47]. In the study of TF colocalization, a block defined as the maximum overlap of TF binding sites is discretized as either bound or not bound for each TF and used as input to SOM. Each neuron in SOM is associated with weights the same dimension as the input and a position in the lower-dimensional map space. The weight for each neuron and its position (and also of its neighbors’) is updated through unsupervised learning process that decays over time. The end result is a map comprised of neurons with distinct colocalization patterns at different genomic loci. Neurons with similar colocalization patterns will be in close proximity within the map space. This lower-dimensional map allows for clustering of regions which have the same TFs bound and therefore also identifies regions bound by high number of TFs (hotspot regions).

116

CHAPTER 7 INTEGRATIVE ANALYSIS OF EPIGENOMICS DATA

PREDICTION OF CHROMATIN AND TRANSCRIPTION BINDING SITES DIRECTLY FROM DNA SEQUENCES USING DEEP LEARNING Genetic variations associated with disease often appear in noncoding parts of the genome. Accumulating evidence has shown that these variants occur in regulatory elements affecting chromosome structure and the modulation of genes [48]. Unfortunately, functional characterization of these noncoding variants remains a challenge. Classical machine learning methods such as hidden Markov model (ChromHMM) and random forest (RFECS) described in previous sections have been used to identify and annotate these regulatory regions [13,15]. However, these methods do not operate directly on the DNA sequences and therefore require preprocessing features (in non-deep learning approach, typically binding regions were first associated with presence or absence of a binding factor and then used as input). Recent advances in artificial neural network lead to the development of state-of-the-art networks with large number of hidden layers known as deep networks. A deep network takes raw data at the lower (input) layer and automatically learns more abstract features that are needed for classification (or regression) through multiple levels of representation learning. These abstract features were obtained by combining and performing nonlinear transformation of outputs from previous layers [49]. These layers of features, which are automatically learned from data, can discover highly complicated structures in high-dimensional data. Deep learning has been successfully applied to image analysis and speech recognition with record breaking results [50,51]. Zhou et al. developed DeepSEA which uses a deep learning architecture called convolutional neural network (CNN) to learn DNA regulatory sequences (motifs) for predicting protein binding sites [10]. Rather than preprocessing relevant features, the network learns regulatory motifs directly from DNA sequences, enabling prediction of binding sites with single nucleotide sensitivity. The model learnt can then be used to improve noncoding variant interpretation. Zhou et al. start by training the model using binding site information from 919 features (125 DNase, 690 TFs, and 104 histones in multiple cell types) from ENCODE and Roadmap Epigenomics projects [2,7]. Each training sample is a 1000-bp sequence centered on 200-bp bin of binding sites. Each 1000-bp sequence (with at least one binding event) is paired with a labeled vector of 919 features (classes). A class is labeled 1 if more than half of 200-bp bin overlap with the peak regions, or 0 otherwise. As a deep learning model capable of learning directly from DNA sequences, the 1000-bp sequence is represented by a 1000  4 binary matrix corresponding to the four nucleotides (Fig. 7.4). Data sets were split into nonoverlapping training, validation, and testing sets. Validation sets were used to select hyperparameters for the network. DeepSEA consists of sequential, alternating convolution and pooling layers which allow learning the sequence motifs at increasing spatial scales. These layers are followed by a fully connected layer which integrates information from the full-length sequences and a final sigmoid layer that outputs the probability of each input sequence belonging to each individual class (Fig. 7.4). The initial convolution layer can be considered as learning a set of position weight matrices (PWMs) [52] that minimizes prediction error. A PWM represents the common sequences (motif) where the target protein is bound. This layer also computes PWM scores with a 1 step moving window and thus provides invariance to small sequence shifts. Higher level layers receive input from lower layers and therefore learn more complex patterns from larger spatial ranges. DeepSEA shares learned predictive regulatory sequences across all classes. This means, for example, a regulatory sequence that

DISCUSSION

117

FIGURE 7.4 Illustration of DeepSEA algorithm. DNA sequences in each peak region are represented as a 1000  4 binary matrix corresponding to the four nucleotides and used as input to DeepSEA. The initial convolution layer can be thought of learning a set of optimal position weight matrices that represents the local binding motif. Higher level layers receive input from lower layers and therefore learn more complex motifs at increasing spatial scales. These layers are followed by fully connected layers which integrate motif information from the full-length sequences and a final layer that outputs the class prediction.

is effective for classification of a specific TF can be used to predict the bindings of another physically interacting protein. On the holdout test sequences, DeepSEA was able to perform classification with high accuracy (median AUC ranges from 0.856 to 0.958). In order to identify functional SNPs, the DeepSEA model was first used to compute 919 class prediction for every pair of sequences carrying either the reference or the alternative allele. DeepSEA predictions for each SNP and evolutionary conservation scores were used to train boosted logistic regression classifiers. This functional prediction on SNPs was shown to outperform previous methods [10].

DISCUSSION As technology advances and more data become available, integrative analysis will become critical as it provides a comprehensive view on the complexity of diseases’ molecular mechanisms. The workflow to learn functional relationships from molecular data involves the following general steps: data cleaning, preprocessing, features extraction, model fitting or clustering, and finally evaluation of the results. Due to the high dimensionality of the data and typically limited number of samples, integrative analysis requires novel statistical and computational solutions. Many methods have been proposed but there is no single method that is universally applicable. Generally, these methods are applied to either make a prediction or to discover the underlying relationship between the variables of interest. For example, a researcher may be interested to understand what factors regulate a gene regulation and how they influence this modulation or simply to predict whether or not a gene is going to be up-/downregulated. There are trade-offs associated with these goals. Optimizing prediction accuracy often comes from decreased interpretability of the model. Model-based methods (e.g., logistics regression, mixture model), although simple and usually fast computationally, require assumptions about the underlying relationships between input and output variables. Model-free algorithms do not make

118

CHAPTER 7 INTEGRATIVE ANALYSIS OF EPIGENOMICS DATA

strong assumptions about the underlying function; however, they are slower to compute and are typically more difficult to interpret. In both types of algorithms, selecting the most informative features are critical for performance and this is essentially the bottleneck, especially for high-dimensional omics data. State-of-the-art deep learning algorithm enables automatic features extraction that is scalable to high-dimensional data, making it a very attractive algorithm for application in genomics. However, deep learning requires a large amount of well-annotated data for the system to successfully learn features automatically from the data. Early applications in genomics have been to study the variation between regions, segmenting the genome to obtain large data sets for training. All of these methods can be a powerful complement to each other. We expect further improvement will help provide deeper understanding and eventually help solve important biological problems.

ACKNOWLEDGMENTS The authors would like to thank Kirsten M. Johnson for critical reading of the manuscript and Trisha D. Shah for creating Fig. 7.1.

REFERENCES [1] Voss TC, Hager GL. Dynamic regulation of transcriptional states by chromatin and transcription factors. Nat Rev Genet 2013;15:69e81. [2] Roadmap Epigenomics Consortium, Kundaje A, Meuleman W, Ernst J, Bilenky M, Yen A, et al. Integrative analysis of 111 reference human epigenomes. Nature 2015;518:317e30. [3] Feng J, Liu T, Qin B, Zhang Y, Liu XS. Identifying ChIP-seq enrichment using MACS. Nat Protoc 2012;7: 1728e40. [4] Nakato R, Shirahige K. Recent advances in ChIP-seq analysis: from quality management to whole-genome annotation. Brief Bioinform 2017;18:279e90. [5] Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 2014. https://doi.org/10.1101/002832. [6] Conesa A, Madrigal P, Tarazona S, Gomez-Cabrero D, Cervera A, McPherson A, et al. A survey of best practices for RNA-seq data analysis. Genome Biol 2016;17:13. [7] ENCODE Project Consortium. An integrated encyclopedia of DNA elements in the human genome. Nature 2012;489:57e74. [8] Cancer Genome Atlas Research Network, Weinstein JN, Collisson EA, Mills GB, Shaw KRM, Ozenberger BA, et al. The cancer genome Atlas pan-cancer analysis project. Nat Genet 2013;45:1113e20. [9] Bernstein BE, Stamatoyannopoulos JA, Costello JF, Ren B, Milosavljevic A, Meissner A, et al. The NIH Roadmap epigenomics mapping consortium. Nat Biotechnol 2010;28:1045e8. [10] Zhou J, Troyanskaya OG. Predicting effects of noncoding variants with deep learning-based sequence model. Nat Methods 2015;12:931e4. [11] Xie D, Boyle AP, Wu L, Zhai J, Kawli T, Snyder M. Dynamic trans-acting factor colocalization in human cells. Cell 2013;155:713e24. [12] Klein H-U, Scha¨fer M, Porse BT, Hasemann MS, Ickstadt K, Dugas M. Integrative analysis of histone ChIP-seq and transcription data using Bayesian mixture models. Bioinformatics 2014;30:1154e62. [13] Rajagopal N, Xie W, Li Y, Wagner U, Wang W, Stamatoyannopoulos J, et al. RFECS: a random-forest based algorithm for enhancer identification from chromatin state. PLoS Comput Biol 2013;9:e1002968. http:// enhancer.ucsd.edu/renlab/RFECS_enhancer_prediction/.

REFERENCES

119

[14] Taslim C, Lin S, Huang K, Huang TH-M. Integrative genome-wide chromatin signature analysis using finite mixture models. BMC Genomics 2012;13:S3. [15] Ernst J, Kellis M. ChromHMM: automating chromatin-state discovery and characterization. Nat Methods 2012;9:215e6. http://compbio.mit.edu/ChromHMM/. [16] Tomazou EM, Sheffield NC, Schmidl C, Schuster M, Scho¨negger A, Datlinger P, et al. Epigenome mapping reveals distinct modes of gene regulation and widespread enhancer reprogramming by the oncogenic fusion protein EWS-FLI1. Cell Rep 2015;10:1082e95. [17] Taslim C, Chen Z, Huang K, Huang TH-M, Wang Q, Lin S. Integrated analysis identifies a class of androgenresponsive genes regulated by short combinatorial long-range mechanism facilitated by CTCF. Nucleic Acids Res 2012;40:4754e64. [18] Goh WWB, Wang W, Wong L. Why batch effects Matter in omics data, and how to avoid them. Trends Biotechnol 2017;35:498e507. [19] Landt SG, Marinov GK, Kundaje A, Kheradpour P, Pauli F, Batzoglou S, et al. ChIP-seq guidelines and practices of the ENCODE and modENCODE consortia. Genome Res 2012;22:1813e31. [20] Bailey T, Krajewski P, Ladunga I, Lefebvre C, Li Q, Liu T, et al. Practical guidelines for the comprehensive analysis of ChIP-seq data. PLoS Comput Biol 2013;9:e1003326. [21] Li Q, Brown JB, Huang H, Bickel PJ. Measuring reproducibility of high-throughput experiments. Ann Appl Stat 2011;5:1752e79. [22] Marinov GK, Kundaje A, Park PJ, Wold BJ. Large-scale quality analysis of published ChIP-seq data. G3 2014;4:209e23. [23] Sheng Q, Vickers K, Zhao S, Wang J, Samuels DC, Koues O, et al. Multi-perspective quality control of Illumina RNA sequencing data analysis. Brief Funct Genomics 2017;16:194e204. [24] DeLuca DS, Levin JZ, Sivachenko A, Fennell T, Nazaire M-D, Williams C, et al. RNA-SeQC: RNA-seq metrics for quality control and process optimization. Bioinformatics 2012;28:1530e2. [25] Wang L, Wang S, Li W. RSeQC: quality control of RNA-seq experiments. Bioinformatics 2012;28:2184e5. [26] Wei LK, Au A. Computational epigenetics. In: Handbook of Epigenetics. 2nd ed. 2017. p. 167e90. [27] Meng C, Zeleznik OA, Thallinger GG, Kuster B, Gholami AM, Culhane AC. Dimension reduction techniques for the integrative analysis of multi-omics data. Brief Bioinform 2016;17:628e41. [28] Rivera CM, Ren B. Mapping human epigenomes. Cell 2013;155:39e55. [29] Farnham PJ. Insights from genomic profiling of transcription factors. Nat Rev Genet 2009;10:605e16. [30] Dekker J, Marti-Renom MA, Mirny LA. Exploring the three-dimensional organization of genomes: interpreting chromatin interaction data. Nat Rev Genet 2013;14:390e403. [31] Bernstein BE, Meissner A, Lander ES. The mammalian epigenome. Cell 2007;128:669e81. [32] Kouzarides T. Chromatin modifications and their function. Cell 2007;128:693e705. [33] Portela A, Esteller M. Epigenetic modifications and human disease. Nat Biotechnol 2010;28:1057e68. [34] Raftery AE. Hypothesis testing and model selection. In: Gilks WR, Richardson S, Spiegelhalter DJ, editors. Markov Chain Monte Carlo in Practice. Chapman and Hall; 1996. p. 163e88. [35] Audet C, Dennis JE. Analysis of generalized pattern searches. SIAM J Optim 2002;13:889e903. [36] Kullback S, Leibler RA. On information and sufficiency. Ann Math Stat 1951;22:79e86. [37] Ishwaran H, James LF. Gibbs sampling methods for stick-breaking priors. J Am Stat Assoc 2001;96:161e73. [38] Ishwaran H. Markov chain Monte Carlo in approximate Dirichlet and beta two-parameter process hierarchical models. Biometrika 2000;87:371e90. [39] Wierer M, Mann M. Proteomics to study DNA-bound and chromatin-associated gene regulatory complexes. Hum Mol Genet 2016;25:R106e14. [40] Bowman GD, Poirier MG. Post-translational modifications of histones that influence nucleosome dynamics. Chem Rev 2015;115:2274e95.

120

CHAPTER 7 INTEGRATIVE ANALYSIS OF EPIGENOMICS DATA

[41] Wang Z, Zang C, Rosenfeld JA, Schones DE, Barski A, Cuddapah S, et al. Combinatorial patterns of histone acetylations and methylations in the human genome. Nat Genet 2008;40:897e903. [42] Hoffman MM, Ernst J, Wilder SP, Kundaje A, Harris RS, Libbrecht M, et al. Integrative annotation of chromatin elements from ENCODE data. Nucleic Acids Res 2013;41:827e41. [43] Bulger M, Groudine M. Functional and mechanistic diversity of distal transcription enhancers. Cell 2011; 144:327e39. [44] Heintzman ND, Hon GC, Hawkins RD, Kheradpour P, Stark A, Harp LF, et al. Histone modifications at human enhancers reflect global cell-type-specific gene expression. Nature 2009;459:108e12. [45] Bonn S, Zinzen RP, Girardot C, Gustafson EH, Perez-Gonzalez A, Delhomme N, et al. Tissue-specific analysis of chromatin state identifies temporal signatures of enhancer activity during embryonic development. Nat Genet 2012;44:148e56. [46] Angermueller C, Pa¨rnamaa T, Parts L, Stegle O. Deep learning for computational biology. Mol Syst Biol 2016;12:878. [47] Kohonen T. Self-organized formation of topologically correct feature maps. Biol Cybern 1982;43:59e69. [48] Zhang F, Lupski JR. Non-coding genetic variants in human disease. Hum Mol Genet 2015;24:R102e10. [49] LeCun Y, Bengio Y, Hinton G. Deep learning. Nature 2015;521:436e44. [50] Mohamed A-R, Sainath TN, Dahl G, Ramabhadran B, Hinton GE, Picheny MA. Deep Belief Networks using discriminative features for phone recognition. In: 2011 IEEE International Conference on Acoustics, speech and signal processing (ICASSP); 2011. https://doi.org/10.1109/icassp.2011.5947494. [51] Krizhevsky A, Sutskever I, Hinton GE. ImageNet classification with deep convolutional neural networks. Commun ACM 2017;60:84e90. [52] Sinha S. On counting position weight matrix matches in a sequence, with application to discriminative motif finding. Bioinformatics 2006;22:e454e63.