Эффективное преобразование последовательностей ACGT в k-мерные частотные векторы

В этой статье я представлю вам быстрый метод преобразования последовательностей ДНК в векторы. Мы рассмотрим метод и код, так как этот шаг не должен замедлять выполнение очень важного последующего анализа.

Последовательности ДНК представляют собой длинные цепочки алфавита ACGT. Поскольку корпус алфавита ACGT увеличивается экспоненциально как k⁴ с количеством букв k, обычно последовательности ДНК представлены в виде векторов маленьких k (3, 4, 5, 7, 9). Векторизация также является важным шагом в передаче последовательностей ДНК переменной длины в модели машинного обучения. Это связано с тем, что модели машинного обучения для обработки полагаются на векторы фиксированной длины. Это стало очень популярным с увеличением доступности чтения Nanopore и PacBio. От эффективности методики будет зависеть время, необходимое для векторизации. В разумном наборе данных мы ожидаем увидеть не менее 1 миллиона чтений. Это может еще больше увеличиться с помощью секвенирования метагеномных образцов на большой глубине, где роль таких векторов играет важную роль в таких сценариях, как биннинг.

Размер вектора

Существует k⁴ возможных различных k-мер, однако на самом деле они ложатся почти наполовину, если мы объединяем обратные дополнения. Однако всегда разумно вычислять их на лету, поскольку нечетные k-мерки и k-mers четного размера может вести себя немного иначе. Идея состоит в том, чтобы подсчитать k-мер и записать их в структуру данных hash-map. Таким образом, мы можем записывать количество k-mer за одно чтение в хэш-карту без изменения размера хэш-карты (Изменение размеров карты C ++ будет дорогостоящей задачей, поскольку она использует сбалансированное дерево. Использование неупорядоченной карты может сэкономить проблемы, но порядок размещения может привести к изменению порядка, которого нам следует избегать). Однако мы можем легко вычислить набор k-мер следующим образом. Здесь мы используем оптимизированный вариант дополнения почтений.

Буквы кодирования и обратные дополнения

A = 1000001
C = 1000011
G = 1000111
T = 1010100

У нас есть 4 буквы, которые нужно закодировать в двоичном формате для повышения эффективности. Для этого нам понадобятся 2 бита. Мы видим, что биты, выделенные жирным шрифтом, напоминают целые числа 0, 1, 3 и 2 соответственно. Поэтому мы используем операцию X ›› 1 & 3 для преобразования букв в целые числа. Обратите внимание, что мы выполняем сдвиг на 1 бит и выполняем побитовую операцию И, чтобы сохранить только два бита. Мы используем следующий код, чтобы получить обратное дополнение.

Я не буду объяснять код построчно для простоты. Однако в основном происходит то, что мы меняем местами биты блоками по 2 бита и вычисляем дополнения с помощью битовых операций. Для сокращения кода используются шестнадцатеричные значения. Этот подход может обрабатывать только k-мер до k = 15, что намного больше, чем обычно используется для векторизации.

Создание векторного заполнителя

Как обсуждалось ранее, нам нужна хеш-карта всех соответствующих k-мер. Это значительно ускоряет процесс, оставляя нам операции O (1). Давайте посмотрим, как будет выглядеть код.

Обратите внимание, что мы выбираем только минимальное значение k-mer из всех возможных k-мер.

Вычисление векторов

Теперь, когда мы рассмотрели, как можно вычислить все возможные k-меры, обратные дополнения и кодирование, давайте посмотрим, как мы можем получить векторы. Во-первых, давайте посмотрим на код, который преобразует одну строку в вектор. Пояснение будет следовать за кодом.

Параметры функции

  1. Ссылка на строку
  2. Размер k
  3. Нормализовать или нет (при необходимости, для будущих вычислений, подобных TF-IDF)
  4. Набор всех k-мер для быстрого копирования хеш-карты

Реализация функции

Сначала мы копируем набор k-мер в локальную переменную и инициализируем хеш-карту. Затем мы инициализируем вектор для хранения наших окончательных значений вектора, которые равны размеру хеш-карты. В первом цикле for мы повторяем строку один раз и выполняем битовую операцию для записи k-mer. Обратите внимание, что мы пропускаем буквы, не относящиеся к ACGT. Мы сохраняем переменную длины, чтобы пропускать такие окна, и учитываем значение val, когда записанная длина равна ожидаемому k. Когда размер достигнут, мы обновляем хеш-карту и продолжаем то же самое.

Во втором цикле мы нормализуем (если normalize = true) или просто записываем векторы в статистика. Теперь у нас есть вектор для поставляемой строки.

Обратите внимание, что в этой функции не выполняется операция копирования строки. Все вычисления являются математическими / побитовыми, что значительно ускоряет работу. Он может обрабатывать миллионы операций чтения за минуты или секунды, если доступны средства многопроцессорной обработки. Это также может быть легко расширено, чтобы вписаться в область программирования GPU (позвольте мне пока оставить этот момент).

Полная реализация

Надеюсь, у вас было хорошее представление о том, как мы можем эффективно вычислять векторы с помощью C ++. Однако, чтобы упростить задачу, я приложил полный код, который вы можете использовать или изменить в соответствии с вашими требованиями (код поддерживает FASTQ, т. Е. Считывает номера строк 2, 6, 10, и т. Д.).

Компиляция

g++ vectorize.cpp -o vectorize -fopenmp --std=c++17

Бег

./vectorize reads.fastq K-SIZE THREADS BATCH_SIZE > vectors.txt

Полная программа

Я надеюсь, что этот код поможет вам улучшить вашу исследовательскую работу и эффективность экспериментов.

Спасибо за прочтение. Ваше здоровье! 😊