logo标签列表

如何高效压缩DNA字符串

c++algorithmcompressiondna-sequencelossless-compression
5

DNA字符串可以包含任意长度的5个字母（A、T、G、C、N）的任意组合。
有没有一种有效的方法来压缩由这5个字母组成的DNA字符串？除了考虑每个字母3个比特位外，我们能否使用更少的比特位进行高效的压缩和检索？有人能提供一个有效的压缩和检索的伪代码吗？

4
只是一个想法……你可以使用7位来编码由3个字母组成的序列，以便浪费更少的位数。 - 463035818_is_not_a_number
3
为什么不使用现成的压缩解决方案？ - Sorin
2
@RoeeGavirel 我刚刚尝试了一下，对于什么序列长度 log2( 5^n) 最接近下一个更大的整数，n=3 似乎是一个不错的选择。 - 463035818_is_not_a_number
1
这需要具体和高效吗？像7zip和gzip这样的工具可以将其压缩到难以改进的大小，但它们需要比专门的解决方案更多的时间。 - Daniel
1
以高效的方式检索 - “随机访问”或“从开头顺序访问”（请参见Roee Gavirel的答案中的(b)）？对于随机访问，5³略低于2⁷：使用7位编码每次三个碱基（N代表什么？）（请参见user463035818的评论）。 - greybeard
显示剩余11条评论
5个回答
7
如果你愿意 (a) 为每个字符拥有不同的位数，且 (b) 总是从开头而非中间读取，那么你可以使用如下代码：
  • A - 00
  • T - 01
  • G - 10
  • C - 110
  • N - 111
从左至右读取时，你只能以一种方式将一串比特拆分成字符。每次读取两个比特，如果它们是 "11"，则需要再读取一个比特才能确定该字符。 这基于哈夫曼编码算法 注意：
我对 DNA 不是很了解，但如果字符的概率不相等（即每个字符出现的概率不是 20%），则应将最短的代码分配给概率更高的字符。
我喜欢这个。理论上，每个符号需要大约log₂(5) ≈ 2.32位。这使用了平均每个符号12/5 = 2.4位。让我想起了UTF编码。 - wally
2
同样由于哈夫曼编码算法，您可以将最小序列（如11 10或01）设置为最常用的序列。因此，总位数需要更小。 - E141
@KonradRudolph - 实际上没有必要“尽可能平衡事物”。即使p(A)=p(C)，它们的编码长度也可以不同。将具有最低概率的符号N赋予最长的单词将创建最短的输出。 - Roee Gavirel
1
@RoeeGavirel 对，如果它们具有相同的频率。我的观点是，我们通常不知道这一点，但相反，我们知道p(C) + p(G) ≈ p(A) + p(T)。 - Konrad Rudolph
1
你也可以仅计算每个符号，然后在最终流中添加它们各自的编码前缀。该前缀只需9-10位，这将节省更多的空间。 - Cássio Renan
显示剩余3条评论
3

您有5个唯一的值，因此需要一个基于5的编码（例如：A=0，T=1，G=2，C=3，N=4）。

在32位中，您可以装下log5(232) = 13个基于5的值。

在64位中，您可以装下log5(264) = 27个基于5的值。

编码过程如下：

uint8_t *input = /* base-5 encoded DNA letters */;
uint64_t packed = 0;
for (int i = 0; i < 27; ++i) {
    packed = packed * 5 + *input++;
}

解码：

uint8_t *output = /* allocate buffer */;
uint64_t packed = /* next encoded chunk */;
for (int i = 0; i < 27; ++i) {
    *output++ = packed % 5;
    packed /= 5;
}
只要结果适合64位，您的代码就可以工作。这很少见;-) - Konrad Rudolph
1
@KonradRudolph只需将输入分成27个值的块并重复，每次迭代输出64位。 - rustyx
这会在边界处浪费比特。 - Konrad Rudolph
1
它每64位浪费不到1/2个比特。我们发现在许多情况下，它比哈夫曼编码更好，特别是如果输入字符的概率不平衡。使用一个好的大整数库（例如GMP），可以将块大小扩展到512位或更多，以进一步减少浪费。 - rustyx
1
我认为它浪费了一点，大约是64 - log₂(5) * 27 ≈每64位1.3个比特。或者你可以说它浪费了64 / log₂(5) - 27 ≈每64位0.56个符号。不过仍然相当不错。 - wally
1
有很多压缩方法，但主要问题是你想要压缩什么数据？ 1. 从测序机（fastq）获得的原始未对齐序列数据 2. 对齐数据（sam / bam / cram） 3. 参考基因组
你应该重新排序读取，将来自相近基因组位置的读取放在一起。例如，这将使通常的gzip压缩效果提高3倍。有很多方法可以做到这一点。例如，您可以将fastq对齐到bam，然后导出回fastq。使用后缀树/数组查找相似的reads，这是大多数比对器的工作方式（需要大量内存）。使用最小化器-非常快速，低内存解决方案，但不适合具有许多错误的长读取。良好的结果来自debruijn图构建，用于此目的（也称为de-novo比对）。 像霍夫曼/算术这样的统计编码会压缩到1/3（可以将霍夫曼流传递给二进制算术编码器以获得另外20％）。
这里最好的结果是基于参考的压缩-只需存储参考和对齐读取之间的差异。
这里几乎无法做任何事情。使用统计编码，您可以获得每个核苷酸2-3位。
0
我们可以结合Roee Gavirel的想法和以下方法来得到更紧凑的结果。由于Roee的想法仍然规定我们的五个字符中有两个要映射到一个3位字，所以序列的某些部分中至少有一个字符不出现但是一个3位字出现了，这些部分可以映射到2位字，从而减少我们的最终结果。
切换映射的条件是存在这样一个部分：至少有一个五个字符中的字符不出现，而一个3位字出现的次数比我们的节前缀长度多两倍的位数。如果我们按照可能的字符（例如按字母顺序）排序，给定三位表示特定缺失字符（如果有多个缺失，则选择第一个）或没有缺失，我们可以立即为其他四个字符分配一致的2位映射。
我们的前缀有两个想法：
（1）
- 3位：缺失的字符（如果没有，则使用Roee的编码进行节的编码）； - x位：表示节长度的常数位数。对于最大长度为65000的节，我们可以分配x = 16。
为了证明前缀的使用，我们需要一个部分，在其中五个字符中没有出现一个，并且一个3位字出现39次或更多。
（2）

  • 3位：缺失的字符（如果没有，则使用Roee的编码进行该部分）；

  • x位：前缀下一部分中的位数-取决于最长部分可能有多少个字符。 x = 6意味着最大部分长度可以是2^(2^6)! 不太可能。对于长度最大的部分为65000，我们可以分配x = 4;

  • 前缀上一部分指示的位数，表示当前部分的长度。

在上面的例子中，我们的前缀长度可以在11到23位之间变化，这意味着为了证明其使用，我们需要一个部分，在其中五个字符中没有出现一个，并且一个3位字出现在23到47次或更多。
0
坦率地说，我会从一些版本的Lempel-Ziv压缩开始（这是一类包括通用的gzip压缩格式的压缩算法）。请注意，有些评论说通用压缩算法在原始基因数据上不起作用，但它们的有效性取决于数据如何呈现给它们。
请注意，大多数通用压缩程序（如gzip）以每字节为基础检查其输入。这意味着在3比特/碱基处“预压缩”基因组数据是适得其反的；相反，在通过通用压缩器之前，应将未压缩的基因组数据格式化为每个碱基一个字节。Ascii“AGTCN”编码应该很好，只要您不通过包括空格、换行或大小写变化来添加噪音。
Lempel-Ziv压缩方法通过识别其输入中的重复子字符串，然后通过参考先前的数据对它们进行编码。我希望这类方法可以在适当呈现的基因组数据上表现出良好的效果。一种更具基因组特异性的压缩方法可能会有所改进，但除非我不知道基因组编码上存在某些强烈的非局部约束，否则我不会期望有显著的改进。
如果数据没有明显的重复模式，那么这种压缩方式是否也能减少每个字符的占用空间呢？ - גלעד ברקן
这可能取决于压缩方法的细节，但如果文件足够长，我预计它的表现可能与一阶霍夫曼编码器差不多。 - comingstorm
1
请注意，如果数据确实没有任何重复模式，那么根据定义，您只能压缩一阶统计数据（即单个符号频率）。在这种情况下，您无法比一阶哈夫曼编码器更好地完成任务... - comingstorm
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接