Алгоритм сжатия ДНК-последовательностей для хранения геномов на сервере с очень быстрой распаковкой
Встала необходимость держать ряд геномов на веб-сервере (различного рода утилиты для биоинформатиков, генетиков и биологов/медиков, в целом). А это, как оказывается, стало недешево :) Объёмы большие, доступ нужен файловый (так проще для скриптов), возникает явное желание как-то запаковать так, чтобы они занимали мало места, но при этом очень быстро и без большой нагрузки на CPU распаковывались. Стандартные варианты, издавна используемые в биоинформатике (потоковые на базе gzip), не очень хороши: жмут в 2.5 - 4 раза (в принципе, неплохо), но большие затраты CPU и памяти.
Геном изначально - это текстовый файл с малым алфавитом. Пока используем придуманную ещё в стародавние времена идею: поскольку алфавит состоит всего лишь из 5 букв определений для нуклеотидов (A, G, T, C, N), то используем полубайтовое кодирование. Конкретно:
А 0001
C 0010
G 0100
T 1000
N 1111
Преобразуем текст в обычный файл. Размер файла сокращается в два раза (чуть больше, чем в два за счёт удаления символа возврата каретки), операции сравнения можно проводить побитово, что увеличивает производительность скриптов.
Хочется большего. Навскидку, идеи:
- Четвертьбайтовое кодирование, но тогда пропадает N (это обозначение любого нуклеотида, т.е. A или G или C или T) и поиск подстрок становится дольше
- Замена последовательности одинаковых нуклеотидов, например, AAAAAAAAA, на что-то вида ЧИСЛО;СИМВОЛ. Минусы: разворачивать в память, дополнительные коды парсинга (хотя 0000 никто не запрещал)
Плюс, эпидемия COVID-19 (РНК-вирус, а не ДНК) однозначно добавила необходимость хранения РНК, т.е. добавляется пятый символ - U вместо Т (понятно, что можно просто помечать, что это РНК, а не ДНК, но тогда теряется универсальность).
Гиганты, вроде Qiagen, используют, судя по всему, полубайтовое кодирование. Но на то они и гиганты, могут себе позволить большие затраты.
Может, мы что-то не видим, лежащее на поверхности, и есть простой способ сэкономить пространство не в 2, а, скажем, в 3 раза?
В качестве примера - подборка вариаций вируса иммунодефицита человека (не содержит N), формат FASTA (смотреть как тестовый файл, 39 МБ) и 17-я хромосома человека (текст, фактически одна строка, 79 МБ). Доступны на Яндекс-диске: https://disk.yandex.ru/d/4E-_ed6QT6AVug
Ответы (1 шт):
Если считать что все пять букв встречаются равновероятно, то теоретический предел плотности упаковки равен log5256 (≅ 3.4454). Практические комбинации для блоков из различного количества байтов, упорядоченные от самой плотной упаковки к самой неплотной:
байты нуклеотиды плотность % от макс. плотности 7 24 3.4286 99.5112 5 17 3.4000 98.6819 8 27 3.3750 97.9563 6 20 3.3333 96.7470 3 10 3.3333 96.7470 4 13 3.2500 94.3283 2 6 3.0000 87.0723 1 3 3.0000 87.0723
Верхняя граница - восемь байт, потому что распаковка/запаковка требуют целой арифметики. Современный компьютер обрабатывает за раз 64 бита - восемь байт.
Последние две строки таблицы: в один байт можно запаковать три нуклеотида. Потери - 13% от теоретического максимума.
Первая строка: 24 нуклеотида в семи байтах с плотностью очень близкой к максимальной. Потери около половины процента.
Упаковка/распаковка не сложная и очень быстрая. Произвольный доступ тоже возможен. Даже поиск по упакованному файлу можно сделать, но это сложнее.
P.S. Для алфавита из шести символов теоретический максимум плотности упаковки 3.0948. Таблица очень простая - самая плотная упаковка - три символа на байт. Около трёх процентов потерь:
байты нуклеотиды плотность % от макс. плотности 8 24 3.0000 96.9361 7 21 3.0000 96.9361 6 18 3.0000 96.9361 5 15 3.0000 96.9361 4 12 3.0000 96.9361 3 9 3.0000 96.9361 2 6 3.0000 96.9361 1 3 3.0000 96.9361