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

【问题标题】：Python: working with bits. Coding nucleotides with zeros and onesPython：使用位。用 0 和 1 编码核苷酸
【发布时间】：2015-03-08 07:43:07
【问题描述】：

我想在 Python 中使用位编码对核苷酸“A”、“G”、“C”和“T”进行编码。例如：

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

为了构建一个包含 k-mer 的巨大字典，例如：

dic = { 'ATGACTGACT':231, 'AAATGACGGAC':500 ... }

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

我不确定这是否可以做到，如果可以，我如何使用 Python 来做到这一点

提前致谢！

已编辑：对不起，我没有正确解释自己。

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

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

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

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

目前我正在使用 itertools 构建该字典。

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

我猜另一个问题是每个 k-mer 都需要转换为该位编码才能访问字典

【问题讨论】：

不清楚：您的 dic 示例中仍然有完整和位版本的 k-mer，那么内存减少在哪里？另一个重要的问题 - 之后你打算用这些位做什么？部分搜索？使用指标？
00001100 会编码什么？亚特？助教？ AATA？
对于可变长度字符串，您可以使用 bitstring 之类的模块（特别是用作 dict 键的 ConstBitStream）。我不知道这实际上会如何影响性能/内存使用。顺便问一下，你能成为prematurely optimizing吗？我不是说你是；问题是，持续 4 倍的内存改进是否值得额外的复杂性？现在是否应该进行更改？
为了放大@Lack，除非你已经证明存在内存不足，否则不要针对预期的问题进行优化。此外，您应该探索Biopython，而不是重新发明轮子。使用 Python 的部分理由是大量的第三方库可能已经满足您的需求。
我不认为我过早优化。实际上，我的程序已经在运行，并且我使用了 Biopython。问题是处理人类基因组需要 7 个小时，现在我正在尝试另一种方法。

标签： python bit bioinformatics

【解决方案1】：

您可以将 kmers 转换为二进制文件，但正如 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))

典型输出：

{'AATGC': (54, 5), 'ATGC': (54, 4), 'ATGACTGACT': (215883, 10)}
(54, 5) AATGC
(54, 4) ATGC
(215883, 10) ATGACTGACT

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

【讨论】：

【解决方案2】：

这完全符合您的要求

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 在他们的 cmets 中提出的反对意见确实很重要，我认为您应该认真重新考虑您的方法。

【讨论】：

【解决方案3】：

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

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

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

此代码为输入中的第一个字符提供最低位位置，因此"AT" 变为0b0100，或4，而"TA" 变为0b0001 或1。如果您希望编码将第一个字母视为最重要的字母，请在生成器表达式中使用 enumerate(reversed(kmer)) 而不是 enumerate(kmer)。

正如其他人评论的那样，这些整数仅在给定长度k 内是唯一的。如果不同长度的字符串仅在尾随As 的数量上不同（例如"ATG"、"ATGA"、"ATGAA"、"ATGAAA" 等，则将给出相同的整数作为编码36)。

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

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

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

一种选择是使用dict 类的某些方法，以避免将结果集中的所有值初始化为0。如果您将增量代码更改为 d[key] = d.get(key, 0) + 1，则无论 key 是否已在字典中，它都将起作用。

另一种选择是使用collections.Counter 而不是常规的dict。 Counter 类专为计算输入序列中的项目实例而设计，这似乎正是您正在做的事情。它认为它尚未看到的任何键的值都为0。

【讨论】：