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

Question

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

所需的映射：

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

字符串 = 'aa'; 映射到索引 -> 0

字符串 = 'tt'; 映射到索引 -> 0

字符串 = 'at'; 映射到索引 -> 1

通过手工，我能够计算出具有正向和反向补码折叠的所有 mers 的数组的长度为 10，并且特定索引将代表以下 mers： AA、AT、AC、AG、TA、TC、TG、CC、CG、GC

我很难想出一个通用算法来了解更大尺寸 k 的可能 mer 数量。 count 数组中应该分配多少个单元格？

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

第一个函数将获取 mer 字符串并返回与 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，有两种可能性：要么它只有一个反向补码，要么是它自己的反向补码（“回文”mer）。所以如果有p回文k-mers 那么我们知道数组大小应该是p + (4**k - p)/2。

• 对于 k 奇数，没有回文 mers，因为中间核苷酸不能是它自己的补充。所以数组的大小应该是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。