"""
Using the motif module of biopython
Note that the python tutorial does not contain information about the
Motif package. Documentation is available in the directory
Doc/cookbook/motif of the biopython distribution.
"""
from Bio import Motif
# Let's create an instance of the E2F1 motif (downloaded from the
# jaspar database):
motif=Motif.read(open("MA0024.1.pfm"), "jaspar-pfm")
# the format method displays the motif in a variety of formats:
print motif.format('transfac')
# You can also get a "logo" image that represents the motif:
motif.weblogo('e2f1.png')
# The length method does what we would expect:
print len(motif)
# The jaspar format for a motif specifies a matrix of counts: how many times
# each nucleotide occurs in each position of the motif.
# This information is stored in the 'counts' instance variable of a motif.
# Let's see what design choice the biopython folks did for storing this
# information:
motif.counts
# The matrix of counts is converted into a matrix of frequencies.
# This is called a position weight matrix (PWM).
# The PWM associated with a motif can be accessed as:
motif.pwm()
# Notice that the pwd assigns a non-zero frequency to every nucleotide at every
# position, even if none have been observed. Those are the result of what are
# called "pseuodo-counts" which will be explained next in the context of scoring
# a sequence with respect to the motif.
# The consensus sequence is the mostly like sequence represented by the motif.
# It can be generated by:
motif.consensus()
# Let's address the question of how well does a sequence matches the motif.
# We'll do that by assigning a score that is a sum of lod-odds ratios.
# A motif instance has a variable associated with it called _log_oddds:
motif._log_odds
# The maximum score that can be achieved for the motif
motif.max_score()
# By default the background is assumed to be uniform:
motif.background
def compute_background(file_name) :
from Bio import SeqIO
background = {'A':0.0, 'C':0.0, 'G':0.0, 'T':0.0}
count = 0.0
for seq_record in SeqIO.parse(file_name, "fasta") :
seq = seq_record.seq
for nucl in seq :
if nucl in background :
background[nucl] += 1
count += 1
for nucl in background :
background[nucl] /= count
return background
# Let's see that the max score changes when we modify the background:
motif.max_score()
motif.background = compute_background('human_promoter500.fasta')
motif._log_odds_is_current = False
motif.max_score()
# Let's look for matches of the motif in DNA sequences
from Bio.Seq import Seq
seq = Seq("AATTTCGCGCAA")
hits = motif.search_pwm(seq, threshold = 0.7 * max_score)
# Note that the search_pwm method returns an iterator.
# This iterator returns tuples that contain the position and score of all the hits:
for pos, score in hits :
print pos, score
# Let's now search for motif hits in the promoter regions of human genes
from Bio import Motif
motif=Motif.read(open("MA0024.1.pfm"), "jaspar-pfm")
max_score = motif.max_score()
from Bio import SeqIO
for seq_record in SeqIO.parse('human_promoter500.fasta', "fasta") :
seq = seq_record.seq
id = seq_record.id
search_results = motif.search_pwm(seq, threshold=0.7 * max_score)
for pos, score in search_results :
print id, pos, score, '+'
# We need to understand better what it means when there is a hit
# on the negative strand
from Bio import Motif
motif=Motif.read(open("MA0024.1.pfm"), "jaspar-pfm")
max_score = motif.max_score()
from Bio import SeqIO
for seq_record in SeqIO.parse('human_promoter500.fasta', "fasta") :
seq = seq_record.seq
id = seq_record.id
rc = seq.reverse_complement()
search_results = motif.search_pwm(seq, threshold=0.7 * max_score)
for pos, score in search_results :
if pos > 0 :
print id, pos, score, '+'
print seq[pos:pos + len(motif)]
else :
print id, pos, score, '-'
print seq[-pos : -pos + len(motif)]
search_results = motif.search_pwm(rc, threshold=0.7 * max_score, both=False)
for pos, score in search_results :
print id, pos, score, '-'
print rc[pos:pos + len(motif)]