Local genomic features predict the distinct and overlapping binding patterns of the bHLH-Zip family oncoproteins MITF and MYC-MAX

Abstract

MITF and MYC are well-known oncoproteins and members of the basic helix-loop-helix leucine zipper (bHLH-Zip) family of transcription factors (TFs) recognizing hexamer E-box motifs. MITF and MYC not only share the core binding motif, but are also the two most highly expressed bHLH-Zip transcription factors in melanocytes, raising the possibility that they may compete for the same binding sites in select oncogenic targets. Mechanisms determining the distinct and potentially overlapping binding modes of these critical oncoproteins remain uncharacterized. We introduce computational predictive models using local sequence features, including a boosted convolutional decision tree framework, to distinguish MITF versus MYC-MAX binding sites with up to 80% accuracy genomewide. Select E-box locations that can be bound by both MITF and MYC-MAX form a separate class of MITF binding sites characterized by differential sequence content in the flanking region, diminished interaction with SOX10, higher evolutionary conservation, and less tissue-specific chromatin organization.

Keywords: MITF; MYC-MAX; boosted convolutional decision tree; cobinding factors; machine learning.

Figure 1:

Motifs of (a) MITF, (b) MYC-MAX, and (c) SOX10 inferred from their respective ChIP-seq data (Supplementary Methods).

Figure 2:

(a) Density of SOX10 motifs around MITF-bound and MYC-MAX-bound E-boxes. SOX10 motif shows a strong co-localization 30–150 bps from MITF-bound E-boxes, but does not co-localize with MYC-MAX E-boxes. (b) ChIP-seq read density of SOX10 enrichment around MITF ChIP-seq peaks and MAX ChIP-seq peaks.

Figure 3:

(a) Distribution of average GC nucleotide frequency in 100 bp windows centered at E-boxes bound by both MITF and MYC and those bound by MITF exclusively. (b) Distribution of average CpG di-nucleotide frequency in 100 bp windows around E-boxes bound by both MITF and MYC and those bound by MITF exclusively. (c) Distinguishing characteristics of the overlapping E-box class, including H3K27ac, H3K4me3 histone modifications, gene expression response to siMITF knockdown, SOX10 co-localization, and E-box evolutionary conservation. Binomial test p-values for testing the difference between the two subclasses of E-boxes are indicated (Supplementary Methods).

Figure 4:

(a) Area under the receiver-operating characteristic (ROC) curve for three different models trained to classify between MITF- and MYC-MAX-bound sequences. From left to right: BDT model with bHLH motifs removed; BDT model with full set of motifs; BCDT model. (b) Heatmap of the partial dependence slopes, ordered by slope, of the eight features with greatest importance for the BDT models with bHLH motifs removed (left) and using the full set of motif features (right). Positive (negative) slope values indicate a positive (negative) association between the presence of a particular motif and the model’s prediction that a sequence is bound by MITF (MYC-MAX). The relevant TRANSFAC IDs are as follows: (LEF1:M00805), (YY1:M00793), (TP53:M00761), (KLF12:M00468), (E2F-1:M00428), (TFE:M01029), (SREBP-1:M00220), (MYC: M00799).

Figure 5:

(a) Density plot showing a positive correlation (Pearson r = 0.70) between percent GC content and the probability of a sequence being MYC-MAX-bound according to our BCDT model trained to classify between MITF- and MYC-MAX-bound sequences. (b) Normalized histogram of GC content percentage for MITF-bound, MYC-MAX-bound, and random DNase I hypersensitive sequences, demonstrating relative GC enrichment in MYC-MAX-bound sequences.