Validation of skeletal muscle cis-regulatory module predictions reveals nucleotide composition bias in functional enhancers

Abstract

We performed a genome-wide scan for muscle-specific cis-regulatory modules (CRMs) using three computational prediction programs. Based on the predictions, 339 candidate CRMs were tested in cell culture with NIH3T3 fibroblasts and C2C12 myoblasts for capacity to direct selective reporter gene expression to differentiated C2C12 myotubes. A subset of 19 CRMs validated as functional in the assay. The rate of predictive success reveals striking limitations of computational regulatory sequence analysis methods for CRM discovery. Motif-based methods performed no better than predictions based only on sequence conservation. Analysis of the properties of the functional sequences relative to inactive sequences identifies nucleotide sequence composition can be an important characteristic to incorporate in future methods for improved predictive specificity. Muscle-related TFBSs predicted within the functional sequences display greater sequence conservation than non-TFBS flanking regions. Comparison with recent MyoD and histone modification ChIP-Seq data supports the validity of the functional regions.

Figure 1. Selection of clones for differential expression analysis.

The selection is divided into 2 phases, where the clones selected for Phase 2 are a subset of all clones tested in Phase 1. Phase 1 and Phase 2 samples are from different plasmid preparations.

Figure 2. Dinucleotide frequencies in responding regions vs. non-responding regions.

a) Muscle Regulatory Regions. Muscle Validated = 19 validated muscle regions in this study. Muscle Reference = 28 muscle reference regions from literature. Muscle Non-Responders = all regions that were tested in this study and found not to drive gene expression. b) Pleiades Curated Regulatory Regions. Pleiades Curated All = 1341 curated regulatory regions from all species. Pleiades Curated Human = 631 curated regulatory regions in humans only. Non-Responders = all regions that were tested and found not to drive gene expression.

Figure 3. Examples of positive regions.

Muscle TFBS hits (threshold of 80%) and the phastCons conservation profile for the region are shown as well. When square brackets are shown, they indicate the original CRM prediction. a) Positive sequence from the muscle set. The muscle-specific TFBS are located in regions of high sequence conservation. b) Positive sequence from the non-muscle. This sequence showed the most consistent increase in reporter expression, with all 12 replicates determined as significantly up-regulated in muscle. c) Positive sequence from the background set. Despite the clear cluster of muscle-specific TFBS located in the region of high sequence conservation, none of the CRM prediction tools could classify this as a muscle CRM.

Figure 4. Phylogenetic depth analysis of TFBSs in responding and non-responding regions using phyloP (46-way, hg19).

The x-axis grouping indicates the species used to calculate the phyloP scores. TFBSs were searched using all vertebrate profiles from the JASPAR CORE collection using the threshold of 0.8. TFs with at least 2-fold increase in phyloP scores (46wayAll) in the TFBS positions over the non-TFBS positions in the responding regions were identified. The score ratios for these TFs were compared among the three region sets. a) Average phyloP scores for the predicted TFBS positions and non-TFBS positions in each region set. b) Ratios of phyloP scores for the TFBS positions and non-TFBS positions in each region set.

Figure 5. Histone modifications in the responding and non-responding regions.

Proportion of the regions that overlap with ChIP-Seq peaks from Asp et al. are displayed. (MB = Myoblasts, MT = Myotubes).