Discovery of cell-type specific regulatory elements in the human genome using differential chromatin modification analysis

Abstract

Chromatin modifications have been comprehensively illustrated to play important roles in gene regulation and cell diversity in recent years. Given the rapid accumulation of genome-wide chromatin modification maps across multiple cell types, there is an urgent need for computational methods to analyze multiple maps to reveal combinatorial modification patterns and define functional DNA elements, especially those are specific to cell types or tissues. In this current study, we developed a computational method using differential chromatin modification analysis (dCMA) to identify cell-type-specific genomic regions with distinctive chromatin modifications. We then apply this method to a public data set with modification profiles of nine marks for nine cell types to evaluate its effectiveness. We found cell-type-specific elements unique to each cell type investigated. These unique features show significant cell-type-specific biological relevance and tend to be located within functional regulatory elements. These results demonstrate the power of a differential comparative epigenomic strategy in deciphering the human genome and characterizing cell specificity.

Figure 1.

Illustration of the framework used in identifying CSREs. (A) The data profiles of nine cell types as characterized by nine marks and one control. The raw ChIP-seq reads were mapped to 200 bp bins and the signals were binarized using a Poisson null model. (B) For each bin, the dissimilarity of the resulting binary vectors between two different cell types was measured using hamming distance. For each cell type, the DMS of a bin was the summation of pairwise hamming distance (sPHD) computed between it and other cell types. (C) The DMS profile of each cell type was normalized across the genome. Then, each column of the matrix was multiplied by the corresponding Z-scores to consider the variance in the column. (D) Wavelets smoothing strategy was adopted to smooth the resulting differential profile of each cell type. CSREs were extracted by selecting suitable height and length parameters. Statistical significance (P-value) of each CSRE was calculated by a non-parametric test.

Figure 2.

Relationship between CSREs and various genomic features. (A) The distribution of CSREs in six different genomic regions, including promoter, 5′ UTR, 3′ UTR, exon, intron and intergenic regions. (B) The fold enrichments of the CSREs in the six different genomic regions. (C) Box plot of the distance between the intergenic CSREs and the nearest TSSs, compared with those of randomly generated ones. For each intergenic CSRE, the random one was an arbitrarily selected genomic element from the same chromosome with the same length. Then, the distances between the random regions to their nearest TSS were computed. (D) The normalized proportion of CSREs in each chromosome (1–22, X) in all cell types. (The bar of chromosome 22 in K562 was truncated to 0.03 for visualization, and its real number is 0.056) (E) Bar plot of the CV (defined as the ratio of the standard deviation to the mean) of normalized proportion of CSREs in each cell type.

Figure. 3.

Functional relevance and cell-type specificity of CSREs. (A) The proportion of CSREs belonging to single, two or more cell types. (B) Overlaps of CSREs between each pair of cell types. The values in the diagonal correspond to the number of identified CSREs in nine cell types and the value in row i column j records the number of CSREs in cell type i overlapped by those in cell type j. (C) CSRE neighboring genes tended to be more significantly connected than was expected, with indicated by (asterisk). (D) CSRE neighboring genes are more likely to be non-housekeeping genes with the exception of those in K562 and GM12878. Dashed vertical line on the left side represent the P-value threshold (). (E) CSRE neighboring genes show distinct functional enrichments highly relevant to the corresponding cell type contexts. We chose the top five enriched GO terms in each cell type, and was used to generate the heat map. (F) Enriched phenotypes of SNPs located in the CSREs in three cell types.

Figure 4.

Transcription levels of CSRE neighboring genes. (A) Boxplots comparison of ‘CSRE neighboring’ and ‘Other’ genes. The significance of the difference was calculated using two-sample Wilcoxon tests with . The ‘CSRE neighboring’ gene group had higher expression levels, with exception of the H1 ES cell line. (B) Bar plot of the mean intensity of H3K27me3, a repressive chromatin modification, among all CSREs. (C) The comparison of transcription diversity of all CSRE neighboring genes and remaining genes. For each gene, the CV was used to measure the diversity of expression levels across the nine cell types examined. Boxplot comparison of the CV was shown, and the statistical significance was evaluated using two-sample Wilcoxon test.

Figure 5.

Relationship of CSREs to DHSs and EP300 binding sites. The CSREs overlap (A) DHSs in eight cell types and (B) EP300 ChIP peaks, in three cell types, significantly more than randomly simulated ones. The histograms showing the overlap distributions are calculated based on random simulated CSREs, and the red curves are plotted based on kernel density estimate. The black arrows indicate the true number of overlaps. Statistical significance is measured by one-sample Wilcoxon test with e–16 for all cases.

Figure 6.

Illustration of two distinctive modification patterns revealed by CSREs: the BCR-ABL fusion gene and the INSIG1 gene. (A) The left two plots show the binary modification profiles of the two component genes that make up the BCR-ABL fusion gene in all nine cell types investigated. BCR and ABL gene are covered by one and five CSRE, respectively. Five red bottom lines indicate those five adjacent CSREs covering ABL gene. The top right two plots show the corresponding detailed modification patterns in K562, where the Philadelphia translocation results in the fusion gene. (B) The average intensity comparison of three marks before and after the BCR gene TSS. (C) Binary modification patterns of the CSRE encompassing the promoter region of the INSIG1 gene. (D) The dramatic loss of five active marks including H3K4me1/2/3, H3K27ac and H3K9ac in HepG2 (note that the heights of corresponding bars are almost zero). The modification intensity in the ‘others’ group was calculated based on the combination profiles of the other eight cell types.