Integrative analysis of haplotype-resolved epigenomes across human tissues

Abstract

Allelic differences between the two homologous chromosomes can affect the propensity of inheritance in humans; however, the extent of such differences in the human genome has yet to be fully explored. Here we delineate allelic chromatin modifications and transcriptomes among a broad set of human tissues, enabled by a chromosome-spanning haplotype reconstruction strategy. The resulting large collection of haplotype-resolved epigenomic maps reveals extensive allelic biases in both chromatin state and transcription, which show considerable variation across tissues and between individuals, and allow us to investigate cis-regulatory relationships between genes and their control sequences. Analyses of histone modification maps also uncover intriguing characteristics of cis-regulatory elements and tissue-restricted activities of repetitive elements. The rich data sets described here will enhance our understanding of the mechanisms by which cis-regulatory elements control gene expression programs.

Extended data Figure 1. Active enhancers cluster along developmental lineages

a) Pie charts showing fractions of tissue-restricted and non-tissue-restricted strong enhancers and promoters. b) Hierarchical clustering with optimal leaf ordering based on all H3K27ac marked highly active enhancers. Four major clusters are represented: early embryonic cell-types (blue), a large set of meso/endoderm-derived tissues (dark green), a set consisting of ectoderm-derived brain tissues (red) and a small cluster of mesoderm cell lines (purple), which bridged the early embryonic lineages with the somatic tissues. It is worth noting that although TRO did not fall within any clusters, it shared the highest degree of similarity with the early embryonic cell lines. On a subsequent level, two clusters are seen separating endoderm-derived tissues (gray line) and mesoderm-derived tissues (green line). Heart tissues are denoted by yellow asterisk. c) Clustering of tissues by promoters histone acetylation status shows grouping of tissues that are of similar types but are less evident in germ-layer divisions than clustering of enhancers.

Extended data Figure 2. Tissue-restricted enhancers are enriched for TF motifs important for cell identity and/or function

Significantly enriched motifs (p-value<10e-10) across all 28 tissues are divided into 29 clusters (method described in Supplementary Information). An overall p-value is generated for the enrichment of each tissue for each cluster. The figure illustrates –log(p-value) of a) pancreas b) anterior caudate and c)liver-restricted enhancer motif enrichment for the various clusters. For ease of visualization, any cluster with p-values greater than 0.05 is denoted 0. Red highlighted text refers to a subset of motif for TFs with literature support (See Supplementary Information) to have function in a) the pancreas, b) the brain and c) the liver.

Extended data Figure 3. Endogenous retroviruses (ERVs) are enriched for active cis-regulatory element marks in a tissue-restricted fashion

a) A clustered heatmap showing the H3K27ac enrichment (RPKM) of all mappable elements of the 3 classes of ERVs. b) Distribution of the Shannon-entropy of H3K27ac across enhancers, promoters and 3 classes of ERVs is shown as a density curve, demonstrating that H3K27ac enrichment of ERVs are more tissue-restricted than promoters and slightly less than enhancers. c) Boxplots illustrating the H3K27ac enrichment of 127 mappable members of the HERV-H subfamily across all tissue/cell-types. The enrichment in H1 hESCs is significantly higher than all other cell/tissues-types (p-value<1.4e-9, Wilcoxon test). d) UCSC genome browser snapshots showing example of an HERV-H element harboring H1-restricted active promoter marks, corresponding RNA-seq signal and H3K36me3 enrichment. It is note worthy that this particular element has been annotated in Refseq as the ES cell Related Gene (ESRG), a human-specific long non-coding RNA gene. e) UCSC genome browser snapshots showing example of a LTR12C element harboring TRO-restricted active enhancer chromatin marks. f) A matrix illustrating the average H3K27ac enrichment for subfamilies of class I ERVs across all cell- and tissue-types. LTR12C subfamily (green arrow) shows enrichment of H3K27ac across many distinct cell-types and tissues.

Extended data Figure 4. cREDS are enriched with dynamic histone mark signatures in different tissues and have putative cis-regulatory functions

a) Heatmaps showing the enrichment (RPKM) of the H3K27ac, H3K4me3 and H3K4me1 at MES-restricted enhancers (n=650), which are predicted as promoters in other tissues, across all 28 samples. The red box highlights the histone modifications in MES cells. b) A schematic of the pGL3-enhancer vector used in these luciferase-reporter assays (top) and the activity of 10 selected cREDS with promoter signatures and a negative control region cloned 5’ to the reporter gene after transfection into H1 hESCs (bottom). Luciferase activity of each region is normalized to the average activity of the negative controls. c) A schematic of the pGL3-promoter vector used in these luciferase-reporter assays (top) and the activity of 11 selected cREDS with enhancer signatures and a negative control region cloned 3’ to the reporter gene after transfection into H1 hESCs (bottom). Luciferase activity of each region is normalized to averaged activity of negative control regions. Error bars reflect standard deviation between 3 technical replicates

Extended data Figure 5. VISTA validated enhancers also possess dynamic histone modification signatures across tissues

Example screen shots of VISTA validated enhancers and the patterns of activity in vivo are displayed along with histone modification patterns in representative tissues (adapted from VISTA enhancer browser).

Extended data Figure 6. cREDS show enrichment of CAGE signal and putative enhancer functions in zebrafish reporter assay

a) UCSC genome browser screen shots show the 2 cREDS elements (Grey shading) harboring enhancer and promoter signatures in distinct tissue types. When compared to CAGE datasets from the FANTOM5 project, these elements show substantial overlap with transcript signals (red and blue signals indicate CAGE signal on the forward and reverse orientation, respectively). b) Selected cREDS (same elements as above) with enhancer marks in left ventricle shows heart-restricted enhancer activity, as indicated by GFP expression, in 3 days post-fertilization (3 dpf) zebrafish larvae. In parallel pT2MX negative control did not show any GFP expression. White arrow indicates location of the 3dpf zebrafish heart. For enhancer 1, 13 out of the 38 surviving embryos showed similar patterns. For enhancer 2, 18 out of the 35 surviving embryos showed similar patterns. None of the 30 surviving embryos, injected with the control vector, showed any appreciable GFP signal in the heart. (Scale bar = 50 μm)

Extended data Figure 7. Identification of widespread allelically expressed genes

a) Fraction of genes with allelically biased expression in each sample. Y-axis indicates number of samples and x-axis indicates fraction of allelically biased genes among informative genes (more than 10 SNP-containing short reads). b) Distribution of fold change of allelically biased genes between P1 and P2 alleles. c) The occurrence of allelically biased imprinted and other genes is shown. X-axis refers to the number of individual donors where corresponding allelically expressed genes are commonly detected. d) A density plot showing the fraction of sample-restricted genes with allelically biased expression (grey). Three tissue samples were randomly selected and, sample-restricted allelically expressed genes were defined, which includes random variance effect. The random selection was repeated 10,000 times. Shaded blue box indicates the range of fractions of individual-restricted allele biased genes in all analyzed tissues-types (n=10). The fraction of sample-restricted allelically biased genes is lower than individual-restricted allele biased genes in Figure 2e. e) Fold change of allele biased gene expression between two alleles are shown as scatter plot. X-axis is for the fold changes in one randomly selected tissue in each donor and y-axis is for the fold changes in all other remaining tissues in the corresponding donor. Allelic bias in one tissue is highly correlated with allelic bias in other tissues in the same individual. f) A histogram illustrates the proportions of allelically expressed genes in donor 2 (left) and 3 (right) defined in various numbers of tissues. The fraction of all testable genes or allelically expressed genes (y-axis) is calculated for the number of tissues where they are called as active (x-axis). The results indicate that the majority of allelically biased genes, as oppose to testable genes, are restricted to 1 or 2 tissue samples. KS-test was performed between allele biased genes and testable genes (p-value < 2.2e-16).

Extended data Figure 8. Allele biased chromatin states

a) Boxplots illustrating haplotype-resolved ChIP-seq signal enrichment on the two alleles at promoter regions. The P1 or P2 allele biased promoter regions were defined by H3K27ac signals and then H3K4m1, H3K4me3, and H3K9me3 signals were presented for the corresponding promoter regions. All chromatin states are consistent according to the allele biased H3K27ac patterns. KS-test was performed for p value calculation. b) Allelically biased enhancers were tested in thymus from donor 1, pancreas from donor 2 and 3. H3K27ac enrichment was tested by allele-specific ChIP-qPCR. Two control enhancers were included and showed to have no allelic biases in thymus or pancreas from donor 2 (top right and bottom left, respectively).

Extended data Figure 9. Putative enhancers with identical genotypes in different individuals exhibit similar biases in histone acetylation

a) Scatter plots of P1 allele biased enhancer activities for pairwise comparison of allele biased enhancers in donor 1 (n=85) and donor 2 (n=4,427). X- and y-axis indicate P1 allele bias. b) Scatter plot of reference allele bias of enhancer activities for pairwise comparison of allele biased enhancers in all tissues from all three donors and lymphoblastoid datasets obtained from a previous study (n=309).

Extended data Figure 10. Analyses of concordant allelically biased gene-enhancer pairs

a) Frequency of allelically expressed genes according to the distance between concordantly allele biased enhancer-gene pairs. Blue bars represent data obtained from whole chromosome-spanning haplotype blocks while green bars represent data obtained from simulated 300kb haplotype blocks. 56% of enhancer-gene pairs are more than 300kb apart. b) Accumulation curve showing fraction of allelically biased genes that have at least one concordantly allelic enhancer within a given distance (x-axis). Up to 83% of allelically expressed genes are within 2Mb of a concordantly biased allelic enhancer. c) The frequency of allele biased enhancers in donor 1, 2, and 3. Y-axis indicates fraction of enhancers and x-axis indicates frequency of allelically biased enhancers. KS-test was performed between allele biased enhancers and testable enhancers. d) Bar plots presenting the number of enhancers overlapping with DHS-QTLs and H3K27ac-QTLs for allelic enhancers, testable enhancers, and random control regions (*** - p value <10e-5).

Figure 1. Epigenome profiles of tissues reveal cREDS with dynamic histone modification signatures

a) Schematic of the cell/tissue-types profiled and their progression along developmental lineages. Samples include embryonic stem cells (H1), early embryonic lineages (mesendoderm cells(MES), neural progenitor cells (NPC), trophoblast-like cells (TRO) and mesenchymal stem cells (MSC)) and somatic primary tissues, representative of all three germ layers (Ectoderm: hippocampus (HIP), anterior caudate (AC), cingulate gyrus (CG), inferior temporal lobe (ITL) and mid-frontal lobe (MFL); Endoderm: lung (LG), small bowel (SB), thymus (TH), sigmoid colon (SG), pancreas (PA), liver (LIV) and IMR-90 fibroblasts; Mesoderm: duodenum smooth muscle (DUO), spleen (SX), psoas (PO), gastric tissue (GA), right heart ventricle (RV), right heart atrium (RA), left heart ventricle (LV), aorta (AO), ovary (OV) and adrenal gland (AD)). b) Heatmaps show H3K27ac, H3K4me3 and H3K4me1 enrichment (RPKM) at predicted lung enhancers (n=1,321), which are defined as promoters in other tissues, across all 28 samples. Red box highlights the signatures in lung. c) A UCSC genome browser snapshot of a region on chromosome 20, showing the chromatin states of a cREDS element (gray shading) predicted as a promoter in psoas and an enhancer in lung. d) A boxplot of RNA-seq signals (RPKM) overlapping ±1kb of cREDS enhancers, cREDS promoters, non-cREDS control enhancers and non-cREDS control promoters. (*** indicates p-value<10e-142, Wilcoxon test) e) RNA-seq and chromatin states of a cREDS element (gray shading) is shown for a region on chromosome 17 in H1 and IMR-90. Arrow indicates an alternate exon incorporated in IMR-90.

Figure 2. Widespread, individual-specific allelic bias in gene expression

a) Genome browser snapshots illustrate completeness and resolution of haplotypes resolved in donor 4. Y-axis indicates the number of variants within 100kb windows. The density of all (blue), phased (orange) and unphased (grey) variants across chromosome 1 are shown. b) Proportion of genes with allelically biased expression among informative genes and the number of tissue samples derived from each donor (n_tissue) are described. c) Boxplot illustrates occurrence of imprinted and other allelically biased genes across samples. (*** - p-value<9.9e-5, KS-test) d) Including only tissues with 2 or 3 equivalent samples derived from distinct donors (n_tissue=10), genes with allelic imbalances were defined as common between individuals (consistent bias among same tissue-type from multiple donors) or as individual-restricted. Random control represents average from randomly selected samples (10,000 iterations). e) Fold change of gene expressions between alleles in AD from donor 2 (x-axis) is compared to all other tissues from donor 2 (y-axis). f) A histogram illustrates the proportions of allelically expressed genes in donor 1 defined in various numbers of tissues. The fraction of all testable genes or allelically expressed genes (y-axis) is calculated for the number of tissues where they are identified as active (x-axis)(p value<2.2e-16, KS-test).

Figure 3. Characterization of allele bias in chromatin states at cis-regulatory elements

a) Boxplots present haplotype-resolved ChIP-seq reads at promoter or gene bodies (H3K27ac: n=744, p-value=10e-14, H3K4me1: n=32 p-value=0.035, H3K4me3: n=177, p-value=0.0047, H3K27me3: n=12, p-value=0.43, H3K9me3: n=27, p-value=0.13 and H3K36me3: n=291, p-value=4.3e-6, KS-test). b) Allelically biased gene expression of IQCE is concordant with chromatin marks at the promoter (grey) and gene body. c) Proportion of allelic (n=11,714) and non-allelic (n=89, 599) among all informative enhancers (n=101,313) across 20 tissues. d) A snapshot showing a SNP (rs138143205) with H3K27ac bias towards the G allele in both LV donors (Left). Bar chart illustrates the number of H3K27ac reads corresponding to the P1 versus P2 alleles in both donors (Right) (*** - p-value<10e-19, binomial test). e) Scatter plots show strong correlation of the P1 allele bias of enhancer activities among two different tissue-types from donor 3 (n=4,427) and f) among the P1 allele bias in donor 3 (x-axis) and the allele bias of corresponding genotypes in donor 1 or 2 (y-axis) at allelic enhancer in the same tissue-type (n=447).

Figure 4. Allelic histone acetylation at enhancers is associated with allelically biased gene expression

a) Average distance of allelic (5% FDR) and non-allelic enhancer to the closest allelically expressed gene is significantly different (n=3,829 *** - p-value<2.2e-16, KS-test). b) Genome browser snapshots show an allelic enhancer within the intron of the allelically expressed A4GALT gene (P1- red, P2 – blue) on chromosome 22 across 3 samples. c) Density plot presents the fraction of concordant allelic bias between allelically expressed genes and allelic enhancers in terms of distance. The allelic enhancer-gene pairs were defined with FDR cutoff values of 5% (n=14,082)(black), 1% (n=6,057)(blue) and 0.1% (n=2,362)(yellow). Permutated control of a set of enhancer-gene pairs was included (n=14,082)(grey). Distance between allelically biased enhancer-gene pairs and fraction of concordant allelic bias are denoted by x- and y-axes, respectively (p-value<2.2e-16, KS-test). d) Fractions of tissue-restricted enhancer-gene pairs (y-axis) that show concordant (blue) or discordant (orange) allelic biases in the same tissue, are presented across a range of Pearson correlation coefficients (x-axis) (p-value < 2.2e-16, KS-test, random permutated control concordant pairs = 50%). (e) Overlap between eQTLs and allelic enhancers, testable enhancers or random control regions are shown. Error bars represent standard deviations. Testable enhancers and random control regions were generated 10,000 times with the same numbers as allelic enhancers (*** - p-value<10e-5).

Figure 5. Motif disruption by genetic variants is concordant with allelic H3K27ac biases at enhancers

a) Differential GABPA binding motif scores between two alleles (P1-P2 motif scores) in LV is correlated with the proportion of H3K27ac reads corresponding to the P1 allele (top). Values range from negative to positive, indicating P1 and P2 motif disruption, respectively. An example on chromosome 12 illustrates P1, with a motif preserving C allele, has higher H3K27ac enrichment and the P2, with the motif disrupting T allele, has little H3K27ac enrichment (bottom). b) Three examples (FLI1 in SX, SPDEF in SG, and TEAD in AO) of motif-disrupted enhancers demonstrate allelic biased activities. The variant location and genotypes of P1 and P2 alleles are marked in motif sequence. c) All possible motif disrupted enhancers-gene pairs within the indicated distance window are defined with concordant allelic bias (blue, gene-enhancer pairs with biases towards the same allele) or discordant allelic bias (red, gene-enhancer pairs with biases towards different alleles). Only TH, LV and AO were considered due to the availability of Hi-C data. Short-range pairs are defined if any allelically expressed genes are located <20kb away. (*** - p-value<2.5e-5, binomial test).