Ок. Вся программа:
import random
from scipy.stats import chisquare
size = 160000
nuc = 'ATGC'
def rand_seq(nuc, size):
seq = ''.join([random.choice(nuc) for N in xrange(size)])
return seq
def combine(length):
strings = [N for N in 'ATGC']
L=1
while L < length:
strings = [string+N for N in 'ATGC' for string in strings] #С каждым циклом удлиняем на одну букву, создавая четыре слова из одного
L= L+1
return strings
def words(seq, word_length):
words = [N*word_length for N in 'ATGC'] #Учесть все одинаковые, которые иначе никак не создадутся
for word in combine(word_length):
word = ''.join(map(str, word))
words.append(word) #Добавляем слово
words.append(word[::-1]) #Добавляем перевернутое слово, чтобы все учесть
return words
def get_word_dict(seq):
word_dict = {}
for length in range(1, 6):
for word in words(seq, length):
word_dict[word] = seq.count(word) #Сам подсчет
return word_dict
#Собственно, начинаем использовать вышенаписанные функции
seq = rand_seq(nuc, size)
word_dict = get_word_dict(seq)
for key in sorted(word_dict.iterkeys()):
word_dict[key] = (len(key), word_dict[key]) #В значение ключа записываем его длину для удобства, и рассчитанное число.
length = 2
numbers = []
labels = []
for key in sorted(word_dict.iterkeys()):
if len(key) == length:
numbers.append(word_dict[key][1])
labels.append(key)
named = zip(numbers, labels)
named = sorted(named)
chisq = chisquare(numbers)
print 'P-value for uniform ditribution', chisq[1]
pylab.bar([i for i in range(len(numbers))], numbers)
pylab.xticks([i+0.4 for i in range(len(numbers))], labels, rotation='vertical')
pylab.show()
Идея в том, чтобы сгенерировать все варианты слов из этих 4х букв определенной длины, потом посчитать, сколько каждое из них встречается.