Functional characterization of motif sequences under purifying selection.

3.3 MB·English

by  Chen, De-Hua

Published online 8 January 2013 Nucleic Acids Research, 2013, Vol. 41, No. 4 2105–2120
doi:10.1093/nar/gks1456
Functional characterization of motif sequences
under purifying selection
De-Hua Chen1, Andrew Ying-Fei Chang2, Ben-Yang Liao2 and Chen-Hsiang Yeang1,*
1Institute of Statistical Science, Academia Sinica, Taipei, Taiwan, ROC and 2Division of Biostatistics and
Bioinformatics, Institute of Population Health Sciences, National Health Research Institutes, Zhunan, Miaoli
County, Taiwan, ROC
Received September 18, 2012; Revised and Accepted December 13, 2012
ABSTRACT promoters retain significantly coherent expression
profiles, and those genes are over-represented in
Diverse lifeforms are drivenbytheevolutionofgene
the functional classes involved in gene regulation.
regulatory programs including changes in regulator
The validation results reveal the dependencies
proteins and cis-regulatory elements. Alterations of
between natural selection and functions of cis-regu-
cis-regulatory elements are likely to dominate the
latory elements and shed light on the evolution of
evolution of the gene regulatory networks, as they
gene regulatory networks.
are subjected to smaller selective constraints
compared with proteins and hence may evolve
quickly to adapt the environment. Prior studies on INTRODUCTION
cis-regulatory element evolution focus primarily on
Diverse life forms are largely driven by conservation and
sequence substitutions of known transcription variations of the gene regulatory circuits. Recent progress
factor-binding motifs. However, evolutionary models in high-throughput technologies such as next-generation
for the dynamics of motif occurrence are relatively sequencingplatformsandDNAmicroarraysenablesbiolo-
rare,andcomprehensivecharacterizationoftheevo- gists to map the regulatory networks and investigate their
lution of all possible motif sequences has not been evolution across multiple species. For instance, studies in
pursued. In the present study, we propose an algo- evolutionary developmental biology (EvoDevo) compared
the gene regulatory networks for animal development and
rithmtoestimatethestrengthofpurifyingselectionof
discovered conserved cores responsible for body plan for-
a motif sequence based on an evolutionary model
mation and variable modules modifying species-specific
capturing the birth and death of motif occurrences
phenotypessuchastheshapesoflimbsorwings[e.g.,(1,2)].
on promoters. We term this measure as the ‘evolu-
One remarkable feature from the gene regulatory
tionary retention coefficient’, as it is related yet
networks of multiple species is the conservation of their
distinctfromthe canonical definitionof selectionco- constituent proteins (3). Most proteins possess multiple
efficient in population genetics. Using this algorithm, functions(pleiotropic),hencearesubjectedtotightselective
weestimateandreporttheevolutionaryretentionco- constraints.Alterationsonproteinsequences(e.g.,changes
efficients of all possible 10-nucleotide sequences ontheDNA-bindingdomainofatranscriptionfactor)may
from the aligned promoter sequences of 27748. affect many partners (e.g., changes on the bindings of all
orthologous gene families in 34 mammalian species. targetsofatranscriptionfactor),thusarelikelytobedele-
terious. In contrast, alterations on cis-regulatory elements
Intriguingly, the evolutionary retention coefficients
have local effects and thus enable the systems to evolve in
of motifs are intimately associated with their func-
an incremental fashion. Consequently, evolution of non
tionalrelevance.Top-rankingmotifs(sortedbyevolu-
protein-coding regions in general and cis-regulatory
tionary retention coefficients) are significantly
elements in particular plays a critical role in the evolution
enrichedwithtranscriptionfactor-binding sequences
of the gene regulatory systems.
according to the curated knowledge from the Early studies of cis-regulatory element evolution focus
TRANSFAC database and the ChIP-seq data onidentificationofconservedtranscriptionfactor-binding
generated from the ENCODE Consortium. Moreover, motifs(4)anddetectionofconservedregionsongenepro-
genes harbouring high-scoring motifs on their moters (5). Sequence conservation alone, however, does
*To whom correspondence should be addressed. Tel:[phone]; Fax:[phone]; Email: [email protected]
(cid:2)TheAuthor(s)2013.PublishedbyOxfordUniversityPress.
ThisisanOpenAccessarticledistributedunderthetermsoftheCreativeCommonsAttributionNon-CommercialLicense(http://creativecommons.org/licenses/
by-nc/3.0/),whichpermitsunrestrictednon-commercialuse,distribution,andreproductioninanymedium,providedtheoriginalworkisproperlycited.
2106 NucleicAcidsResearch,2013,Vol.41,No.4
notsufficetoaccountfortheevolutionofgeneregulatory are highly enriched with the processes of transcriptional
systems.Comparisonofknowncis-regulatoryelementson regulation.Theresultsprovideacomprehensivepictureof
closelyrelatedspeciesindicateshighratesofturnoverand the evolution of cis-regulatory elements.
divergence (6–10). These changes may yield gains and
losses of cis-regulatory elements (19), modify the regula-
toryprograms(1,2)orareaccompanied bycompensatory MATERIALS AND METHODS
mutationstomaintainstableregulatoryprograms(11,12).
Data sources
Likeprotein-coding regions, evolutionofcis-regulatory
elements is driven by a variety of mechanisms including Aligned 5kb upstream sequences of 27748 orthologous
sequence substitutions (18), gene duplications (13), gene families from 34 mammalian species were extracted
tandem repeat insertions and deletions (14). Cis-regula- from the UCSC Genome Browser (21). Supplementary
tory elements are added or deleted on the promoters/en- Table S1 and Figure S1 report the names and the phylo-
hancers according to these mechanisms. Ideally, a genetic tree of the selected species.
complete model for the evolution of cis-regulatory Tovalidate thefunctionalrelevance ofthehigh-scoring
elements should be based on the models of all individual motifs, we downloaded external datasets from the follow-
mechanismsformolecularevolution.Inpractice,mechan- ing sources: the consensus motifs of transcription factor-
ismsotherthansequencesubstitutionsarehardtoformal- binding sequences from the TRANSFAC database (22),
ize.Consequently,themajorityofquantitativemodelsfor 407ChIP-seqdatafilesfromtheENCODEdatabase(23),
cis-regulatory element evolution are derived from DNA microarray data of human and mouse tissue gene
sequence substitution processes. Several studies use simu- expressions (24) and RNA-seq data of human tissue gene
lationstoexaminetheeffectsofsequencemutationsonthe expressions (25), the annotations and member genes of
rates for cis-regulatory element evolution [e.g., (15,16)]. 3201 Gene Ontology (GO) categories (27) and pathway
Others startwithsequence substitution modelsin popula- information from three databases (28–30).
tion genetics and attempt to identify the cis-regulatory
elements under selection [e.g., (17–19)]. Despite the Quantifying the strength of natural selection
fruitful outcomes generated from these studies, they of motif sequences
suffer from two major limitations. First, they focus pri-
We define a motif as a collection of sequences with the
marily on the deviation of observed sequences from a
same length. Over time motifs are created, annihilated or
known regulatory element (e.g., a transcription factor-
maintainedinaspecifiedregion(e.g.,agenepromoter)by
binding motif) rather than the changes of regulatory
sequence substitutions of the constituting nucleotides.
element occurrence on promoters. Alterations on motif
Motifs undergoing purifying selection would possess
counts can be more critical for gene regulation than
slower rates of annihilation than those without selective
specific sequence variations, as the former modulate the
constraints. Accordingly, we quantify the strength of
number of transcription factors bound on promoters.
natural selection of a motif by comparing the empirical
Second, all the current studies only examine a collection
distribution of its occurrences over multiple species with
of known transcription factor-binding motifs. Complete
the one generated by a neutral evolutionary model. The
characterization of the evolution of all possible motif se-
evolutionary model of motif occurrences and the algo-
quencesofafixedlengthislacking.Thischaracterization,
rithm of evaluating the evolutionary retention coefficients
however, is critical for discovering new regulatory
of motifs are described below.
elements and comprehending their evolution on genomes.
Previously, we proposed an evolutionary model and an A Poisson process model of sequence substitution
algorithmtoquantifythestrengthofnaturalselectionofa We adopt the simplest model of sequence substitution
motif sequence (20). The evolution of motif occurrence assuming all nucleotides at all positions and across all
was formulated as a birth–death process, whereas the lineagestransitionwithanequalrate(31).Inaninfinitesi-
rates of motif additions and deletions were derived from maltimeintervaldt,thenucleotidesequenceofaposition
substitutions of their constituent sequences. The evolu- transitions to another base with probability (cid:2)dt. n ðtÞ
s
tionary retention coefficient of a motif was defined as a denotes the cumulative number of sequence changes at
penalty to slow down motif death, and the evolutionary time t. The transitions within the time interval [t, t+dt]
retention coefficient value maximizing the log likelihood is as follows
ofthedatawasestimated.Inthepresentstudy,weextend
this model and evaluate the evolutionary retention coeffi- Pðnsðt+dtÞ¼ðN+1ÞjnsðtÞ¼NÞ ¼(cid:2)dt: ð1Þ
cientsofallthe410=104857610-nucleotidesequenceson Pðnsðt+dtÞ¼NjnsðtÞ¼NÞ ¼1(cid:2)(cid:2)dt:
thepromotersof27748orthologousgenefamiliesfrom34
and n ðtÞ has a Poisson distribution
mammalian species. Intriguingly, evolutionary retention s
coefficients of the 10-mer sequences are significantly ð(cid:2)tÞN
associated with the tendency of transcription factor- PðnsðtÞ¼Njnsð0Þ¼0Þ¼ N! e(cid:2)(cid:2)t: ð2Þ
binding events and expression coherence of the genes
harbouring the motifs. By examining the annotations of The maximum likelihood estimate of (cid:2) is simply the total
the top-ranking motifs, we find many of them match the numberofsequencechangesdividedbythetotallengthof
GC-rich binding sequences of the transcription factors. the time interval considered. In this work we estimated (cid:2)
Furthermore, genes harbouring the top-ranking motifs from the aligned 5kb promoter sequences of the 27748
NucleicAcidsResearch,2013,Vol.41,No.4 2107
gene families over the 34 mammalian species. For each We now extend the analysis to the entire promoter of
position of the aligned promoters in each gene family, length l (cid:5)l . Suppose motif occurrence at time t is
s m
we observed the sequences in the terminal nodes (the 34 nðtÞ¼n and the n occurring motifs are not overlapped.
extant species) of the species tree and inferred the se- Each motif instantiation can be annihilated with a rate
quences in the internal nodes (ancestral species) by a (cid:2)l r . Hence the ‘death rate’ on the entire promoter is
m 10
dynamic programming algorithm (32). We then counted the rate on an l -mer window multiplied by n:
m
the total number of sequence changes along all branches
Pðnðt+dtÞ¼n(cid:2)1jnðtÞ¼nÞ¼(cid:2)l r ndt: ð5Þ
of the species tree for all positions and all gene families m 10
andthetotallengthsofthetimeintervals,andcalculated(cid:2) There are l (cid:2)l n positions unoccupied by motif se-
s m
accordingly. From the empirical data, (cid:2)=0.8371. quences, and the maximum number of (possibly
overlapped) l -mer windows is l (cid:2)l n(cid:2)l +1. Each of
A birth–death model for the evolution of motif occurrences m s m m
these l -mer windows can generate a new motif. Hence
A motif M(cid:3)BlmðB(cid:4)fA,G,TgÞ is defined as a collection m
the ‘birth rate’ on the entire promoter is approximately
ofnucleotidesequencesoffixedlengthl .Wefirstconsider
m the rate on an l -mer window multiplied by
the sequence evolution of l consecutive positions. There m
m l (cid:2)l n(cid:2)l +1:
are 4lm possible sequences that can occur in this lm-mer s m m
window, and each sequence s2Blm can be labelled as Pðnðt+dtÞ¼n+1jnðtÞ¼nÞ¼(cid:2)l r ðl (cid:2)l n(cid:2)l +1Þdt:
m 01 s m m
either a member of the motif ðs2MÞ or not ðs2=MÞ.
ð6Þ
These sequences comprise an undirected graph G=(V,
E), where a node (cid:3)2V denotes an l -mer sequence and Equations(6)and(5)specifyabirth–deathprocess(34)of
m
an edge e=(v , v ) denotes a sequence pair v and v dif- motifoccurrenceonapromoteroflengthl .Thedistribu-
1 2 1 2 s
ferent at one position. The evolution of l -mer sequences tion P ðtÞ(cid:4)PðnðtÞ¼nÞ of motif occurrences over time
m n
can be viewed as a Markov random walk on G. In an in- can be expressed as a system of differential-difference
finitesimaltimeinterval,asequencecanonlytransitiontoa equations:
neighboring node in G and the rate of transition is (cid:2)l .
A motif M constitutes a subset of nodes in G (bmlack dPd0tðtÞ¼(cid:5)ð1ÞP1ðtÞ(cid:2)(cid:2)ð0ÞP0ðtÞ:
nodesintheleftdiagramofFigure1),whiletheremaining dPnðtÞ¼(cid:2)ðn(cid:2)1ÞP ðtÞ+(cid:5)ðn+1ÞP ðtÞ(cid:2)ð(cid:2)ðnÞ+(cid:5)ðnÞÞP ðtÞ:
dt n(cid:2)1 n+1 n
nodes are non-motif sequences (white nodes in the left (cid:2)ðnÞ¼(cid:2)l r ðl (cid:2)l n(cid:2)l +1Þ:
m 01 s m m
diagram of Figure 1). We are interested in the transition (cid:5)ðnÞ¼(cid:2)l r n:
m 10
ratefromnon-motifsequencestomotifsequencesandvice
ð7Þ
versa. With a simplifying approximation, we characterize
these transitions with two numbers: r as the fraction of ThesystemisillustratedbytherightdiagramofFigure1.
01
all non-motif ! motif transitions among all transitions The aforementioned model assumes that sequences
from non-motifs, and r as the fraction of all motif ! randomly drift and henceforth no selective pressure is
10
non-motif transitions among all transitions from motifs. exerted on the evolution of motif occurrence. In
P , P , P and P denote the background frequencies of contrast, purifying selection should penalize decrements
A C G T
thefournucleotidesobtainedfromallpromotersofthe34 of motif occurrence. Therefore, the evolutionary model
species. r and r are calculated by the following of motif occurrence under purifying selection largely re-
01 10
formulas: sembles the model for neutral evolution (Equation 7)
P except for a modification of the motif death rate:
r01 ¼PfPðv1,vf2ðÞv21E,v:v21Þ2=2ME:vg12=!Mðvg1,!vð2vÞ1(cid:4),ðvv22Þ2MÞ: ð3Þ (cid:5)0ðnÞ¼(cid:5)ðsnÞ: ð8Þ
r10 ¼ fPðv1,v2Þ2E:v12Mg!ðv1,v2Þ(cid:4)ðv22=MÞ: Themotifdeathrate(cid:5)0ðnÞunderselectionslowsdownthe
fðv1,v2Þ2E:v12Mg!ðv1,v2Þ neutral motif death rate (cid:5)ðnÞ by a factor s>1. We term s
·
Where d( ) is an indicator function and o((cid:3) , (cid:3) ) is the as the evolutionary retention coefficient of a motif.
1 2
nucleotide background probability of (cid:3) at the position Notably, this definition is related yet distinct from the
2
where (cid:3) and (cid:3) differ. For instance, o(AGGC, canonical definition of selection coefficient in population
1 2
AGTC)=P .o((cid:3) ,(cid:3) )’srescaletheweightsoftransitions genetics (35). In population genetics, the selection coeffi-
T 1 2
according to the frequencies of the destination sequences. cient is the decline of the relative fitness of a selectively
For instance, transitions to GC-rich sequences are more disadvantageousgenotypecomparedwiththatofaselect-
likely to occur on mammalian promoters, as they are ivelyfavouredgenotype.Inasufficientlylargepopulation,
over-represented in the CpG islands (33). the selectively advantageous genotype will appear with a
n(t)denotesthenumberofmotifoccurrenceattimet.In higher frequency than that of a genotype without selec-
an l -mer window, nðtÞ2f0,1g as the sequence is either a tion.Inthisregard,boththecanonicalselectioncoefficient
m
motif or not. The transitions of nðtÞ hence conform with and evolutionary retention coefficient aim for capturing
the following equations thestrengthofpurifyingselectionfromtheobservedgeno-
types. However, despite the common goals shared by the
Pðnðt+1Þ¼1jnðtÞ¼0Þ¼(cid:2)l r dt:
m 01 ð4Þ two measures, the evolutionary retention coefficient is
Pðnðt+dtÞ¼0jnðtÞ¼1Þ¼(cid:2)lmr10dt: distinct from the canonical selection coefficient in two
2108 NucleicAcidsResearch,2013,Vol.41,No.4
Figure 1. Left:Asequencespaceoffixedlengthasagraph.Anodedenotesasequence,andanedgedenotestwosequencesdifferingatoneposition.
Black nodes are members of a motif and white nodes are non-motifs. Dotted edges denote transitions between motifs and non-motifs. Solid edges
denote transitions within motifs and non-motifs. Right: The state transition diagram of a birth–death model. State n denotes the count of motif
occurrence on a promoter. lðnÞ and mðnÞ denote the birth and death rates emanating from state n.
aspects. First, the evolutionary retention coefficient differential equations. To estimate s that maximizes the
bypasses the abstract notion of relative fitness and log likelihood in Equation (9), we used a binary search
directly tackles the consequences of purifying selection— to find the optimal s over the interval ½0,20(cid:7).
elevation of motif occurrence frequencies. Second, the ca-
nonicalselectioncoefficientexaminestheallelefrequencies
Comparison of selection coefficients and evolutionary
of a single site, whereas the evolutionary retention coeffi-
retention coefficients in simulated data
cient is inferred from the frequencies of motif occurrence
in a consecutive region of the genome. We will further Toelucidate therelation betweenselectioncoefficientsand
clarifytherelationbetweenthesetwoscoresinsimulation evolutionary retention coefficients, we simulated haploid
studies. sequence evolution with varying selection coefficients and
compared the evolutionary retention coefficients estimated
Estimating the evolutionary retention coefficients fromtheobserveddatawiththegivenselectioncoefficients.
from empirical data Givenapromotersequenceoffixedlength(30nucleotides)
We estimated the evolutionary retention coefficient of a andamotif(10nucleotides),wedefinetherelativefitnessof
lm-mermotiffromthealigned5kbpromotersequencesof the promoter as f(cid:4)1(cid:2)ð2(cid:2)minðk,2ÞÞ(cid:6), where k is the
34 mammalian species with the following procedures. number of motif occurrence on the promoter and (cid:6) the
First, we divided a promoter into multiple segments of selectioncoefficient.Thesequencecontaining(cid:8)2motifin-
fixed length ls ¼30l. Segments with >10% gaps in any stancespossesses the highestfitness,whereas the sequences
species were discarded. This partition reduces the containing1and0motifinstancepossessintermediateand
number of valid terms in Equation (7), hence greatly low fitness, respectively.
simplifies subsequent estimation. We simulated promoter evolution according to both
Second, we treated humans as the reference species and sequence substitutions of individual positions and purify-
assumedthatalterationsofmotifcountsfromthereference ing selection dictated by motif occurrence. One promoter
to another species followed the birth–death process. t sequence was randomly generated in the first generation.
denotes the distance between humans and another species In each of the following generations, each sequence
xinthephylogenetictree,n0 andn1 themotifcountsinthe produced 10 progenies with a Poisson mutation rate of
segments of humans and species x, respectively. For each 0.02 per position. Among the progenies from the same
combinationoft(orspeciesx),n0 andn1,wethencounted cohort, 100 of them were selected with probabilities pro-
fðt,n0,n1Þ,thetotalnumberofsegmentswithn0andn1motif portional to their relative fitness, and the remaining se-
instances in the counterparts of humans and species x. quences were eliminated. This process of sequence
Third, the log likelihood of motif occurrences can be substitutions and selection lasted for 100 generations.
expressed as For each selection coefficient, we generated randomly 20
XXX motif and initial promoter sequences and simulated their
L¼ fðt,n n ÞlogPðnðtÞ¼n jnð0Þ¼n Þ+C: ð9Þ
0 1 1 0 evolution separately. Finally, we repeated simulations for
t n0 n1 the following selection coefficient values: 0.01, 0.05, 0.1,
where C is a constant and PðnðtÞ¼n jnð0Þ¼n Þ denotes 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.49.
1 0
the conditional probability derived from the birth–death In each simulation, the coalescent tree of the 100
model under selection by 6. applying the death rate of observed sequences was recorded. We chose the
Equation (8) in Equation (7). Given the relatively short sequence closest to other observed sequences as the refer-
length of segments, we only considered motif occurrences ence and evaluated their phylogenetic distances according
up to 3 and restricted the terms in Equation (7) to n(cid:6)3 tothestructuresandbranchlengthsofthecoalescenttree.
accordingly. For each fixed value of evolutionary reten- The evolutionary retention coefficient of each designated
tion coefficient s, we solved PðnðtÞ¼n jnð0Þ¼n Þ motif can be consequently inferred from the simulated
1 0
numerically by the finite difference method for ordinary sequences.
NucleicAcidsResearch,2013,Vol.41,No.4 2109
An exhaustive evaluation of evolutionary retention ENCODE consortium (23). 407 ChIP-seq data files were
coefficients of all 10-mers on mammalian promoters extractedfromtheENCODEwebsite(http://hgdownload.
cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeHai
Using the aforementioned algorithm, we estimated the
bTfbs/). Each filereports the sequences offragments con-
evolutionary retention coefficients of all 410 ¼1048576
taining the binding sites of one protein in one cell type.
10-mer sequences on the aligned 5kb promoter sequences
For each motif, we constructed a simple null model
of34mammalianspecies.Toacceleratecomputations,we
assuming its occurrence on the entire human genome
ran the estimation procedures on two PC cluster systems
followed a Poisson process. The rate of motif occurrence
simultaneously.Eightjobswereassignedinparalleltothe
per position is N, where N denotes the number of motif
HP DL360 G7 servers containing dual Intel Xeon E5520 L
occurrenceintheentiregenomeandLdenotesthegenome
CPUs with 2.27GHz and 24GB main memory, and 10
length.SupposeinaChIP-seqfile,thetotallengthoffrag-
jobs were assigned in parallel to the HP BL460C servers
ments is l and the number of motif occurrence is n. Then
containing Intel Xeon CPUs with 3.16GHz and 16GB
the Poisson rate of motif occurrence in the designated
main memory. The total running time was 6048 hours.
fragmentsis(cid:7)¼NlandtheP-valueformotifenrichmentis
Notably, although the theoretical framework we present L
can handle more general motifs (i.e., a collection of nu- XN (cid:7)m
cleotide sequences), in this work we only investigate the P¼ e(cid:2)(cid:7): ð11Þ
m!
single sequence motifs (i.e., occurrences of a particular m¼n
10-mer sequence), as their evolutionary retention coeffi-
cients can be exhaustively calculated with limited
Coherence of expression profiles in human
computing resources.
and mouse genes
Third, we showed that human and mouse genes
Functional validation of motif sequences under harbouring high-scoring motif sequences tended to have
selective pressure coherent expression profiles compared with genes
harbouring low-scoring motif sequences. We used both
We incurred the following four tests to validate the
oligonucleotide microarray data (24) and RNA-seq data
functional relevance of motif sequences under selection.
(25) in defining expression levels and co-expression of
Enrichment of TRANSFAC motifs human and mouse genes. For the microarray data, we
First, we demonstrated that evolutionary retention coeffi- obtained the expression information of human genes and
cients were correlated with enrichment of transcription mouse genes from the Gene Atlas V2 dataset (http://
factor-binding motifs extracted from TRANSFAC (22). symatlas.gnf.org/SymAtlas/). This dataset comprises
A non-parametric statistical test was used to evaluate oligonucleotide microarray data in 63 human and 58
the enrichment of transcription factor-binding motifs in mouse normal tissues sampled from animal bodies. We
high-scoring sequences. In brief, all the 104857610-mer assigned the expression data from probe sets to corres-
sequences were sorted by their evolutionary retention co- ponding Ensembl genes following (37,38). The expression
efficientsinadescendingorder.168397ofthesesequences levels of a gene in a specific tissue were averaged among
matchedcompletelyorpartiallywithtranscriptionfactor- replicates. For the RNA-seq data, we obtained that of 11
bindingmotifsinTRANSFAC.WedefinedF ðxÞoverthe human tissues from GEO Series GSE13652 from the
1
normalized rank x(cid:4)rank2½0,1(cid:7) resembling the cumula- University of Toronto (25) (brain/liver/muscle/cerebral
410
tive distribution function (CDF) of TRANSFAC motifs cortex) and GSE12946 from MIT (26) (adipose/breast/
over the sorted 10-mer sequences. colon/heart/lung/lymph node/testes). The raw 32-mer
RNA-seq sequence reads were mapped to human
# (cid:2)TRANSFAC motifs in sequences1!½x(cid:9)410(cid:7)(cid:3) genome (Ensembl version v56), and RNA-seq-based
F ðxÞ¼
1 # TRANSFAC motifs in sequences gene expression levels were calculated according to
(39,40).
ð10Þ
The expression profile divergence between two genes in
F ðxÞ should have a high area under the curve if thehumanormousegenomewasdefinedby1(cid:2)r,wherer
1
TRANSFAC motifs are enriched in the top-ranking se- isthePearson’scorrelationcoefficient ofexpression levels
quences. In contrast, the null hypothesis stipulates that across the tissues. In the present study, we specifically
TRANSFAC motifs are evenly distributed along the examined co-expression of genes that are not paralogs or
sorted sequences and the corresponding CDF is genes located on different chromosomes. The chromo-
F ðxÞ¼x. The maximum deviation between F ðxÞ and somal coordinates and annotations of paralogous rela-
0 1
F ðxÞ gives rise to a statistic of the Kolmogorov– tionships of human and mouse genes based on Ensembl
0
Smirnov test. This method is similar to the Gene Set v62 were obtained through BioMart (http://www.biomart
Enrichment Analysis (GSEA) (36). .org/).
Enrichment of protein-binding sites from the Functional enrichment of genes harbouring
ENCODE data high-scoring motifs
Second,wedemonstratedthatthetop-rankingmotifswere Fourth, we showed that human genes harbouring
enrichedintheprotein-bindingsitesofthehumangenome high-scoring motif sequences were enriched with certain
reported from the ChIP-seq data generated by the functional classes. Four sources pertaining to functional
2110 NucleicAcidsResearch,2013,Vol.41,No.4
information of human genes were extracted: the GO instancesonhumanpromoters.Thespecies(indicesinthe
categories (27), the curated pathway databases of horizontal axis) are sorted by their phylogenetic distances
Reactome (28), Biocarta (29) and the NCI-Nature cur- to humans. All (both high-scoring and control) motifs
ations (30). For a given motif, we extracted the genes havehighlevelconservationbetweenhumansandanthro-
harbouring the motif on their promoters and calculated poid primates (chimpanzees, gorillas, orangutans,
hyper-geometric P-values of enrichment for each GO macaques, indices 2–5) and marmosets (Callithrix
category and pathway. The enriched functional classes jacchus, index 6), and the conditional probabilities drop
for both top-ranking motifs (evolutionary retention coef- abruptly beyond marmosets. For instance, the median
ficient (cid:8)3:0, 231 motifs) and control motifs (231 motifs conditional probabilities between humans and chimpan-
surrounding the median of the sorted list) were reported. zees(index2),orangutans(index4)andmarmosets(index
6) are 0.861, 0.697 and 0.320 respectively, whereas the
median conditional probability between humans and the
RESULTS
Philippine tarsiers (Tarsius syrichta, index 7) drops below
Evolutionary retention coefficients are correlated with 0.074. However, the top-ranking motifs retain consider-
selection coefficients in simulated data ably higher level conservation than control motifs in all
selected mammals. For instance, the median conditional
We first demonstrate the resemblance between canonical
probabilities of the top-ranking motifs between humans
selection coefficients and motif evolutionary retention co-
and guinea pigs (Cavia porcellus, index 14), horses
efficients with simulated data. For each of 11
(Equus caballus, index 21) and African elephants
pre-determined selection coefficient values, we simulated
(Loxodonta africana, index 28) are 0.074, 0.138 and
the evolution of 20 phylogenies, where each phylogeny
0.070 respectively, whereas those of the control motifs
constituted 100 generations and 100 observed promoter
are 0.030, 0.072 and 0.031 respectively.
sequences in their leaves. We then inferred the motif evo-
Notably, in Figure 3 all but one of the top 23 motifs
lutionaryretentioncoefficientfromtheobservedpromoter
have relatively low level conservation compared with the
sequences of each simulated phylogeny. The left part of
remaining top-ranking motifs. By examining those se-
Figure 2 shows the (pre-determined) canonical selection
quences (Table 1), we found they all belonged to the
coefficientsandthe(inferred)motifevolutionaryretention
Alu-J repeat elements (41). They exhibit background
coefficients of the 220 phylogeny instances. Overall, the
level conservation between humans and most other
two scores exhibit a positive correlation coefficient
species but undergo a massive number of insertions on
(r=0.678). Because direct evaluation of relative fitness
the promoters of gray mouse lemurs (Microcebus
in a population is often challenging, the high correlation
murinus, index 8) and small-eared galagos (Otolemur
between the two quantities indicates that the motif reten-
garnettii, index 9) (Supplementary Table S3). These inser-
tion coefficient is a reasonable measure for the selective
tions violate the Poisson process model of sequence sub-
strength of motifs.
stitutions, increase the counts of fðt,n ,n Þ where n >0,
0 1 1
and therefore elevate the evolutionary retention
Summary of evolutionary retention coefficients of
10-mer motif sequences coefficients.
Weevaluatedtheevolutionaryretentioncoefficientsofall
High-scoring motifs are enriched with transcription
410 ¼1048576 10-mer sequences among the 5kb pro-
factor-binding sites
moters of 27748 gene families in 34 mammalian species.
TherightpartofFigure2showstheempiricaldistribution Transcription factor-binding sites likely accommodate
of evolutionary retention coefficients, and Supplementary some motif sequences under selective pressure. We con-
Table S2 reports the sorted evolutionary retention coeffi- firmed the dependency of transcription factor-binding
cients of these sequences. As expected, the majority of sites and evolutionary retention coefficients with two
sequences possess low evolutionary retention coefficients: external datasets. First, we verified that sequences with
the median value is 0.508 and the scores of about 80% of higher evolutionary retention coefficients tended to
the sequences (834552 of 1048576) are below 1.0. We match the transcription factor-binding motifs reported in
considered the first 231 (0.022%) sequences with evolu- the TRANSFAC database (22). A simple check on se-
tionary retention coefficients (cid:8)3:0 as the top-ranking quences sorted by evolutionary retention coefficients
motif sequences and employed further analysis to these provides obvious evidence: 77 of the top 231 10-mer se-
sequences. quences match TRANSFAC motifs, whereas only 35 of
Motifs with high evolutionary retention coefficients the 231 sequences in the middle and 22 of the 231 se-
have slower death rates, thus should exhibit high level quences in the bottom of the sorted list match
conservation on promoters. For each motif, we define a TRANSFAC motifs. Supplementary Table S4 reports
conservationmeasureastheprobabilityofitspresenceon the TRANSFAC match on top-ranking, middle and
the promoter of a mammalian species, conditioned on its bottom control motifs.
presence on the orthologous promoter of humans. In addition to observations on small subsets of se-
Figure 3 displays the conditional probabilities of motif quences, we also quantified this dependency on the
presence of the 231 top-ranking motifs (top panel) and entire sorted list. Denote X a random variable indicating
231 control motifs (bottom panel) with evolutionary re- the match of sorted sequences with TRANSFAC motifs.
tention coefficients near the global median and with (cid:8)50 PðX¼xÞ indicates the probability that a sequence with
NucleicAcidsResearch,2013,Vol.41,No.4 2111
Figure 2. Left:Thescatterplotofcanonicalselectioncoefficientsandevolutionaryretentioncoefficientsonsimulateddata.Eachpointdenotesthe
scores obtained from 100 simulated sequences derived from one common ancestor over 100 generations. Right: Empirical distribution of selection
coefficients among the 410¼1048576 10-mer sequences. The probabilities are displayed in a log scale.
50
100
s)
otif
m
ol 150
ntr
o
c
2:
6
4
−200
2
3
2
s,
otif
m
p 250
o
1: t
3
2
−
1
x (300
e
d
n
otif i
m
350
400
450
5 10 15 20 25 30
species index
Homo sapiens Callithrix jacchus Mus musculus Oryctolagus cuniculus Equus caballus Erinaceus europaeus Dasypus novemcinctus
Figure 3. Conservationofmotifoccurrencebetweenhumansandanotherspecies[P(motifoccursinaspeciesjmotifoccursinhumans)]forthetop
231motifsand231controlmotifsfromthemiddleoftherankedlist.Thehorizontalaxisdenotesthespeciesindexwithanincreasingdistancefrom
humans(sameasthespeciesorderinSupplementaryTableS1).Theverticalaxisdenotesthemotifindexfromhighselectioncoefficients(top)tolow
selection coefficients (bottom). The top-ranking and control motifs are separated by a white line. Colours in the heat map denote the levels of
conditional probabilities between 0 (black) and 1 (bright red).
2112 NucleicAcidsResearch,2013,Vol.41,No.4
d)
e
u
n
AP2 (conti
Annotation Alu-JAlu-JAlu-JAlu-JAlu-JAlu-JAlu-JEGR1,DEAF1NRF1RFX1SP1,EGR1,DP1,PAX5,EGR1,SP1,AP2 MITFAP4,E12,MITF,E47 RFX1 aMEF2AP4
83649433906 8 5645182 94 322544253 8382627263109737
97700841612 1 1181975 98 427311005 7141605460447488
18764143655 2 5492983 53 091287241 4785759332032772
eff 0025793639077339446777 40 20620956430754 0987 785145159377655736 16038367524233250896877671635649
o 79540975321 9 8776665 54 444433333 3322222221111111
C 6.5.5.5.5.4.4.4.4.4.4. 3. 3.3.3.3.3.3.3. 3.3. 3.3.3.3.3.3.3.3.3. 3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.
Sequence CAGCAACCTCAGCTCACAGCCAACCTCAAATCACAGCAACCACAGCAACCGGTTGCTGTGTTGAGGTTGCGCCGCCGCCGGCGCCTGCGCGTTGCCATGGGGGGCGGGGC GGGGGCGGGG CAAGTTCTTTTCTTCTTTGATATTGATTCTAGGAGGAGGATGCTGCTGCTGGTTTGTATTCAGCTCACAG TTTAATCAAAGTTTTAATTA CCTCAGTTTCATTTGCATTTCTTATAAATTCCTAGCTCACGGAGCTGCTGCCAGCTGTGGCATTTCCTGCAGCAGCTGCTCAGCTGCTGC TTTTCATTTACCTCAAACTCCCAGCTCCAGGGTTGCTATGATTTCCTGTTGTTAATTTTATATTTATGAAGAAATGCAAAAAACTGAGGCTTAATTAATACATTAATAAATATTTTTATAGAGCAGCTGCCTTCCATTTTCTCATTTAATTTTTCTATGC
k
Ran 3691215182124273033 36 39424548515457 6063 666972757881848790 939699102105108111114117120123126129132135138
2 7 7
1 4 4
E E E
1, 2 7, 2, 2,
P P 4 1 1
Y1D A E E E
Annotation Alu-JAlu-JAlu-JAlu-JAlu-JAlu-JAlu-JAlu-JEGR1 NFY Alu-J NF1AP4 SP1E2F,RB,YSP1,EGR1,PAX5AP2E4BP4 EGR1,SP1, RAR-alpha MITF,AP4, MITF,AP4, AP4,MITF,
44608493665 5 7148886 96 220756567 9565069044758270
29524874564 0 6892069 19 247399263 4001475124053945
94497121062 7 1785991 65 396495366 3928509779907568
eff 5640655477602516401515 85 87951962482585 2392 816048279478686236 22039774564636281596898271666252
o 01741987432 0 8776665 54 444433333 3322222221111111
C 7.6.5.5.5.4.4.4.4.4.4. 4. 3.3.3.3.3.3.3. 3.3. 3.3.3.3.3.3.3.3.3. 3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.
entioncoefficients Sequence TGAGGTTGCTGCAACCTCAAAGGTTGCTGTGAGGTTGCTGGTTGCTGTGAGCTCACAGCAAACCTCAAACGTTTGAGGTTCGCCGCCGCCTCCAGCTGTGTCTGATTGGC CTGTGAGCTA TGATTGGCTGAGCTGCTGCTTTTTTCATTATTTAATGTTTCCCCTCCCCCGCCGCCATCTGCCCCGCCCC TTTAAATAAAAAGTTCTTTG GGGCGGGGCCTTATCTTGATCATTTGTTTTTAGTTTAATTGGAAACTGAGGAGGAGGAGGGGAGGAGGAGTGTCTATTTCTTTGTATTTA CTGCAGCTGCGCCTCAGTTTTTATTTCAAAGTGTCTATTTCAGCTGCTCCCTGCTGCTGGGAAACTGAGGGCAGCTGCTGAAATATATTGTTTTTCTTCATTTTCATTAATGCAGCTGCTACCTCCAACTCTATTTATAAGGCCCCGCCCTTTATTTCAA
et
r k
nary Ran 2581114172023262932 35 38414447505356 5962 656871747780838689 929598101104107110113116119122125128131134137
o
uti
ol
v
e
of A
5
motifsinterms Annotation Alu-JAlu-JYY1,E2F,RBAlu-JAlu-JAlu-JAlu-JAlu-JAlu-JEGR1,DEAF1Alu-J Alu-J NF1 E12 NFY AP4,E12,MITF,E47 NFYNRF1 NFY MITF,AP4,E12,E47 E12 NF-AT1FOXM1,STAT E12,HEB
g
kin 0938984025142128268 17 75126154540341 5742 763922773101801410 53463216841428023650556351116494
op-ran Coeff 7.79066.27915.89815.45995.20654.99054.85334.72274.48654.31844.2322 4.0896 3.91773.81793.72553.68263.65333.63173.6012 3.54933.4961 3.48163.47063.45193.44213.40663.38333.37463.36473.3536 3.33613.30523.30103.28253.26263.25033.24053.23033.21833.19813.19053.18553.17503.16753.16353.1555
t
Annotationsof1. Sequence AGCAACCTCAACAGCAACCTCCGCCATCTTTTGCTGTGAGACCTCAAACTTTTGAGGTTGTGCTGTGAGCCTCACAGCAAGCTGTGAGCTCCGCCGCCGCTAGCTCACAG AGTTTGAGGT CTGATTGGCTCCCCCCCCCCCTGCTGCTGCTCTTGCTCAGAGCCAATCAGTTGATTCTTAGCTGCTGCTG AAACTGGTTTGCAGCAGCTG TTAATTTCAATTTATTGTAGCCCCGCCCCCCCAATCAGCGGCGCATGCGCATTTATTGTACTCTGATTGGTTTTAAATGCCAGCAGCTGC CTTTTGTTGTCTCAGTTTCCCAAATATTTGCTATTTTTAGGCAGCAGCAGCAGCTGTGGTAATTATTTTGTCATTTCCTCTACATTTCCCTGTTTTAATTTTGTATTTATCTAGCTCACAGCTGCTGCTTCCAGCAGGTGAATATATTGAATTCTTAGAA
e k
Tabl Ran 1471013161922252831 34 37404346495255 5861 646770737679828588 919497100103106109112115118121124127130133136
NucleicAcidsResearch,2013,Vol.41,No.4 2113
2
F
R
N
B,
E
R
C
ation AP1,
not F6,
n T
A A
5327605918057307968938632233145
9124237122318986478852653018820
9765896274649163371366728598419
eff 43393124221916050088817876716967605749454237302723201610070301
o 1111111110000000000000000000000
C 3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.
uence TAATTTCATATTGATTCAGCAGCTATATTTGTATCCATTTTTCTTCTTATTTCCTGTGAGCTAGGCCTAGCTTTCCAAATTCCTCTTTTATTTTGATTAACTTCTTCCATTTATTTAATTTTACATTTTTTGTTTTTCTATGCAATTGATTCCTTCTTTGTTCAAATTTTATTTATCAACTCCCCAGCAGGTATTTTCTTACTGTTTATGATGTCATTAATTCTCTTCAAGCAATTTCTAATTTTATT
Seq TTTATGAAAAACTCTGCATTTTATATATTCAAGATTATCTTACATACCGAAATGTTTTTTAC
k
Ran 141144147150153156159162165168171174177180183186189192195198201204207210213216219222225228231
a
ph F2
al 3
nnotation FX1 U.1BX1A,HNF3 KX6,HNF-1 U.1 P4,E12,E47 OU3F1,POUTAT5A
A R PP N P A PS
1594124996042708435705429115720
5513595401352786361974091400174
4366527051675769865522469608943
eff 44423425232018150391817977736968625850474338332723211712070302
o 1111111110000000000000000000000
C 3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.
uence CCTAGCTCATTTCCTTTTTTCTTCTGCTATGGCTGGTTTGTTCAAACTATTAGCATTTCATCTTATGTTTCACCAGCTGTTCCTCTTTTAAATAAATAATTTCATTTTCATTGCAGCAGCATTATTTTTATTTGTTGAGCTGCTTCCTGTTTCTTATTTTGTTTTTCTGCAGCTGGTTTTAAATACATCATCTTTAAATGACATTTCCTCATTATTTCTGAGCCCTGCAGCTTTTATTGATCAAATAA
Seq AGTTGTGTAACTTCTTAATTCTTAATATCAAAAAAGCTTGCTCACATGATCTTCTCTCTATT
k
Ran 140143146149152155158161164167170173176179182185188191194197200203206209212215218221224227230
3
0
X
O
F
2,
A
X
O F1
Annotation IPF1FOXA1,F POU2F1 RFX1 NRF1,EL NFY
606212419324272580269895873660
6352495811001493865853715116654
3627617636706832728680543537263
eff 46433626232218150494828177747069646054474439333024221813080502
o 1111111110000000000000000000000
C 3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.3.
Continued1. Sequence ATTTCCTCTGATTTTCATTTTTCTGCTACTTTATTATCTTTAATTAATAATGTTTACTTACTCTTTTGTTCTATTAATTTTTTGTGATTTTGTTTTTCTTTTCCTGCTCTCTTGAATTTTAAATGCAAATAGTTTAATTAAATTTCCTCTTCTATTAATTTTATTTTGAATGTGACCTTGTTGATTCTGGTTATAAATAAATTAACTTTTCCATGGCAACATCATCATCACATCTGTAAAGTAATTATTTCTTCCTGGAGCTCTAATTACTGTTTCATTTTCAGCCTAGCTTTTTACATTTTCTGATTGG
e k
Tabl Ran 139142145148151154157160163166169172175178181184187190193196199202205208211214217220223226229
2114 NucleicAcidsResearch,2013,Vol.41,No.4
normalized rank xð0(cid:6)x(cid:6)1Þ matches TRANSFAC proteins such as SP1 (42), AP2 (43), NRF1 (44) and
motifs. If evolutionary retention coefficients are E2F1 (44). Reciprocally, the ChIP-seq files of seven
uncorrelated with the presence of TRANSFAC motifs, proteins contain at least 10 enriched motif sequences:
then all TRANSFAC motifs should be evenly distributed ELF1, GABP, YY1, ERG1, RAD21, POL2 and PAX5.
along the normalized ranks, and X follows a uniform dis- Some of these proteins (such as SP1, AP2, POL2, YY1,
tribution. Therefore, enrichment of TRANSFAC motifs RAD21) ubiquitously regulate many genes, thus their
onhigh-scoringsequencesisquantifiedbythedeviationof binding motifs yield high evolutionary retention
the empirical distribution of X from a uniform coefficients.
distribution. Four motif–protein combinations enriched in
Figure 4 plots the empirical CDF of X (F (x) in ENCODE files correspond to exact match with the
1
Equation 10) and the CDF of a uniform distribution TRANSFAC data. Motifs 33 (GGGGCGGGGC) and
(F (x)). F (x) lies above F (x) for all x2[0,1], indicating 36 (GGGGGCGGGG) match the SP1-binding motif in
0 1 0
that sequences with high evolutionary retention coeffi- TRANSFAC. Motif 7 (CCGCCATCTT) matches the
cients are more likely to match TRANSFAC motifs YY1-binding motif, and motif 170 (CTTCCTCTTT)
than those with low evolutionary retention coefficients. matches the PU.1-binding motif.
The P-value of the Kolmogorov–Smirnov test <10(cid:2)325.
Second, we showed that the top-ranking motifs were Genes harbouring high-scoring motifs tend to retain
enriched in the protein-binding DNA fragments reported functional coherence
from the ENCODE data (23). We downloaded 407 files
Genes sharing the same protein-binding sequences on
from the ENCODE website. Each file reports the
protein-binding DNA fragments generated by one their promoters are likely co-regulated by the same tran-
ChIP-seq experiment with a specified antibody and cell scription factors. Consequently, we expect genes
type. The 407 files cover 59 proteins (transcription harbouring high-scoring motifs to possess functional co-
factors, RNA polymerase II, nucleosome-binding herence. We validated this prediction with two tests using
proteins, etc.), where multiple ChIP-seq experiments externaldata.First,usingthemethoddescribedin(39,40),
with distinct cell types and replicates were undertaken weevaluatedthedivergenceofexpressionprofilesofgenes
foreachprotein.Foreachmotif,wequantifiedthesignifi- from two human expression datasets and one mouse ex-
cance of its enrichment in an ENCODE file with a null pressiondataset.Thedistributionofexpressiondivergence
model assuming that its occurrence followed a Poisson on genes harbouring the top 5000 motifs was compared
process with a rate (cid:7)¼Nl, where N denoted the number with the distribution on the genes harbouring the bottom
L
of motif occurrence in the entire genome, L the genome 5000 motifs. Intriguingly, genes harbouring the top 5000
length and l the total fragment length in the file. motifs have consistently lower expression divergence than
We evaluated the enrichment P-values for the genes harbouring bottom 5000 motifs across all three
top-ranking motifs in each ENCODE file. About 12% datasets. The Wilcox test P-values of the deviation
of the motif-file combinations (11063 of 94017) exhibit between the two gene sets are significant across the three
significant enrichment (P(cid:6)10(cid:2)20). To reduce errors datasets: 2.047 (cid:10) 10(cid:2)14, 2.414 (cid:10) 10(cid:2)4 and 1.670 (cid:10) 10(cid:2)12
generated by individual ChIP-seq experiments, we respectively. Furthermore, by ruling out the two
grouped the results of the same proteins together and confounding factors for co-expression—co-localization
countedthefractionsoffilesineachgroupwithsignificant of genes on the same chromosomes and paralogous
enrichment. There are 182 motif-protein combinations genes sharing the same ancestry—the deviation of
with at least 5 ENCODE files and significant enrichment expression divergence between genes harbouring top
P-values ((cid:6)10(cid:2)20) in at least 80% of the constituent 5000 and bottom 5000 motifs remains pronounced.
ENCODE files. The number of enriched motif–protein Table 3 reports the significance of the deviation of
combinations drops considerably in control motifs. expression divergence between gene pairs harbouring top
Among the 231 control motifs in the middle of the 5000andbottom5000motifs.Thedeviationofexpression
sorted list, 80 motif–protein combinations retain the divergence suggests that genes harbouring motifs of high
same level of enrichment. Furthermore, among the 231 evolutionary retention coefficients tend to retain
control motifs in the bottom of the sorted list, only 38 functional coherence.
motif–protein combinations retain the same level of en- Second,weextractedthehumangenesharbouringeach
richment. Intriguingly, both TRANSFAC and of the top 231 motif sequences and assessed their over-
ENCODE data indicate that levels of enrichment shrink representations in 3201 GO categories and 889 pathways
byhalffromthetoptothemiddleandfromthemiddleto from three sources. Supplementary Table S5 reports the
the bottom of the sorted list. functional categories and pathways significantly enriched
Table2showsthe182motif–proteincombinationswith (hyper-geometricP(cid:6)0:001)witheachtop-rankingmotif.
significant enrichment. Six motifs are enriched in the There are 45 motif-functional class pairs with significant
ChIP-seq data of at least 10 proteins. These motifs are enrichment.Incontrast,thereareonly9significantmotif-
heavily biased toward GC-rich sequences: GCGCCTGC functional classpairs amongthe231control motifs inthe
GC(index27),GGGGCGGGGC(index33),GCCCCGC middle of the sorted list.
CCC(index56),GGGCGGGGCC(index65),GCGCAT By examining the functional enrichment results in
GCGC (index 76) and GGCCCCGCCC (index 134). The Supplementary Table S5, we found that many top-
GC-rich sequences match the binding motifs of several ranking motifs were highly enriched in functional classes