Алгоритм срыва вперед и назад дополнения последовательности ДНК в python?

Question

У меня есть коллекция файлов fasta, содержащих много фрагментов последовательности ДНК. Я пытаюсь подсчитать суммарные вхождения k-mers, которые можно найти в каждом файле. Хорошей частью подсчета k-mers является то, что можно создать один массив размером 4 ** k, где k - размер используемого k-mer. Файлы последовательности, которые я обрабатываю, представляют собой короткие последовательности считывания из машин для секвенирования нового поколения, и поэтому предполагать, что все чтения из 5 '-> 3' не могут быть выполнены. Лучший способ решить эту проблему - сопоставить все k-mers, наблюдаемые с одним индексом для последовательных и обратных последовательностей.

Желаемая отображение:

с к = 2 & исходного индекса для массива равен 0

строка = «аа»; карты для указания -> 0

string = 'tt'; map to index -> 0

string = 'at'; карты индексировать -> 1

По руке я был в состоянии понять, что массив для всех меров с распадом прямых и обратным комплементы будет иметь длину 10, а конкретные показатели будут представлять следующие меры: А.А. , AT, AC, AG, TA, TC, TG, CC, CG, GC

У меня возникли проблемы с рассмотрением обобщенного алгоритма, чтобы узнать количество возможных мер для больших размеров k. Сколько ячеек должно быть выделено в массиве count?

В моем существующем коде я придумал эти три функции для обработки фрагментов, создания обратного дополнения и сопоставления mer (или обратного дополнения) с индексом.

Эта первая функция примет строку mer и возвращает индекс, который относится к mer в массиве размером 4 ** k.

def mer_index_finder(my_string, mer_size): 
    # my_string = my_string.lower() 
    char_value = {} 
    char_value["a"] = 0 
    char_value["t"] = 1 
    char_value["c"] = 2 
    char_value["g"] = 3 
    i = 0 
    j = 0 
    base_four_string = "" 

    while(i < mer_size): 
     base_four_string += str(char_value[my_string[i]]) 
     i += 1 

    index = int(base_four_string, 4) 

    return index 

Эта функция обрабатывает все фрагменты ДНК и отображает счетчики с индексом в массиве размера 4 ** к

def get_mer_count(mer_size, file_fragments, slidingSize): 
    mer_counts = {} 
    for fragment in file_fragments: 
     j = 0 
     max_j = len(fragment) - mer_size 
     while(j < max_j): 
      mer_frag = fragment[j:j+mer_size] 
      mer_frag = mer_frag.lower() 
      if("n" not in mer_frag): 
       try: 
        mer_counts[mer_frag] += 1 
       except: 
        mer_counts[mer_frag] = 1 
      j += slidingSize 

    myNSV = [0] * (4**mer_size) 
    for mer in mer_counts.keys(): 
     mer_index = mer_index_finder(mer, mer_size) 
     # examples showing how to collapse, 
     # without shrinking the array 
     # rev_mer = make_complment_mer(mer) 
     # print rev_mer 
     # rev_index = mer_index_finder(rev_mer, mer_size) 
     # min_index = min(mer_index, rev_index) 
     # print mer_index,"\t",rev_index,"\t",min_index 
     # myNSV[min_index] += mer_counts[mer] 
     myNSV[mer_index] = mer_counts[mer] 

    return myNSV[:] 

Наконец эта функция будет принимать Меры и произвести обратный комплемент :

def make_complment_mer(mer_string): 
    nu_mer = "" 
    compliment_map = {"a" : "t", "c" : "g", "t" : "a", "g" : "c"} 
    for base in mer_string: 
     nu_mer += compliment_map[base] 
    nu_mer = nu_mer[::-1] 
    return nu_mer[:] 

Похоже, что должно быть очевидным способом, чтобы всегда знать, сколько ячеек массива должны иметь при сворачивании в прямом и обратном дополняет вместе, и есть экзамен в литературе, и некоторые пакеты, показывающие это, были выполнены; однако я не нахожу, где в исходном коде они могут генерировать эти вычисления.

Вторая часть этого вопроса: как вы узнаете, является ли mer прямым или обратным дополнением, не проверяя оба?

Пример:

(вперед)

AAGATCACGG

(дополнение)

TTCTAGTGCC

(обратный комплемент)

CCGTGATCTT

В моем приведенном выше коде я беру нижний из двух индексов, но кажется, что должен быть способ понять это, не найдя индекс для каждого мера дважды: один раз вперед и один раз в качестве обратного дополнения.

TL; DR Каков размер массива, если форматированные и обратные дополнения сопоставлены с одним и тем же индексом?

Edit: Чтобы определить размер массива, используя ответ, который я модифицировал get_mer_count(), чтобы включить в него следующие строки, чтобы создать размер индекса:

array_size = (4 ** mer_size)/2 
if mer_size % 2 == 0: 
    array_size += 2**(mer_size - 1) 

myNSV = [0] * array_size 

Answer 1

Для каждого k -mer, есть два возможности: либо он имеет ровно одно обратное дополнение, либо это его собственный обратный комплимент («палиндромный» мер). Поэтому, если есть p palindromic k -mers, тогда мы знаем, что размер массива должен быть p + (4**k - p)/2.

• Для k нечетных, нет палиндромных меров, так как средний нуклеотид не может быть своим собственным комплиментом. Поэтому массив должен иметь размер 4**k/2.

• Для k уже тогда k = 2*j для некоторых j. Мер является палиндромным тогда и только тогда, когда его первая половина - обратный комплимент второй половины. Есть 4**j возможных первых половинок, поэтому есть p = 4**j = 2**k palindromic k -mers. Таким образом, используя нашу формулу выше, массив должен иметь размер p + (4**k - p)/2 = 2**k + (4**k - 2**k)/2.