Nature Neuroscience Review

RNA-seq: Initial Methods of Data Processing and Annotation

The raw data produced by RNA-seq (see Figure 1 for initial pipelines of RNA-seq data analysis) is—for each biological sample—tens to hundreds of million short sequences (called reads, typically 50-100 bp) that correspond to random fragments of expressed RNAs present in the original tissue. The first step in analyzing such data is to assess the quality of these reads, which greatly influences downstream bioinformatics outputs. For that purpose, fastqc (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) can be used. Fastqc is a lightweight, highly efficient and low profile (i.e., requires reduced memory and yields less excessive outputs) program that simply requires raw sequencing reads in fastq format for initial quality control (QC) assessments. Fastqc, however, does not address RNA-seq specific questions, such as exonic vs. intronic alignments, transcript detection rates, etc., making QC determinations for subsequent downstream analyses difficult. To do so, RNA-SeQC²⁶ is often used after fastqc, allowing investigators to examine numerous additional parameters associated with RNA-seq sample quality including yield, alignment and duplication rates, GC bias, contaminating rRNA content, regions of alignment (exon vs. intron vs. intragenic), continuity of coverage, 3′/5′ biases and counts of detectable transcripts, among others. Although RNA-SeQC is a comprehensive program that addresses most QC issues, its usage can be cumbersome requiring a great deal of computational power. In addition, ngs.plot²⁷ can be used in connection with RNA-SeQC or fastqc to generate gene body plots for RNA-seq data, thereby providing an intuitive visualization tool for investigators to examine overall coverage patterns of sequencing reads from transcription start sites (TSSs) to transcription end sites (TESs). This tool allows for the identification of common sequencing abnormalities, such as strong 3’ biases. Although underreported, such QC data are extremely important to interpretations of RNA expression in brain, as differences observed in transcript abundance between control and experimental conditions, although important, are often small in magnitude, owing to a variety of unique challenges encountered when working with neural tissues (Box 1). Since low quality sequencing results in decreased signal to noise, as well as potential biases, such small differences may be inadvertently masked or amplified, thereby leading to high false positive and negative rates.

Following initial QC assessments, short sequence reads are mapped to appropriate transcriptomes using splicing-aware alignment tools. Such mapping requires both a reference genome sequence and a description of transcripts (e.g., a format such as GFF) rendering alignment results highly specific to a particular version of the transcriptome, which one can refer to using a database name, as well as a release number (e.g., Ensembl 67). When working with a species in which a comprehensive reference genome is not available, de novo transcriptome assembly can be performed using tools such as Cufflinks²⁸ and Trinity²⁹. A popular choice for RNA-seq data alignment is Tophat³⁰, which wraps around Bowtie as the basic alignment tool but provides additional workflows to accomplish tasks specific to RNA-seq data analysis (e.g., splicing detection). Since spliced alignment is critical to RNA-seq analysis, it has attracted much research effort in recent years. For example, a new alignment program, referred to as STAR³¹, has been gaining significant market share and attention. In a recent study comparing 11 programs for RNA-seq alignment³², STAR achieved impressive performance across numerous benchmarks including alignment yield, basewise accuracy, mismatch and gap placement, exon junction discovery and suitability for transcript reconstruction. STAR is an ultra-fast, sensitive and precise tool that reduces alignment time from 1-2 days to a few hours for RNA-seq samples containing ~100 million reads, based on 4-core workstation processing capabilities. STAR is limited, however, in that its memory footprint can easily reach 32 GB or more, which is prohibitive for many common servers.

Analysis of Transcripts, Splice Variants and Non-coding RNAs

The next step in data analysis is to infer expression levels of individual transcripts from aligned reads. Simply put, if tissue samples derived from control conditions generate on average 1,000 reads for a given RNA, and samples from an experimental condition generate 2,000 reads, one can conclude a 2-fold induction of that gene’s expression. Such analyses are rarely so simple, because of the expression of multiple transcripts from a single gene and because such transcripts share a majority of their sequences, often differing by only a few or dozens of base pairs. Cufflinks³³ solves this problem through the use of a statistical learning approach that assigns a fraction of the total read count to each transcript based on a maximum likelihood principle. In addition, Cufflinks attempts to quantify splicing events by grouping transcripts into TSS groups. A TSS group is a collection of transcripts that share TSSs. Such grouping is based on the rationale that alternative splicing events, such as exon skipping, will only happen between transcripts sharing the same TSS. Two types of splicing events are thereby characterized: 1) alternative splicing, which is defined between different transcripts that are within the same TSS group; and 2) alternative promoter usage, which is defined between transcripts with different TSSs. It should be noted that Cufflinks only reports these events at the level of TSS groups, as well as entire transcripts, and does not provide detailed information regarding which exons demonstrate alternative splicing. In order to obtain more detailed splicing information, Mixture of Isoforms (MISO) model³⁴ can be used. MISO is based on a different design from that of Cufflinks and works on a predefined set of splicing events, such as exon skipping, intron retention, mutual exclusion and alternative 3’ UTR. MISO uses a Bayesian approach, which allows researchers to combine prior knowledge with obtained data for more reliable inferences, to derive relative abundance levels of spliced transcripts. Such detailed information provides an advantage when performing integrative analyses between alternative splicing and other forms of epigenomic regulation, such as the contribution of histone modification states and transcription factor binding events to regulation of alternative splicing. Additional information on alternative transcripts can be obtained from sequencing nuclear RNA as opposed to total cellular RNA. Nuclear RNA contains much larger amounts of sequenced introns, which can provide invaluable information about splicing mechanisms³⁵. RNA-seq analysis of brain tissue has revealed an order of magnitude more alternative transcript production compared to that inferred from older microarray and related technologies¹⁶.

RNA-seq also enables the detection and quantification of several types of non-coding RNAs, which are proving to play crucial roles in biological regulation. A subset of long non-coding RNAs (lncRNAs)—similar to protein-coding genes—contain polyadenylated (polyA) tails, which allow for their detection by RNA-seq regardless of the RNA purification procedure used (e.g., ribozero or polyA selection). Identification of non-polyA lncRNAs, however, can only be accomplished through purification procedures preserving total RNA (i.e., polyA+ and polyA−). Since lncRNAs exist in high abundance in mammals, the Ensembl database has recently incorporated a large number of lncRNAs into its gene collection. Thus, predefined lncRNAs can now be analyzed alongside of protein-coding genes in the same sample. If investigators wish to predict a large number of novel, and still un-annotated, lncRNAs, however, this can be accomplished by utilizing histone modification state data to define candidate regions. To do so, ChIP (chromatin immunoprecipitation)-seq data for euchromatic H3K4me3 (trimethylated Lys4 of histone H3) and H3K36me3 need to first be obtained, as described in the next section. Using the intersection of H3K4me3 and H3K36me3 (so called ‘K4-K36 domains,’ as defined through “peak calling”—see below), RNA-seq reads at these domains can be extracted using programs such as bedtools³⁶ to estimate lncRNA abundance levels.

MicroRNA (miRNA) sequencing experiments, which must be run separately from standard RNA-seq experiments, investigate mature miRNAs that are typically 22 nucleotides in length. Since most high-throughput sequencing machines produce reads that are significantly longer than mature miRNAs, the short reads obtained from miR sequencing often contain of portion of adapter sequences. The first step of miRNA analysis is thus to remove adapter sequences using tools such as fastx (http://hannonlab.cshl.edu/fastx_toolkit/). As an alternative, the Cutadapt program³⁷ can be used for these purposes. Cutadapt is useful in that it supports a larger range of sequencing platforms, including color-space reads, which allow basepairs to be encoded in color to reduce sequencing error rates (https://www.biostars.org/p/43855/), in comparison to fastx. Following adapter removal, short reads can be mapped to the reference genome similar to that of data obtained from long RNA sequencing. miRBase³⁸ provides a comprehensive collection of miRNA sequences. Bedtools can then be used to extract read counts for each of the annotated miRs obtained from alignment. If novel prediction of miRNAs from a given sample is desired, the miRDeep³⁹ program can be used. Many additional families of small RNAs (Table 2) can be identified and analyzed with similar approaches^40-46,38,47.

Differential Analysis: Approaches, Advantages and Disadvantages

Differential analysis refers to the process of identifying differences in RNA expression levels of individual genes, or of individual splice variants of a single gene, between control and experimental samples. This is not straightforward, as there are numerous tools available, each of which is associated with high rates of false positive or false negative discovery and generates very different lists of regulated genes when applied to the same sequencing data. Generation of ideal differential analytical tools is therefore a focus of great interest in the field.

The first step in this process often involves summarizing gene- or gene variant-level read counts using a popular python program called htseq (http://www-huber.embl.de/users/anders/HTSeq/). According to individual needs, bedtools can also be used for gene count summarization. All read counts across genes and samples are then imported into a data matrix so that each row represents a gene, and each column represents a sample. This data matrix serves as the input for downstream differential analyses. Numerous differential analysis tools have been developed in recent years. Two popular choices are DESeq⁴⁸ and edgeR⁴⁹, with both methods based on negative binomial testing, which provide an exact test (generalization of the Poisson distribution model) that is ideal for modeling biological variances of read count data. Variance estimates are often problematic for RNA-seq datasets, as sample sizes from animal models are typically small with an N=3 (three biological replicates) used per condition in most experiments. Therefore, numerous statistical methods are used to exploit the relationship between mean- and variance-related information obtained from neighboring genes to stabilize variance estimations. DESeq uses local regression to model mean-variance relationships, while edgeR uses Bayesian methods to “borrow” information from neighboring genes. Both methods generate satisfactory results, and their lists are often consistent. Recently, a new method called voom⁵⁰ has been developed, which does not depend on negative binomial testing and instead models mean-variance relationships on log-transformed read counts, thereby assigning a precise ‘weight’ for each gene. These weights are then entered into the limma⁵¹ empirical Bayes analysis pipeline. This approach provides access to a large collection of analysis tools developed originally for microarrays, which makes “voom” a particularly attractive option. Furthermore, DESeq has been criticized for its extremely conservative outputs, and its authors have acknowledged that their method yields high false negative rates⁵². Recently, the statistical power of the program has been improved with the introduction of DESeq2 (http://www.bioconductor.org/packages/release/bioc/html/DESeq2.html). Nevertheless, experience with RNA-seq suggests that a larger number of biological replicates are needed to derive the most reliable data^53,54.

Cuffdiff, another popular tool for RNA-seq differential analysis, and part of the cufflinks pipeline, is a natural choice for investigators using cufflinks for initial identification of transcripts and splice variants. However, we, along with numerous other groups, have found cuffdiff to generate a high degree of false positives. The authors of cufflinks claim to have improved the reliability of detection in cuffdiff2⁵⁵, but a recent comparison pointed out that cuffdiff2 may be too conservative⁵⁶. Another reason for avoiding cuffdiff is its restrictive workflow. Cuffdiff only accepts short read alignments and performs read counts, GC content adjustments, 3’ bias corrections and differential analyses in one integrated workflow; based on a standard four-core workstation, it would take 1-2 days to complete analysis for a 20 sample (10 vs. 10) dataset and may crash for datasets with >80 samples. Using htseq for gene counts in connection with limma or DESeq greatly reduces the analysis time to <1 hour.

Performing all of the steps of RNA-seq analysis, including quality assessment, alignment, gene count summarization and differential analysis, can be tedious if investigators have hundreds/thousands of samples (e.g., from studies of clinical populations), as well as dozens of comparisons, to run and examine. Therefore, a new Python-based computational pipeline has been developed SPEctRA (scalable PipEline for RNA-seq Analysis) (https://github.com/shenlab-sinai/complete-seq/) to perform all these tasks with a single command accepting various parameter settings as a configuration file. This pipeline is designed to function within different computational environments. For example, when performing analyses under a portable batch system (PBS) cluster, computer software that allocates computational tasks, SPEctRA will generate batch scripts, automatically submit them and wait for them to finish. A new feature of this pipeline, which will allow it to run in Amazon cloud where both CPU power and memory can be configured on demand for each run, is currently being developed.

ChIP-seq: Initial Methods of Data Processing and Annotation

As with RNA-seq, the raw data obtained from ChIP-seq are tens to hundreds of million short reads per sample that correspond to genomic regions bound to the DNA-binding protein subjected to ChIP (e.g., a modified histone, transcription factor, etc.)¹⁹. Quality control and sequence alignment are the first steps in analyzing such ChIP-seq data. Any of several short read aligners, such as bowtie⁵⁷ or BWA⁵⁸, can be used to align reads to a reference genome and assess the enrichment of the DNA-binding protein of interest. Alignments are exported as BAM⁵⁹ files for further analysis. Sequencing quality can then be assessed using fastqc, as described above. If ChIP-seq samples are generated from sonication-based fragmentation methods, then duplicated alignments need to be removed, a procedure easily accomplished with samtools⁵⁹ or picard (http://picard.sourceforge.net/). In doing so, PCR duplicates are identified as reads aligned to the same genomic location or on the same strand of fragmented DNA. During sonication, genomic DNA is theoretically fragmented at random, thereby rendering it unlikely that two unique reads will align to the same location or on the same strand of DNA. Since micrococcal nuclease (MNase) digestion, as used in “native” ChIP, preferentially digests genomic sequences containing specific nucleotide sequences (i.e., the rate of MNase cleavage 5’ of A or T nucleotides occurs at a 30 times greater rate in comparison to cleavage at G or C nucleotides), the probability of obtaining fragments from the same location or strand is increased. Therefore, removal of these reads with native ChIP samples, unless excessively high (thresholds can be determined using DANPOS⁶⁰, as described below), is not typically performed. Removal of repeat reads can be problematic when analyzing repeat regions of the genome, which are now known to be important in biological regulation, thus framing a challenge for current bioinformatics innovations. The number of unique reads per sample is an important criterion for the quality of ChIP-seq data. Since antibodies play a key role in the success of any ChIP-seq experiment¹⁹, it is vitally important to determine the efficiency and specificity of the antibody being used. For this purpose, the phantomPeak⁶¹ tool can be employed, as suggested by the Encyclopedia of DNA Elements (ENCODE) consortium⁶², to examine the distribution of cross-correlations between represented strands of DNA to determine peak enrichments (see next section).

It is also helpful to visualize ChIP-seq enrichment profiles to ensure quality and diagnose any problems that may exist. Two approaches—localized inspection and global visualization—are used. For local inspection, the IGV genome browser⁶³, which runs on all platforms without uploading data, is popular. With the tools provided by the IGV, BAM file are first converted to TDF files and then loaded onto the browser to display coverage between two chromosomal coordinates. For global visualization, ngs.plot (https://code.google.com/p/ngsplot/)²⁷, which allows for inspection of both average and ‘laid out’ coverages as curves or heatmaps, respectively, at various functional genomic regions, such as TSS, TES, gene body, CpG islands, enhancers, exons and DNase I hypersensitive sites (DHSs). ngs.plot, which is simple to use and supports numerous genomes, can accommodate large alignment files with relatively small memory footprints. A ChIP-seq quality assessment pipeline has been developed (https://github.com/nyshao/chip-seq_preprocess) to perform all the above steps with a single command. This pipeline uses a similar design to the RNA-seq pipeline described previously. When dealing with large sample numbers and multiple comparisons, this pipeline saves significant amounts of time and effort.

Another important quality control consideration during the initial analysis/alignment of ChIP-seq data requires an in depth understanding of the differences between ‘uniquely mapped’ vs. ‘unique reads/tags,’ two terms that are often used interchangeably in the field but are distinctly defined. ‘Uniquely mapped’ reads are identified by all alignment software programs and exclude sequences that align to multiple genomic locations. These excluded sequences likely represent repetitive regions of the genome, or non-repetitive genomic loci that are extremely similar in sequence (although the latter becomes increasingly unlikely with greater read lengths). ‘Unique reads/tags,’ however, refer to reads remaining after PCR de-duplication (i.e., the non-redundant fraction or NRF) using tools such as samtools. PCR de-duplication essentially excludes reads (all copies but one) that align to the same genomic location. These non-unique reads are often PCR duplicates resulting from over-amplification during library processing, most likely due to low starting material or poor antibody efficiency. These reads are thus often removed to avoid PCR amplification bias. However, arguments now exist^64,65 to suggest that these duplicative reads may provide increased dynamic range to ChIP signals. For example, if one were to imagine two 100 bp regions of the genome (A and B), where factor X binds strongly to region A and weakly to region B, and if all duplicated reads mapping to these regions were removed (assuming a sequencing read length of 100 bp), the resulting effect would be two reads for each of these regions (one mapping to each strand of DNA), an output that would prevent the investigator from distinguishing the differential binding strengths of factor X to these regions. Therefore, although removing duplicates is a highly conservative measure to prevent PCR amplification bias, it is also likely that this process similarly removes real signals that might be informative for determination of binding interactions. Having said this, without excluding these duplicates, one runs the risk of generating high levels of false positive findings due to increased signal. In our experience, keeping two reads for any identified duplicates (instead of only one) dramatically increases the sensitivity of downstream data analysis and genome browser viewing with IGV, while retaining low false positive discovery rates.

Peak calling—identifying a genomic region that displays a significant level of a DNA-binding protein above background—is an important task in analyzing ChIP-seq data. Currently, dozens of peak calling tools exist, and these methods can generally be separated into two categories: sonication- and MNase-digestion based. For sonication-based peaks, Model-based Analysis of ChIP-Seq (MACS)⁶⁶ and Hypergeometric Optimization of Motif EnRichment (HOMER)⁶⁷ are used. In our experience, both methods work very well for punctuated peaks (e.g., H3K4me3), however, if peaks are broad and diffusive (e.g., total histone H3 or H3K9me2), detection can be challenging, and HOMER is generally recommended, as systematic evaluations of broad peak detection are often difficult. For MNase-digested peaks, DANPOS⁶⁰, which detects not only basal enrichment but also various nucleosomal events (e.g., nucleosomal positioning and occupancy and “fuzziness” between the two), can be used. After peaks are detected, it is useful to determine their location within the genome, such as genes, gene deserts or pericentromeres. A regional analysis tool, which is part of the diffReps⁶⁸ package—more recently extended as a standalone program (https://github.com/shenlabsinai/region_analysis)—has been developed to perform this task. This program features single command usage and can assign peaks to one of eight distinct categories including proximal promoter (+/− 250 bp of a TSS), promoter 1k (1 kb upstream of a TSS), promoter 3k, gene body, gene desert, pericentromere, subtelomere and other intergenic loci. Regional annotation is a very useful feature, and other alternatives exist to perform this function including the ChIPpeakAnno⁶⁹ package. After peaks are annotated, enrichment of specific chromatin marks is easily visualized for their genomic distribution using piecharts or plots.

Differential Analysis: Approaches, Advantages and Disadvantages

Differential analysis of ChIP-seq data aims to identify genomic loci or broader regions that display significant changes in enrichment between control and experimental conditions. Although it is natural for investigators to consider differential analyses as an extension of peak calling, in reality the former cannot be measured accurately using standard peak calling methods. Currently, one differential analysis tool⁷⁰ based on peak calling exists, however, the numerous challenges associated with this approach make it disadvantageous. Differences in antibody efficiency, sequencing biases and manual handling of samples (i.e., sequencing library preparation), can produce very different signal-to-noise ratios between samples, which confounds peak calling. Insufficient sequencing depth (i.e., too few reads per sample) can also result in specific genomic regions having different coverage across different samples and complicate differential analyses. For example, consider the following situation: a 3 kb peak in biological replicate #1 under condition A; two 1 kb peaks with 1 kb gap in biological replicate #2 under condition A; a 2.5 kb peak shifted to the left with lower enrichment in biological replicate 1 under condition B; and a 2 kb peak shifted to the right with higher enrichment in biological replicate 2 under condition B. Such heterogeneity in peak calling, which is common in analyses of brain (see Box 1) makes comparisons difficult.

Further exaggerating this problem is the fact that different peak calling methods, or the same peak calling method with different parameter settings, can yield very different peaks using the same ChIP-seq dataset. Biological replicates are thus essential in neuroepigenomics research to enhance statistical power and increase precision. However, most peak calling methods to date were not designed with biological replicates in mind. To address these challenges, one can use a sliding window-based strategy to ensure that the entire genome is scanned and scored continuously, and that significant regions can be extracted for further analyses. diffReps⁶⁸, which was developed to address this need, is a Perl-based program that features single command usage and has been applied to several projects. diffReps scans genomes with a predefined window size and performs one of four distinct statistical tests, identifies samples passing pre-defined cutoffs, merges replicates, performs multiple testing corrections and reports results. Using benchmark standards, diffReps is highly sensitive and efficiently controls for false positives (see Figure 2 for examples of ChIP-seq analysis outputs). It is also important to assess the reproducibility among biological replicates as an additional quality control measurement during data processing. To do so, programs such as corrgram⁷¹ can be used, which generates Pearson’s correlation coefficients between ChIP-seq signals derived from multiple signals. Alternatively, irreproducible discovery rates (IDRs) can be derived following peak calling to determine the number of ‘bound’ regions observed between replicates, as described⁷².

Important questions for ChIP-seq experiments are the read depth required to appropriately make assumptions based on aforementioned analyses, and to determine if experimental sequencing depths are sufficient for the questions being asked. In general, researchers can follow the guidelines established by the ENCODE consortium⁷², which indicate a minimum of 10 or 20 million uniquely mapped reads for factors displaying punctuated or broad peaks, respectively. One can also perform saturation analyses for each individual factor to assess the sufficiency of sequencing depth. In essence, such analyses repeatedly perform peak calling on a series of samples from the original ChIP-seq library using increasing sampling rates, followed by plotting the number of peaks assigned vs. the sampling rate to identify plateaus indicative of sequencing depth saturation. One caveat for such analyses, however, is that narrow peak binding factors, such as transcription factors, will typically display increased numbers of binding sites/peaks with increasing numbers of reads obtained. This is due to the fact that such factors rarely saturate at the number of reads that are practical for saturation analysis.

Another challenge in analyzing ChIP-seq data involves appropriately annotating chromatin states while examining combinations of chromatin modification patterns in a given sample. To achieve this, an automated computational system, referred to as ChromHMM⁷³, has been developed to integrate ChIP-seq information from multiple histone modifications, chromatin factors, transcription factors, etc. to accurately assess combinatorial and spatial patterns of marks/binding factors in biological samples. Such analyses provide a novel computational tool for ‘learning’ chromatin states, characterizing biological functions and achieving genome-wide visualizations of such annotations. Although numerous methods now exist for analyzing ChIP-seq data, independent biological validation remains an essential step to confirm the accuracy of any data derived from genome-wide approaches.

Overlaying ChIP-seq and RNA-seq Data

In an effort to understand the transcriptional and epigenomic mechanisms underlying RNA expression, a major goal of current research is to merge ChIP-seq and RNA-seq datasets. In fact, if one analyzes a sufficiently large number of epigenomic endpoints (numerous histone modifications, chromatin remodeling factors, transcription factors and other chromatin regulatory proteins) it should in theory be possible to identify “chromatin signatures” that reflect specific modes of transcriptional regulation, such as gene activation or repression as well as gene priming or desensitization.

Accomplishing this task remains challenging given the vast amounts of data (terabytes) involved. One initial approach is to first reduce observed ChIP-seq events to individual genes. This thereby reduces the problem to essentially comparing two gene lists—one from differential analysis of ChIP-seq data and the other from differential analysis of RNA-seq data. GeneOverlap (http://www.bioconductor.org/packages/devel/bioc/html/GeneOverlap.html), a bioconductor package for testing and visualizing gene overlaps, can be used for this purpose. However, such analyses oversimplify the problem. First, a given histone modification can show opposite changes across the span of a given gene. Also, all histone modifications regulate gene expression in a highly cooperative but complex fashion. For example, the histone modification, H3K4me3, is highly enriched at active gene promoters and is often used as an indication of transcriptional initiation. However, H3K4me3 levels at some genes can increase, in concert with other DNA-binding proteins (e.g., RNA polymerase II [Pol 2] serine 5 phosphorylation), in response to a specific stimulus without increasing the expression of its cognate gene, a phenomenon thought to indicate gene “poising.” And in a much smaller number of examples, H3K4me3 has been linked to gene repression when it associates with ING2 (inhibitor of growth family-2), a chromatin regulatory protein⁷⁴. On the other hand, the histone mark, H3K9me3, which typically enriches within the gene bodies of silenced genes along with numerous other corepressor proteins, may exhibit increased enrichment with distinct binding partners to promote alternative splicing⁷⁵.

Previously, numerous groups^76,77 have attempted to use histone modifications to predict gene expression at baseline states in cultured cells and have achieved good performance. These groups were able to generate predicted expression levels correlating with real expression levels with coefficients as high as 0.8. Their approach worked by separating gene bodies into dozens of bins, where the relative enrichment of histone marks was determined. In this case, enrichment values serve as observed data and corresponding gene expression levels serve as target variables. Such training data were then fed into machine-based learning methods, such as support vector machines or generalized linear models, to learn how the behavior of histone modifications relates to gene expression. This approach, however, has not yet been applied to brain tissue at rest and, in particular, to data examining the relationship between stimulus-induced changes in epigenomic endpoints and gene expression, perhaps the most essential question facing molecular neuroscientists interested in examining the role of epigenomic mechanisms in mediating the impact of environmental manipulations on the brain.

We recently developed a bioinformatics approach that was successfully applied to studying the actions of repeated cocaine exposure in the mouse nucleus accumbens¹⁶, a key brain reward region important in addiction. We first correlated ChIP-seq data for six histone modifications and for total Pol 2 binding with changes in gene expression as determined by RNA-seq. While we were able to demonstrate >50 distinct chromatin signatures that correlate with altered RNA expression, the correlations were relatively weak and not deterministic. This suggests that a far larger number of chromatin endpoints are needed to make such predictions; indeed, current research is revealing large numbers of previously unappreciated histone modifications (see below). Additional considerations unique to the analysis of brain tissue are also likely involved (see Box 1).

We used a similar approach to determine whether repeated cocaine-induced changes in alternative splicing in nucleus accumbens correlate with changes in chromatin modifications⁵⁰. We focused on genomic regions most closely associated with splicing, such as variant exons, alternative donors, alternative acceptors and their neighboring intronic regions, and examined numerous chromatin modifications at these regions. We applied this information extraction procedure to all known transcripts of the mouse genome and calculated the log2 fold changes between repeated cocaine and saline and loaded the values into a data matrix so that each row represents a transcript and each column represents a histone mark-gene region combination (e.g., a H3K4me3 peak at the proximal promoter of that gene). We also used cufflinks to determine the expression changes of these transcripts. K-means clustering was applied on the chromatin modification matrix to identify clusters of transcripts that show similar patterns. These clusters were then correlated with transcripts’ expression changes to determine whether there is any significant overlap. We were able to identify 29 such clusters or chromatin modification signatures. Further analysis of these signatures using motif finding revealed several important splicing factors and transcription factors deduced to be important in their regulation. One of the splicing factors, A2BP1 (Rbfox1/Fox-1), was validated by showing that knocking out A2BP1 in adult nucleus accumbens blocks cocaine’s regulation of several of the genes identified in this genome-wide analysis. Moreover, A2BP1 knockout was shown to attenuate behavioral responses to cocaine¹⁶. These findings illustrate how the overlay of ChIP-seq and RNA-seq data can reveal fundamentally new insight into the biological basis of a complex brain disorder.

Overlaying DNA Methylation with Gene Expression Analyses

Methylation of cytosine bases in DNA (5-methylcytosine or 5mC) has historically been viewed as an important mode of gene repression. However, many alternative forms of DNA methylation have been demonstrated in recent years, most importantly, 5-hydroxycytosine (5hmC), which is enriched in brain and associated with gene activation⁷⁸. Despite the implication of DNA methylation in numerous neuropsychiatric phenomena, virtually all studies to date have focused on individual candidate genes, with genome-wide explorations of brain sorely lacking.

Several methods, described in the companion review¹⁹, are used to obtain a genome-wide map of 5mc and 5hmc, but doing so at single nucleotide resolution remains a challenge due to the sequencing costs involved. These methods include genome-wide bisulfite sequencing (BS-seq)—including oxidized BS-seq (oxBS-seq) to distinguish 5mc and 5hmc, RRBS (reduced representation bisulfite sequencing) which focuses on CG rich regions, and meDIP-seq which involves immunoprecipitation of genomic DNA fragments with an antibody directed against 5mC or 5hmC followed by deep sequencing. The human “methylome” can be obtained from a chip-based method, however, available chips do not distinguish between 5mC and 5hmC.

The first step in analyzing BS-seq data is to trim the adaptors and low quality 3’ ends. Trim Galore! (http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/), which wraps around Cutadapt and fastqc to provide extra functionality for RRBS libraries, can be used for this purpose; it reports the sequencing QC by calling fastqc after trimming. Then the reads are aligned to the reference genome. Generally, the aligners need a preprocessing step to tolerate the conversion of C to T. Two strategies of aligners are available. One is a wild-card aligner, such as BSMAP (https://code.google.com/p/bsmap/)⁷⁹, as all Cs in a reference genome are replaced by a wild-card letter (e.g., “Y”), with both C and T possibly aligned to Y. Another approach is a three-letter aligner, such as Bismark (http://www.bioinformatics.babraham.ac.uk/projects/bismark/)⁸⁰; here, all Cs in reads and a reference genome are converted to Ts, whereas Gs are converted to As equivalent to C-to-A conversions on the reverse strand, and the aligners process alignments in only three-letter alphabets (A, G and T). Wild-card aligners may reach better genomic coverage while potentially introducing bias to increase DNA methylation levels, compared with three-letter aligners⁸¹. After the alignment, the quantification of absolute DNA methylation can be achieved by counting the alignments, referred to as “methylation calling”.

The next step is to detect sites of differential methylated regions (DMRs) in experimental vs. control conditions. Basic analyses like t-tests or Wilcoxon rank tests can be applied; details about more advanced methods have been reviewed recently⁸¹. oxbs-sequencing-qc (https://code.google.com/p/oxbs-sequencing-qc) is a pipeline based on Bismark, including trimming of adaptors and low quality, alignment and methylation calling. Users need to implement the detection of DMRs by other analytical tools. MOABS (http://code.google.com/p/moabs/)⁸² is another pipeline, developed by the authors of BSMAP. It calls and accepts the alignment results from BSMAP or Bismark, and implements the detection of DMRs based on a Beta-Binomial hierarchical model.

As the final step of data analysis, DMR events detected by BS-seq or oxBS-seq are annotated by region_analysis or ChIPpeakAnno, and then integrated with differential expressed genes detected by RNA-seq, or with ChIP-seq characterizations of other epigenomic endpoints, as described in the ChIP-seq section above. Approaches similar to ChIP-seq are used to analyze sequencing data obtained from meDIP-seq experiments. Reads are first evaluated for QC and then aligned to a reference genome using the various tools outlined above, with the number of aligned reads used for methylation calling⁸³. In addition to canonical CpG methylation, highly conserved non-CpG methylation (mCH, where H is either A, T or C) may also play an important role in the central nervous system. It was recently reported that non-CpG methylation accumulates in neurons, but not glia, of the cerebral cortex during development¹³.

Mass Spectrometry in Epigenomics Research

Mass spectrometry (MS) has become a powerful tool for the analyses of DNA methylation and histone modifications¹⁰⁰ with some newer applications to small non-coding RNAs¹⁰¹. However, much of this has been at the global chromatin level and not at specific genes¹⁰². Identification and quantification of histone PTMs from a variety of cells and tissues have been the focus of several studies, and typically involve examining the protein by Bottom Up, Middle Down or Top Down MS¹⁰².

Bottom Up MS refers to digestion of the protein into small (<3 kDa) size peptides, and analysis of those peptide typically by nanoLC-MS/MS. These experiments are fairly high-throughput and easier to perform than Middle or Top Down MS. The Bottom Up MS approach has been coupled with many diverse platforms and quantification schemes (e.g., stable isotope labeling approaches)¹⁰³ to compare the histone PTM profiles from various cellular sources and states. This Bottom Up approach has been applied to the global analysis of core histone PTMs and variants from whole brain of adult C57BL/6 mice¹⁰⁴. This large-scale proteomics analysis identified more than 10,000 peptides creating a dataset that containing 146 modification sites on 1475 peptides in various combinatorial patterns on short peptides. Among these histone peptides, 58 novel sites of modification were discovered. Bottom Up MS in fact has led the way to identify novel marks on histones, such as crotonylation, succinylation, malonylation and 2-hydroxyisobutyrylation of lysine residues^105-107 and serine/threonine O-acetylation¹⁰⁸.

Middle Down MS involves interrogation of polypeptides over the 3 kDa molecular weight range. These MS experiments allow for long-range spanning PTMs to be identified on the same peptide to determine a long distance combinatorial code, and has mainly been focused on the N-terminal tails of the histone proteins. Currently, it has been determined by this approach that there exists >200 combinatorially modified N-terminal forms of histone H3, with about half as many N-terminal forms existing on histone H4¹⁰⁹. These numbers are far lower than those that are theoretically possible, seemingly indicating that much specificity exists for creating multiply modified histone forms. Kelleher and co-workers used this approach to identify histone PTM patterns and variants from whole brain, cerebral cortex, cerebellum and hypothalamus from 7-8 week old Sprague–Dawley rats¹¹⁰. Interesting patterns of modifications were found, including more silencing modifications such as H3K9me3 in the cerebellum.

Analysis of intact histone proteins (Top Down MS) is a challenging task, but allows for the characterization of combinations of modifications across the entire protein sequence. Top Down MS has been used to identify over 700 unique histone forms across all the core histones¹¹¹ but still remains the most technically challenging MS experiment to perform.

Data Processing for Histone PTM Analyses

As can be inferred from the above overview, the data analysis of histones is very complicated due to the many different types of modifications that can be detected, the distinct combinations of these PTMs that can be found and the low level stoichiometry of many histone PTMs, which can often give rise to false positive assignments. Currently, there are numerous computational approaches for both the qualitative and quantitative analysis of histone proteomics datasets. The general workflow for proteomics data is that peptides are resolved by liquid chromatography, typically electrosprayed ionized into the mass spectrometer and full MS spectra of the precursor ions collected. The precursor peptide ions are then selected and fragmented to obtain tandem mass spectra (MS/MS spectra) by various approaches including high energy C-trap dissociation (HCD), collision induced dissociation (CID) or electron transfer dissociation (ETD). Peptides and associated modifications are identified by searching the MS/MS spectra against an organism specific database.

As histone PTM patterns are challenging, several in-house developed algorithms to deal with these datasets have been developed. (e.g., MILP or PTMap)^112-114. However, several other database searching engines are commercially available for analyzing histone proteomics data with the most common being Mascot¹¹⁵ or SEQUEST¹¹⁶ but with many others emerging with different specificities or advantages^117-119. Searching for a large number of PTMs in one search drastically increases the number of candidate peptides, slowing computational searching speed, increasing false positive IDs and even identifying less spectra. Therefore, it has been shown with histone data to best search for classes of modifications one at a time to reduce false positive misassignments. Additionally and equally important to these analyses for histone PTM proteomics datasets is the use of a site localization scoring algorithm¹²⁰. As histone PTMs occur on peptides with multiple acceptor amino acids (e.g., several lysine residues within a single peptide), it is necessary to either manually confirm the fragment ions or use a localization scoring algorithm to confidently determine on which residue the modification site occurs.

Combining Proteomic with Genomic and Epigenomic Analyses

Most histone proteomics data have been generated from global chromatin extractions, and hence the genomic loci involved are usually lost. The strength in being able to combine genomics-epigenomics level information with proteomics data is that an unbiased assessment of histone PTMs could be performed and mapped back to genomic loci. Several approaches based on tagging of proteins or using antibodies to isolate nucleosomes have been performed^121-123 and have resulted in a proteomic portrait of the genomic landscape. For example, histone PTMs and other proteins associated with MSL (male-specific lethal) associated proteins in Drosophila S2 cells were recently mapped by quantitative proteomics¹²³. A ChIP-seq and proteomics study also found that Brd (bromodomain) proteins are found on active Hox genes in cultured HEK cells, with distinct patterns of histone H4 acetylation¹²⁴. However, both of these studies map histone PTMs, not to a specific locus, but rather over a distribution of enriched genes. Experiments that isolate DNA specific sequences have still greater power to establish protein profiles and histone PTM patterns at a distinct gene. One of the first attempts at this was by Kingston and co-workers using PiCh (Proteomics of isolated Chromatin segments) to isolate telomeric repeats¹²⁵. Other approaches such as ChAP–MS and iChIP have also all been shown to be very useful for isolation of genomic material for mass spectrometry interrogation^126,127, however, it remains to be seen if they will be highly adopted.