psp-gen

Usage:

psp-gen -pos <primary sequences> -neg <control sequences> [options]

Description

psp-gen is used to allow MEME to perform discriminative motif discovery—to find motifs overrepresented in one set of sequences compared to in another set. It takes two files as input—the sequence file to be input to MEME, (the "primary" file) and a "control" sequence file of sequences believed not to contain the same motifs as in the "primary" file. psp-gen creates a file for use by MEME that encapsulates information about probable discriminative motifs. psp-gen records its chosen motif width in the file, and MEME is able to adjust the data when it tries different motif widths.

Input

Primary Sequences File

A file containing FASTA formatted sequences which are to be used as the primary set in PSP calculation.

Control Sequences File

A file containing FASTA formatted sequences which are to be used as the control set in PSP calculation.

Output

A FASTA-like PSP format is written to standard output, it contains a prior for every position of every sequence in the primary set.

Options

Option Parameter Description Default Behaviour
General Options
-h Print a short usage message and exit. Run as normal.
-minwmin w The minimum width to use with selecting the "best" width for PSPs. The minimum width is set to 4.
-maxwmax w The maximum width to use with selecting the "best" width for PSPs. The maximum width is set to 20.
-dna  Use the DNA alphabet. The DNA alphabet is assumed.
-protein  Use the protein alphabet. The DNA alphabet is assumed.
-rna  Use the RNA alphabet. The DNA alphabet is assumed.
-alphfile Use the alphabet defined in the file. The DNA alphabet is assumed.
-triples Use spaced triples instead of whole-word matches. Recommended for protein. Whole-word matches are used.
-fixedstart When using the -triples option, only consider triples starting at the start of the site. Allow triples to start anywhere in a width w site.
-equivLetter Groups Any sequence of letter that appears together should be matched as equal. Separate letter groups using "-", e.g. -equiv "IVL-HKR" means treat all occurrences of I, V or L as the same, and all occurrences of H, K or R as the same. Note that in alphabets where '-' is an allowed symbol separate groups may be specified by repeating the -equiv option for each group. No letter groups are used.
-revcomp Both strands are considered when calculating PSPs for complementable alphabets. Only the given strand is considered.
-scaleminmin score The lowest score value after scaling. The lowest score is set to 0.1, unless the -scalemax option is set in which case the lowest score is max score - 1.
-scalemaxmax score The highest score value after scaling. The highest score is set to 0.9, unless the -scalemin option is set in which case the highest score is min score + 1.
-maxrange Choose the width with the biggest difference between minimum and maximum scores (before scaling). Choose the width with the biggest maximum score (before scaling).
-raw Output scores instead of priors. The program output PSPs.
-reportscores Send to standard error primary and control file names, min and max scores and min and max widths, tab-separated. Do not report scores to standard error.
-verbose Send detailed stats to standard error, reporting frequency of each score. No extra information is reported.

Examples

In each example, the sequence file is the same as the postive set used for psp-gen. While there is no check that you have done this because MEME runs as a separate program (other than that the sequence names need to match, as do the number of sites), using a primary set different to the sequence set MEME searches is unlikely to be useful.

Generate a discriminative prior for a DNA sequence set (pos-dna.fasta) likely to contain TF binding sites, given another DNA sequence set unlikely to contain TF binding sites (neg-dna.fasta):
psp-gen -pos pos-dna.fasta -neg neg-dna.fasta > dna.psp

Use these position-specific priors in MEME, searching both strands:
meme pos-dna.fasta -psp dna.psp -oc meme_out_dna -revcomp

Generate a discriminative prior for a protein sequence set (pos-prot.fasta) associated with a specific function, given another protein sequence set unlikely to relate to that function (neg-prot.fasta):
psp-gen -pos pos-prot.fasta -neg neg-prot.fasta -protein -triples -maxrange > prot.psp

Use these position-specific priors in MEME:
meme pos-prot.fasta -psp prot.psp -oc meme_out_prot

View usage information:
psp-gen -h

Detailed description

The psp-gen program generates a position-specific prior (PSP) in a data file, which can be used as an additional input to MEME to bias the search to sites more likely to contain a motif (or motifs) of interest. Without using a PSP, MEME assigns the same prior probability to all sites. A PSP is generated by psp-gen using two sequence sets, each in FASTA format:

The basic principle of a discriminative prior is to start from the question:

What fraction of all instances of a w-wide subsequence (or w-mer) across both sequence sets occur in the primary set?

The intuition is that any w-mer that is common in the primary set but not the control set could form part of motif of interest, based on the way the two sets are chosen.

PSP calculation starts from the following equation based on the D-prior described in the Narlikar et al. RECOMB 2007 paper "Nucleosome Occupancy Information Improves de novo Motif Discovery". We count the number of instances Xi,j(w) of the subsequence w wide starting at position j in sequence Xi of the primary set, and similarly count Yi,j(w) in the control set, then calculate for each site D i , j = X i , j ( w ) X i , j ( w ) + Y i , j ( w ) . As pseudocounts, we add 1 to the enumerator, and 1 + |X|/|Y| to the denominator. The purpose of the pseudocounts is to avoid giving a high score to infrequent subsequences that are not (or rarely) found in the control set. Once we have scored all sites for a given value of w, we scale scores to a range that can be set on the command line, with the default range [0.1,0.9].

We repeat scoring for the whole range of motif widths under consideration, and choose a "best" width using one of two methods: either we choose the value of w that maximizes the score over all sequences, or the value of w that maximizes the difference between the maximum and minimum score over all sequences. We also allow an option of using spaced triples instead of whole words to match subsequences. See Variations below for more details on triples and choosing the "best" width.

Once we have scored every available site and chosen a specific value of w, we convert the scores to probabilities, based on the proportionalities Pr ( Z i = 0 ) 1   and   Pr ( Z i = j ) D i , j 1 - D i , j , where Zi = j means that the site for a motif in sequence Xi is found at location j (with j = 0 signifying there is no site in the given sequence), and Di,j is the score for sequence Xi at the site starting at location j, numbered from 1. We convert these proportionalities to a PSP by normalizing them to add up to 1 over each sequence.

Running with MEME

Since MEME runs as a separate program, there is no way to check that the primary set is the same as the sequence set given to MEME. There are however checks that every name in the PSP file matches a name in the MEME sequence set, and that the number of sites in the PSP file for a particular sequence name matches the number of sites in the sequence file seen by MEME. While it does not make sense to use a completely different sequence set for PSP generation and for running MEME (with the resulting PSP file), MEME can use a PSP file that does not provide probabilities for every sequence. Any sequences not in the PSP file are given a uniform prior probability. MEME documentation includes more detail on how PSPs are used.

Variations

For each variation, "default" means that you do not need to specify a command-line switch to turn on that setting.

Performance scalability

When you use exact-word (w-mer) search, psp-gen run time scales linearly in each of |X|, |Y| and maxw-minw. Formally, time is O((|X| + |Y|)(wmax - wmin)). For triples, psp-gen scales linearly in each of |X| and |Y|, and quadratically in wmax, so run time is O((|X| + |Y|)(wmax)2).

Extra memory in both cases scales linearly with |X| + |Y|.

Since triples have a higher-order term in time complexity, we predict run time scalability in more detail than for exact word match. Adding in constants and a linear term, time complexity for triples is approximated by (where W = wmax-2, the number of variations in width used by the algorithm):

t = (k1|X| + k2|Y|)W2 + k3(|X| + |Y|)W

On a real machine (iMac, 3.06GHz Intel Core 2 Duo, 6MB L2 cache, 4GB RAM) the following is a good predictor of the upper bound on run time for likely values of |X|, |Y| and range of widths. Scale any results according to the speed of the actual machine on which you intend to run psp-gen. The formula reflects runs on randomly generated data; real data will likely result in faster runs because random data defeats some optimizations in psp-gen.

t=(4.37x10-6|X|+ 2.46x10-06|Y|)W2 + 6.25x10-06(|X| + |Y|)W s

For example, for |X| = |Y| = 60,000 and wmax = 50 (i.e., W = 48), estimated run time = 980s, the same in this case as measured run time. The following command line reproduces this run:

      time psp-gen -neg pos60k.fasta -pos neg60k.fasta -protein -minw 3 -maxw 50 -triples > priors60kmaxw50.psp
      

Known problems