在 python 中折叠 DNA 序列的正向和反向补码的算法？

Question

我有一组包含许多 DNA 序列片段的 fasta 文件。我正在尝试计算可以在每个文件中找到的 k-mers 的总出现次数。计算 k-mers 的好处是可以创建一个大小为 4**k 的数组，其中 k 是正在使用的 k-mers 的大小。我正在处理的序列文件是来自新一代测序仪的短读取序列，因此无法假设读取全部来自 5' -> 3' 端。解决该问题的最佳方法是将观察到的所有 k-mers 映射到正向和反向补体序列的单个索引。

所需映射：

k = 2 & 数组的起始索引为 0

string = ‘aa’; 映射到索引 -> 0

string = ‘tt’; 映射到索引 -> 0

字符串=‘在’；映射到索引 -> 1

通过手算，我能够计算出所有正向和反向补码崩溃的 mer 的数组长度为 10，具体索引将代表以下 mer： AA、AT、AC、AG、TA、TC、TG、CC、CG、GC

我很难想出一个通用算法来了解较大 k 的可能 mer 数。应在计数数组中分配多少个单元格？

在我现有的代码中，我提出了这三个函数来处理片段、生成反向补码并将 mer（或反向补码）映射到索引。

第一个函数将获取 mer 字符串和 return 与 4**k 大小数组中的 mer 相关的索引。

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

此函数处理所有 DNA 片段并将计数映射到大小为 4**k 的数组中的索引

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[:]

最后这个函数会取一个mer并生成反向补码：

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

在我上面的代码中，我取了两个索引中较低的一个，但似乎应该有一种方法可以解决这个问题，而不必为每个 mer 找到索引两次：一次正向和一次反向补充.

TL;DR 如果正向和反向补码映射到同一个索引，数组的大小是多少？

编辑：为了使用答案找到数组大小，我修改了 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，有两种可能性：要么恰好有一个反向互补，要么就是它自己的反向互补（一个"palindromic" mer）。所以如果有 p 回文 k-mers 那么我们知道数组大小应该是 p + (4**k - p)/2.

对于k奇数，没有回文聚体，因为中间的核苷酸不能是它自己的互补。所以数组的大小应该是 4**k / 2.
对于 k 甚至 k = 2*j 对于某些 j。一个 mer 是回文的当且仅当它的前半部分是其后半部分的反向恭维。有 4**j 个可能的前半部分，所以有 p = 4**j = 2**k 个回文 k-mers。因此使用我们上面的公式，数组的大小应该是 p + (4**k - p)/2 = 2**k + (4**k - 2**k)/2.

在 python 中折叠 DNA 序列的正向和反向补码的算法？

Algorithm to collapse forward and reverse complement of a DNA sequence in python?

python

algorithm

bioinformatics