article: In mammals, numerous autosomal genes are subject to mitotically stable monoallelic expression (MAE), including genes that play critical roles in a variety of human diseases. Due to challenges posed by the clonal nature of MAE, very little is known about its regulation; in particular, no molecular features have been specifically linked to MAE. In this study, we report an approach that distinguishes MAE genes in human cells with great accuracy: a chromatin signature consisting of chromatin marks associated with active transcription (H3K36me3) and silencing (H3K27me3) simultaneously occurring in the gene body. The MAE signature is present in ∼20% of ubiquitously expressed genes and over 30% of tissue-specific genes across cell types. Notably, it is enriched among key developmental genes that have bivalent chromatin structure in pluripotent cells. Our results open a new approach to the study of MAE that is independent of polymorphisms, and suggest that MAE is linked to cell differentiation.
  A variety of genetic and epigenetic factors affect the relative expression levels of the two copies of each given gene in diploid cells. In addition to cis- and trans-regulatory variation (Gilad et al. , 2008), there are at least three major kinds of non-Mendelian phenomena that control allele-specific expression in mammals. One is the X chromosome inactivation (Lyon, 1961): in female embryos, around the time of implantation, about half of the cells choose to inactivate the maternal X, and the rest inactivate the paternal X, affecting most of the X-linked genes (Carrel and Willard, 2005; Berletch et al. , 2010; Yang et al. , 2010). Another is imprinting: genes such as IGF2 and H19 are expressed either from one allele, either paternal or maternal (Glaser et al. , 2006). Finally, a significant fraction of mammalian autosomal genes are subject to monoallelic expression (MAE), which reflects a mitotically stable allele-specific expression with different allelic states in clonal lineages. MAE is observed in olfactory receptor genes (Chess et al. , 1994), as well as genes coding for immunoglobulins and some cytokines (Pernis et al. , 1965; Bix and Locksley, 1998; Holländer et al. , 1998). Using genome-wide analyses of allele-specific expression, we and others have added a surprisingly large number of the autosomal genes in human and mouse to the MAE class (Gimelbrant et al. , 2007; Jeffries et al. , 2012; Zwemer et al. , 2012; Li et al. , 2012b), including genes implicated in a number of human diseases, such as Alzheimer’s disease (APP) (Bertram and Tanzi, 2012) and cancer (DAPK1) (Raval et al. , 2007). MAE affects about 10% of ∼4000 tested genes in human lymphoblastoid cells (LCLs) and about 15% of more than 1300 assessed genes in analogous mouse cells (Gimelbrant et al. , 2007; Zwemer et al. , 2012). Our growing appreciation of the prevalence of MAE only underscores how little we know about its biology. The only existing, large-scale sets of data are collected in clonal lymphocyte cell lines in vitro. The limited number of analyzed clones is insufficient to generate a complete catalog of MAE genes in that cell type, and little is known about the prevalence of MAE in other cell types. Virtually nothing is understood about the establishment of MAE during development and differentiation. Mechanistically, allelic choice has been linked to changes in chromatin states in some special cases: imprinting (Wen et al. , 2008), olfactory receptor gene choice (Magklara et al. , 2011), and immunoglobulin-kappa gene rearrangement (Farago et al. , 2012). In contrast, for hundreds of other autosomal MAE genes, no molecular features have been associated with establishment and maintenance of allelic choice. Similarly, there is no general understanding of MAE’s function. A major technical bottleneck in addressing these questions is the clonal nature of MAE (Figure1—figure supplement 1). Like X inactivation, MAE is masked in polyclonal samples, and obtaining monoclonal cell populations is challenging for most tissue types, particularly so in vivo. Moreover, genome-wide methods are limited by the availability of polymorphisms. In this study, we report a fundamentally new approach to the detection of monoallelic expression. In contrast to other methods, it does not require any allele-specific information, instead relying on a specific chromatin pattern as a proxy for MAE. We use this approach to address questions about MAE’s prevalence, development, and function.
  Histone modifications, in their diversity, present rich combinatorial possibilities for controlling gene transcription (Barski et al. , 2007; Ernst and Kellis, 2010; Filion et al. , 2010). They therefore offer a constrained yet rich set of data for systematic analysis. To identify histone modifications that might be specific to the MAE genes, we compared chromatin marks associated with known MAE genes against those for known biallelically expressed (BAE) genes. Previous observations in human and mouse cells suggested that a gene could show MAE in one cell type but not another (Gimelbrant et al. , 2005; Gimelbrant and Chess, 2006; Gimelbrant et al. , 2007). We therefore hypothesized that if there were a correspondence between MAE state and chromatin modifications, it would be more pronounced when comparing data from the same cell type, while not necessarily from exactly the same cells. Since the largest available sets of known human MAE and BAE genes were identified in human lymphoblasts (Gimelbrant et al. , 2007) (see also Dataset S1 in Dryad, Nag et al. , 2013), we used these sets to compare with histone modification data for GM12878 lymphoblastoid cells deposited by the ENCODE project (Dunham et al. , 2012). We focused on the eight marks that were investigated in a broad variety of cell types: H3K27me3 (histone H3 Lys-27 trimethylation), H3K36me3, H3K4me2, H4K20me3, H3K27ac (histone H3 Lys-27 acetylation), H3K4me1, H3K4me3, H3K9ac (Figure1—figure supplement 2A). To enable our analysis, we reduced complex patterns of histone modifications to two simple features capturing signal intensity in two distinct spatial domains: the proximal promoter signal (for any given modification: ChIP-Seq signal intensity integrated over the 2. 5 kb region upstream of transcription start) and the signal integrated over the whole gene body (green and red areas in Figure 1A; see ‘Materials and methods’, Figure 1—figure supplement 2). We then set out to analyze in a systematic way whether some combination of these measured signals can reliably distinguish known MAE genes and known biallelic genes. 10. 7554/eLife. 01256. 003Figure 1. Genes with monoallelic expression have specific chromatin signature within the gene body. (A) Assessment of histone modifications. The mapped ChIP-Seq signals for the listed modifications were derived from the total signal over the gene-body or promoter region: shown is the gene body signal for the two most informative chromatin marks H3K36me3 (green) and H3K27me3 (red). EBF1 gene was shown to be MAE, ABCC1 was shown to be biallelic in lymphoblastoid cells (Gimelbrant et al. , 2007). ChIP-Seq data in GM12878 lymphoblasts were generated by the ENCODE project. Graphics adapted from UCSC genome browser (http: //genome. ucsc. edu/; Meyer et al. , 2013). Height of the signal tracks was set 0–8. (B) High confidence MAE (blue) and biallelic (gold) autosomal genes in the training set are separated by the gene body signal for H3K27me3 and H3K36me3 in GM12878 cells. Light blue area illustrates partitioning of this space by the most optimal classifier (DT2F). Solid line demarcates external border of ‘Neutral’ setting; dotted line shows more restrictive ‘Precision’ setting and is a graphical representation of the boundary identified by an alternating decision tree (DTree), which was the best-performing machine learning method applied to the features after feature selection. Of 270 high confidence MAE genes, 268 had data for both H3K27me3 and H3K36me3. Of these, 204 (76%) are within predicted MAE region. (C) Distribution of all autosomal RefSeq genes in GM12878 cells according to gene body signal for H3K27me3 and H3K36me3. Genes are color-mapped according to their expression level in GM12878 cells, from lowly expressed in red to highly expressed in yellow. Silent transcripts (RPKM <= 0. 1) are shown in gray. Solid and dotted lines as in 1B. (D) Fraction of predicted MAE genes as a function of gene expression level. Left vertical axis: absolute number of predicted MAE (blue) and non-MAE genes (gold) per expression level bin. Right axis: fraction of predicted MAE genes (red circles) per same bin. Expression bins are 0. 1 log10 units of RPKM in GM12878 cells. (E) Genome distribution of predicted MAE and biallelic genes and their expression level. Shown is chromosome 19; other autosomes are similar. Blue—genes predicted as MAE; gold—genes predicted as biallelic. Position along the chromosome corresponds to transcription start site of the gene; marker length reflects gene expression level in GM12878 cells. Only genes with RPKM > 1 are shown. DOI: http: //dx. doi. org/10. 7554/eLife. 01256. 00310. 7554/eLife. 01256. 004Figure 1—figure supplement 1. Chromatin signature of monoallelic expression allows its detection in monoclonal and polyclonal samples. Detection of MAE by expression bias is not possible in polyclonal cell populations as both paternal and maternal transcripts are present, making expression appear biallelic. H3K36me3 is indicated by green circles and H3K27me3 is indicated by red circles. DOI: http: //dx. doi. org/10. 7554/eLife. 01256. 00410. 7554/eLife. 01256. 005Figure 1—figure supplement 2. Building and performance of chromatin feature classifiers. (A) The mapped ChIP-Seq signals for the listed modifications were derived from the total signal over the gene-body (green) or 2. 5 kb promoter region (red). EBF1 gene was shown to be MAE in lymphoblastoid cells (Gimelbrant et al. , 2007). ChIP-Seq data in GM12878 lymphoblasts were generated by the ENCODE project. Graphics adapted from UCSC genome browser (http: //genome. ucsc. edu/; Meyer et al. , 2013). Height of the signal tracks was set 0–8. (B) Comparison of precision and recall of different classifier types when using distinct sets of chromatin features. False positive (FP) and false negative (FN) calls for training set of MAE and BAE genes are shown as function of the increasing cost of false positive errors. Classifiers shown: DT–Decision Tree; NB–Naïve Bayes. Feature sets: ‘7 features’–gene body signal for H3K27me3 and H3K36me3; and promoter signal for H3K27me3, H3K36me3, H3K4me2, H4K20me3, and H3K27ac; ‘2 features (also called DT2F) ’—only gene body signals for H3K27me3 and H3K36me3. Neutral and Precision settings were chosen, respectively, for best recall, and for the optimal combination of recall and precision. (C) Comparison of the 2-feature (GeneBody) and 7-feature (GenePromoterAndBody) classifiers. Similarity of precision and recall values suggests that the two chromatin marks, H3K27me3 and K3K36me3, account for most of the difference between MAE and BAE genes. DOI: http: //dx. doi. org/10. 7554/eLife. 01256. 00510. 7554/eLife. 01256. 006Figure 1—figure supplement 3. Distribution of various promoter and/or normalized gene body signal combinations in GM12878 cells in our training set. High confidence MAE (blue) and biallelic (gold) autosomal genes in the training set do form clusters in some of these cases but fail to achieve as clear a separation of MAE and biallelic genes as does normalized H3K27me3 and H3K36me3 gene body signal combination (Figure 1B). Data are shown for a few representative cases. DOI: http: //dx. doi. org/10. 7554/eLife. 01256. 006 Using all available features, we wanted to identify parameters that would offer the optimal tradeoff between low false positive rate and the largest possible number of correctly identified MAE genes. After systematic exploration detailed in Figure 1—figure supplement 2B, we chose two relative penalty tradeoff settings: ‘Neutral’ (1: 1) and ‘Precision’ (8: 1). The best-performing classifier, Decision Tree (DT), identified at the more relaxed ‘Neutral’ stringency setting almost 80% of the known MAE genes with 20% false positive rate. At a more stringent, precision-optimized setting, the DT classifier identified only 10% false positives, maintaining ∼60% of true positives (Figure1—figure supplement 2B). Strikingly, performance of the DT classifier with just two gene-body features, H3K27me3 (associated with gene silencing) and H3K36me3 (associated with gene transcription), was as good as the performance of any full-featured model (Figure1—figure supplement 2C). Other feature combinations also had some discriminatory power, but not as significant (Figure 1—figure supplement 3). This suggests that H3K27me3 and H3K36me3, taken together, account for most of the distinction between MAE and BAE genes. Figure 1B also illustrates the partition of the phase space defined by these two signals by the most optimal classifier (two-feature Decision Tree; DT2F; see ‘Materials and methods’) at both high and medium stringency settings (Precision and Neutral). Performance of the DT2F classifier was at least as good as the performance of any classifier, including any other tested feature combination (Figure 1—figure supplement 2). We concluded that known autosomal MAE and BAE genes in this study show consistent, characteristic differences with respect to the gene body ChIP-seq signal for H3K27me3 and H3K36me3. Co-occurrence of these two signals in the same gene body is thus a specific chromatin signature of known MAE genes. Importantly, MAE identified in clones is reflected by the chromatin signature in a polyclonal sample. To broadly test this chromatin signature as a predictor of the MAE state, we identified novel genes with the MAE chromatin pattern, and assessed their allelic bias in monoclonal cell lines. Specifically, we calculated gene-body signal for H3K27me3 and H3K36me3 for all autosomal RefSeq genes based on the ENCODE ChIP-seq data for GM12878 lymphoblastoid cells (Figure 1C). In these cells, the classifier predicted as monoallelic 1315 (13%) of 10,322 autosomal genes with moderate or higher expression (RPKM > 1), and 15 genes (3%) of 450 highly expressed genes (with RPKM > 100) (Figure 1D; detailed in Dataset S1 in Dryad, Nag et al. , 2013). As a group, the predicted MAE genes were spread throughout all autosomes. They were interspersed with biallelic genes and showed a variety of expression levels (Figure 1E). This is consistent with MAE genes identified earlier (Gimelbrant et al. , 2007; Zwemer et al. , 2012). To test for transcriptional allelic bias in predicted MAE genes, we used two complementary approaches. For a broad, genome-wide analysis, we used mRNA-seq in two independent (with different direction of X-inactivation) monoclonal cell lines, DF1 and DF2. We derived these lines from GM12878 line (see ‘Materials and methods’), because both its own genome and its parental genomes were fully sequenced by the 1000 Genomes Project (Consortium, 2010). We generated the RNA-seq data from the DF1 and DF2 clones and analyzed them using a custom analysis pipeline (see ‘Materials and methods’). Of the 48 X-linked genes with SNPs covered by 10 or more reads, 43 (∼90%) showed positive evidence of clone-specific allelic bias; the rest were inconclusive; positive evidence for equivalence was not detected in any X-linked genes (Dataset S2 in Dryad, Nag et al. , 2013). By contrast, of the 3270 autosomal genes with comparable coverage, 1167 (35%) showed positive evidence of equivalent expression of both alleles and no evidence of bias in either clone; 269 (8%) had allelic bias, while the rest were inconclusive. Altogether, there were 5001 autosomal genes with at least one SNP covered at any level (Dataset S2 in Dryad, Nag et al. , 2013). Of these, 1021 genes were predicted to be MAE by the DT2F classifier at the Neutral setting and 236 genes at the Precision-optimized setting. For a quantitative bias comparison, we used all expressed genes. For additional control, we also used equally sized sets of genes with matched levels of expression, randomly chosen from these remaining genes. Both control groups showed quite small mean allelic bias, about 60: 40. By contrast, genes predicted MAE by the Neutral DT2F classifier showed about 75: 25 bias, and the genes predicted by the classifier at Precision setting had mean bias of 90: 10 (Figure 2B). To estimate probability of error, we sampled 10 sets of 40 predicted genes and 10 equally-sized control sets, matched by expression; the comparison of the mean bias in the 10 sampled sets showed highly significant difference (p<9e−05; non-paired t test). Thus the predicted MAE genes had significantly higher bias than the control genes. In subsequent analysis, we used the neutral classifier setting in order to maximize the number of candidate MAE genes and scrutinize predictive properties of this less stringent setting. 10. 7554/eLife. 01256. 007Figure 2. Prediction testing with RNA-Seq. (A) Representative examples of allelic counts in data from two clones (DF1 and DF2) derived from GM12878 cells. Shown are total maternal (Mat; ‘pink’) or paternal (Pat; ‘blue’) counts for X-linked genes and autosomal monoallelic genes illustrating that the direction of allelic bias is clone-specific. (B) Mean allelic bias in different groups of genes in DF1 and DF2 clones as assessed by the RNA-Seq analysis. ‘50’ corresponds to perfect balance between alleles; ‘100’ to perfectly monoallelic expression. ‘Precision’ and ‘Neutral’—all informative expressed genes predicted as MAE using corresponding settings of the DT2F classifier; ‘All’—all informative expressed genes from GM12878 cells; ‘RPKM matched’—predicted biallelic genes, matched by the expression level to the predicted MAE genes (shown are the mean and standard deviation for 10 permuted sets of genes). (C) Definitions of allelic bias and lack of bias. Unbiased genes (gold) pass equivalence test (with 2: 1 boundaries in either direction; equivalence area is light yellow); biased (blue) pass binomial test; genes that pass neither statistical test are called inconclusive (gray). See Figure 2—figure supplement 1 for X-chromosome analysis results according to this scheme. (D) Fraction of genes showing allelic bias in DF1 and DF2 clones as assessed by RNA-seq. Biased genes were called in at least one clone based on FDR-corrected binomial test and displayed at least 2: 1 bias. Unbiased genes were called based on passing the equivalence test in at least on clone and not passing the bias test in the other clone. (E) Allele-specific analysis of RNA-Seq data from DF1 and DF2 clones. Experimentally determined allelic states of autosomal genes. Predicted monoallelic and biallelic status is based on the neutral DT2F classifier. Assignments of genes as biased, unbiased or inconclusive (indeterminate in both clones). Color-coding as in panel C. DOI: http: //dx. doi. org/10. 7554/eLife. 01256. 00710. 7554/eLife. 01256. 008Figure 2—figure supplement 1. Allelic bias calling on X-Chromosome. Application of the bias definition to X-linked genes in the DF1 and DF2 clones. Each rectangle corresponds to one gene; only informative genes are shown; genes are in the correct order on the chromosome but positions are not to scale. Blue—biased with paternal expression; magenta—biased with maternal expression (as defined in Figure 2C). DOI: http: //dx. doi. org/10. 7554/eLife. 01256. 008 Next, we used the RNA-Seq data to categorize predicted and control genes as biased, unbiased, or indeterminate (Figure 2C). Biased expression was identified based on FDR-corrected binomial testing and allelic skewing of at least 2: 1 (see ‘Materials and methods’). Importantly, rejection of the bias hypothesis by this test does not automatically mean the gene could be called unbiased. Therefore, we used equivalence testing (Limentani et al. , 2005), with equivalence boundaries corresponding to the two-fold imbalance; genes that failed both tests were called indeterminate. Genes predicted by the DT2F neutral classifier were enriched for genes with positively identified allelic bias; the precision classifier setting, as expected, yielded still better enrichment but fewer positively identified genes (Figure 2D, E). This RNA-Seq approach confirms MAE predictions on a whole-transcriptome level, but it has significant limitations. Insufficient coverage depth leaves an overwhelming majority of genes as ‘indeterminate’ (Figure 2E). This results in underestimation of both the true positive and the true negative rates. Furthermore, a large majority of known MAE genes (about 85%) show biallelic expression in some clonal lineages (Gimelbrant et al. , 2007; Zwemer et al. , 2012). This is highly important when considering any validation experiments, since even exhaustive analysis of just two independent clones would miss monoallelic expression in many such genes that would happen to show biallelic expression in the two assessed clones. To validate MAE predictions more conclusively, we measured allelic bias in a greater number of independent clones. To simultaneously increase both coverage depth and the number of biological samples, we designed a targeted extra-deep RNA-Seq assay (allele-specific targeted sequencing; AST-Seq) that allowed us to precisely quantify allele-specific expression of a subset of genes in an increased number of clones (see Figure 3A). 10. 7554/eLife. 01256. 009Figure 3. Prediction testing with allele-specific targeted sequencing (AST-Seq). (A) Schematic representation of the deep barcoding approach. As an illustration, analysis of three clones with no multiplexing is shown, each with a different allelic bias at a SNP of interest. Random-primed cDNA or genomic DNA are used as templates for PCR1, using gene-specific primers with universal tails. The next step associated universal amplicon tails in each sample with two barcodes (PCR2); this allows for barcoding a large number of samples with limited number of secondary primers. For a given sample, all amplicons share the same two barcodes. Barcoded amplicons from all samples are pooled, and sequencing adaptors attached. After sequencing and deconvolving by barcode, allelic hits are counted. (B) Representative allelic counts using the AST-seq. Allelic bias was assessed in two clonal lines, DF1 and DF2, derived from GM12878 and four clones, H2, H7, H14, and H16, from GM13130 (‘H0’) cells. Target SNPs were chosen to be informative in both cell lines. Genomic DNA (gray) was used as a control for allelic bias introduced in amplification; only unbiased assays were pursued. Shown are representative assays for X-linked genes (as control), and examples of genes predicted MAE or biallelic based on the chromatin signature in GM12878. Pink: expression bias towards reference (Ref) allele; blue: expression bias towards alternative (Alt) allele; gold: unbiased expression; no color: counts below threshold—data ignored. Note that, as expected, genes with clone-specific MAE could be biallelic in some clones. (C) Summary of the AST-seq analysis for all tested genes in six clonal samples. Biased (blue) and unbiased (gold) expression as defined in Figure 2C. DOI: http: //dx. doi. org/10. 7554/eLife. 01256. 009 To assess both false negative and false positive rates for predictions by the DT2F classifier, we chose a set of predicted, unconfirmed MAE genes expressed in both assessed clones, and a comparable random set of predicted biallelic genes (see ‘Method note 2’). Previously, we had derived and characterized several independent clones from GM13130 lymphoblastoid cells (Gimelbrant et al. , 2007). Starting with four of these clones and the two clones from GM12878, we selected SNPs heterozygous in both genotypes. To control for possible genotyping errors and amplification bias, we used genomic DNA from the same cells. After removing SNPs that did not pass the equivalence test in the gDNA (cf. Figure 2C), we had SNPs in 17 predicted MAE genes and 28 predicted biallelic genes. As templates, we used DNA and DNase-treated nuclear RNA from original cell lines and the clones (Gimelbrant et al. , 2007); as a positive control for expression bias, we included X-linked genes Dataset S3 in Dryad (Nag et al. , 2013). We performed two sequencing experiments, with overlapping samples, obtaining 1. 6 M reads on average per sample per run. Figure 3B shows a representative subset of all classes of assessed genes and associated allelic counts (complete data in Dataset S3 in Dryad, Nag et al. , 2013). While alleles were equally represented in genomic DNA, the X-linked genes (e. g. , PIR and XIAP) showed strong clone-specific bias; genes predicted as biallelic (e. g. , CEP110 and HDLBP) showed no significant bias. Most autosomal genes predicted as monoallelic (e. g. , FRMD6 and IGF1R) showed the characteristic pattern of clone-specific monoallelic expression: strong bias in some clones, and biallelic expression or opposite bias in others. The complete results of this experiment (see Dataset S3 in Dryad, Nag et al. , 2013) are summarized in Figure 3C. Of the 17 predicted autosomal MAE genes, 13 exhibited strong allelic bias in at least some of the assessed clones, with 12 of those showing different direction of bias between individual clones (false positive rate of about 24%). Note that the false positive rate should be treated as an overestimate, since two of the unconfirmed genes failed in the GM13030-derived clones, while two others failed altogether. In addition, the number of clones we assayed, while higher than in the first experiment, is not sufficient to establish at high level of certainty that these are truly biallelic genes. Conversely, of the 28 predicted biallelic genes, none showed significant allelic bias in expression in any of the clones. We can thus conclude that the MAE chromatin signature is a specific and sensitive predictor of clone-specific MAE status. We asked how the allelic distribution of the active and repressive chromatin marks in clonal cell lines relates to the transcriptional allelic bias. A multiplexed padlock probe approach (Zhang et al. , 2009) enabled us to assess allelic bias in heterozygous exonic SNPs in two clones with GM12878 genotype, and four clones from GM13130 cells. After removing assays that failed equivalence test in gDNA in all samples, we had 482 SNPs (see Dataset S4 in Dryad, Nag et al. , 2013). We used this approach to assess allelic bias in H3K27me3 and H3K36me3 ChIP samples simultaneously with cDNA from the same cells, as well as ChIP input and genomic DNA controls. Figure 4A shows representative examples of allelic counts for all classes of expression bias using a small number of SNPs in the two GM12878 clones (complete dataset in Dataset S4 in Dryad, Nag et al. , 2013). Allelic expression bias is evident in an imprinted gene SNRPN, on the X chromosome (SLC25A43 and XIAP), and on autosomal MAE genes. Biased expression is accompanied by an asymmetric distribution of the two chromatin marks: H3K36me3 is associated with the higher transcribed allele, while H3K27me3 with the fully or partially silenced allele. 10. 7554/eLife. 01256. 010Figure 4. Correlation of allelic bias in expression with bias in chromatin marks. (A) Representative examples of allelic counts in SNPs assessed with multiplexed targeted sequencing using padlocked probes. Apart from shown clones, additional clones from GM13130 individual were assessed (see ‘text’). Shown as control are an imprinted gene SNRPN and X-linked genes. Other genes were predicted MAE or biallelic based on their gene-body chromatin signature in GM12878 cells. Measurement for each SNP is shown as read counts for the reference and alternative (Ref/Alt) alleles as designated in dbSNP. Color-coding as in Figure 3B. Analysis summarized in this figure is based on 482 SNPs within 458 genes. (B) Correlation of allelic bias in H3K27me3 with allelic bias in cDNA. All informative SNPs were put into one of three bins according to their cDNA allelic bias: unbiased (Figure 2C); significantly biased towards reference allele; significantly biased towards alternative allele. For each of these groups, allelic bias for the ChIP sample from the same clone was assessed and analyzed using Kruskal–Wallis non-parametric ANOVA test. (C) Same as B, for H3K36me3. (D) Fraction of SNPs in predicted MAE and biallelic genes, showing allelic bias of 2: 1. Bias calls were made by binomial testing in cDNA, H3K36me3 and H3K27me3 ChIP for SNPs, using data from the padlock probe experiments. DOI: http: //dx. doi. org/10. 7554/eLife. 01256. 010 We used the complete dataset Dataset S4 in Dryad (Nag et al. , 2013) to evaluate global relationships between expression bias and chromatin allelic bias. In order to pool data from two individuals, one of which (GM13130) lacked complete genotypes for parents, we assessed SNP bias as reference and alternative alleles (rather than maternal or paternal bias). SNPs in cDNA were assigned to one of three bins: reference allele bias; no bias; and alternative allele bias. For these groups, allelic bias in H3K27me3 (Figure4B) and H3K36me3 (Figure 4C) was determined. In unbiased loci, both H3K27me3 and H3K36me3 were equally represented. In contrast, preferential expression of an allele was associated with elevated levels of H3K36me3 and decreased levels of H3K27me3 on that allele. Both effects were highly significant (p<2 × 10e−9). Genes predicted to have MAE were about fourfold over-represented among genes where SNPs showed significant bias (Figure 4D). SNPs with skewed H3K27me3 and H3K36me3 distribution were highly enriched in the genes predicted as MAE (p<10e−6 and p=0. 01, respectively; two-tailed Fisher’s exact test). This suggests that the asymmetric distribution of the histone modifications is to a large extent due to the genes that have the chromatin signature of monoallelic expression. Using RNA–DNA FISH, we have previously shown (Gimelbrant et al. , 2007) that individual MAE genes identified in lymphoblasts also show monoallelic expression in vivo, in peripheral blood mononuclear cells (PBMCs). We used the trained DT2F classifier to analyze available ChIP-seq data from PBMCs (Bernstein et al. , 2010). In this and the following analyses we applied the classifier to all genes with evidence of transcription (excluding only genes that were not expressed; with RPKM<0. 1), to achieve the most comprehensive possible coverage (‘Method note 1’). The overall distribution of genes in the DT2F phase space was very similar in cells in vitro and in vivo, with two principal clusters at high H3K36me3 and high H3K27me3 (Figure 5A). Predicted genes in PBMC largely overlapped with our predictions for LCLs: out of 2057 genes showing the MAE chromatin signature in PBMC, and called in GM12878 cells, 83% (1712) were also predicted MAE there. This is particularly striking since the experimental data for these two related cell types were collected in two different laboratories. This suggests that the genes with MAE signature in lymphoblasts are quite similar to those in the related ex vivo cells. 10. 7554/eLife. 01256. 011Figure 5. MAE chromatin signature in multiple cell types. (A) Comparison of H3K27me3 and H3K36me3 ChIP-Seq gene body signal distribution for the autosomal genes in GM12878 cells (red) and in primary peripheral blood monocytes (PBMC; blue). Silent genes (RPKM < 0. 1) are excluded in either case. Both datasets were collected by ENCODE; PBMC data: GSE16368. Note that the dots are made more transparent than in Figure 1 to make clear the overall shape of the distribution. (B) Overlap of predicted MAE genes in GM12878 and PBMC (silent genes are excluded). (C) Tissue-specific distribution of MAE genes. Overlap between predicted MAE genes in three cell types as labeled. Within dotted circle: genes expressed in all three lines (and MAE in at least one). Outside dotted circle: MAE genes showing cell type-specific expression (predicted MAE and expressed in that cell type, but silent in at least one of the two other cell lines). (D) Similarity of predicted MAE gene sets in seven cell types: GM12878—lymphoblasts, K562—acute myelocytic leukemia, H1ESC—human embryonic stem cells, HSMM—human skeletal muscle myocytes, HUVEC—human vascular epithelium, HMEC—human mammary epithelium, HCC1954—breast cancer. Similarity assessed according to the Jaccard similarity measure. In the heat map, darker gray corresponds to higher similarity. ChIP-Seq and RNA-Seq data sources: see Dataset S1 and S2 in Dryad (Nag et al. , 2013). (E) Cumulative number of predicted MAE genes as a function of the number of cell lines assessed. Counted are only genes expressed in all analyzed cell lines. Order of addition of cell lines was sampled by permutation, shown are mean values ± standard deviation. (F) Cumulative number of all genes and predicted MAE genes as a function of the number of cell lines assessed. Counted are all genes with any evidence of expression in at least one cell type. (G) Gene Ontology (GO) analysis of genes predicted MAE in indicated cell types. Most over- and under-represented categories for GM12878 cells are also shown for other cell types (in each cell line, predicted MAE genes are compared to all expressed genes). Horizontal axis: −log10 (p), after Benjamini–Hochberg correction. Gray lines correspond to −log10 (p) values as noted. DOI: http: //dx. doi. org/10. 7554/eLife. 01256. 01110. 7554/eLife. 01256. 012Figure 5—figure supplement 1. Overall distributions of the normalized H3K27me3 and H3K36me3 gene body signal in the analyzed cell types. Comparison of H3K27me3 and H3K36me3 ChIP-Seq gene body signal distribution for the autosomal genes in a number of human cell lines (as noted). Each dot represents an autosomal gene. Silent genes (RPKM < 0. 1) are excluded. Light blue area illustrates partitioning of this space by the most optimal classifier (DT2F). DOI: http: //dx. doi. org/10. 7554/eLife. 01256. 012 We then asked if that similarity is in fact due to biological relatedness of LCLs and PBMCs; alternatively, it could be due to the signature being uniform in all cell types. We applied our analysis to all autosomal genes expressed in several cell types analyzed by ENCODE: vascular and mammary epithelial cells (HUVEC and HMEC), skeletal muscle (HSMM), embryonic stem cells (H1ESC), a leukemia line K562 and breast cancer line HCC1954, drawing on deposited data to analyze with our model: ChIP-Seq of H3K27me3, H3K36me3, Input control, and RNA-Seq data (Figure 5C–D). In the DT2F space, genes formed two major clusters in all assessed cell types (Fig. 5—figure supplement 1). Although the precise positions of these clusters varied somewhat, the overall distribution remained consistent. Importantly, the DT2F classifier trained on the GM12878 data consistently covered the lower cluster of genes in different assessed cell types. In contrast to the similarity between PBMC and GM12878 data collected by different labs, cells of different types, showed pronounced differences in the set of genes classified by the DT2F as MAE, even though they were assessed in the same laboratory. This, along with the similarity in the general structure of the feature space (Fig 5—figure supplement 1) suggests that the biological properties of the cell type are captured by the model to a greater extent than other potential sources of variation and provides indirect support for the use of our method in other cell types. Figure 5C illustrates the prevalence of tissue-specific differences by comparing the predicted MAE genes in three cell lines. Genes predicted as MAE in only one of the three cell types are typically not expressed in all the three cell types. Strikingly, in each cell type between 1/3 and half of the predicted MAE genes showed cell type-specific expression. (In Figure 5C, compare 1184 predicted MAE genes expressed in HMEC but silent in one or both other cell lines, with 495 genes that are MAE only in HMEC but expressed also in HSMM and HUVEC cells, or with 395 genes that are expressed and predicted MAE in all three cell types.) We interpret these observations to mean that the H3K36me3–H3K27me3 gene body signature is a general feature of monoallelic expression in multiple types of human cells, while a particular set of affected genes can be cell type-specific. It should be noted that the correspondence of the chromatin signature to monoallelic expression has only been systematically validated in clonal LCL lines (see Figure 2 and Figure 3). Analysis in nonclonal cell populations (such as the ENCODE cell lines or isolated ex vivo cells) has been mostly restricted to single-cell techniques such as the fluorescent in situ hybridization approach which we have used in PBMCs (Gimelbrant et al. , 2007). Until a more general study is completed, it remains formally possible that the chromatin signature does not reflect MAE in some or many cell types. The analysis in the rest of this section should be understood as being performed on genes that have a particular chromatin signature, which at least in one cell type strongly correlates with monoallelic expression. The greater the biological similarity between cell types, the greater the overlap in MAE genes (Figure 5D). Overall, cell lines fell into two groups: GM12878, K562 and H1ESC had a larger percentage of transcribed genes with MAE chromatin signature (about 25%), while the rest had a lower percentage (about 15%). A large number of predicted MAE genes in ES cells suggest that the allelic choice may occur fairly early in development. When assessing overlap between predicted MAE genes in different cell types, several patterns became evident. Only 61 genes were predicted as MAE in all of the seven analyzed cell types. Among commonly expressed genes (∼10,000 genes with at least some evidence of expression in all analyzed cell lines), many are predicted as MAE only in some lines, implying tissue-specific MAE. Among ubiquitously expressed genes within this group, addition of more cell lines rapidly led to a plateau in the cumulative number of MAE genes, making up about 20% of the genes in this set (Figure 5E). In contrast, when all genes are considered (including cell type-specific ones), the number of genes expressed in at least one cell line and the cumulative number of predicted MAE genes keep rising, since a large fraction of the cell type-specific genes are predicted as MAE (Figure 5F). With seven cell lines, the total number of genes with evidence of expression in at least one cell line (RPKM>0. 1) was 18,248, among which the total number of genes with the chromatin signature of MAE in at least one cell line was 8716 (48%). Extrapolation of this trend indicates that, with a large number of cell types, about half of all genes would be predicted to have MAE in one or more cell type. Qualitatively similar results remain if a higher expression threshold is applied, such as RPKM>1 (39% at seven cell lines). Taken together, these analyses lead to two unexpected conclusions: (a) among widely expressed genes, about 20% have a chromatin pattern characteristic for MAE in at least one tissue; (b) taking genes with tissue-specific expression as a group, we estimate that between 39% and 48% have the MAE pattern. The genes with MAE chromatin signature share some pronounced features as a group. Among predicted MAE genes, genes coding for cell surface proteins and those involved in multicellular developmental processes were heavily over-represented. Conversely, housekeeping genes and genes specific to intracellular organelles were highly under-represented. When the same gene ontology categories were assessed in other analyzed cells lines, they were similarly over- or under-represented at high significance levels (Figure 5G). Since the overlap in actual genes is relatively modest, this suggests that in different cell types MAE affects different genes involved in the same type of activity. A particularly important group consists of bivalent genes, which are also associated with overlapping active and inactive histone modifications. In embryonic stem cells, bivalent genes are transcriptionally silent, and their promoters carry overlapping peaks of H3K4me3 and H3K27me3 (Mikkelsen et al. , 2007). These genes play a crucial role in determining cell fate during differentiation: upon reaching a point of lineage commitment, this ‘poised’ chromatin allows these genes to rapidly resolve into either an active or inactive state. To our surprise, we found that a remarkably large fraction of the known bivalent genes (>80%) acquire the MAE chromatin signature in at least one of the differentiated cell types, much higher even when compared to other genes that are also silent in ES cells (Figure 6A). Our observations indicate that upon activation, these master regulator genes preferentially resolve into a state with the characteristic MAE chromatin signature. 10. 7554/eLife. 01256. 013Figure 6. MAE chromatin signature in bivalent genes. (A) Overrepresentation of bivalent genes among predicted MAE genes. Predicted: predicted MAE in at least one cell line and not silent (RPKM > 0. 1); not predicted: not predicted MAE in any cell line. Groups of genes: (top) not bivalent, silent in hESC (RPKM < 0. 1); (middle) not bivalent, not silent in hESCs; (bottom) reported bivalent in human ES cells. (B) A speculative model of MAE establishment in bivalent genes. Genes with bivalent/poised chromatin in promoters are silent in undifferentiated stem cells; two alleles have symmetric distribution of active and inactive histone marks. When such gene is activated upon reaching a point of cell fate determination, either one of the alleles (or both) can become transcriptionally active. The initial choice is stochastic, but it is stable in the clonal progeny. Asymmetric histone modifications in the gene body reflect activity of the alleles. DOI: http: //dx. doi. org/10. 7554/eLife. 01256. 013
  Using a systematic machine learning approach, we identified and then experimentally validated a specific and sensitive signature for MAE: a gene-body overlap between chromatin marks for active transcription (H3K36me3) and gene silencing (H3K27me3). Interestingly, the promoter signal carried very little additional information relative to the gene-body features (Figure1—figure supplement 3). This could be due to some combination of biological reasons with much lower noise associated with gene body features due to their much greater length. It should also be noted that we focused on promoters and gene bodies because they are unequivocally associated with particular genes. It might be informative to study the involvement of other regulatory elements, such as enhancers. Upon applying this approach to a variety of cell types, we estimate that up to 20% of ubiquitously expressed genes and more than one-third of tissue-specific genes showed the chromatin signature of MAE. Since detection of the MAE signature does not rely on SNPs or on sample clonality, we expect that this approach will make feasible the analysis of MAE in primary samples, including systematic comparisons of normal and diseased tissue. Our findings fill three major gaps in our understanding of MAE. First, we have positively identified a molecular correlate of MAE with strong predictive power. We mapped MAE signature genome-wide in multiple cell types. Finally, we uncovered an unexpected relationship between MAE and known drivers of differentiation, bivalent genes, implying a functional role of MAE in this process. The overlap of MAE and bivalent genes poses a potential mechanistic puzzle. In the promoters of bivalent genes, both active and inactive modifications are simultaneously present on both alleles (Bernstein et al. , 2006). Moreover, direct evidence showed individual nucleosomes carrying H3K4me3 or H3K36me3 along with H3K27me3, albeit on the opposite histone H3 tails (Voigt et al. , 2012). At the same time, asymmetric distribution of H3K27me3 and H3K36me3 within the gene body of MAE genes is consistent with a very simple asymmetric model: H3K27me3 is known to be associated with transcriptionally silenced chromatin, and H3K36me3 is associated with gene bodies of the actively transcribed genes (Black et al. , 2012). This leads us to speculate that transition of the bivalent genes from transcriptionally silent state with ‘poised’ chromatin into an active state often happens independently for the two alleles. In that scenario, either allele can become stably inactive, with gene body chromatin enriched with H3K27me3, while the other becomes active and enriched with H3K36me3 (Figure 6B). Once the poised state has resolved into active or inactive state, it locks in that state; with some probability, both alleles become active, and the biallelic state is locked. While the molecular details of this process remain unclear, this model makes testable prediction that for bivalent genes, MAE is established during lineage commitment events, such as ES cell differentiation. Our finding of a specific combination of molecular markers associated with autosomal MAE creates opportunities for a deeper understanding of mechanisms underlying this epigenetic phenomenon. The presence of H3K27me3 immediately implies that histone methyltransferase EZH2, and possibly other components of the PRC2 complex (Margueron and Reinberg, 2011), could be involved in MAE maintenance. Our experimental observations suggest that H3K27me3 and H3K36me3 marks are asymmetrically distributed between the alleles of genes with MAE. Our observations are consistent with the idea that a uniform molecular mechanism might be involved in MAE maintenance genome-wide. Our study strongly implies that very similar sets of genes are subject to MAE in different individuals: the presence of the H3K36me3/H3K27me3 chromatin pattern in one individual corresponds to MAE in the same cell type in another individual. A parsimonious hypothesis is that the capacity for MAE is genetically controlled by regulatory DNA sequences, presumably in cis to the affected MAE genes. Consistent with this notion, significant inter-species conservation exists among genes subject to MAE in human and mouse (Zwemer et al. , 2012). Epigenetic regulation of monoallelic expression can profoundly affect the relationship between genotype and molecular phenotype (Pereira et al. , 2003; Zuo et al. , 2007). One intriguing possibility is raised by a recent observation that genomically identical clones can be dramatically different in terms of tumorigenic properties and drug resistance (Kreso et al. , 2012). The knowledge of the chromatin signature of MAE will help in exploring the relationship between MAE and functional cellular states in normal development and disease.
  Alignment of RNA-seq data was done using Bowtie2 (Langmead and Salzberg, 2012), using paired end alignment, seed length = 28, and max seed mismatches = 2. RPKM was calculated using cufflinks v. 2, with a gtf file for Homo sapiens, UCSC, hg19 obtained from http: //cufflinks. cbcb. umd. edu/igenomes. html. For comparison, cufflinks was run using the multiple aligned files option simultaneously on all aligned files (bam format) of the ENCODE cell lines. The HMEC and HCC1954 were subsequently run together using the same option separately from the rest and quantile normalization to GM12878 reference was performed on the RPKM values. In-house analysis pipeline in perl and awk was used to generate an SNP-masked reference and to obtain mapped read counts for each SNP. The reference was derived from the hg19 genome by removing non-transcribed regions using gtf annotation, and masking SNP loci. Long-distance bias effects (McDaniell et al. , 2010) were also removed using in silico sequencing simulation. The reads covering SNP loci were tallied into maternal and paternal ‘hits’. In-house Matlab analysis pipeline was used to calculate binomial p-value with FDR-correction and perform equivalence testing for each gene. Allelic bias was statistically identified from the resulting SNP allelic counts. As a first step, counts for multiple SNPs spanning the same gene were tallied, since SNPs were never less than one read length apart and each hit was therefore considered the result of an independent Bernoulli trial of a random variable representing gene bias, with complete maternal bias corresponding to parameter value 0, lack of bias to value 0. 5, and complete paternal bias to value 1. In processing AST-seq data, which had very high coverage per SNP (generally >103), we used equivalence testing and rejection of equivalence as the test for BAE and MAE with equivalence boundaries of 2: 1 bias (corresponding to the maximum likelihood estimate parameter value of less than 0. 33 or more than 0. 67). Due to the much lower coverage in the whole genome RNA-seq experiment, we used FDR-corrected binomial testing of the tallied counts for each of the heterozygous genes for which hits were obtained, along with a threshold of 2: 1 bias (corresponding to the maximum likelihood estimate parameter value of less than 0. 33 or more than 0. 67) to filter out statistically significant but biologically weak bias. To call unbiased (biallelic) genes in each clone, we used equivalence testing, with equivalence boundaries corresponding to the above 2: 1 bias. While the use of 2: 1 bias as threshold is conventional (Zwemer et al. , 2012; Li et al. , 2012a), the results remained qualitatively robust with other thresholds (e. g. , 3: 1) used. This was due to the fact that a disproportionate fraction of MAE genes showed extreme biases. In the analysis of data from the multiplex capture of common coding SNPs using padlock probes (see below), we restricted the analysis to SNPs with 100 reads or more.
summary: Understanding how genes are activated and silenced is one of the central challenges in modern biology. These processes underpin the development of a fertilized egg into a complex organism, and they can also lead to life-threatening diseases when they go wrong. There are two copies of each gene in a human cell, a maternal copy and a paternal copy, and it is thought that both copies are usually regulated together. However, there are exceptions to this rule: for certain genes only the maternal copy is expressed as a protein in some cells, whereas the paternal copy is expressed in other cells. This form of gene regulation, which is called monoallelic expression, can result in neighboring cells heading down very different paths. In extreme cases, depending on the differences between the two copies of the gene, cells that express one copy may function normally, while cells where the other copy is activated will start forming tumors. However, despite these potentially grave consequences, and early results which suggested that monoallelic expression affected a large number of human and mouse genes, it has proved to be a major technical challenge to identify these genes in most cell types. Now, Nag, Savova et al. have discovered a molecular signature that can be used to detect monoallelic expression. The signature was found in chromatin, the densely packed structure formed by DNA and proteins inside the cell nucleus. Nag, Savova et al. discovered that the genes that are subject to monoallelic expression are bound with proteins that are modified in two contrasting ways. One modification, which is usually a sign of gene silencing, is prevalent on the inactive copy of the gene, and the other, which often marks active genes, is chiefly present on the active copy. Nag, Savova et al. report that these modifications are found in different sets of genes in different cell types, indicating distinct genome-wide patterns of monoallelic expression. The chromatin signature approach lets them estimate the fraction of human genes that are subject to monoallelic expression. This number is surprisingly high: about 20% of commonly expressed genes and more than one-third of tissue-specific genes. In a particularly intriguing finding, almost all bivalent genes—a subset of genes that are involved in determining the fate of cell during development—are estimated to become monoallelic when they are activated. In addition to these unexpected findings, the chromatin signature approach opens the door to exploring monoallelic expression as a form of gene regulation in all types of cells and, ultimately, to understanding how it is involved in both normal development and in disease.
keywords: chromosomes and gene expression
  developmental biology
title: Chromatin signature of widespread monoallelic expression