在python中摺疊正向和反向互補的DNA序列的算法？

Question

我有一個包含許多DNA序列片段的fasta文件集合。我試圖計算每個文件中可以找到的k-mers的總髮生次數。計數k-mers的好處在於可以創建大小爲4 ** k的單個數組，其中k是所使用的k-mer的大小。我正在處理的序列文件是來自新一代測序機器的短讀序列，因此假定讀數全部來自5' - > 3'末端，這是無法完成的。解決這個問題的最好方法是將觀察到的所有k聚體映射到正向和反向互補序列的單個索引。

期望的映射：

，其中k = 2 &爲數組起始索引爲0

串=「AA」; maps to index - > 0

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

string ='at'; 映射到索引 - > 1

手動我能夠弄清楚，正向和反向互補摺疊的所有mers陣列的長度爲10，具體指數將代表以下mers： AA ，AT，AC，AG，TA，TC，TG，CC，CG，GC

我很難想出一個廣義算法來知道更大的k值可能的mers數。在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的陣列中的**畝

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如果正向和反向補碼映射到相同索引，數組的大小是多少？

編輯： 要使用我修改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聚體，有兩個可能性：或者它只有一個反向補充，或者它是它自己的反向補語（一種「迴文」）。所以如果有p迴文k -mer，那麼我們知道陣列大小應該是p + (4**k - p)/2。

• 對於k奇，有沒有迴文聚體，由於中間核苷酸不能將自己的讚美。所以陣列的尺寸應該是4**k/2。

• 對於k即使然後k = 2*j對於一些j。阿梅爾是迴文，當且僅當其前半部分是其下半部分的反面恭維。有4**j可能的前半部分，因此有p = 4**j = 2**k迴文k -mer。因此，使用我們上面的公式，陣列的尺寸應該是p + (4**k - p)/2 = 2**k + (4**k - 2**k)/2。