Enrichment & classification

Discovery returns many motifs. The two steps below tell you where each one sits relative to the binding site, and which known factor it corresponds to.

Motif enrichment and positional distribution

KDMMap performs the positional enrichment analysis, comparing bound against unbound sequences within a sliding window:

Function: `KDMMap`

Parameter	Meaning
`motifs`	a `kdmotifs` object obtained from `KDMFind`
`pseq`	array of sequences: a set of equal length bound sequences
`nseq`	array of sequences: a set of equal length unbound sequences (sequence length must be the same as `pseq`)
`halfWin`	integer: defines the sliding window size (actual window size is `2*halfWin+1`) used to calculate the signal along the provided sequences
`sig.level`	the significance level at which a motif is considered centrally enriched (default: `0.05`)
`experiment.type`	string: either `"ChIP"` or `"Clip"`, used to define the regions along the sequences which are considered "central" (default `"ChIP"`)
Output	a list containing two elements: `results`, a data frame with the motifs identified as both centrally distributed and enriched in the bound sequences; and `profiles`, a list containing signal and hits along the sequences for each motif, used for plotting purposes (see `plot.motif.map`)

For each sequence, KDMMap slides a window, scores the motif in every window and keeps the maximum-scoring position. A score threshold chosen to maximise classification accuracy between bound and unbound sequences decides which of those maxima count as hits, so each sequence contributes at most one hit.

On the motifs obtained in the previous step and the bound and unbound sequences of our dataset:

centrality<-KDMMap(motifs=motifs,pseq=NRF1$train$pseq,nseq=NRF1$train$nseq,halfWin=15)

The first element of the output contains a data frame with the motifs which resulted centrally enriched:

centrality$results

      motif Start End bin_location       PE        CE motif_auroc    PE_p.value    PE_q.value    CE_p.value    CE_q.value np_in np_out nn_in nn_out

18455  M-105   -40  40            0 2.652382 90.250000    0.983796 9.388931e-106 1.267506e-103  0.000000e+00  0.000000e+00   361     91     4      9
17691   M-59   -38  38            0 2.615022 25.600000    0.841544  7.069717e-70  1.875921e-68 8.632840e-289 6.474630e-288   256     86    10     33
180751  M-92   -39  39            0 2.583608 32.500000    0.766384  2.325509e-53  1.569719e-51 4.258871e-233 2.613398e-232   195     62     6     16
26375   M-67   -64  64            0 1.653126 16.066667    0.805500  2.830310e-29  9.097424e-29 3.173155e-222 2.519859e-221   241     63    15     22
17303  M-121   -37  37            0 2.043038 10.000000    0.806176  2.787635e-27  3.763307e-25 3.362378e-125 6.877591e-125   180    136    18     61
28995   M-57   -74  74            0 1.444295  5.000000    0.809280  3.713062e-22  1.193484e-21 9.553374e-131 4.447260e-130   280     70    56     43
20295   M-98   -45  45            0 1.981893  7.375000    0.715312  3.960973e-28  5.347314e-26 2.593876e-101 2.059843e-100   177     87    24     65
22371  M-122   -51  51            0 1.535141  4.725490    0.821372  4.002203e-17  4.502478e-16 8.999370e-103 2.825384e-102   241    169    51    113
20651   M-53   -46  46            0 1.735484  4.772727    0.766004  2.912104e-22  3.573945e-21  1.210498e-87  1.907666e-86   210    140    44    111
28251   M-77   -71  71            0 1.380198  3.865672    0.791848  5.221001e-15  2.135864e-14  7.175777e-93  1.405411e-91   259     94    67     67
16515    M-1   -35  35            0 1.742592  4.558824    0.747740  9.129409e-15  6.616866e-13  3.852132e-59  1.368521e-58   155    182    34    114
18075   M-15   -39  39            0 1.592845  4.073171    0.742032  3.087111e-12  8.416485e-11  1.632904e-57  4.322393e-57   167    190    41    131
22035    M-5   -50  50            0 1.344868  2.756757    0.722620  8.038073e-08  8.347230e-07  4.429307e-45  8.542236e-45   204    200    74    132
27203   M-95   -67  67            0 1.270818  2.821918    0.723272  5.873790e-07  2.832006e-06  3.303985e-47  1.205508e-46   206    117    73     81
220351  M-26   -50  50            0 1.506649  2.719298    0.680516  1.349700e-10  1.656450e-09  7.107858e-32  1.573051e-31   155    119    57     95
14055   M-70   -29  29            0 1.560351  3.000000    0.737544  7.270253e-09  6.404828e-07  5.197548e-33  8.158942e-33   141    271    47    171
13631   M-13   -28  28            0 1.651754  3.027027    0.629112  9.001952e-09  6.076318e-07  4.726331e-26  1.595137e-24   112    208    37    171
25511   M-17   -61  61            0 1.238490  2.562500    0.725428  2.036558e-05  1.255477e-04  1.521125e-40  2.976113e-40   205    157    80    107
26091   M-38   -63  63            0 1.261607  2.595238    0.759484  1.415064e-06  1.289711e-04  1.022365e-44  2.706261e-44   218    148    84     76
14891   M-83   -31  31            0 1.686179  3.263158    0.528352  5.293780e-06  8.627250e-05  6.911120e-16  4.541193e-14    62     95    19     92
12771    M-6   -26  26            0 1.805369  2.944444    0.531032  4.537321e-06  1.566227e-04  4.996895e-12  2.872774e-11    53     96    18     80

The columns of centrality$results:

Column	Meaning
`motif`	name of the motif
`Start, End`	the boundaries of the region in which the motif shows higher values of PE and CE
`bin_location`	the center of the region
`PE`	Positional Enrichment: the fraction of bound-sequence hits that fall inside the region, divided by the fraction expected under a uniform positional distribution — PE = (H_r,pos / H_s,pos) / p₀, with p₀ = L_r / L_s, the region length over the sequence length
`CE`	Class Enrichment: the fraction of bound sequences with a hit in the region, divided by the same fraction in the unbound sequences — CE = (H_r,pos / N_pos) / p₀, with p₀ = H_r,neg / N_neg, where N_pos and N_neg are the numbers of bound and unbound sequences
`motif_auroc`	the AUROC of the motif's maximum score per sequence, discriminating bound from unbound sequences. It is reported as -(1 - AUROC) when the motif scores higher in the unbound sequences, so its absolute value is always ≥ 0.5 and its sign gives the direction
`PE_p.value, PE_q.value`	significance of PE: an upper-tail binomial test p-value, and its multiple-testing-corrected q-value
`CE_p.value, CE_q.value`	significance of CE: an upper-tail binomial test p-value with H_r,pos successes out of N_pos trials, and its multiple-testing-corrected q-value
`np_in, np_out`	number of hits inside and outside the region for bound sequences: `np_in` is H_r,pos, and `np_in + np_out` is H_s,pos
`nn_in, nn_out`	number of hits inside and outside the region for unbound sequences: `nn_in` is H_r,neg

PE and CE are as defined in the manuscript (Eq. 14-15).

Centrality profiles

The second element contains the scores and the number of hits for all the motifs along all the sequences: it can be used to visualize the centrality profiles as follows, for the most enriched motif:

plot.motif.map(centrality,"M-105")

Two panels for motif M-105. Left: the signal along the sequence position peaks sharply around position zero for the positive (bound) sequences and stays flat and low for the negatives. Right: hit counts cluster tightly around the centre in the positives. — `M-105` is **centrally enriched**: both signal and hits concentrate around the peak summit in the bound sequences, while the unbound ones stay flat.

And, for a motif which is not central but is still enriched with respect to the unbound sequences:

plot.motif.map(centrality,"M-19")

Two panels for motif M-19. Left: the signal for positive sequences runs consistently above the negatives across the whole window, but with no peak at the centre. Right: hits are spread evenly across all positions in both groups. — `M-19` is **enriched but not central**: the bound sequences carry more signal than the unbound ones everywhere, with no concentration at the summit.

Motifs classification

KDMMatch annotates the discovered motifs by information-content-aware matching against a reference set already embedded in KDM space:

Function: `KDMMatch`

Parameter	Meaning
`query.motifs`	a `kdmotifs` object obtained from `KDMFind` containing the motifs we want to classify
`target.motifs`	a `kdmotifs` object containing the motifs against which we want to classify the query
`target.parameters`	an object obtained through the helper function `buildKDMReference` which contains parameters for the linear models used to correct for target motifs IC content
`best.hit`	logical: if `TRUE` the function reports the best hit for each query motif. Otherwise the function reports all the matches for each motif
`threshold`	the significance level at which matches are reported (default: `0.05`)
Output	a data frame containing the matching info

The columns of that data frame:

Column	Meaning
`query`	the name of the query motif
`target`	the name of the matching target motif
`p.value`	the p-value for the matching
`q.value`	the Bonferroni corrected p-value (depends on the number of target motifs)
`similarity`	the similarity between the query and target motifs
`rank`	for each query motif, the similarity rank with the targets (1 the most similar; if `best.hit=TRUE` this is always 1)

Building a reference

To classify our motifs we need an annotated reference collection of KDM. To build such a collection the helper function buildKDMReference can be used:

Function: `buildKDMReference`

Parameter	Meaning
`motifs`	either a meme file name (see the MEME motif format) containing PWMs definitions, or an object of class `kdmotifs`
`info`	a data frame with the number of rows corresponding to the number of motifs. It must contain a `name` column with the names to assign to the motifs. It can contain other columns of your choice (i.e. the protein to which each motif is associated, or the experiment the motif comes from). If left unassigned the motif names will be `REF-1`, `REF-2` and so on
`dictionary`	the dictionary used to transform PWMs into KDMs
Output	a list containing `motifs`, the `kdmotifs` object representing the motifs; `info`, the passed info data frame; and `parameters`, an object which contains parameters for the linear models used to correct for target motifs IC content

For instance we can produce a reference for the HOCOMOCO v13 in vivo motifs collection as follows:

PWMREF<-buildKDMReference("H13INVIVO_meme_format.meme",dictionary=dictionaryTF)

In this example though, we use the reference we assembled for the paper from the ENCODE experiments. It can be loaded from the TF_ENCODE_KDM_REF.Rdata file in the examples folder:

load("TF_ENCODE_KDM_REF.Rdata")

Now we can classify the motifs we obtained before with:

classification<-KDMMatch(query.motifs=motifs,target.motifs=TF_ENCODE_KDM_REF$motifs,best.hit=TRUE,target.parameters=TF_ENCODE_KDM_REF$parameters,threshold=0.05)

obtaining:

head(classification)

  query       target       p.value       q.value similarity rank
1   M-1 TF_REF_03129  6.935578e-61  4.427673e-57  0.4082271    1
2   M-2 TF_REF_02873  3.535196e-05  1.504579e-02  0.2306516    1
3   M-3 TF_REF_05828 1.021659e-113 6.522269e-110  0.7008772    1
4   M-4 TF_REF_05436  2.420400e-20  1.545183e-16  0.4073220    1
5   M-5 TF_REF_00166  3.785907e-83  2.416923e-79  0.6171039    1
6   M-6 TF_REF_01372  5.509329e-09  1.172385e-05  0.2199218    1
. . .

It is possible to match the obtained target motifs to the corresponding information:

head(TF_ENCODE_KDM_REF$info[match(classification$target,TF_ENCODE_KDM_REF$info$name),])

             name target  cell  experiment
3129 TF_REF_03129   NRF1 HepG2 ENCSR853ADA
2873 TF_REF_02873  KDM4B HepG2 ENCSR280SCF
5828 TF_REF_05828    MAZ MCF-7 ENCSR288IJC
5436 TF_REF_05436 ZNF319  K562 ENCSR231PDA
166  TF_REF_00166   E2F6  A549 ENCSR000BTC
1372 TF_REF_01372   NRF1    H1 ENCSR000ECC

KDM — source on GitHub · preprint