.
There are multiple pages of motifs ordered by decreasing score, and you can either go sequentially or
jump to another page number, or search for a specific motif name. By clicking an individual motif cluster
you can see detailed information, including the position frequency matrix (PFM) and
elementary motifs aligned to the cluster motif.
Fig. 3. Motif output results dialog
Fig. 4. List of motifs (portion of the page)
The bottom portion of the page with output results (Fig.3) can be used to save motifs (PFMs) for future
use. You need to specify file name and, if needed, some additional desctiption. The computer provides a
suggested file name that you can change.
5. Opening a file with motifs
In the My CisFinder web page (Fig. 1), click the button "Open motif file" to explore saved motif files
including preloaded public databases of motifs (HOCOMOCO, X-TFBS, HOMER, JASPAR, and SwissRegulon). Use
the pull-down list in the dialog box to select the motif file and click the button "Open motif file". The
menu of motif-related analysis (Fig 5) includes displaying the list of motifs (as shown in Fig. 4),
comparison of motif similarity with motifs in other files, motif clustering, improving motifs by resampling,
and searching for motif occurence in DNA sequences.
Fig. 5. Menu for exploring motifs.
6. Comparing motifs in two motif files
Any two files with motifs can be compared with each other by similarity (i.e., correlation of
PWMs). Such comparison is helpful for annotating newly generated motifs when they are compared to
annotated public databases of motifs (HOCOMOCO, X-TFBS, HOMER, JASPAR, and SwissRegulon).
In the Exploring motifs menu (Fig. 5), click the button "Compare motifs" to open a dialog box
where you can select another file with motifs and select the correlation treshold. Option "best"
means to find the best matching motif from the second file for each motif in the open file.
Click the "Continue" button to start the
analysis. The output screen aligned motifs in the first file (highlighted) with motifs of the second file.
7. Clustering DNA motifs
Motifs in any motif file can be clustered by similarity (correlation of
PWMs). In the Exploring motifs menu (Fig. 5), click the button "Cluster motifs" to open a dialog box
for selecting the correlation threshold. Click "Continue" to poceed. The
results of analysis are similar to the results of finding motifs,
8. Improving motifs by resampling
Motif PFM can be fine-tuned using the resampling method, which is an iterative
procedure where a PFM is used to search the sequence and select matching locations, and then,
the PFM is updated based on nucleotide frequency distribution at selected locations.
In the Exploring motifs menu (Fig. 5), click the button "Improve motifs" to open a dialog box
for selecting input data files.
Because resampling is a computationally-intensive task, the number of motifs is limited to 50.
You need to select the "test" sequence file and, optionally, the "control"
sequence file.Then click the button "Improve motifs" to open the page with parameters, most of which
are the same as for motif identification (see above). Three resampling methods are available:
(1) regression (default), (2) simple resampling, and (3) subtraction. Simple resampling
does not use the control sequence, whereas regression and subtraction use the control
sequence (see details in Algorithms). Select the number of iterations
and sub-iterations. The difference between them is that sub-iterations
using the same set of identified sites and only the rank of matching to motif is modified,
whereas iterations start with finding a new set of sites that match the motif.
Sub-iterations are much faster than real iterations and they can increase the speed
of the overall procedure.
In the paramenters page click the button "Continue" to start the analysis. The results of analysis
are similar to the results of finding motifs, except that there are no elementary motifs.
9. Searching motifs in a sequence file
Prediction of TF binding sites is based on finding locations in the genome sequence that match
the given binding motif represented by a PFM. In the Exploring motifs menu (Fig. 5),
click the button "Search motifs in sequence" to open a dialog box for selecting the
sequence file. Then click the button "Search motifs in sequence" to
open the parameter page. There are two ways to control the stringency of
motif matching. The first method is to control the number of false positive matches
per 10 Kb of sequence length. The default is 5 false positive matches per 10 Kb of
random sequence. The second method is to select "Use existing thresholds" in the
"Number of false positives" pull-down menu. In this case, the program will use
pre-defined match thresholds uploaded together with the
motif file. Thresholds should be listed in the column called "Threshold"
(see file formats), and a match is estimated as a sum of elements
of the PWM (log-transformed PFM, log10). If thresholds are pre-defined, then the
stringency can be adjusted by increasing/decreasing these thresholds by a certain
constant parameter "Add this value to the score threshold").
The program removes overlapping and redundant motifs as follows. If two instances
of the same motif overlap, then only one of them (with a higher match score) is retained.
For example, palindromic motifs that match in both dirsctions are counted only once.
In addition, if two different motifs overlap by more than 75% length, then only one
is shown in the output. The option of removing overlapping and redundant motifs can
be turned off; in this case all matches will be listed in the output.
To analyze the probability distribution of binding sites (see below) it is necessary
to specify the sequence interval for counting motifs (100 bp is a default).
Click "Continue" to poceed. The results of search are described below; they can be
saved for future use.
10. Exploring motif search results
In "My CisFinder" web page (Fig. 1), click the button "Open search motifs" to open a dialog box for
selecting the output file of search results. Then click the button "Open motif search" to
process the results and open a new web page as shown in a Fig. 6 example. At the top of the page there
are links to three tables in text format: Search results (large file), motif frequency distribution
(small file), and motif abundance table (small file).
Fig. 6.Exploring results of motif search
Information on individual
motifs and individual sequence segments (e.g., ChIP-seq peaks) is accessible by clicking buttons "Show all
motifs" and "Show all sequences", respectively. Clicking on these buttons opens lists
of all motifs and all sequences, and each motif name or sequence name is a link to detailed information
about the motif or sequence. The detailed information for a motif includes the frequency distribution
graph, as shown in Fig. 7. In this example, motifs were searched in 2000-bp long sequences centered at
binding sites of POU5F1 identified by ChIP-seq. The peak in the middle (at 1000 bp) corresponds
to the binding site of POU5F1. Motif OCT-3N-SOX on the left, which is an alternative motif with a 3-bp spacer
between OCT and SOX binding motifs, has a much lower frequency than the main motif OCT-SOX shown on the right.
The central peak of OCT-3N-SOX is detectable but its height is close to the background. In contrast,
the main motif OCT-SOT has a very high central peak dominating over the low background in the flanks.
Fig. 7. Frequency distribution of composite motifs OCT-SOX and OCT-3N-SOX in 2-Kb
genomic regions centered at locations identified with ChIP-seq for transcription factor POU5F1
(data from Chen et al. 2008). Magenta line shows all sites matching to motifs
whereas blue, green, and gray lines show 50, 25, and 12.5% sites with the highest
match score.
The abundance table shows the number of motif matches in each sequence. It can be used
for classifying sequences based on the composition of binding motifs.
Click on the button "Show all motifs" to get a list of motifs with the total
count of matches for each motif. Click on the button "Show all sequences" to get
a list of sequences with the total count of matches for each sequence. If
genomic coordinates and/or attributes were uploaded in association with the
sequence file, then the table of sequences will include the coordinates and attributes
(e.g., closest genes). It is also possible to search sequences by name, number
of motif matches, or genes. If you click
on the sequence name, you will get a page with detailed information about the sequence and
motif matches. All motif matches are listed, and their coordinates are shown. If
coordinates were uploaded, then
hyperlinks will appear which lead to corresponding locations in genome the gene
browser: UCSC genome browser.
From this page, you can visualize the location of motifs by
selecting motifs (up to 3 motifs can be selected) and then clicking the button
"Show sequence". The sequence with highlighted binding motifs will appear.
11. File upload and file formats
All files used in CisFinder are in plain text format and use tabs for delimiting
fields within one line. There is a total of 3 main types of files in the CisFinder: sequences, motifs,
and motif search results. Motifs can be submitted as PFMs and as ;attern sequences,
which are converted to the PFM form. These 3 types of files
may have 2 optional lines of annotation that start with the keyword "Parameters:"
and "Headers:." Parameters may specify the origin of the file, such as algorithm
parameters for derivation, whereas headers specify the columns of tab-delimited lines.
The menu for uploading new data files is shown in Fig. 8.
Fig. 8. Menu for new file upload.
11.1. Sequence files are in fasta format: the sequence name is preceded by ">" sign
and the following lines contain the sequence. Maximum sequence length = 25Mb. Minimum
recommended length = 5Kb; the program will process smaller files but has a low
chance of finding a motif.
Example:
>PET070528
AACCCAAAGTATGATATGCTATGATAGATAACCAAAAGGTAATATTATGAAATTTTTATCAACTATAATTATATAACTTG
AAACTGTTTCCTAAATCCGCCCTAGAGCTTACACAAAGCTGAGGGAAGTTTGCTGGAAAGTTCAGGCTGAGTGGGATGTT
>PET070400
TACTATTGGCGCTTCAATCAGTATTCGTCTTTTATAATACAATAATGCTATTTTGGATAAGTAAGTTTCTATTCAAGGAC
ACGTGTGGGCAGCTGTAACACTAATAATGTCCCATAAATAAGCGAGCAGAGCACATACTGCTGAGACAGACATGTAAGAA
............................................
To extract sequences from genome, make a tab-delimited text file with genome coordinates
ither in the bed format: (chr start end name) or alternative format (name chr strand length start).
In thr latter case, strand is either "+" or "-".
In the "Upload new data file" page of the menu, select "Genome locations".
Alternatively you can paste the contents of the file
into the top area. Select file type = Genome locations. Specify the name of the file to be
generated in the "File name:" field. This file will have fasta format, hence it's name
should end with ".fa" extension. File description is optional.
To extract a subset of sequences
go to the "Upload New Data" menu (Fig.8), open "File type" pull-down menu and select
"Subset of sequences". In the "Related to" pull-down menu, select the file from which
sequences are going to be extracted, and paste a list of sequence names to be extracted
(or select a file with these sequence names). Specify a new file name
in the "File name:" field, add description, and push the "Upload" button. If the
original sequence has associated files (coordinates, attributes) then this
information is transferred automatically to the new file.
To extract regions from each sequence of already uploaded file, select file_type option
"Sub-regions of sequences", fill out th "File name" field, and enter the coordinates of
the region to be extracted (positions are counted from zero) into the paste box at the top.
For example, coordinates "0-500" mean to extract the first 500 bp of the sequence.
Sub-sequences can be used for
finding motifs over-represented in a specific region. For example, if you have extracted 2000-bp
long sequences centered at TF binding sites identified with ChIP-seq, then to extract the central
200-bp regions of these sequences, specify the region: "900-1100". Motif frequency in this
region can be compared to control flanking regions that are from 500 to 1000 bp away from
binding sites. To extract flanking regions, use coordinates "0-500,1500-2000". In this case,
two regions will be extracted from each sequence. If the
original sequence has associated files (e.g., attributes), then this
information will be transferred automatically to the new file.
11.2. Motif files are submitted either in the MEME format or as follows. The first line starts with a ">" sign
followed by motif name. The same line may contain additional tab-delimited fields:
| Pattern | (=consensus),
|
| PatternRev | (=reverse consensus),
|
| Threshold | (=threshold score),
|
| Nmembers | (=number of member motifs in a cluster),
|
| Freq | (=number of motif matches),
|
| Ratio | (=enrichment ratio of motif matches),
|
| Info | (=information content of motif PFM),
|
| Score | (=motif score used for ordering),
|
| FDR | (=False Discovery Rate),
|
Here is an example, where an empty line is placed at the end of each
matrix.:
| Parameters: | matchThresh=0.8500 | nucleotideOrder=A,C,G,T
|
| Headers: | Name | Pattern | PatternRev | Freq | Ratio | Info | Score | p | FDR | Palindrome | Nmembers
|
| >SOX9 | CCWTTGTT | KAACAAWG | 12813 | 4.5857 | 9.407 | 176.0950 | 0.0000 | 0.0000 | 0 | 151
|
| 0 | 2 | 72 | 13 | 13
|
| 1 | 3 | 72 | 2 | 23
|
| 2 | 46 | 1 | 2 | 51
|
| 3 | 1 | 2 | 1 | 96
|
| 4 | 1 | 1 | 1 | 96
|
| 5 | 2 | 4 | 93 | 1
|
| 6 | 3 | 2 | 2 | 93
|
| 7 | 6 | 26 | 13 | 55
|
|
|
| >OCT | HATGCWAA | ATTWGCAT | 404 | 3.1375 | 10.055 | 18.8350 | 0.0000 | 0.0000 | 0 | 4
|
| 0 | 93 | 3 | 3 | 2
|
| 1 | 1 | 1 | 2 | 96
|
| 2 | 1 | 1 | 97 | 1
|
| 3 | 3 | 76 | 10 | 11
|
| 4 | 61 | 3 | 1 | 36
|
| 5 | 80 | 2 | 17 | 2
|
| 6 | 96 | 1 | 2 | 1
|
| 7 | 3 | 15 | 14 | 69
|
Motifs can be also uploaded as a pattern (consensus), then the file is formatted as follows:
| Parameters: | matchThresh=0.8500 | nucleotideOrder=A,C,G,T
|
| Headers: | Name | Pattern | TFgroup | Info
|
| >MIT_001NRF1 | RCGCANGCGY | NRF1 | 16
|
| >MIT_002MYC | CACGTG | MYC | 12
|
| >MIT_003ELK | SCGGAAGY | ELK | 14
|
| >MIT_005NFY | GATTGGY | NFY | 13
|
| >MIT_009ATF | TGAYRTCA | ATF | 14
|
| >MIT_010YY1 | GCCATNTTG | YY1 | 16
|
11.3. Search results are formatted as a tab-delimited text file with the following
columns: MotifName, SeqName (sequence name), Strand, Len (motif length), Start (starting position counted from 0),
Score (matching score), Sequence (matching sequence).
The "Parameters:" line should specify the name of the sequence (file_fasta) file that was searched
and the motif name that contain PFM (file_motif). Here is an example of search results:
Parameters: file_motif=public-motif_pluripotent file_fasta=public-P300_binding.fa
Headers: MotifName SeqName Strand Len Start Score Sequence
STAT3 Chen2008P300000002-900-1100 - 9 119 4.3932 TTCCCGGAA
TEF Chen2008P300000005-900-1100 - 8 87 3.3398 AGGAATGC
NANOG Chen2008P300000004-900-1100 + 9 116 3.5202 CCACTTCCT
KLF Chen2008P300000004-900-1100 + 8 96 3.7154 CTCCACCC
KLF4 Chen2008P300000004-900-1100 + 9 109 4.0565 GCCACACCC
SOX9 Chen2008P300000003-900-1100 - 8 60 3.7513 CCATTGTT
11.4. Sequence attributes files have multiple columns (tab-delimited). The first column
has sequence names. Column names are listed in the "Headers" line (which is a starting line).
Use the header "Symbols" for gene symbols, and "Distances" for distance to gene transcription start site
(TSS). Several gene symbols can be separated by commas; in this case, distances are also
comma-separated. Distances are negative if a sequence is located upstream of TSS and
positive if it is downstream of TSS. Here is an example of attribute file:
Headers: Name Symbols Distances
Chen2008P300000001 Sulf1 156
Chen2008P300000006 Tmem131 -46164
Chen2008P300000010 Pecr,Tmem169 22150,-22285
Chen2008P300000015 Armc9,Htr2b 33225,-76052
To upload a file with attributes, select "Sequence attributes" in the "File type" pull-down
menu (Fig. 7); select the sequence (fasta) file with which genome attributes file is
associated (both should have matching sequence names) in the "Related to:" pull-down
menu. Do not fill the "File name:" box because the attribute file will be automatically
renamed as the name of sequence file plus ".attr" ending. Do not fill in a description either.
12. Algorithms
12.1. Estimating position frequency matrix (PFM) from 8-mer word counts
The proposed algorithm is based on estimating PFM directly from 8-mer word counts in the
test and control sequences. We will start with a numerical example to make the algorithm
clear, and then will proceed to the formal justification of the method.
Consider a specific 8-mer word W (e.g., W = "ATGCAAAT") which has T(W) = 200 matches
(instances) in the set of test sequences and C(W) = 50 matches in the set of control
sequences. For simplicity, we will not consider the distribution of word instances in
individual sequences within the test set and control set and will count only the total
number of instances as if all sequences in a set are concatenated. However, the CisFinder
has an option to count only one match of each word per sequence. The total length of all
test sequences and all control sequences is assumed to be the same (3 Mb, in this example).
For word W, we define a nucleotide substitution matrix [Wpi] which contains words that are
derived from W by placing a nucleotide i in position p (Fig.9).
Fig. 9. Nucleotide substitution matrix derived from word W = "ATGCAAAT".
The frequency of each
word from the nucleotide substitution matrix counted in the same target sequence makes the
frequency substitution matrix (Fig. 10, 11). For convenience, we will use brief notations
Tpi = T(Wpi) and Cpi = C(Wpi) for the frequency of word Wpi
matches in the test and control sequence sets (elements of frequency substitution matrixes):
Fig. 10. Frequency substitution matrixes for the test and control sequence sets.
The method for estimating PFMs is
|  | | (e1) |
where jpi is the estimate of PFM element,
and Tpi and Cpi are the counts
of word Wpi, in the test and control sequences, respectively. Because word
counts are random variables, they may appear smaller in the test sequences than in control
sequences by chance, yielding a negative PFM element. To avoid this, negative differences
are replaced with zero and then normalized (Fig.11).
|  | | (e2) |
If test and control sequence sets have different total lengths,
then the number of word counts in the control sequences is adjusted by total sequence
length.
Fig. 11. Estimation of PFM in 3 steps: subtraction of frequency substitution matrixes,
replacing negative values by zeroes, and normalization.
This method is justified from the following model. Let us assume that a TF binds to a set
of locations in the genome where corresponding DNA sequences can be aligned together.
Using this alignment, we can estimate the frequency, fpi, of each nucleotide i
in each position of aligned sequences, p, which is the element of the PFM. We
further assume that binding strength of the TF is additive with no interaction
between positions. This simplification is justified by the fact that all existing
databases use PFMs to describe TF motifs, and this strategy works reasonably well.
A successful ChIP-seq experiment generates a set of genome locations which is
enriched in TF binding sites. The test set of sequences typically consists of
200-bp long segments centered at projected binding sites, whereas the control
set of sequences can be taken away from binding locations (e.g., 500 bp away).
In example shown in Fig. 12, there are 2 control sequence fragments near each projected
binding site. It does not matter that control sequences are longer than test sequences
because the frequency of words in the control set of sequences, Cpi, is
automatically adjusted to sequence length. However, if word frequencies are counted only
once per individual sequence, then sequence length should be the same in the test
and control data sets.
Fig. 12. Example of selecting test and control sequences.
Consider a word W with a sequence of nucleotides that corresponds to the maximum values
of the PFM at each position. This word is then used to generate frequency substitution
matrixes [Tpi] and [Cpi] for the test and control sets of
sequences, respectively. Each instance of word Wpi in the test or control sequences can
either correspond to a true binding site of the TF (we call it functional) or not
(non-functional). Factors determining functionality of different instances of the
same DNA word are largely unknown and may include sequence context and chromatin status.
Because the probability of TF binding is proportional to PFM elements at each position
(assumption of additive contribution of each position to TF binding), the number of
functional instances, FT(Wpi), of word Wpi in the test sequences
is proportional to fpi:
|  | | (e3) |
The total number of instance of word Wpi in test sequences
equals the sum of functional, FT(Wpi), and non-functional,
NT(Wpi), instances:
|  | | (e4) |
Similarly, the total number of instance of word Wpi in control sequences
equals the sum of functional, FC(Wpi), and non-functional,
NC(Wpi), instances. Although the functional instances are enriched in the test
sequences compared to control sequences, some functional instances may be present
in the control sequences because of possible false negatives in ChIP data. Because
non-functional instances of the word are not affected by the ChIP procedure, their
counts are equal in the test and control sets of sequences:
NT(Wpi) = NC(Wpi). Then, the numerator in (e1) is
|  | (e5)
|
Because functional instances of word W are over-represented in the test set of
sequences compared to control, the final sum in (e5) is always positive and the
difference (Tpi - Cpi) is proportional to fpi. Thus,
the equation (e1) gives a true estimate of fpi in the PFM. This reasoning remains
correct if the word W is shorter than the full binding motif or includes a gap. However,
the word should be long enough to capture the informative portion of the motif so
that it remains strongly over-represented in the set of test sequences compared to control.
Because the PFM is estimated as a difference between word counts in the test and control
sets of sequences (e1,e2); thus, the variance of PFM elements is equal to the sum of
variances of word counts in the test and control sequences. The variance of word counts
is very close to the mean, which is expected from the Poisson distribution (we checked
it using pseudo-random sequences generated with the 3rd order Markov process). For
example, if word counts are 120 in the test set of sequences and 40 in the control
set (i.e., 3-fold over-representation), then the relative error (accuracy) is equal
to sqrt(120+40)/(120-40) = 0.158.
12.2. Implementation of the method for PFM estimation
A successful ChIP-seq experiment generates a set of genome locations that are enriched in TF
binding sites. For a test sequence set, we usually extract 200 bp sequence segments centered
at a peak of projected TF binding sites. For a control sequence set, we usually extract
500 bp sequence segments starting from nucleotide positions 400 bp away from both ends of
200 bp test sequence segments. (However, the CisFinder allows users to choose different
sequence lengths). The CisFinder identifies binding motifs of TFs using direct counts
of all possible 8-mer words with and without gaps in both test and control sequence sets
(Fig. 13). This word length was selected experimentally based on the observation that
longer words have too few matches in target sequences, whereas shorter words may fail to
capture the most informative portion of the motif and show lower rates of over-representation.
(Note: the command-line version of CisFinder allows the use of 6- and 10-mer words).
Word counts are stored in the array of integers. Although there are many different ways
to insert gaps in the 8-mer words, we consider only 8 specific patterns of gap insertions
(Fig. 11). We found that this limited set of gap insertion effectively helps to capture
composite motifs with multiple functional elements. For example, search for a word
"ATGCAAAT" with a 2 bp gap in the middle is equivalent to the search for word "ATGCNNAAAT".
PFM is then estimated for each word based on >1.5 fold (default threshold) enrichment in
the test sequences compared to the control sequences using (e2). The fold enrichment
criterion is optional (it can be set to 1).
Fig. 13. Words with and without gaps used for searching of over-represented motifs.
Over-representation of word counts in the test sequences, T, compared to the control
sequences, C, is then evaluated using a z-score which is estimated based on the hypergeometric
probability distribution. Let us first consider the case where only one instance of each
word is counted per each sequence of equal (or approximately equal) length. The proportion
of sequences, q, with a given word in the set of test sequences is compared with the
proportion of sequences, p, with the same word in the combined set of test and control
sequences (if the null-hypothesis is true then test and control sets of sequences can
be combined) with z-score:
| z = (q - p)/sqrt[p*(1-p)*(N-n)/(N-1)/n],
|
where n is the number of test sequences, and N is the
number of combined test and control sequences. If multiple instances of each word are
counted per each sequence, then the method is modified as follows. The set of test
sequences with the total length T is split into T/m segments of length m, where m is the
actual length of the word including gaps. Each instance of the word is then associated
with the segment where it starts. Because overlapping instances of the same word are
counted as one instance, there are not more than one instance associated with the same
segment. Similarly, the set of control sequences with the total length C is split into
C/m segments of length m. We use the same equation for the z-score (see above) where q
is the proportion of test segments with the word, p is the proportion of test and
control segments with the word, n = T/m, and N = (T + C)/m. Although occurrences of
word instances in adjacent segments may be weakly correlated, the hypergeometric
distribution gives a reasonable approximation of the z-score.
To fill the gaps and extend the length of PFMs, the test and control sequences are searched
again for the exact match to each word with z > 1.643 (to satisfy the condition of
p < 0.05 for one-tail z-test). Each match of the word in the test (or control) sequence
is then examined for nucleotides in the gaps and flanking sequences (2 bp at each side)
that are not included in the word. In this way, we can count nucleotide frequencies in
gaps and flanking regions and estimate the PFM for these positions using (e2) (Fig. 14).
The program is also designed to trim flanking sequences if they are not informative (if
the ratio of maximum frequency to minimum frequency is <3). To increase the information
content of PFMs, we use the contrasting procedure: the median of minimum PFM values at
each position is subtracted from all PFM values; negative values are then replaced by
zero; and the PFM is re-normalized.
Fig. 14. Extending the position frequency matrix (PFM) over
flanks and gaps with equation (e2)
The frequency distribution of nucleotides in the flanking regions and in gaps of a
certain word may differ substantially between the test and control sets of sequences,
and this difference can be used to increase the statistical power for identification
of significant motifs. Frequency distributions of nucleotides (counted for each
nucleotide and each flanking/gap position) in the test and control sequences were
compared using the G-test (Sokal and Rohlf 2001). Assuming that this test is independent
from the z-test for over-representation of word counts (see above), we combined
p-values from these tests using Fisher's method (Hess and Iyer, 2007. BMC Genomics, 8:96).
Then, False Discovery Rate (FDR) is estimated using method of Benjamini and
Hochberg (1995) to account for simultaneous testing of multiple hypotheses. The program
generates at least 100 top-score motifs, and additional motifs are limited to those
that satisfy the criterion of FDR < 0.05. Flanking positions are trimmed if the ratio
of maximum frequency to minimum frequency is <3.
12.3. Clustering of motifs
Motifs are clustered based on similarity (Fig. 15) and/or co-occurrence. Various methods
were proposed to measure the similarity of PFMs including Bayesian models. Here we use
a simpler method and measure similarity by Pearson correlation between elements of the
corresponding position-weight matrixes (PWMs) for all overlapping positions, where PWM is
derived from PFM by log-transformation: : xij = log(pij/qj),
xij is the weight of nucleotide j in position i, pij is the probability to
find nucleotide j in position i, and qj is the background frequency of
nucleotide j. For simplicity, here we assume equal background
frequencies (qj = 0.25); zero probabilities are avoided by adding pseudocount =1
to nucleotide counts in the PFM. Offset and orientation of motifs is selected based
on maximum correlation, restricted to the minimum overlap of 6 bp and maximum
overhang of 2 bp. Because correlation is estimated for a minimum of 6 overlapping
positions, there are at least 24 points (6 positions x 4 nucleotides) for estimating
correlation. Thus, even low correlation is significant (e.g., r=0.5; d.f.=22; p<0.05).
The default correlation threshold in CisFinder (r=0.7) is always significant, however the
user can adjust the correlation threshold to increase or decrease the size of clusters.
We use single-linkage clustering, and then each cluster is checked for homogeneity. If the
cluster is not homogeneous, it is separated into sub-clusters using the second round of
clustering. Sub-clustering is done iteratively starting from a pair of seed motifs by
sequential addition of most similar motifs and re-estimating the combined PFM for the
sub-cluster. Each pair of motifs is characterized by the score = r*m1*m2, where r is
the correlation between PWMs, and m1 and m2 are the number of linked members for motif
1 and motif 2, respectively. Then the pair with the highest score is selected as a seed
for the sub-cluster. This procedure is different from the original single-linkage
clustering because motifs are added to the sub-cluster based on the similarity to
the combined PFM of all motifs that are already included into the subcluster, whereas
original clustering is based on the similarity between individual (non-combined) motifs.
Motifs are added until no motif within the cluster can be added to the
sub-cluster using the given threshold of similarity. If all elements in the cluster
appear to be in the same sub-cluster, then the cluster is considered homogeneous.
Otherwise, the elements of the sub-cluster are removed from the cluster, and the same
algorithm is applied to the remaining elements.
Fig. 15.Clustering of motifs based on similarity of PFMs
The advantage of clustering PFMs compared to clustering words (as in RSAT) is that PFMs
contain more information than words alone.
Words differ qualitatively (the nucleotide either matches or mismatches) whereas PFMs
differ quantitatively (i.e., the probability of each nucleotide correlates between two PFMs).
Motifs within the same cluster are then arranged using the hierarchical clustering with
cluster flips to place similar motifs near each other (Eisen et al. 1998). Then, the PFM for the entire
cluster is estimated as the weighted average of member PFMs using local information
content at each position p (estimated for the position p and its two neighboring positions
at both sides) multiplied by motif abundance as a weight. Finally, a sequence logo (Schneider and Stephens 1990) is
generated from the PFM (Fig. 16).
Co-occurrence of word instances in the test sequence is
used as an alternative criterion for clustering. This method is generally less accurate
because of the limited number of word pairs in the sequence. However, it works better
for clustering PFMs with high level of self-similarity (after position shift by 1-4 bp),
because their relative position cannot be uniquely identified based on correlation.
The use of co-occurrence clustering method is optional and can be used either globally
or only for PFMs with high level of self-similarity.
Fig. 16.Sequence logo generated from combined PFM
12.4. Search for motifs that match to a PFM
Finding DNA sites that match a specific PFM in a given sequence is computationally intensive
if a matching score is estimated sequentially at each position of the sequence as in
MatInspector or MATCH. We propose a faster method to identify DNA sites, which
is based on a lookup table. For each motif represented by a PFM, we selected the most
informative stretch of 8 nucleotides, which is used as a core. Then a lookup table
is generated which specifies all PFMs from the list, whose cores match sufficiently
well to each possible 8-mer word. The length of 8 bp for the core is selected because
the number of all 8-mers is small enough to keep the lookup table in the computer memory,
and 8-mer words are enough specific to be linked with only a few PFMs that match them.
A match score is defined as a log likelihood that a specific sequence matches a matrix;
it is equal to the sum of those elements of the PWM (log-transformed PFM) which
correspond to nucleotides at each position of the sequence. The match score for the
full matrix, Tfull = T8 + Tresid, where T8
is the match score for 8-mer core and Tresid
is the match score for the residual of the matrix. The program finds the threshold
value R8 for the match score T8, which ensures that the match
score for the full matrix
exceeds the given threshold Rfull with probability 0.999 if T8 > R8:
| R8 = Rfull - F-1(0.999) | | (e6) |
where F is the cumulative probability distribution of Tresid when matching to a
random sequence. The value of F-1(0.999) is estimated by Monte-Carlo simulation.
A PFM is included into the lookup table for the 8-mer core word if T8 > R8. The query
sequence is scanned sequentially, and for each position only those matrices are tested
that are in the lookup table for the specific 8-mer word that starts at this position.
Although this method may miss up to 0.1% of matching sites, we consider it a reasonable
tradeoff for the increase in computation speed by several orders. Based on (e6), these
missed sites always have a poor match to the core motif. Although the match score of
missed sites formally exceeds the threshold, the quality of these sites is low from the
biological point of view because of the poor match to the core motif.
12.5. Resampling algorithms
Three algorithms are available for resampling: simple resampling, regression, and
difference-based. In all algorithms, the motif matching score is estimated as a sum of
elements of the PWM that correspond to the sequence at each position. Match threshold
is adjusted according to the expected number of false positives per 10,000 bp in
control sequence.
Simple resampling: in each iteration, find matching motifs with score ≥s and
then set the PFM equal to the frequency of nucleotides in matching motifs.
Regression: in each iteration, find matching motifs with score ≥s and
estimate the frequency of nucleotides in matching motifs in the test and control
sequences. Do linear regression of log-transformed nucleotide frequencies in the test
sequence versus corresponding log-transformed nucleotide frequencies in the control
sequence. Then the PFM is adjusted to bring individual points closer to the
regression line. The idea behind the method is to equalize the stringency of matches
between motif positions and individual nucleotides.
Difference-based: in each iteration, find matching motifs with score ≥s and
estimate the frequency of nucleotides in matching motifs in the test and control
sequences, then set the PFM equal to the difference between nucleotide frequency
in the test file and nucleotide frequency in the control file (replace negatives
with zero).
13. Enable scripted windows in the browser
CisFinder often opens new windows for data processing rather than continue
in the same window. Thus, you may need to enable pop-up windows to
use it successfully.
Terminology
- ChIP, ChIP-chip, ChIP-seq
- Chromatin immunoprecipitation
= proteins are crosslinked to DNA to which
they are bound; DNA is fragmented by sonication; antibody attached to beads
selectively bind to a protein of interest, and beads are magnetically pulled
purifying the protein; DNA is amplified and then either hybridized to a
tiling microarray with oligos from promoter regions (ChIP-chip), or sequenced
(ChIP-seq). DNA sequences are compared to the genome and locations of protein
binding are identified.
False Discovery Rate (FDR)
- FDR is used instead of p-values for simultaneous testing of multiple hypotheses.
It is defined as the expected proportion of false positives among tests (motifs
in this case) that are considered significant. FDR is intermediate in stringency
between p-values and Bonferroni correction. FDR becomes more stringent as the number of true
positives decreases. Here we used the method of Benjamini and Hochberg (1995) for
estimating FDR.
