Python：使用位。用0和1编码核苷酸

Question

我想用Python中的位编码来编码核苷酸“A”、“G”、“C”和“T”。例如：

'A' = 00
'G' = 01
'C' = 10
'T' = 11

为了建立一个包含k-mers的巨大dict，比如：

^{pr2}$

我认为这可以减少dict所需的内存量，因为“ATGC”需要4个字节，但同一个单词需要8个位编码。在

我不确定这是否可以实现，如果可以，我如何使用Python实现呢

提前谢谢！在

编辑：很抱歉我没有正确地解释自己。在

我想要的是遍历一个由'ATGC'和一个大小为k的滑动窗口组成的序列，并计算每个k-mer在该序列中出现的次数。例如：

'ATGAATGAA' # with a sliding window of 5 would be
 dic = { 'ATGAA':2, 'TGAAT':1, 'GAATG':1, 'AATGA':1, }

因为在开始读取序列之前，我想用大小为k的“AGTC”的所有可能的组合来构造dict，以便以每个k-mer作为键访问dict，并将其值求和1，所以我想知道dict上的k-mer是否可以用位编码存储。或多或少：

dic = {1011001010: 3, 0000110011: 666, ...  etc }

目前我正在用itertools构建dict。在

# k-mers of size 8
{''.join(x):0 for x in itertools.product('ATGC', repeat=8)}

我想另一个问题是，为了访问dict，每个k-mer都需要转换成该位编码

Answer 1

这完全符合你的要求

In [11]: d={'A':'00','G':'01','C':'10','T':'11'}

In [12]: int('0B'+''.join([d[c] for c in 'ATGACTGACT']),2)
Out[12]: 215883

In [13]: int('0B'+''.join([d[c] for c in 'ATGACTGACT'[::-1]]),2)
Out[13]: 925212

In [14]: 

但是pmod和Ignacio Vazquez Abrams在他们的评论中提出的反对意见是非常重要的，我认为你应该认真地重新考虑你的方法。在

Answer 2

正如@gbofi的回答所示，将k-mer转换为0和{}之间的整数非常简单。另一种方法是进行数学编码：

def kmer_to_int(kmer):
    return sum(4**i * "ATGC".index(x) for i, x in enumerate(kmer))

我没有测试过这是否比构建一个二进制字符串然后将其转换为int更快

此代码为输入中的第一个字符指定最低位位置，因此"AT"变成{}，或者{}和{}变成{}或{}。如果希望编码将第一个字母视为最重要的，请在生成器表达式中使用enumerate(reversed(kmer))而不是{}。在

正如其他人评论的那样，这些整数只对给定长度k唯一。如果长度不同的字符串只在尾随的A个数上有所不同（例如"ATG"，"ATGA"，"ATGAA"，"ATGAAA"，等等，所有这些都编码为36）。在

至于你更广泛的目标是在一个更大的序列中计算特定k-mers的出现，我不确定你是否会看到用这种方式编码k-mers的优势。这些好处可能取决于数据集的详细信息。在

整数键的一个优点是它们允许您使用列表而不是字典来保存计数。您可以用lst = [0] * 4**k构建一个适当的列表，然后用lst[kmer_to_int(kmer)] += 1增加一个值。在词条数量相同的情况下，列表的开销确实比字典低，但我不确定差异是否会大到足以提供帮助。在

如果数据是稀疏分布的（也就是说，许多4**k个可能的k-mer序列从未出现在输入中），使用列表可能仍然会浪费大量内存，因为列表总是4**k个元素。更好的方法可能是使用其他一些方法来简化稀疏数据的dict代码。在

一种选择是对dict类的某些方法进行修改，以避免将结果集中的所有值初始化为0。如果您将增量代码改为执行d[key] = d.get(key, 0) + 1，那么无论{}是否已经在字典中，它都可以工作。在

另一个选择是使用collections.Counter，而不是常规的dict。Counter类是专门为计算输入序列中的项的实例而设计的，这似乎正是您所做的。它认为它还没有看到的任何键的值是0。在

Answer 3

您可以将kmer转换为二进制，但是正如Ignacio指出的那样，您仍然需要知道它们的长度，所以您可能还需要存储这些数据。所以，对于很长的序列，这仍然可以节省内存空间。在

下面是一些示例代码，它获取序列，对它们进行编码并再次解码：

encoding_map = {'A': 0, 'G': 1, 'C': 2, 'T': 3}
decoding_lst = ['A', 'G', 'C', 'T']


def encode(k):
    code = 0
    for ch in k:
        code *= 4
        code += encoding_map[ch]
    return code, len(k)


def decode(enc):
    code, length = enc
    ret = ''
    for _ in range(length):
        index = code & 3
        code >>= 2
        ret = decoding_lst[index] + ret
    return ret


kmers = ['ATGACTGACT', 'ATGC', 'AATGC']

kmerdict = {k: encode(k) for k in kmers}

print(kmerdict)

for key, enc in kmerdict.items():
    print(enc, decode(enc))

典型输出：

^{pr2}$

顺便说一句，不管序列有多长，Python应该能够处理编码和解码，因为整数可以扩展到足以容纳数字的位。在