[Biopython] codon usage

Anirban Bhattachariya anbhat at utu.fi
Fri Jan 15 16:07:51 UTC 2010


Hi,

I found this script on "http://www.pasteur.fr/recherche/unites/sis/formation/python/exercises/seqrandom_count_codons_plot.py" which is supposed to count codon usage and plot them in a bar plot but this not working since some of the modules used in the script does not exist anymore.

How can modify the script to make it usable or is there a better way to do that?

here is the code:

import Bio.Fasta
from sys import *
from string import *
from dna import codons
from mutateseq import mutateseq

file = argv[1]
handle = open(file)
it = Bio.Fasta.Iterator(handle, Bio.Fasta.SequenceParser())
count = {}
count_random = {}
seq = it.next()
while seq:
	for codon in codons(seq.seq.tostring()):
		if count.has_key(codon):
			count[codon] += 1
		else:
			count[codon] = 0
	mutableseq = seq.seq.tomutable()
	mutateseq(mutableseq,span=1000,p=0.1)
	for codon in codons(mutableseq.tostring()):
		if count_random.has_key(codon):
			count_random[codon] += 1
		else:
			count_random[codon] = 0

	seq = it.next()

handle.close()

#--------------------------------------------------------
# bar charts of codons frequencies
#   - for legibility, 2 charts are built
#   - both random and normal frequencies are dsplayed

from tkplot import *
from Numeric import *

def codon_sort(a,b):
	if a < b:
		return -1
	elif a > b:
		return 1
	else:
		return 0

for codon in count.keys():
	if not count_random.has_key(codon):
		count_random[codon] = 0

for codon in count_random.keys():
	if not count.has_key(codon):
		count[codon] = 0

labels=count.keys()
labels.sort(codon_sort)

w1=window(plot_title='Count codons',width=1000)
y=array(count.values())[:len(count)/2]
x=arange(len(y)+1)
w1.bar(y,x,label=labels[:len(count)/2])

w2=window(plot_title='Count codons(2)',width=1000)
y=array(count.values())[(len(count)/2)+1:]
x=arange(len(y)+1)
w2.bar(y,x,label=labels[(len(count)/2)+1:])

y=array(count_random.values())[:len(count_random)/2]
x=arange(len(y)+1)
w1.bar(y,x,label=labels[:len(count_random)/2])

y=array(count_random.values())[(len(count_random)/2)+1:]
x=arange(len(y)+1)
w2.bar(y,x,label=labels[(len(count_random)/2)+1:])




Regards,
Anirban




More information about the Biopython mailing list