ngs.plot: Quick mining and visualization of next-generation sequencing data by integrating genomic databases

ngs.plot workflow and algorithms

The workflow of ngs.plot is depicted in Figure 1. Initially, ngs.plot searches through its database to find the genomic coordinates for the desired regions and uses them to query the alignment files of an NGS dataset. It then calculates the coverage vectors for each query region based on the retrieved alignments. It finally performs normalization and transformation on the coverage and generates two plots. One plot is an average profile that is generated from the mean of all regions. This plot provides the overall pattern at the regions of interest. The other plot is a heatmap that shows the enrichment of each region across the genome using color gradients. The heatmap can provide three-dimensional details (enrichment, region, and position) of the NGS samples under study.

A user can specify the plotting regions using a genome name, such as “mm9” and a region name, such as “genebody”. Further options are provided to choose a particular type of region, otherwise the default is used. For example, exons are classified into “canonical” (default), “variant”, “promoter”, etc.; CGIs are classified into “ProximalPromoter” (default), “Promoter1k”, “Promoter3k”, etc.; gene lists can be provided to create subsets of the regions. For convenience, we have provided the gene names/IDs in both RefSeq [23] and Ensembl [24] format. To be more flexible, a user can also use a BED (https://genome.ucsc.edu/FAQ/FAQformat.html#format1) file for custom regions. A BED file is a simple TAB-delimited text file that is often used to describe genomic regions. This is particularly useful if a user performed peak calling for a transcription factor and would like to know what is happening at or around the peaks.

The alignment files must be in BAM [25] format, which is now used widely for short read alignments. ngs.plot conforms to the SAM specification [25] of BAM files and can work with any short read aligner. A BAM file is compressed and indexed for efficient retrieval. In ngs.plot, the “physical coverage” instead of the “read coverage” is calculated for both ChIP-seq and RNA-seq. This is achieved by extending each alignment to the expected DNA fragment length according to user input. The coverage data are then subjected to two steps of normalization. In the first step, the coverage vectors are normalized to be equal length and this can be achieved through two algorithms. The default algorithm is spline fit where a cubic spline is fit through all data points and values are taken at equal intervals. The alternative algorithm is binning where the coverage vector is separated into equal intervals and the average value for each interval is calculated. This first step of length normalization allows regions of variable sizes to be equalized and is particularly useful for genebody, CGI, and custom regions. In the second step, the vectors are normalized against the corresponding library size – i.e., the total read count (only the reads that pass quality filters are counted) for an NGS sample to generate the so called Reads Per Million mapped reads (RPM) values. The RPM values allow two NGS samples to be compared regardless of differences in sequencing depth.

We have implemented many functions to manipulate the visual outputs of an ngs.plot run, as follows:

RNA-seq mode

ngs.plot can accurately calculate coverage for RNA-seq (Figure 2A). RNA-seq experiments are unique because the short reads are derived from messenger RNAs and other expressed RNAs, many of which result from exon splicing. The ngs.plot database contains the exon coordinates for each transcript so that the coverage vectors for exons are concatenated to simulate RNA splicing in silico.

Coverage extraction

Coverage extraction is at the core of the ngs.plot workflow. This process often consumes a lot of computational resources because of the large size of genomes (e.g., the human genome has approximately 3 billion nucleotides) and because alignment files are also very large (on the order of tens of GB). In the history of ngs.plot, we first used a strategy called “run-length encoding” (RLE) to represent genomic coverage vectors. RLE uses a very simple approach so that consecutive and repetitive values are represented by the value and number of repeats. For example,

Genomic annotation databases

We developed a genome crawler that fetches various genomic annotations from public databases, and processes and saves them into R binary tables (Figure 2C). R binary tables are very easy to create and their columns are indexed by R internally. This helps to avoid setting up local databases, which turns out to be a convenience for users. Currently, we considered Ensembl [24], UCSC [29], and ENCODE [11] [see Additional file 1: Table S2], and will incorporate more public databases in the future. Ensembl and UCSC provide classic genomic features such as genes, transcripts, exons, and CGIs, while ENCODE provides more recent epigenomic features such as enhancers and DHSs. Because these databases host genomic information at different servers that are setup by separate groups of people, there is no uniformity in constructing the URL for a specific genome. Sometimes, a large database (such as Ensembl) may store different classes of species, such as animal and plant, using slightly different naming schemes. To address this issue, we used JSON format to manually create configuration files for each naming scheme so that an automated pipeline can pull data from different sources. New naming schemes can be handled by simply adding JSON configuration files. The files that are downloaded by the genome crawler include Gene Transfer Format (GTF), Gene Prediction (GP), BED, and MySQL database inquiries, each of which is processed by a separate program module. The RefSeq annotations downloaded from UCSC are in GP format, which can be converted into GTF files using the “genePredToGtf” utility from UCSC. The GTF files are parsed by custom scripts to generate uniformly formatted text files that are further converted into R binary tables. The gene annotations are used to derive gene deserts. Locations about heterochromatic regions such as centromeres and telomeres are downloaded from UCSC and are used to derive pericentromeres and subtelomeres. All the gene annotations, gene deserts, pericentromeres and subtelomeres are used to build a genome package for the “region analysis” utility (https://github.com/shenlab-sinai/region_analysis) on the fly, which is used to perform location-based classifications on CGIs and DHSs. In total, more than 60 million functional elements have been incorporated into ngs.plot’s database so far (Table 1). Additional genomes can be added at any time as needed. The functional elements for each genome are packed into a compressed archive file that can be installed on demand by a user. A Python script (named ngsplotdb.py) is provided to manage the locally installed genomes. In the following, we describe each type of functional element and how they are processed.

Galaxy integration

ngs.plot command interface features simple and easy usage. This allows users to blend ngs.plot with other bioinformatic and Unix tools seamlessly. However, the command interface may be intimidating to wet lab biologists. Therefore, we developed a plug-in so that ngs.plot can be integrated into Galaxy [22] – a very popular web-based bioinformatics platform, which allows users to build their own point-and-click workflows using various tools. The plug-in features an easy-to-use graphic interface that can typically generate a figure in 3-4 steps. We have created a wiki-page to demonstrate such an example: https://code.google.com/p/ngsplot/wiki/webngsplot. Currently, this plug-in requires a locally installed Galaxy instance and is not available on the main Galaxy server.

Website and community involvement

ngs.plot’s hosting website provides manuals, source code, installation files, and links to many other resources. The source code is tracked by Google’s git server and is open for public contributions. To facilitate users in using ngs.plot, we have created nine wiki-pages so far and will keep adding new ones. Issue tracking is used for users and developers to report bugs and make suggestions. As this manuscript is being written, users from all over the world have downloaded ngs.plot for hundreds of times. We have also created an online discussion group for users to ask questions and help one another. So far, there are 51 active members who have contributed to 69 topics. We also use this opportunity to collect opinions from users so that we can improve the program further.

NGS data processing

The NGS data used in this manuscript were obtained from the Sequence Read Archive (SRA, http://www.ncbi.nlm.nih.gov/Traces/sra). The accession numbers and references of the datasets are listed in Table S1 [see Additional file 1]. ChIP-seq data were aligned to the reference genome by Bowtie [39]. Peak calling was accomplished by use of MACS [40] using default parameters. RNA-seq data were analyzed by the Tuxedo Suite [41]. The differential chromatin modification sites were detected by diffReps [42] using default parameters and the FDR cutoff was set as 0.1.

Benchmarking the performance of ngs.plot

To benchmark different coverage extraction methods, we used a ChIP-seq dataset that we previously published [43]. H3K9me2 is a histone mark that displays dispersive enrichment patterns and is often associated with gene silencing [44]. The ChIP-seq samples were derived from a mouse brain region (nucleus accumbens) where two biological conditions were assessed: chronic morphine and chronic saline administration. For each condition, three biological replicates were analyzed. We merged and sorted the alignment BED files for the three biological replicates under saline conditions and used BEDTools [45] to create a large BAM file that contains nearly 250 million alignments. From this file, 10, 20, 40, 80, and 160 million alignments were randomly sampled to create a series of BAM files that increase in alignment size exponentially. Different methods were used to extract coverage vectors for the TSS ± 5 Kb regions of all protein coding genes (~20,000). A number of metrics such as run time, memory usage, and file size were measured for different alignment sizes. All tests were performed on a Linux workstation with two 2.4 GHz CPU cores and sufficient memory.

At first, coverage needs to be pre-calculated for Tabix, bigwig, and RLE. This takes a long time to complete and the run time is strongly associated with the alignment size (Figure 3A). It takes samtools around 1,000 s to calculate the coverage for a 10 million read BAM file and more than 5,000 s for a 160 million read BAM file. RLE is much faster but involves a more rapid increase in time than samtools: it takes 80 s for a 10 million read BAM file and more than 800 s for a 160 million read BAM file. This is because RLE tries to load all alignments into memory and then performs calculations in a batch while samtools does the calculations by reading alignments in a stream. After coverage calculations, Tabix and bigWig also require the coverage files to be indexed. The indexing is more than 10 times faster than coverage calculation and shows strong association with the alignment file size (Figure 3A). Tabix is faster than bigWig: this is most likely because bigWig uses more than one index for different zoom levels [28].

Memory usage is a big problem for RLE. Even for the 10 million read BAM file, it uses 6 GB to finish the run, while for the 160 million read BAM file, it uses 75 GB (Figure 3B). In contrast, the memory footprint for Tabix indexing is very small: it uses ~50-60 MB for all BAM files. bigWig uses more memory for indexing than Tabix but is still reasonably small: at 160 million alignments, it uses 2 GB to finish the run (Figure 3B).

File size is another important metric. Both Tabix and bigWig create large coverage files that strongly associate with alignment file size (Figure 3C): at 160 million alignments, the Tabix coverage file is 1.2 GB while the bigWig coverage file is 1GB. As a comparison, RLE files are three times smaller: at 160 million alignments, the RLE file is 330 MB. Both Tabix and BAM have very small index file sizes. For BAM, the index remains around 6 MB for all alignment sizes while Tabix index is three times smaller. For bigWig, the index is an integral part of the format and its size is unknown to us.

By grouping regions into chunks we can save resource in index loading. This strategy worked well in our tests (Figure 3D). Based on a 10 million alignment file, it took the BAM method 1,200 s to load all TSS ± 5 Kb regions into memory for chunk size of 1. For a chunk size of 10, the time was reduced to less than 200 s – a six fold reduction. The time was further reduced to 88 s for a chunk size of 100. The other two methods – Tabix and bigWig – enjoyed similar degrees of time reduction by use of region grouping. It should be noted that Tabix used much less time than BAM at small chunk sizes. This is expected since the Tabix index is much smaller than the BAM index (Figure 3C). bigWig used the longest time among the three methods at small chunk sizes (1-100), suggesting its index is larger than the other two.

In our tests, Tabix was implemented with the Rsamtools [46] package and bigWig was implemented with the rtracklayer [47] package. Note that Tabix is a generic index program for text entries. After the texts are loaded into memory, they must be converted into binary representation of numerical numbers. The rtracklayer package, however, will unfortunately merge and sort the query regions before coverage vectors are retrieved. This means that the loaded coverage vectors are mixed and must be distinguished between the query regions for them to be useful for our purposes. All of the above operations require a significant amount of computational resources. At chunk size of 100, it took bigWig >3,500 s and Tabix >2,500 s to finish the operations (Figure 3E). In comparison, it took BAM only 64 s to calculate the coverage vectors on-the-fly. In the end, we abandoned support for Tabix and bigWig for this reason. A future goal of the field is to re-write the extraction functions in Rsamtools and rtracklayer extensively in order to optimize the retrieval time. Once that is done, we can add support for these two file formats.

Finally, we tested the entire process of coverage extraction and calculation for both BAM and RLE for different alignment sizes with regard to time and memory usage. The BAM method was tested with chunk size of 100 that is the default value for ngs.plot. BAM functioned superiorly compared to RLE on both metrics at all alignment sizes (Figure 3F). BAM’s run time only slightly increased from 143 s to 165 s for 10 and 160 million alignments; and its memory usage remained stable: less than 1 GB for all alignment sizes. In contrast, RLE used 4.1 GB RAM at 10 million alignments and increased to 15.9 GB RAM at 160 million alignments. RLE’s run time was also significant: at 160 million alignments, it took >1,200 s to finish – seven times longer than BAM.

In summary, the BAM strategy we chose in ngs.plot is a versatile, low profile approach that works robustly even with very large alignment files. This approach was introduced in ngs.plot v1.64 and has remained the approach of choice ever since.

Analysis of Tet1 and 5hmC ChIP-seq data in the differentiation of P19.6 cells

An easy-to-use and flexible visualization method of NGS data is necessary for computational biologists to formulate and validate hypotheses quickly. To demonstrate the power of ngs.plot, we used ChIP-seq data [see Additional file 1: Table S1] to study the relationship between Tet1 (ten eleven translocation protein-1), a methylcytosine dioxygenase, and 5-hydroxymethycytosine (5hmC) in the differentiation of mouse embryonal carcinoma P19.6 cells. P19.6 cells can be differentiated into neurons or glia by exposure to retinoic acid (RA) [48], and are widely used in research on stem cell differentiation. Tet family proteins play important roles in the conversion of 5-methylcytosine (5mC) into 5hmC, 5-formylcytosine (5fC), and 5-carboxymethylcytosine (5caC) in DNA [49–51], and are important regulators in the maintenance and differentiation of embryonic stem cells (ESC) [52, 53]. 5fC and 5caC are present low abundance in mammalian genomes [54] and are difficult to be detected by ChIP-seq. Therefore, we focus on 5hmC in this study. 5hmC is known to be enriched at TSSs, exons, CGIs, and enhancers [55–57]. The distributions of Tet proteins and 5hmC across the genome roughly overlap, while Tet1 and Tet3 prefer CpG enriched regions [50, 58]. This preference is at least partially due to their CXXC domains [58].

First, we used ngs.plot to investigate the enrichment profiles of Tet1 and 5hmC in P19.6 cells at different genomic regions, including genebodies, CGIs, exons, and enhancers (Figure 4A & Additional file 1: Figure S2). As the genebody plot (Figure 4A) shows, Tet1 is most enriched at TSSs but generally depleted at genebodies. CGI plots (Figure 4A & Additional file 1: Figure S2A) indicate that Tet1 is enriched at all kinds of CGIs at similar levels (~0.5-0.9 RPM) and demonstrates a clear drop of enrichment at flanking regions (±3 Kb). This suggests that the CXXC domain of the Tet1 protein highly prefers CpG abundant regions. In addition, Tet1 shows some enrichment at exons as well as enhancers (Figure 4A & Additional file 1: S2B) but the enrichment levels are weaker than that of CGIs, with enhancers being the weakest. As we expected, the enrichment patterns of 5hmC are highly similar to those of Tet1 (Figure 4A & Additional file 1: Figure S2), indicating concordance between the two marks. All of the above plots can be generated by ngs.plot with only one command for each. The user only needs to input into ngs.plot which regions and samples to examine and the size of the flanking regions.

5hmC plays an important role in stem cell differentiation, where its conversion from 5mC is mediated by Tet1 [49, 51, 55, 59]. The activities of enhancers are known to be specific to differentiated cell types and are often marked by the dynamics of 5hmC [60]. Here we illustrate the role that Tet1 plays in the conversion of 5hmC by studying the differential sites of Tet1 between control and RA-induced P19.6 cells (Figure 4B & Additional file 1: Figure S2). diffReps is a powerful program to detect differential chromatin modification sites using ChIP-seq data [42]. We used diffReps to find 7,735 (Increased: 3,762, Decreased: 3,973) Tet1 differential sites in total. To restrict the analysis to enhancers, we used H3K27ac as a mark for active enhancers [33]. Peak calling using MACS was performed in both control and RA-induced P19.6 cells and the two peak lists were combined to obtain 135,280 H3K27ac enriched sites (excluding the TSS ± 3 Kb regions). The peak list was used to filter the Tet1 differential sites that are not in enhancer regions. After filtering, we obtained 507 increased and 1,875 decreased enhancer-specific Tet1 sites induced by RA, whose genomic coordinates are then converted into two separate BED files. ngs.plot was applied on each BED file to plot the enrichment of both Tet1 and 5hmC (Figure 4B). It can be seen clearly that the trends of 5hmC dynamics follow those of Tet1 dynamics, with an overall consistency ratio of 82% (Tet1 increased sites: 74%; decreased sites: 84%). Their log fold changes are also weakly correlated (Pearson’s r = 0.46, Spearman’s ρ = 0.32, both with P < 2.2E-16). This is a vivid example illustrating how different computational tools can be used to identify biologically meaningful genomic regions and then feed them into the ngs.plot program for visualization.

Integrative analysis of poised and active promoters in ESC

Integrative analysis using genomic sequence information and multiple NGS samples is essential to investigate gene transcription and epigenomic regulation. ngs.plot’s ability to graph both ChIP-seq and RNA-seq samples allows a user to quickly establish correlations between different epigenomic marks and associated gene expression levels. Here, we demonstrate this feature of ngs.plot by use of multiple ChIP-seq samples, including several histone marks (H3K4me3, H3K27ac, and H3K27me3) and transcription factors (Suz12, Oct4, and Tet1), and an RNA-seq sample, from mouse ESCs (mESCs) [see Additional file 1: Table S1]. H3K4me3 is a promoter-enriched histone mark that is generally associated with transcriptional activation [33]. H3K27ac is an activation mark that locates at both promoters and enhancers [33]. The enrichment of both H3K4me3 and H3K27ac provides a signature of CpG-related promoters [61]. H3K27me3 is catalyzed by the Polycomb group proteins and is implicated in the silencing of genes [62]. The enrichment of both H3K4me3 and H3K27me3 marks the so-called “bivalent” domains that are prevalent in ESCs. They maintain the silencing or low expression of many genes in ESCs, which are poised for activation in differentiated cell types [33, 63]. Suz12 is a subunit of the Polycomb repressive complex 2 (PRC2) – a transcriptional repressor that catalyzes H3K27me3 [64]. Oct4, also known as POU5F1, is a critical transcription factor in the self-renewal of ESCs [65].

We divided all promoters (TSS ± 3 Kb) of the coding regions of genes into two groups, namely, Polycomb-targeted (PT, n = 5,132) and non-Polycomb-targeted (nPT, n = 19,013), based on the presence or absence of H3K27me3 peaks. To reveal the relationship between genomic sequences and epigenomic regulation, we sorted all promoters within each group based on their CG di-nucleotide percentages (CGP) and entered the gene lists into ngs.plot’s configuration files. We ran ngs.plot with its ranking algorithm set to “none” so that it used the input order. We also used a DNA input sample to pair with each epigenomic mark so that ngs.plot’s bam-pair functionality plots log fold changes. The use of the input sample is to counteract various biases introduced in ChIP-seq experiments [66]. All of the ChIP-seq samples within each group were then plotted with one command by use of the configuration file (Figure 5). We also plotted the RNA-seq sample using the same gene list with another command using the “RNA-seq” mode (Figure 5).

Figure 5 shows that the PT group has lower gene expression levels than the nPT group, indicating that genes containing the H3K27me3 mark are suppressed. Conversely, the activation mark H3K27ac shows lower enrichment in the PT group. As previously reported [67], H3K27ac is mutually exclusive with H3K27me3. However, another activation mark, H3K4me3, appears to be enriched in both groups. H3K4me3’s enrichment in the PT group demonstrates the prevalent existence of the “bivalent” domain in mESCs. The heatmaps of Figure 5 also indicate that there are strong correlations between certain epigenomic marks as well as with gene expression. To quantitatively measure these correlations, we used the plotCorrGram.r script included in the ngs.plot package to calculate and visually demonstrate all pairwise correlations between the samples. All the correlation coefficients and p-values are presented in [Additional file 2]. The corrgram is presented in [Additional file 3]. As expected, H3K27me3 and Suz12 show very strong correlations in both groups (PT: r = 0.89, nPT: r = 0.94, both with P = 0). Interestingly, CGPs show a moderate correlation with gene expression in the nPT group (r = 0.69, ρ = 0.73, both with P = 0), but this correlation is significantly decreased in the PT group (r = 0.20, ρ = 0.24, both with P = 0). As we mentioned above, Tet1 has a preference for CG-rich regions due to its CXXC domain. A moderate correlation is observed between Tet1 and CGPs in the PT group (r = 0.53, ρ = 0.54, both with P = 0), while a weak correlation is observed in the nPT group (r = 0.30, ρ = 0.37, both with P = 0). Tet1 also shows a moderate correlation with Oct4 in both groups (PT: r = 0.59, 0.57, both with P = 0; nPT: r = 0.76, ρ = 0.71, both with P = 0). It has been reported that Tet1 can replace the role of Oct4 in inducible pluripotent stem cell (iPSC) reprogramming, a process that is implicated in the regulatory circuit of ESCs [68]. This example demonstrates the ngs.plot’s capability to quickly correlate multiple epigenomic marks with other genomic features and with gene expression and creates figures that are publication-ready. A user can use these figures to gain biological insights into their NGS data and even generate novel hypotheses.

Examination of RNA-seq 3’ bias

The RNA-seq mode of ngs.plot can perform exon splicing in silico and this functionality can be exploited during RNA-seq quality control. For instance, in studies of human postmortem brain tissue, a major problem is that the RNA samples are often severely and variably degraded, as measured by the RNA integrity number (RIN) [69]. An RNA sample with low RIN is often associated with strong 3’ bias, which can impair the ability to otherwise assess the sample’s mRNA quantity. To demonstrate this, we analyzed an in-house RNA-seq dataset (unpublished) from human postmortem brain tissue obtained from two individuals with schizophrenia: one sample has an acceptable RIN (=7.8) and the other sample has a very low RIN (=3). The figure (Figure 6) generated by ngs.plot shows that the sample with low RIN is clearly biased towards 3’ in comparison to the sample with high RIN. A plot like this provides a visual inspection of the read coverage of RNA-seq samples and can help an investigator derive useful information from suboptimal tissue, while guiding decisions regarding whether a sample should be discarded or not.

Project name

ngs.plot.

Project home page

https://code.google.com/p/ngsplot/.

Operating system

Platform independent.

Programming language

R and Python.

Other requirements

R package doMC; Bioconductor package BSgenome, Rsamtools and ShortRead.

License

GNU GPL3.

Any restrictions to use by non-academics

Contact Lisa Placanica (lisa.placanica@mssm.edu) or the technology transfer office of Mount Sinai.