我有一组fasta文件，里面有很多dna的序列片段。我试图计算每个文件中k-mers的总出现次数。计算k-mers的好处是可以创建一个4**k大小的数组，其中k是所用k-mer的大小。我正在处理的序列文件是来自新一代测序机的短读序列，因此假设所有的读取都来自5'->3'端，就无法完成。解决这一问题的最佳方法是将观察到的所有k-mer映射到一个单一的正、反补体序列索引。
所需映射：
当k=2时，数组的起始索引为0
字符串='aa'；
映射到索引->0
字符串='tt'；
映射到索引->0
字符串='在'；
映射到索引->1_
通过手工，我能够计算出所有mer的数组（前补和后补的折叠）的长度为10，并且特定的索引将表示以下mer：
AA、AT、AC、AG、TA、TC、TG、CC、CG、GC
我很难想出一个通用算法来知道更大尺寸k的可能mer的数量。在count数组中应该分配多少个单元？
在我现有的代码中，我提出了这三个函数来处理片段，生成反向补码，并将MER（或反向补码）映射到索引。
第一个函数将获取mer字符串并返回与4k大小数组中的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是正补语还是反补语而不检查两者？
例子：
（向前）
阿加塔卡格
（补充）
TTCTAGTGCC公司
（反补）
关贸总协定
在上面的代码中，我取了两个索引中较低的一个，但似乎应该有一种方法来解决这个问题，而不必为每个mer找到两次索引：一次向前，一次作为反向补码。
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
最佳答案
对于每个k-mer，有两种可能：要么只有一个反向补语，要么是它自己的反向补语（一个“回文”mer）因此，如果存在p回文k-mers，那么我们知道数组大小应该是p + (4**k - p)/2。
奇怪的是，没有回文mers，因为中间的核苷酸不能是它自己的补充所以数组的大小应该小于cc>。
对于k甚至对于4**k / 2对于某些k当且仅当mer的前半部分与后半部分相反时，mer才是回文。可能存在前半部分，因此存在回文因此，使用上面的公式，数组的大小应该小于cc>。
关于python - 在python中折叠DNA序列的正向和反向互补的算法？，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/40952719/