我们怎样才能有效地压缩DNA串

Question

DNA字符串可以是任何长度，包括5个字母(A、T、G、C、N)的任何组合
压缩包含5个字母(A、T、G、C、N)的字母表DNA串的有效方法是什么。我们是否可以使用较少的比特数有效地压缩和检索，而不是考虑每个字母表3个比特？有人能提出一个有效压缩和检索的伪代码吗？

Answer 1

如果您愿意(a)每个字符都有不同的比特大小，并且(b)您总是从开始读取，而不是从中间读取，则可以。然后，你可以有一个类似的代码：

• A-00
• T-01
• G-10
• C-110
• N-111

从左到右阅读只能用一种方式将比特流拆分为字符。你一次读取2个比特，如果它们是"11"，你需要再读取一个比特才能知道它是什么字符

这是基于霍夫曼编码算法

注：
我对DNA了解不多，但如果字符的概率不相等(意味着每个字符的概率为20%)。您应该将最短的代码分配给概率较高的代码。

Answer 2

您有5个唯一的值，因此您需要一个基于5的编码(例如a=0，T=1，g=2，C=3，N=4)。

在32位中，您可以拟合log₅(2³²)=13基-5值。

在64位中，您可以拟合log₅(2⁶⁴)=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;
}

Answer 3

有很多方法可以压缩，但主要问题是要压缩什么数据？1.来自测序机的原始未比对测序数据(fastq)2.对齐的数据(sam/bam/class)3.参考基因组

1. 您应该重新排列您的读数，将来自相近基因组位置的读数彼此靠近。例如，这将允许普通的gzip压缩效果提高3倍。有很多方法可以做到这一点。例如，您可以将fastq与bam对齐，然后导出回fastq。使用后缀树/数组查找类似的读取，这是大多数对齐器的工作方式(需要大量内存)。使用最小化器-非常快速、低内存的解决方案，但不适合有很多错误的长时间读取。良好的结果来自debruijn图构造，该构造用于此目的(也称为从头对齐)。像huffman/athmetic这样的统计编码将压缩到1/3(可以将huffman流传递给二进制算术编码器以获得另外20%的增益)
2. 这里的最佳结果来自基于引用的压缩——只需存储引用和对齐读取之间的差异
3. 在这里几乎无能为力。使用统计编码，每个核苷酸可以得到2-3位

Answer 4

坦率地说，我将从一些版本的Lempel-Ziv压缩(一类包括通用gzip压缩格式的压缩算法)开始。我注意到一些评论说，通用压缩算法在原始基因组数据上效果不佳，但其有效性取决于数据如何呈现给他们。

请注意，大多数通用压缩程序(如gzip)以每个字节为单位检查其输入。这意味着以3位/碱基对基因组数据进行"预压缩"会产生反作用；相反，在通过通用压缩器运行之前，应该将未压缩的基因组数据格式化为每基一个字节。Ascii"AGTCN"编码应该很好，只要你不通过包含空格、换行符或大写变体来添加噪音。

Lempel-Ziv压缩方法的工作原理是识别输入中的重复子串，然后参考前面的数据对其进行编码；我认为这类方法应该在适当呈现的基因组数据上做得相当好。一种更具基因组特异性的压缩方法可能会在这方面有所改进，但除非对基因组编码有一些我不知道的强大的非局部约束，否则我预计不会有重大改进。

Answer 5

我们可以使用Roee Gavirel的想法和以下内容的组合来获得更紧密的结果。由于Roee的想法仍然规定将五个字符中的两个映射到一个3位单词，因此序列中至少五个字符之一没有出现但3位单词之一出现的部分可以映射到2位单词，从而减少我们的最终结果。

切换映射的条件是，如果存在一个部分，其中五个字符中至少有一个没有出现，而其中一个3位单词只出现了一次，是我们的部分前缀长度(以位为单位)的两倍多。如果我们对可能的字符进行排序(例如，按字母顺序)，给定三个比特来指示特定的缺失字符(如果有多个缺失字符，我们按顺序选择第一个)或没有缺失字符，那么我们可以立即为其他四个字符分配一致的2比特映射。

我们的前缀有两个想法：

(1)

• 3位：缺失的字符(如果没有，我们对该部分使用Roee编码)；

• x位：表示区段长度的恒定位数。对于最大长度为65000的截面，我们可以指定x=16。

为了证明前缀的使用是合理的，我们需要一个部分，其中五个字符中的一个没有出现，而三位单词中的一位出现了39次或更多。

(2)

• 3位：缺失的字符(如果没有，我们对该部分使用Roee编码)；

• x位：前缀下一部分的位数-取决于最长部分的字符数。x=6意味着最大部分长度可以是2^(2^6)！不大可能发生的对于最大长度为65000的截面，我们可以指定x=4；

• 在前缀的前一部分中指示的比特数，指示当前段长度。

在上面的例子中，我们的前缀长度可能在11到23位之间变化，这意味着为了证明它的使用，我们需要一个部分，其中五个字符中的一个没有出现，而三位单词中的一位出现了23到47次或更多。