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

Question

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

所需的映射：

与k = 2＆数组的起始索引为0？

string ='aa'; 映射到索引-> 0？

string ='tt'; 映射到索引-> 0？

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

通过手工，我能够弄清楚所有具有正向和反向互补折叠的聚体的数组的长度为10，并且特定索引将表示以下聚体：AA，AT，AC，AG，TA，TC，TG，CC，CG，GC

我在考虑通用算法时遇到麻烦，无法知道较大k的可能合并器的数量。在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

（反补）

贸易总协定

在上面的代码中，我采用了两个索引中的较低者，但似乎应该有一种方法可以解决这一问题，而不必为每个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聚体，那么我们知道数组大小应该为p + (4**k - p)/2。

• 对于k奇，有没有回文聚体，由于中间核苷酸不能将自己的赞美。因此，数组应具有size 4**k / 2。

• 对于k即使在当时k = 2*j一些j。当且仅当上半部是其后半部的反向补语时，mer才是回文。有4**j可能的第一半，所以有p = 4**j = 2**k回文k聚体。因此，使用上面的公式，数组应具有size p + (4**k - p)/2 = 2**k + (4**k - 2**k)/2。