Как мы можем эффективно сжать нить ДНК

У вас есть 5 уникальных значений, поэтому вам нужна кодировка base-5 (например, A = 0, T = 1, G = 2, C = 3, N = 4).

В 32 бита вы можете поместить журнал₅(2³²знак равно 13 базовые 5 значений.

В 64 бита можно уместить журнал₅(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;
}

Существует множество способов сжатия, но главный вопрос в том, какие данные вы хотите сжать?
1. Необработанные невыровненные секвенированные данные с секвенатора (fastq)
2. Выровненные данные (sam / bam / cram)
3. Справочные геномы

1. Вы должны изменить порядок чтения, поместив чтения с близких позиций генома близко друг к другу. Например, это позволило бы обычному сжатию GZIP в 3 раза лучше. Есть много способов сделать это. Вы можете, например, выровнять fastq к bam и затем экспортировать обратно в fastq. Используйте Suffix Tree / Array, чтобы найти похожие чтения, как работает большинство выравнивателей (требует много памяти). Используйте минимизаторы — очень быстрое решение с низким объемом памяти, но не подходит для длительного чтения со многими ошибками. Хорошие результаты получены при построении графа debruijn, который используется для этой цели (так называемое выравнивание de-novo).
   Статистическое кодирование, такое как Хаффман / Арифметика, будет сжато до 1/3 (можно передать поток Хаффмана в двоичный арифметический кодер, чтобы получить еще 20%).
2. Наилучшие результаты здесь получаются из сравнений на основе ссылок — просто сохраняйте различия между ссылками и выровненным чтением.
3. Здесь мало что можно сделать. Используя статистическое кодирование, вы можете получить 2-3 бита на нуклеотид.

Честно говоря, я бы начал с некоторой версии сжатия Lempel-Ziv (класс алгоритмов сжатия, который включает в себя универсальные gzip формат сжатия). Я отмечаю, что в некоторых комментариях говорится, что алгоритмы сжатия общего назначения плохо работают с необработанными данными генома, но их эффективность зависит от того, как данные им предоставляются.

Обратите внимание, что большинство программ сжатия общего назначения (например, gzip) проверить их входные данные для каждого байта. Это означает, что «предварительное сжатие» данных генома со скоростью 3 бита / основание контрпродуктивно; Вместо этого вы должны отформатировать несжатые данные генома по одному байту на базу, прежде чем запускать их через универсальный компрессор. Кодирование Ascii «AGTCN» должно быть в порядке, если вы не добавляете шум, включая пробелы, переводы строки или изменения в заглавных буквах.

Методы сжатия Лемпеля-Зива работают, распознавая повторяющиеся подстроки в их входных данных, затем кодируя их со ссылкой на предыдущие данные; Я ожидаю, что этот класс методов должен делать достаточно хорошую работу с должным образом представленными данными генома. Этот метод сжатия, более специфичный для генома, может улучшить это, но если нет какого-то сильного нелокального ограничения на кодирование генома, о котором я не знаю, я бы не ожидал существенного улучшения.

Мы можем использовать комбинацию идеи Роуи Гавиреля и следующего для еще более трудного результата. Поскольку идея Роу все еще предусматривает, что два из наших пяти символов должны быть сопоставлены с 3-битным словом, участки последовательности, в которых хотя бы один из пяти символов не отображается, но одно из 3-битных слов может быть сопоставлено с 2 слова, уменьшающие наш конечный результат.

Условие переключения отображения — если существует раздел, в котором хотя бы один из пяти символов не появляется, а одно из 3-битных слов появляется только один раз более, чем в два раза больше длины префикса раздела в битах. Если мы упорядочим наши возможные символы (например, в алфавитном порядке), учитывая три бита, указывающие на определенный пропущенный символ (если пропущено более одного, мы выбираем первый по порядку) или ни одного пропущенного, мы можем немедленно назначить согласованное 2-битное отображение для остальных четырех персонажей.

Две идеи для наших префиксов:

(1)

• 3 бита: отсутствующий символ (если его нет, мы используем кодировку Роу для раздела);

• x битов: постоянное количество битов, представляющих длину секции. Для секций максимальной длины 65000 мы могли бы назначить x = 16.

Чтобы оправдать использование префикса, нам нужен раздел, в котором один из пяти символов не появляется, а одно из 3-битных слов появляется 39 или более раз.

(2)

• 3 бита: отсутствующий символ (если его нет, мы используем кодировку Роу для раздела);

• x bits: количество битов в следующем разделе префикса — зависит от того, сколько символов может быть в самом длинном разделе. x = 6 подразумевает, что максимальная длина секции может быть 2 ^ (2 ^ 6)! Навряд ли. Для секций максимальной длины 65000 мы могли бы назначить x = 4;

• количество битов указывается в предыдущей части префикса, указывающего текущую длину раздела.

В приведенном выше примере длина нашего префикса может варьироваться, скажем, от 11 до 23 битов, что означает, что для оправдания его использования нам понадобится раздел, в котором один из пяти символов не появляется, а одно из 3-битных слов появляется между От 23 до 47 раз или больше.