Наука. JavaScript. Скорость. Можем ли мы выбрать только один?

К сожалению, идея о том, что эти три концепции являются взаимоисключающими, широко распространена среди многих разработчиков (и ученых). Пакеты биоинформатики написаны биологами, а затем разработчиками, часто на Python и Perl, и часто не оптимизированы для скорости (или даже удобства использования). У JavaScript непростая история, изобилующая жалобами на детали реализации (Брендан Эйх написал всего за десять дней, пошли!). Принято считать, что если вам нужна скорость, вам нужно работать с оборудованием. Как с помощью JavaScript достичь скорости, близкой к металлической?

V8 JavaScript Engine от Google прошел долгий путь. В этой статье я собираюсь объяснить, как я написал и реализовал алгоритм выравнивания вырожденной последовательности ДНК без пропусков в JavaScript (браузер + узел) с использованием ArrayBuffers и побитовых операторов.

Полная библиотека JavaScript, содержащая алгоритм и другие полезные функции, доступна в репозитории github, keitwhor / NtSeq. Прежде чем я углублюсь в мельчайшие детали алгоритма NBEAM (Nucleotide Bitwise Exhaustive Alignment Mapping), я хотел бы прояснить некоторые определения как для биологов, так и для разработчиков.

Во-первых, для биологов: приведенный здесь алгоритм предназначен для сопоставления последовательностей без пропусков и не реализует матрицу замещения. Он оценивает только точное совпадение: 1 за совпадение (вырожденный нуклеотид или нет) и 0 за отсутствие совпадения. Я уверен, что ее можно расширить, чтобы использовать более сложную систему подсчета очков. Он лицензирован Массачусетским технологическим институтом, так что сходите с ума! Просто отдайте должное там, где полагается.

Для разработчиков: немного жаргона!

Нуклеотид:
«буква» ДНК. Обычно A, T, G или C. В РНК U используется вместо T.

Вырожденный нуклеотид:
буква, обозначающая один или несколько возможных нуклеотидов. Есть 15 возможных символов, 16, если вы включаете нет совпадений (представлены как -). Для получения дополнительной информации вы можете ознакомиться с нотацией нуклеиновых кислот в Википедии.

Без пробелов:
Последовательность ДНК не содержит пробелов или нулевых нуклеотидов. Хотя сам алгоритм поддерживает пробелы, он не предпринимает никаких усилий для проверки вариантов последовательностей с пробелами. Он всегда будет считать пропуск в предоставленной последовательности как «несоответствие».

Матрица замещения:
Способ оценки нуклеотидных сравнений на основе вероятности случайной мутации одного нуклеотида в другой тип нуклеотида в течение определенного периода времени. Этот тип оценки здесь не используется, но его нужно уточнить.

Геном:
Вы, наверное, слышали об этом. Строка нуклеотидов, содержащая всю генетическую информацию о конкретном организме.

Эта проблема

Во-первых, давайте посмотрим, что мы хотим сделать. У нас есть две последовательности, назовем их s eqSearch и seqGenome. Проблема, которую мы пытаемся решить, заключается в том, что мы хотим найти A) лучшее соответствие seqSearch в seqGenome и B) все последующие «близкие» совпадения seqSearch в seqGenome, ранжированный по сходству, а также по положению, в котором они совпадают.

Итак, если seqSearch - «ATGC», а seqGenome - «ATGGCATGC», мы ожидаем, что список совпадений будет выглядеть следующим образом (с индексированной позицией 0):

Ранг: 1, Позиция: 5, Матчи: 4/4, Результат: ATGC
Ранг: 2, Позиция: 0, Матчи: 3/4, Результат: ATGG
Ранг: 3, Позиция: 1, Матчи: 2/4, Результат: ТГГК
… и т. Д.

Алгоритм для реализации этого относительно прост:

function search(seqSearch, seqGenome) {
  var sLen = seqSearch.length;
  var gLen = seqGenome.length;
  // Create a Uint32Array to initialize all values to 0
  var map = new Uint32Array(gLen + sLen);
  var curChar;
  var offset;
  for (var j = 0; j < sLen; j++) {
    curChar = seqSearch[j];
    /* As we progress through our seqSearch,
       our offset in our map (the position
       we're matching at) gets changed.
       This is because we're trying to compare
       and aggregate match count based on the
       alignment of seqSearch as compared to
       its first (0-index) position. */
    offset = sLen — j;
    for (var i = 0; i < gLen; i++) {
      if (curChar === seqGenome[i]) {
        ++map[offset + i];
      }
    }
  
  }
  /* Convert map back into regular array,
     so we can use array methods */
  return [].slice.call(map);
};
// Call our search
var alignmentMap = search(seqSearch, seqGenome);
// Map position values
alignmentMap = alignmentMap.map(function(value, index) {
  return {position: index - seqSearch.length, matches: value};
});
// Sort position values
alignmentMap.sort(function(a, b) {
   return a.matches - b.matches;
});

Это очень просто. При рассмотрении агрегирования и сортировки результатов наш вызов map будет выполняться за O (n), а наш вызов sort (с реализацией быстрой сортировки) должен выполнить со средней временной сложностью O (n log n). Однако, поскольку у нас есть вложенный цикл, сравнивающий каждый нуклеотид seqSearch с seqGenome, мы столкнемся с временной сложностью O (n²) для нашего поискового вызова - ура!

Давайте сделаем перерыв прямо здесь - следующий алгоритм не пытается добиться лучшей производительности с точки зрения временной сложности. Он по-прежнему выполняется за O (n²). Однако, зная, что это узкое место таких исчерпывающих сравнений последовательностей, мы можем попытаться оптимизировать этот простой алгоритм поиска, чтобы он был как можно более производительным.

Проблемы первоначальной реализации

Первая проблема с описанной выше реализацией заключается в том, что она не соответствует вырожденным нуклеотидам. Например, если seqSearch содержит «W» (который должен соответствовать «A» или «T»), он точно увеличит счетчик совпадений только в том случае, если seqGenome содержит «W» " также. Не хорошо!

Вторая проблема (и она относится к первой) - это наша условная формулировка. Обычно предсказание ветвления (чтобы узнать больше, прочитайте мой любимый ответ на Stack Overflow за все время) будет пытаться ускорить это для нас:

if (curChar === seqGenome[i]) {
  ++map[offset + i];
}

Если во всей последовательности нет совпадений (скажем, все A против всех T), у вас есть возможность для предсказателя ветвления предположить, что он может пропустить команда, содержащаяся внутри каждый раз, она всегда будет правильной, и вуаля - вы будете работать очень быстро! (Тесты для NtSeq на самом деле показывают этот результат для случаев 0% идентичности между последовательностями. Отлично!)

Однако в случае геномных данных совпадения (более или менее) распределены случайным образом, что означает, что предсказатель ветвления будет работать очень плохо (и всегда проверять условный оператор).

Эти условные проверки требуют времени. Вы выполняете доступ к строке / массиву и сравнение символов для каждого отдельного нуклеотида как в seqSearch, так и в seqGenome.

Что, если бы существовал способ A) адекватно сопоставить вырожденные нуклеотиды и B) сопоставить и подсчитать несколько нуклеотидов в последовательности сразу, а не по одному?

Оказывается, есть способ сделать и то, и другое, и JavaScript отлично с этим справляется!

Сохранение информации о последовательности с использованием двоичных данных

Что, если вместо использования строк для сравнения данных последовательности мы будем использовать целые числа?

Прежде всего, зачем нам это делать?

Начнем с краткого изложения целых чисел. Все целые числа в JavaScript - это 32-битные целые числа со знаком. Это означает, что у них есть 1 знаковый бит и 31 бит, содержащий информацию о числе.

1, представленная в двоичном виде как 32-битное целое число:

00000000 00000000 00000000 00000001

Здесь я разделил байты 32-битного целого числа для удобства чтения. Напоминаем, что каждый раз, когда вы устанавливаете целочисленное значение в JavaScript (например, var a = 1;), вы фактически выделяете 4 байта (32 бита) памяти. .

Ну аккуратно. Оказывается, мы можем использовать это свойство целых чисел для хранения восьми значений вырожденных нуклеотидов в одном целом числе.

Как спросите вы? Помните, что для вырожденных нуклеотидов существует 16 возможных символов (включая «нет совпадений»). Это, как правило, 2⁴ различных возможностей и, таким образом, может быть представлено с использованием только четырех битов данных. Теперь мы можем создать некоторую семантику того, как читать эти четыре бита данных и интерпретировать их как нуклеотиды (или их вырожденные версии).

1000 is A (8)
0100 is T (4)
0010 is G (2)
0001 is C (1)

Это позволяет нам творить чудеса!

Чтобы получить значения вырожденных нуклеотидов:
W = (A | T) = (1000 | 0100) = 1100
N = (A | T | G | C) = (1000 | 0100 | 0010 | 0001) = 1111

Чтобы сопоставить нуклеотиды:
A & A = (1000 & 1000) = 1000
A & W = (1000 & 1100) = 1000
A & T = (1000 & 0100) = 0000

И так далее! Если вы запутались, проверьте Побитовые операторы, прежде чем читать дальше. Вы должны хотя бы понимать ИЛИ (|), И (&), НЕПОДПИСАННЫЙ БИТСДВИГ ВПРАВО (›››) и БИТСДВИГ ВЛЕВО (‹----------------).

Теперь, когда у нас есть все наши потенциальные символы нуклеотидов, хранящиеся в четырех битах данных, мы можем объединить строки из восьми нуклеотидов в одно 32-битное целое число (32/4 = 8). Например, «ATGCATGC» становится 1000 0100 0010 0001 1000 0100 0010 0001.

… Представляем ArrayBuffers!

К счастью, в JavaScript есть действительно простой способ работы с двоичными данными. Это с помощью ArrayBuffer. ArrayBuffer - это просто массив байтов, который может быть сопряжен с 8-битными, 16-битными или 32-битными целыми числами (со знаком или без знака) с использованием Uint8Array, Uint16Array и Uint32Array . Для простоты в этих примерах я не использую целые числа со знаком - буква U перед именами интерфейсов означает беззнаковые.

Помните, что вы можете изменять ArrayBuffers только с помощью типов данных views или Uint8Array (и т. Д.). Создание нового представления данных в ArrayBuffer не изменит базовые данные, а только то, как они представлены. Это означает, что мы можем установить значения ArrayBuffer с помощью Uint8Array и преобразовать его в Uint32Array, просто создав экземпляр нового представления в базовом буфере.

Чтение данных нашей последовательности в массив 32-битных целых чисел теперь можно упростить:

function convertSequenceToArrayOfIntegers(sequenceString) {
  // two nucleotides per byte (4 bits each)
  var sequenceBytes = Math.ceil(sequenceString.length / 2);
  var sequenceBuffer = new ArrayBuffer(sequenceBytes);
  var uint8view = new Uint8Array(sequenceBuffer);
  /*
     lookupTable should be an object mapping each of the 16
     nucleotide characters to their respective 4-bit integer
     value.
     The following would do this for string "AT"
     uint8view[i] starts as 0000 0000 (binary)
     uint8view[i] = (1000) << 4 ("A" in binary bit shift left 4)
     uint8view[i] is now 1000 0000 (binary)
     uint8view[i] |= (0100) ("T" in binary)
     uint8view[i] is 1000 0000 and gets OR (|) with 0100,
       resulting in 1000 0100
     
   */
  for (var i = 0, len = sequenceBuffer.length; i < len; i++) {
    uint8view[i] = lookupTable[sequenceString[i * 2]] << 4;
    uint8view[i] |= lookupTable[sequenceString[i * 2 + 1];
  }
  // sequenceBuffer was modified by uint8view
  // you can now return it as an array of 32-bit ints
  
  return new Uint32Array(sequenceBuffer);
}

Если мы запустим этот метод для обеих входных последовательностей, мы превратим каждую из них в массивы 32-битных целых чисел. Затем мы можем использовать эти массивы для очень быстрого сравнения нуклеотидов.

Сравнение последовательностей с использованием побитовых операторов

А теперь самое интересное! Ранее я упоминал, что мы можем сравнить два 4-битных нуклеотида на совпадение, используя оператор BIT AND (&).

Например, A & A = (1000 & 1000) = 1000.

В этом случае каждый раз, когда вы И два совпадающих нуклеотида, вы получаете последовательность битов, содержащую хотя бы одну 1. Если совпадений нет, последовательность из четырех битов всегда будет 0000.

В поразрядных операторах замечательно то, что (в идеале, работающие на уровне процессора) они занимают только один цикл процессора (ваш процессор, вероятно, работает более 2 миллиарда циклов в секунду) и воздействуют на весь 32-битный регистр сразу. Что это обозначает? Это означает, что мы можем сравнивать восемь вырожденных нуклеотидов одновременно, когда они хранятся как 32-битные целые числа, и это занимает всего долю наносекунды.

«ATGCATGW & ATATWWNN» станет:

1000 0100 0010 0001 1000 0100 0010 1100 &
1000 0100 1000 0100 1100 1100 1111 1111
=======================================
1000 0100 0000 0000 1000 0100 0010 1100

Потрясающие! Теперь, как нам эффективно подсчитать количество совпадений в результирующей битовой строке?

Во-первых, обратите внимание, что: Сравнение двух вырожденных нуклеотидов приведет к получению 4-битной последовательности, содержащей более одного бита. Это означает, что мы не можем просто подсчитать биты для совпадения. (W и W = (1100 и 1100) = 1100).

Однако несовпадающие нуклеотиды всегда приводят к 0000.

Итак, что мы можем сделать, так это собрать флаг соответствия для каждой 4-битной последовательности. Мы делаем это с помощью UNSIGNED RIGHT BIT SHIFT (›››) и BIT OR (|). Мы можем сгруппировать каждый набор из четырех битов в их крайнем правом месте, указывая совпадение или отсутствие совпадения (будь то 1 или 0).

// matches = 1000 0100 0000 0000 1000 0100 0010 1100 (binary)
matches |= matches >>> 1;
/*
  1000 0100 0000 0000 1000 0100 0010 1100 (matches) |
  0100 0010 0000 0000 0100 0010 0001 0110 (matches >>> 1)
  =======================================
  1100 0110 0000 0000 1100 0110 0011 1110
*/
matches |= matches >>> 2;
/*
  1100 0110 0000 0000 1100 0110 0011 1110 (matches) |
  0011 0001 1000 0000 0011 0001 1000 1111 (matches >>> 2)
  =======================================
  1111 0111 1000 0000 1111 0111 1011 1111
*/

Вы заметите, что значение совпадений, хотя и изменено (и не похоже на наше исходное целое число!), На самом деле имеет только восемь важных битов. Это крайние правые биты каждой группы из четырех. Везде, где было совпадение, есть 1, а везде, где совпадения не было, есть 0.

Мы можем удалить ненужные биты с помощью простого & с шестнадцатеричным флагом 0x11111111. (0001 0001 0001 0001 0001 0001 0001 0001).

matches &= 0x11111111;
/*
  1111 0111 1000 0000 1111 0111 1011 1111 (matches) &
  0001 0001 0001 0001 0001 0001 0001 0001
  =======================================
  0001 0001 0000 0000 0001 0001 0001 0001
*/

Красивый! Теперь нужно подсчитать количество 1 -битов в цепочке битов.

Есть несколько способов сделать это. Один - это алгоритм, который выполняется за O (n):

function countBits(n) {
  var count = 0;
  	while (n !== 0) {
  	  count++;
  	  n &= (n - 1);
  	}
  	return count;
}

Также существует алгоритм веса Хэмминга, упомянутый в другом красивом ответе StackOverflow, который выполняется в O (1):

function countBits(n) {
  n = n - ((n >> 1) & 0x55555555);
  n = (n & 0x33333333) + ((n >> 2) & 0x33333333);
  return (((n + (n >> 4)) & 0x0F0F0F0F) * 0x01010101) >> 24;
}

Ага! О (1)! Должно быть лучше! Однако большая проблема, с которой мы сталкиваемся, заключается в том, что наш счетчик бит в нашем фактическом целом числе относительно разрежен.

У нас есть не более восьми бит, которые можно подсчитать. Это означает, что для нашего решения O (n) потребуется максимум около 32 процессорных циклов (условная проверка, приращение, бит И, вычитание - все, умноженное на 8, будет равно 32) и минимум из 0 циклов , тогда как наше решение O (1) Веса Хэмминга всегда занимает около 14 циклов процессора в зависимости от количества операций (приблизительный подсчет , поправьте меня если я ошибаюсь).

Теперь о наших фактических данных - в двух случайно распределенных последовательностях ДНК вы, скорее всего, встретите совпадение только в 25% случаев, а это означает, что наше решение O (n) будет в среднем 8 циклов процессора при столкновении с реальными данными. Это будет лучше, чем вес Хэмминга.

Круто, так что мы воспользуемся этим, не так ли?

Не так быстро! Оказывается, есть лучшее решение O (1), чем Вес Хэмминга. Вы можете использовать таблицу поиска (как массив) для подсчета битов в битовой строке. LookupBits [0] равняется 0, LookupBits [1] - 1, LookupBits [2] - 1, LookupBits [3] - 2… и т. Д. Проблема с этим решением в том, что для этого требуется тонна памяти. Существует 4294967295 различных 32-битных целых чисел без знака, каждое из которых содержит значение в таблице поиска от 0 до 32 (до 3 байтов!), Что означает, что размер таблицы поиска будет 12 ГБ. Фу.

  1. Вероятно, у вас нет 12 ГБ памяти для выделения
  2. Даже если бы вы это сделали, справочная таблица такого размера неуправляема. Промахи кеша замедлят ваш процесс до полной остановки.

Итак, как это можно использовать? Что ж, мы знаем, что на самом деле у нас всего восемь важных бит. Давайте еще раз посмотрим на нашу результирующую битовую строку:

0001 0001 0000 0000 0001 0001 0001 0001

Каждый важный бит содержится в самой правой группе каждых четырех битов. Их тоже всего восемь. Это означает, что на самом деле существует только один байт релевантных данных и 256 возможных комбинаций результатов. Оказывается, мы можем снова использовать битовые сдвиги, чтобы агрегировать наши данные до самого правого байта!

// matches = 0001 0001 0000 0000 0001 0001 0001 0001
matches |= matches >>> 3;
/*
  0001 0001 0000 0000 0001 0001 0001 0001 (matches) |
  0000 0010 0010 0000 0000 0010 0010 0010 (matches >>> 3)
  =======================================
  0001 0011 0010 0000 0001 0011 0011 0011
*/
matches |= matches >>> 6;
/*
  0001 0011 0010 0000 0001 0011 0011 0011 (matches) |
  0000 0000 0100 1100 1000 0000 0100 1100 (matches >>> 6)
  =======================================
  0001 0011 0110 1100 1001 0011 0111 1111
*/

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

0001 0001 0000 0000 0001 0001 0001 0001

стал:

0001 0011 0110 1100 1001 0011 0111 1111

Теперь мы можем сдвинуть 1100 на 12 бит, очистить несущественные биты и выполнить БИТ ИЛИ (|) с 1111 справа:

matches = ((matches >>> 12) & 0xF0) | (matches & 0xF)
/*
  First:
  0000 0000 0000 0001 0011 0110 1100 1001 (matches >>> 12) &
  0000 0000 0000 0000 0000 0000 1111 0000 (0xF0)
  =======================================
  0000 0000 0000 0000 0000 0000 1100 0000
  
  Then:
  0001 0011 0110 1100 1001 0011 0111 1111 (matches) &
  0000 0000 0000 0000 0000 0000 0000 1111 (0xF)
  =======================================
  0000 0000 0000 0000 0000 0000 0000 1111
  Finally:
  0000 0000 0000 0000 0000 0000 1100 0000 |
  0000 0000 0000 0000 0000 0000 0000 1111
  =======================================
  0000 0000 0000 0000 0000 0000 1100 1111
*/

Потрясающие! Теперь у нас есть все данные о совпадениях в крайнем правом байте, и мы можем проверить их в поисковой таблице, содержащей всего 256 записей! (Уф! Вы можете предварительно сгенерировать эту таблицу поиска как Uint8Array, когда ваша программа загружается с использованием веса Хэмминга или нашего примера с разреженным битовым числом выше.)

return bitCountLookup[matches];

Отлично, теперь мы можем подсчитывать совпадения нуклеотидов по восемь сравнений за раз. Все, что осталось сделать, это переписать нашу первоначальную функцию search.

Окончательная реализация

Если вы до сих пор дошли до конца, теперь мы можем вернуться к нашей первой реализации сравнения строк, но заменить нашу условную проверку некоторыми битовыми операциями для запуска полного алгоритма. Я не буду вдаваться в подробное объяснение, но для каждого набора из двух целых чисел мы должны убедиться, что мы «обходим» целые числа друг над другом (смещение битов влево и вправо соответствующим образом), чтобы убедиться, что мы улавливаем потенциальные совпадения в каждой позиции. , и сохраните совпадения относительно смещения.

Следующий код представляет собой полный отрывок из NtSeq.MatchMap в моей Библиотеке NtSeq. Я прокомментировал части, которые могут сбивать с толку. Формулировка и имена переменных немного отличаются от приведенных выше.

MatchMap.prototype.__countMatches = function(int, bitCount) {
  int |= int >>> 1;
  int |= int >>> 2;
  int &= 0x11111111;
  int |= int >>> 3;
  int |= int >>> 6;
  return bitCount[((int >>> 12) & 0xF0) | (int & 0xF)];
};
MatchMap.prototype.__execute = function(queryBuffer, searchSpaceBuffer) {
  
  var queryInts, spaceInts, queryIntsLength, spaceIntsLength, arrLen, mapBuffer, mapArray, A, B, A1, A2, T, cur, pos, move, i, k, adjustNeg, adjustPos, fnCountMatches, bitCount;
  
  queryInts = new Uint32Array(queryBuffer, 4);
  spaceInts = new Uint32Array(searchSpaceBuffer, 4);
  fnCountMatches = this.__countMatches;
  bitCount = __bitCount;
  queryIntsLength = queryInts.length|0;
  spaceIntsLength = spaceInts.length|0;
  arrLen = (queryIntsLength + spaceIntsLength) << 3;
  mapBuffer = new ArrayBuffer(4 * arrLen);
  mapArray = new Uint32Array(mapBuffer);
  for (k = 0|0; k < queryIntsLength; k++) {
    A = queryInts[k];
    // set offset, x << 3 is shorthand for x * 8
    // because there are 8 nucleotides per int
    cur = (queryIntsLength — k) << 3;
    for (i = 0|0; i < spaceIntsLength; i++) {
      // if match without shifting is non-zero, count matches
      (T = A & spaceInts[i]) &&
        (mapArray[(i << 3) + cur] += fnCountMatches(T, bitCount));
    }
    // start walking "A" along searchSpace (seqGenome) by
    // bitshifting left and right, adjust offsets accordingly
    A1 = A >>> 4;
    A2 = A << 4;
    adjustNeg = cur — 1;
    adjustPos = cur + 1;
    // break loop if A1 and A2 have been shifted far enough
    // to zero them both out
    while(A1 || A2) {
      for (i = 0|0; i < spaceIntsLength; i++) {
        B = spaceInts[i];
        pos = (i << 3); // === i * 8
        // if the match result is non-zero, count matches.
        (T = A1 & B) &&
          (mapArray[pos + adjustNeg] += fnCountMatches(T, bitCount));
        (T = A2 & B) &&
          (mapArray[pos + adjustPos] += fnCountMatches(T, bitCount));
      }
      
      // keep "walking" / shifting current integer to each offset
      A1 >>>= 4;
      A2 <<= 4;
      --adjustNeg;
      ++adjustPos;
    }
  }
  
  // return our buffer, we can instantiate a Uint32Array
  // on it if we wish
  
  return mapBuffer;
};

Теперь мы кое-что достигли. ☺

Контрольные точки

Так как же это на самом деле работает? Получается неплохо. Несмотря на такую ​​же временную сложность, он работает в 6–7 раз быстрее, чем «простая» реализация. Вот некоторые данные, полученные 7 февраля 2015 года:

Benchmark         |  naive | search |   naiveScore |  searchScore
------------------------------------------------------------------
1,000,000, 0%     |    9ms |    3ms |    9.00ns/nt |    3.00ns/nt
10,000,000, 0%    |   63ms |    5ms |    6.30ns/nt |    0.50ns/nt
100,000,000, 0%   |  621ms |   60ms |    6.21ns/nt |    0.60ns/nt
1,000,000, 25%    |   15ms |    6ms |   15.00ns/nt |    6.00ns/nt
10,000,000, 25%   |  124ms |   17ms |   12.40ns/nt |    1.70ns/nt
100,000,000, 25%  | 1249ms |  233ms |   12.49ns/nt |    2.33ns/nt
1,000,000, 50%    |   15ms |    2ms |   15.00ns/nt |    2.00ns/nt
10,000,000, 50%   |  131ms |   20ms |   13.10ns/nt |    2.00ns/nt
100,000,000, 50%  | 1305ms |  234ms |   13.05ns/nt |    2.34ns/nt
1,000,000, 100%   |   14ms |    2ms |   14.00ns/nt |    2.00ns/nt
10,000,000, 100%  |  144ms |   18ms |   14.40ns/nt |    1.80ns/nt
100,000,000, 100% | 1471ms |  240ms |   14.71ns/nt |    2.40ns/nt

«Наивный» представляет первую версию алгоритма выравнивания, перечисленную здесь (базовое сравнение строк), а «поиск» представляет полную реализацию с использованием побитовых операторов. Оценка измеряется в наносекундах на нуклеотид, сравнение проводилось на процессоре 2,4 ГГц.

Заголовки испытаний представляют собой общее пространство поиска (длина seqSearch, умноженная на длину seqGenome) и среднюю идентичность последовательностей (совпадение) между двумя последовательностями.

Выводы

Что ж, JavaScript довольно быстро выполняет этот оптимизированный алгоритм битовых операций.

Как уже упоминалось, механизм V8 в node.js в среднем составляет около 5 процессорных циклов на сравнение нуклеотидов (включая доступ, хранение, агрегацию данных), что примерно «близко к металлу», насколько это возможно.

Лично я считаю, что пришло время более серьезно взглянуть на использование JavaScript как «тяжелого» языка. Команда V8 в Google проделала фантастическую работу по повышению производительности своего движка JavaScript. Что касается научных вычислений, потенциал есть. JavaScript чрезвычайно доступен (для работы требуется только браузер и текстовый редактор) и может быть легко представлен молодым ученым и биоинформатикам. Есть также аргументы в пользу масштабирования узловых процессов, но это можно отложить на другой день.

Наконец, описанный здесь алгоритм NBEAM, хотя и реализован на JavaScript, на самом деле не зависит от языка. Пожалуйста, не стесняйтесь изменять его в соответствии с вашими потребностями (мне было бы интересно увидеть реализации GPU!), Но мы всегда ценим кредит. Если кто-то может помочь в проверке новизны NBEAM и помочь опубликовать его, не стесняйтесь обращаться к нам!

Напоминаем, что библиотека NtSeq JavaScript и полная реализация алгоритма NBEAM доступны по адресу https://github.com/keithwhor/NtSeq.

Контакт

Вы можете посетить мой личный веб-сайт по адресу keithwhor.com, подписаться на меня в твиттере @keithwhor, посетить мой github или просто крикнуть привет, если вы случайно увидите меня гуляющим по улицам Сан-Франциско.

Спасибо за прочтение!