This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Wrong Alphabet In Biopython.Motifs When Trying To Pssm.Calculate

I'm trying to find motif occurences in the upstream:

upstream = Seq(mysequence, IUPAC.unambiguous_dna)
pssm.calculate(upstream)

get error:

ValueError: Wrong alphabet! Use only with DNA motifs

if I check alphabets:

>>>pssm.alphabet
IUPACUnambiguousDNA()
>>>upstream.alphabet
IUPACUnambiguousDNA()

The source code where is fails:

  if self.alphabet!=IUPAC.unambiguous_dna:
        raise ValueError("Wrong alphabet! Use only with DNA motifs")

What could I do wrong? What's really the difference between IUPACUnambiguousDNA() and IUPAC.unambiguous_dna?

biopython motif
• 4,783 views • link
updated by Istvan Albert Istvan Albert • written by Jelena_bioinf • Jelena_bioinf

1 answer

The problem is with biopython - it tests instances for equality instead of testing them for having the same class.

IUPACUnambiguousDNA() and IUPAC.unambiguous_dna may actually represent different instances of the same class you can see that by doing something like:

print id(pssm.alphabet)
print id(upstream.alphabet)

These will not test as being equal when checking with the = operator.

When you create the pssm object make sure to use the the IUPAC.unambiguous_dna in the constructor.

• 203 views • link
written by Istvan Albert Istvan Albert

Indeed 

   >>id(upstream.alphabet)
   4414891728
   >>id(pssm.alphabet)
   4428032464

I'm not sure I use alphabet in the constructor properly: 

   instances = []
           for line in alignment:
               instances.append(Seq(line, IUPAC.unambiguous_dna))
           m = motifs.create(instances)
           pwm = m.counts.normalize(pseudocounts=PSEUDOCOUNT)
           background = {'A':0.18,'C':0.32,'G':0.32,'T':0.18}
           pssm = pwm.log_odds(background)

What I also tried to do and it didn't work:

   upstreamList = [Seq(u, pssm.alphabet) for u in upstreamList]

previously it was

   upstreamList = [Seq(u, IUPAC.unambiguous_dna) for u in upstreamList]

The likely reason is that I create pssm in one run and cPickle it and then unpickle when I need to apply. It works fine if I create pssm in the same run as I apply, but constructing matrix takes time, so I'd like to figure out how to pickle it if it's possible at all.

• 0 views • link
written by Jelena_bioinf • Jelena_bioinf

Pickling and unpickling will result in a new instance of the alphabet class, and you'll hit the bug again. You'll have to wait for the next Biopython release where this is fixed, or apply the fix locally. See the mailing list discussion linked to above.

• 0 views • link
written by Peter Peter

Log in to answer this question.

More posts like this
  • pssm.search from Bio.motifs not working with amino acid sequences
    written by stefano • stefano
    I am trying to find a motif out of a list of amino acid sequences.I am using the package `Bio.motifs`. I followed the documentation as …
  • Galaxy "Stitch Maf Blocks" Produced Null Bases
    written by lam • lam
    <p>problem: "stitch MAF blocks" produced NULL bases </p> <p>steps: From <a href='http://genome.ucsc.edu/cgi-bin/das/hg19/dna?segment=chrY:23310464,23331995'>http://genome.ucsc.edu/cgi-bin/das/hg19/dna?segment=chrY:23310464,23331995</a></p> <p>we got the following bases:</p> <pre><code>agtttacccaactcacaggtatttgaggatacaaacccatgactgggcttggctttaaaaattcctatctgaaatttctagt.......... </code></pre> <p>We want to do multiple-species alignment …
  • Assign Pwm In Variable In Biopython
    written by bioLife • bioLife
    <p>I have a PWM (position weight matrix with the site specific frequencies) for a motif, how can I use it in Biopython? In module Motif …
  • Concatenate Sequences With Bioperl
    written by maczerek • maczerek
    <p>Hello, I've got some trivial (I guess) question. I try to write features of a sequence in single output sequences. My code looks like that:</p> …
  • Bioperl: Finding Restriction Enzyme Sites
    written by cindylanzao • cindylanzao
    <pre><code>#!usr/bin/perl -w use Bio::DB::GenBank; use Bio::Restriction::EnzymeCollection; use Bio::Restriction::Analysis; my $gb=new Bio::DB::GenBank(); my $seq=$gb-&gt;get_Seq_by_acc('AF303112'); print $seq-&gt;seq,"\n"; my $all_collection=Bio::Restriction::EnzymeCollection; my $my_enzyme=$all_collection-&gt;get_enzyme('EcoRI'); my $analysis=Bio::Restriction::Analysis-&gt;new(-seq=&gt;$seq); @fragment=$analysis-&gt;fragments($my_enzyme); print "@fragment\n"; my …
  • Biopython Motif To Find Patterns In Sequences
    written by dami.gupta • dami.gupta
    <p>I have a sequence file where for each read, I want to check for instances of a pattern(s). I will demultiplex the file into several …
  • Reverse Complement Of A Sequence Raises Assertionerror: Invalid Alphabet Found
    written by Jelena_bioinf • Jelena_bioinf
    <p>I read a fasta-formatted genome sequence and try to get its reverse complement:</p> <pre><code>genomeSeq = FastaIO.FastaIterator(genomeHandle, IUPACUnambiguousDNA).next() genomeSeq.seq.reverse_complement() </code></pre> <p>but it doesn't work and I …
  • Looking For Tools Tp Draw Motif Logo
    written by Yunfei Li Yunfei Li
    <p>I try to create some motif logos. Now I follow biopython example and write script like following, it works well when dataset is small, but …
  • Fire-Pro Gene Expression Profile
    written by Dhillonv10 • Dhillonv10
    <p>I'm trying to use Fire-pro to find motifs in linear motifs taken from ELM database. However, using the software, the syntax comes as follows: </p> …
  • Import Fasta Sequences To A Motif
    written by RossCampbell RossCampbell
    <p>I need to construct a PWM from every sequence in a fasta file, using biopython. The way I'm trying to do this is to import …