Вы можете использовать умные методы, как предложено в других ответах, но я построю решение, исходя из вашего кода, который почти работает: ваша проблема в том, что каждый раз, когда вы делаете mydict[b] = 1, вы сбрасываете счетчик * От 1002 * до 1.
Минимальное исправление
Вы можете решить эту проблему, протестировав наличие ключа, если нет, создать запись в dict, затем увеличить значение, но в python есть более удобные инструменты.
Минимальным изменением в вашем коде будет использование defaultdict(int) вместо dict. Всякий раз, когда встречается новый ключ, предполагается, что он имеет соответствующее значение по умолчанию для целого числа: 0. Таким образом, вы можете увеличивать значение вместо сброса:
from collections import defaultdict
def cnt(seq):
# instanciate a defaultdict that creates ints when necessary
mydict = defaultdict(int)
if len(seq) % 3 == 0:
a = [x for x in seq]
for i in range(len(seq)//3):
b = ''.join(a[(0+3*i):(3+3*i)])
for base1 in ['A', 'T', 'G', 'C']:
for base2 in ['A', 'T', 'G', 'C']:
for base3 in ['A', 'T', 'G', 'C']:
triplet = base1 + base2 + base3
if b == triplet:
# increment the existing count (or the default 0 value)
mydict[b] += 1
for key in sorted(mydict):
print("%s: %s" % (key, mydict[key]))
else:
print("Error")
Работает как нужно:
cnt("ACTGGCACT")
ACT: 2
GGC: 1
Некоторые возможные улучшения
Теперь давайте попробуем немного улучшить ваш код.
Во-первых, как я писал в комментариях, давайте избежим ненужного преобразования вашей последовательности в список и используем лучшее имя переменной для текущего подсчитанного кодона:
from collections import defaultdict
def cnt(seq):
mydict = defaultdict(int)
if len(seq) % 3 == 0:
a = [x for x in seq]
for i in range(len(seq)//3):
codon = seq[(0+3*i):(3+3*i)]
for base1 in ['A', 'T', 'G', 'C']:
for base2 in ['A', 'T', 'G', 'C']:
for base3 in ['A', 'T', 'G', 'C']:
triplet = base1 + base2 + base3
if codon == triplet:
mydict[codon] += 1
for key in sorted(mydict):
print("%s: %s" % (key, mydict[key]))
else:
print("Error")
Теперь давайте упростим часть вложенного цикла, пробуя все возможные кодоны, заранее сгенерировав набор возможных кодонов:
from collections import defaultdict
from itertools import product
codons = {
"".join((base1, base2, base3))
for (base1, base2, base3) in product("ACGT", "ACGT", "ACGT")}
def cnt(seq):
mydict = defaultdict(int)
if len(seq) % 3 == 0:
a = [x for x in seq]
for i in range(len(seq)//3):
codon = seq[(0+3*i):(3+3*i)]
if codon in codons:
mydict[codon] += 1
for key in sorted(mydict):
print("%s: %s" % (key, mydict[key]))
else:
print("Error")
Теперь ваш код просто игнорирует триплеты, которые не являются действительными кодонами. Возможно, вместо этого вам следует выдать предупреждение:
from collections import defaultdict
from itertools import product
codons = {
"".join((base1, base2, base3))
for (base1, base2, base3) in product("ACGT", "ACGT", "ACGT")}
def cnt(seq):
mydict = defaultdict(int)
if len(seq) % 3 == 0:
a = [x for x in seq]
for i in range(len(seq)//3):
codon = seq[(0+3*i):(3+3*i)]
# We count even invalid triplets
mydict[codon] += 1
# We display counts only for valid triplets
for codon in sorted(codons):
print("%s: %s" % (codon, mydict[codon]))
# We compute the set of invalid triplets:
# the keys that are not codons.
invalid = mydict.keys() - codons
# An empty set has value False in a test.
# We issue a warning if the set is not empty.
if invalid:
print("Warning! There are invalid triplets:")
print(", ".join(sorted(invalid)))
else:
print("Error")
Более причудливое решение
Теперь более изящное решение, использующее cytoolz (вероятно, его необходимо установить, поскольку оно не является частью обычных дистрибутивов Python: pip3 install cytoolz, если вы используете pip):
from collections import Counter
from itertools import product, repeat
from cytoolz import groupby, keymap, partition
# To make strings out of lists of strings
CAT = "".join
# The star "extracts" the elements from the result of repeat,
# so that product has 3 arguments, and not a single one
codons = {CAT(bases) for bases in product(*repeat("ACGT", 3))}
def cnt(seq):
# keymap(CAT, ...) transforms the keys (that are tuples of letters)
# into strings
# if len(seq) is not a multiple of 3, pad="-" will append "-"
# to complete the last triplet (which will be an invalid one)
codon_counts = keymap(CAT, Counter(partition(3, seq, pad="-")))
# separate encountered codons into valids and invalids
codons_by_validity = groupby(codons.__contains__, codon_counts.keys())
# get allows to provide a default value,
# in case one of the categories is not present
valids = codons_by_validity.get(True, [])
invalids = codons_by_validity.get(False, [])
# We display counts only for valid triplets
for codon in sorted(valids):
print("%s: %s" % (codon, codon_counts[codon]))
# We issue a warning if there are invalid codons.
if invalids:
print("Warning! There are invalid triplets:")
print(", ".join(sorted(invalids)))
Надеюсь, это поможет.