我们怎样才能有效地压缩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”，您需要再读一位以了解它是什么字符。

This is based on Huffman Coding Algorithm

注意： 我不太了解DNA，但是如果chars的概率不相等（每个意味着20％）。你应该将最短的代码分配给那些概率较高的代码。

Answer 2

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

在32位中，您可以使log5（232）= 13 base-5值。

在64位中，您可以使用log5（264）= 27个base-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 / cram）3。参考基因组

1. 你应该重新排序你的读数，从相近的基因组位置读取读数。例如，这将允许通常的gzip压缩3倍。有很多方法可以做到这一点。例如，您可以将fastq与bam对齐，然后再将其导出回fastq。使用后缀树/数组来查找类似的读取，大多数对齐方式的工作方式（需要大量内存）。使用最小化器 - 非常快速，低内存的解决方案，但对于具有许多错误的长读取不利。良好的结果来自debruijn图形构造，用于此目的（也称为de-novo alignment）。像霍夫曼/算术这样的统计编码将压缩到1/3（可以将霍夫曼流传递给二进制算术编码器以获得另外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。

为了证明前缀使用的合理性，我们需要一个部分，其中五个字符中的一个不出现，而一个3位字出现39次或更多。

(2)

• 3位：缺少的字符（如果没有，我们使用Roee的编码作为该部分）;
• x位：前缀下一部分的位数 - 取决于最长部分的字符数。 x = 6意味着最大截面长度可能是2 ^（2 ^ 6）！不太可能。对于65000的最大长度部分，我们可以指定x = 4;
• 前缀前一部分中指示的位数，表示当前节长度。

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