[Biopython-dev] [Bug 3135] New: Wrong instance length bug in MEME parser
bugzilla-daemon at portal.open-bio.org
bugzilla-daemon at portal.open-bio.org
Fri Sep 3 09:49:03 EDT 2010
http://bugzilla.open-bio.org/show_bug.cgi?id=3135
Summary: Wrong instance length bug in MEME parser
Product: Biopython
Version: 1.55
Platform: PC
OS/Version: Linux
Status: NEW
Severity: normal
Priority: P2
Component: Main Distribution
AssignedTo: biopython-dev at biopython.org
ReportedBy: johnbaronreid at netscape.net
The MEME parser in biopython 1.55 seems to incorrectly set the length of the
first instance of a motif to 0. Here is an example:
#Sequence, start, length, site
Motif: E-value: 0.000010
seq_3, 213, 0, AGGTGACAGAG
seq_1, 146, 11, AGGTGACAGAG
seq_0, 490, 11, AGGTGACAGAG
seq_0, 83, 11, AGGTGACAGAG
seq_0, 388, 11, AGGAAACAGAG
seq_1, 422, 11, AGGGGACAGAG
seq_1, 79, 11, TGGAGACAGAG
seq_0, 281, 11, TGGGGACAGAG
seq_0, 16, 11, TAGAGACAGAG
seq_1, 228, 11, TTGTGACAGAG
seq_4, 156, 11, AGGGGACAGGG
seq_0, 348, 11, AGGAAAGAGAA
seq_0, 374, 11, AGGAATGAGAG
seq_5, 22, 11, GGGAAACTGAG
seq_3, 486, 11, AAGGGAGTGAG
Here's the code that generated the above:
from Bio.Motif.Parsers.MEME import MEMEParser
import cStringIO
meme_output = cStringIO.StringIO("""
********************************************************************************
MEME - Motif discovery tool
********************************************************************************
MEME version 4.3.0 (Release date: Sat Sep 26 01:51:56 PDT 2009)
For further information on how to interpret these results or to get
a copy of the MEME software please access http://meme.nbcr.net.
This file may be used as input to the MAST algorithm for searching
sequence databases for matches to groups of motifs. MAST is available
for interactive use and downloading at http://meme.nbcr.net.
********************************************************************************
********************************************************************************
REFERENCE
********************************************************************************
If you use this program in your research, please cite:
Timothy L. Bailey and Charles Elkan,
"Fitting a mixture model by expectation maximization to discover
motifs in biopolymers", Proceedings of the Second International
Conference on Intelligent Systems for Molecular Biology, pp. 28-36,
AAAI Press, Menlo Park, California, 1994.
********************************************************************************
********************************************************************************
TRAINING SET
********************************************************************************
DATAFILE= /home/john/Data/Tompa-data-set/Real/hm22r.fasta
ALPHABET= ACGT
Sequence name Weight Length Sequence name Weight Length
------------- ------ ------ ------------- ------ ------
seq_0 1.0000 500 seq_1 1.0000 500
seq_2 1.0000 500 seq_3 1.0000 500
seq_4 1.0000 500 seq_5 1.0000 500
********************************************************************************
********************************************************************************
COMMAND LINE SUMMARY
********************************************************************************
This information can also be useful in the event you wish to report a
problem with the MEME software.
command: meme /home/john/Data/Tompa-data-set/Real/hm22r.fasta -maxsize 1000000
-oc output/run_dataset/Tompa/hm22r/Real -dna -mod anr -revcomp -print_starts
-maxiter 1000 -minw 8 -maxw 20 -minsites 2 -nmotifs 1
model: mod= anr nmotifs= 1 evt= inf
object function= E-value of product of p-values
width: minw= 8 maxw= 20 minic= 0.00
width: wg= 11 ws= 1 endgaps= yes
nsites: minsites= 2 maxsites= 30 wnsites= 0.8
theta: prob= 1 spmap= uni spfuzz= 0.5
global: substring= yes branching= no wbranch= no
em: prior= dirichlet b= 0.01 maxiter= 1000
distance= 1e-05
data: n= 3000 N= 6
strands: + -
sample: seed= 0 seqfrac= 1
Letter frequencies in dataset:
A 0.195 C 0.305 G 0.305 T 0.195
Background letter frequencies (from dataset with add-one prior applied):
A 0.195 C 0.305 G 0.305 T 0.195
********************************************************************************
********************************************************************************
MOTIF 1 width = 11 sites = 15 llr = 159 E-value = 9.8e-006
********************************************************************************
--------------------------------------------------------------------------------
Motif 1 Description
--------------------------------------------------------------------------------
Simplified A 71:439:9:91
pos.-specific C ::::::8::::
probability G 18a37:2:a19
matrix T 31:3:1:1:::
bits 2.4
2.1 *
1.9 * * *
1.6 * * ***
Relative 1.4 * * ****
Entropy 1.2 * * * ****
(15.3 bits) 0.9 *** *******
0.7 ***********
0.5 ***********
0.2 ***********
0.0 -----------
Multilevel AGGAGACAGAG
consensus T TA G
sequence G
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
Motif 1 sites sorted by position p-value
--------------------------------------------------------------------------------
Sequence name Strand Start P-value Site
------------- ------ ----- --------- -----------
seq_3 - 213 4.54e-07 GGCCTTTGGA AGGTGACAGAG
GCGCGGCCAC
seq_1 - 146 4.54e-07 CCCAACAGGA AGGTGACAGAG
GTGGCTCTGG
seq_0 + 490 4.54e-07 AAAACAGCAG AGGTGACAGAG
seq_0 - 83 4.54e-07 CCCAGCAGGA AGGTGACAGAG
GTGGCTCTGG
seq_0 + 388 5.99e-07 ATGAGAGGAG AGGAAACAGAG
CTTCCTGGAC
seq_1 + 422 1.10e-06 ATGAGAGGGG AGGGGACAGAG
GACACCTGAA
seq_1 + 79 1.33e-06 TTGGTGGTAC TGGAGACAGAG
GGCTGGTCCC
seq_0 + 281 3.17e-06 CCTCCCCTGA TGGGGACAGAG
GTCTCATCAG
seq_0 + 16 5.72e-06 CTGGTGACAC TAGAGACAGAG
GGCTGGTCCC
seq_1 - 228 1.18e-05 TTATTTTCCT TTGTGACAGAG
AAACCCAGCA
seq_4 + 156 2.07e-05 TCAAGTCCCA AGGGGACAGGG
AGCAGAAGGG
seq_0 + 348 2.47e-05 GTAGACAGAA AGGAAAGAGAA
AGTAAGGACA
seq_0 + 374 3.14e-05 GGACAAAGGT AGGAATGAGAG
GAGAGGAAAC
seq_5 - 22 4.53e-05 CTCTTGTGTA GGGAAACTGAG
CACGGGGAAC
seq_3 + 486 5.02e-05 CGCCAATGGG AAGGGAGTGAG TGCC
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
Motif 1 block diagrams
--------------------------------------------------------------------------------
SEQUENCE NAME POSITION P-VALUE MOTIF DIAGRAM
------------- ---------------- -------------
seq_3 5e-05 212_[-1]_262_[+1]_4
seq_1 1.2e-05 78_[+1]_56_[-1]_71_[-1]_183_[+1]_68
seq_0 3.2e-06 15_[+1]_56_[-1]_187_[+1]_56_[+1]_
15_[+1]_3_[+1]_91_[+1]
seq_4 2.1e-05 155_[+1]_334
seq_5 4.5e-05 21_[-1]_468
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
Motif 1 in BLOCKS format
--------------------------------------------------------------------------------
BL MOTIF 1 width=11 seqs=15
seq_3 ( 213) AGGTGACAGAG 1
seq_1 ( 146) AGGTGACAGAG 1
seq_0 ( 490) AGGTGACAGAG 1
seq_0 ( 83) AGGTGACAGAG 1
seq_0 ( 388) AGGAAACAGAG 1
seq_1 ( 422) AGGGGACAGAG 1
seq_1 ( 79) TGGAGACAGAG 1
seq_0 ( 281) TGGGGACAGAG 1
seq_0 ( 16) TAGAGACAGAG 1
seq_1 ( 228) TTGTGACAGAG 1
seq_4 ( 156) AGGGGACAGGG 1
seq_0 ( 348) AGGAAAGAGAA 1
seq_0 ( 374) AGGAATGAGAG 1
seq_5 ( 22) GGGAAACTGAG 1
seq_3 ( 486) AAGGGAGTGAG 1
//
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
Motif 1 position-specific scoring matrix
--------------------------------------------------------------------------------
log-odds matrix: alength= 4 w= 11 n= 2940 bayes= 6.7534 E= 9.8e-006
177 -1055 -219 45
-55 -1055 139 -155
-1055 -1055 171 -1055
103 -1055 -19 77
45 -1055 127 -1055
226 -1055 -1055 -155
-1055 139 -61 -1055
215 -1055 -1055 -55
-1055 -1055 171 -1055
226 -1055 -219 -1055
-155 -1055 161 -1055
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
Motif 1 position-specific probability matrix
--------------------------------------------------------------------------------
letter-probability matrix: alength= 4 w= 11 nsites= 15 E= 9.8e-006
0.666667 0.000000 0.066667 0.266667
0.133333 0.000000 0.800000 0.066667
0.000000 0.000000 1.000000 0.000000
0.400000 0.000000 0.266667 0.333333
0.266667 0.000000 0.733333 0.000000
0.933333 0.000000 0.000000 0.066667
0.000000 0.800000 0.200000 0.000000
0.866667 0.000000 0.000000 0.133333
0.000000 0.000000 1.000000 0.000000
0.933333 0.000000 0.066667 0.000000
0.066667 0.000000 0.933333 0.000000
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
Motif 1 regular expression
--------------------------------------------------------------------------------
[AT]GG[ATG][GA]A[CG]AGAG
--------------------------------------------------------------------------------
Time 3.78 secs.
********************************************************************************
********************************************************************************
SUMMARY OF MOTIFS
********************************************************************************
--------------------------------------------------------------------------------
Combined block diagrams: non-overlapping sites with p-value < 0.0001
--------------------------------------------------------------------------------
SEQUENCE NAME COMBINED P-VALUE MOTIF DIAGRAM
------------- ---------------- -------------
seq_0 4.45e-04
15_[+1(5.72e-06)]_56_[-1(4.54e-07)]_187_[+1(3.17e-06)]_56_[+1(2.47e-05)]_15_[+1(3.14e-05)]_3_[+1(5.99e-07)]_91_[+1(4.54e-07)]
seq_1 4.45e-04
78_[+1(1.33e-06)]_56_[-1(4.54e-07)]_71_[-1(1.18e-05)]_183_[+1(1.10e-06)]_68
seq_2 2.03e-01 500
seq_3 4.45e-04
212_[-1(4.54e-07)]_262_[+1(5.02e-05)]_4
seq_4 2.01e-02 155_[+1(2.07e-05)]_334
seq_5 4.34e-02 21_[-1(4.53e-05)]_468
--------------------------------------------------------------------------------
********************************************************************************
********************************************************************************
Stopped because nmotifs = 1 reached.
********************************************************************************
CPU: john-dell
********************************************************************************
""")
parser = MEMEParser()
parsed = parser.parse(meme_output)
print '#Sequence, start, length, site'
for motif in parsed.motifs:
print 'Motif: E-value: %f' % motif.evalue
for instance in motif.instances:
print "%10s, %5d, %5d, %s" % (
instance.sequence_name,
instance.start,
instance.length,
str(instance),
)
#assert instance.length == motif.length
--
Configure bugmail: http://bugzilla.open-bio.org/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are the assignee for the bug, or are watching the assignee.
More information about the Biopython-dev
mailing list