Решение задачи о минимальном размере k-мера, которого нет в геноме, на Python

Задача:
При каком минимальном числе 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()

Leave a Reply