Задача:
При каком минимальном числе k существует нуклеотидная последовательность длины k, не встречающаяся в геноме человека? Сколько таких последовательностей при найденном k?
Подпоследовательности длины k, содержащиеся в биологической последовательности, в биоинформатике называют k-мерами (k-mers). Обычно этот термин используется в контексте вычислительной биологии и анализа последовательностей, в которых k-меры состоят из нуклеотидов. В данной задаче предлагается определить, при какой минимальной длине k-мера существуют k-меры, не встречающиеся в геноме человека.
Для решения использовался геном человека GRCh38, который предварительно был обработан через командную строку так, чтобы файл содержал только заглавные буквы A, T, G или C.
Код на Python, решающий задачу:
import time
from itertools import product
import re
seq = ''
with open('GRCh38_genome_new.fna', 'r') as file:
# write the genome into a string, removing control characters
seq=file.read().replace('\n', '').replace('\r', '')
print('The file has been read!')
for k in range(1, 16):
print('Starting k: ', k)
start = time.time()
match = 0
for word in product('ATGC', repeat=k):
# go through the k-mers that are generated using product
a = re.search(''.join(word), seq)
# looking for k-mer in genome
# move to the next iteration, if there is no match
if a == None:
continue
# add 1 to counter, if there is a match
match += 1
# if there is no matches for al least one k-mer for this k,
# then we have found what we were looking for
if match < 4**k:
print('Min k:', k)
print('The amount of missing k-mers:', 4**k - match)
end = time.time()
print('k =', k, ':', end - start)
break
end = time.time()