[Biopython] codon usage
Anirban Bhattachariya
anbhat at utu.fi
Fri Jan 15 11:07:51 EST 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