Быстрый подсчет типов нуклеотидов в большом количестве последовательностей

В порядке важности:

1. Хороший код для этой задачи будет 100% I / O-связанный. Ваш процессор может считать символы намного быстрее, чем ваш диск может качать их в процессор. Таким образом, первый вопрос ко мне: какова пропускная способность вашего носителя? Какова ваша идеальная пропускная способность RAM и кеша? Это верхние пределы. Если вы нажали на них, нет смысла смотреть дальше на ваш код. Возможно, ваше решение для повышения уже есть.

2. std::map поиски относительно дороги. Да это так O(log(N)), но твой N=5 маленький и постоянный, так что это ничего вам не говорит. Для 5 значений карта должна будет использовать около трех указателей для каждого поиска (не говоря уже о том, насколько это невозможно для предиктора ветвления). Ваш count Решение имеет 5 поисков по карте и 5 обходов каждой строки, в то время как ваше ручное решение имеет поиск по карте за каждый нуклеотид (но только один обход строки).

   Серьезное предложение: используйте локальная переменная для каждого счетчика. Они почти наверняка будут помещены в регистры процессора и, следовательно, по существу бесплатны. Вы никогда не будете загрязнять свой кеш счетчиками таким образом, в отличие от map, unordered_map, vector и т.п.
   Замена абстракции на повторение, как это обычно не очень хорошая идея, но в этом случае довольно немыслимо, что вам когда-либо понадобится значительно больше счетчиков, поэтому масштабируемость не проблема.

3. Рассматривать std::string_view (что потребует другого метода чтения файла) избегать создания копий данных. Вы загружаете все данные в память с диска и затем для каждой последовательности копируете их. Это на самом деле не нужно и (в зависимости от того, насколько умный ваш компилятор) может вас утомить. Тем более, что вы продолжаете добавлять к строке до следующего заголовка (что является более ненужным копированием — вы можете просто считать после каждой строки).

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

Кстати, boost::count это тонкая обертка вокруг std::count по крайней мере в этой версии.

Я думаю, что вы все сделали правильно: написание хорошего и удобочитаемого кода, затем определение его как узкого места в производительности и проверка, сможете ли вы сделать его быстрее (возможно, сделав его немного уродливым).

Причина, по которой он настолько медленный, заключается в том, что у вас есть косвенный доступ все время или 5 сканирований одной и той же строки.

Вам не нужна карта, используйте 5 целых чисел и увеличивайте их отдельно. Тогда это должно быть быстрее, чем boost::count версия, потому что вы не пересекаете строку 5 раз, и она будет быстрее, чем map или unordered_map увеличивается, потому что у вас не будет n косвенных доступов.

так что используйте что-то вроде:

switch(char)
{
case 'A':
++a;
break;
case 'G':
++g;
break;
}
...

Как люди, предложенные в комментариях, попробуйте что-то подобное

enum eNucleotide {
NucleotideA = 0,
NucleotideT,
NucleotideC,
NucleotideG,
NucleotideN,
Size,
};

void countSequence(std::string line)
{
long nucleotide_counts[eNucleotide::Size] = { 0 };

if(line[0] != '>') {
for(int i = 0; i < line.size(); ++i)
{
switch (line[i])
{
case 'A':
++nucleotide_counts[NucleotideA];
break;
case 'T':
++nucleotide_counts[NucleotideT];
break;
case 'C':
++nucleotide_counts[NucleotideC];
break;
case 'G':
++nucleotide_counts[NucleotideC];
break;
case 'N':
++nucleotide_counts[NucleotideN];
break;
default :
/// error condition
break;
}
}

/// print results
std::cout << "A: " << nucleotide_counts[NucleotideA];
std::cout << "T: " << nucleotide_counts[NucleotideT];
std::cout << "C: " << nucleotide_counts[NucleotideC];
std::cout << "G: " << nucleotide_counts[NucleotideG];
std::cout << "N: " << nucleotide_counts[NucleotideN] << std::endl;
}
}

и вызывать эту функцию для каждого содержимого строки. (Не проверял код.)

Если это основная задача, которую вам нужно выполнить, возможно, вас заинтересует решение awk. Различные проблемы с файлами FASTA очень легко решаются с помощью awk:

awk '/^>/ && c { for(i in a) if (i ~ /[A-Z]/) printf i":"a[i]" "; print "" ; delete a }
/^>/ {print; c++; next}
{ for(i=1;i<=length($0);++i) a[substr($0,i,1)]++ }
END{ for(i in a) if (i ~ /[A-Z]/) printf i":"a[i]" "; print "" }' fastafile

Это выводит на ваш пример:

>Header 1
N:1 A:7 C:6 G:8 T:8
>Header 2
A:10 C:10 G:11 T:12

нота: Я знаю, что это не C ++, но часто полезно показать другие средства для достижения той же цели.

Тесты с awk:

• TestFile: http://hgdownload.cse.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz
• разархивированный размер: 2.3G
• всего записей: 5502947
• Всего строк:

Сценарий 0: (время выполнения: слишком долго) Первый упомянутый сценарий крайне медленный. Использовать только на небольших файлах

Сценарий 1: (время выполнения: 484,31 с) Это оптимизированная версия, где мы делаем целевой счет:

/^>/ && f { for(i in c) printf i":"c[i]" "; print "" ; delete c }
/^>/ {print; f++; next}
{   s=$0
c["A"]+=gsub(/[aA]/,"",s)
c["C"]+=gsub(/[cC]/,"",s)
c["G"]+=gsub(/[gG]/,"",s)
c["T"]+=gsub(/[tT]/,"",s)
c["N"]+=gsub(/[nN]/,"",s)
}
END { for(i in c) printf i":"c[i]" "; print "" ; delete c }

Обновление 2: (время выполнения: 416,43 сек) Объедините все подпоследовательности в одну последовательность и сосчитайте только одну:

function count() {
c["A"]+=gsub(/[aA]/,"",s)
c["C"]+=gsub(/[cC]/,"",s)
c["G"]+=gsub(/[gG]/,"",s)
c["T"]+=gsub(/[tT]/,"",s)
c["N"]+=gsub(/[nN]/,"",s)
}
/^>/ && f { count(); for(i in c) printf i":"c[i]" "; print "" ; delete c; string=""}
/^>/ {print; f++; next}
{ string=string $0 }
END { count(); for(i in c) printf i":"c[i]" "; print "" }

Обновление 3: (время выполнения: 396,12 с) Уточните, как awk находит свои записи и поля, и злоупотребляйте этим за один раз.

function count() {
c["A"]+=gsub(/[aA]/,"",string)
c["C"]+=gsub(/[cC]/,"",string)
c["G"]+=gsub(/[gG]/,"",string)
c["T"]+=gsub(/[tT]/,"",string)
c["N"]+=gsub(/[nN]/,"",string)
}
BEGIN{RS="\n>"; FS="\n"}
{
print $1
string=substr($0,length($1)); count()
for(i in c) printf i":"c[i]" "; print ""delete c; string=""}

Обновление 4: (время выполнения: 259,69 сек) Обновите поиск регулярных выражений в gsub, Это создает достойное ускорение:

function count() {
n=length(string);
gsub(/[aA]+/,"",string); m=length(string); c["A"]+=n-m; n=m
gsub(/[cC]+/,"",string); m=length(string); c["C"]+=n-m; n=m
gsub(/[gG]+/,"",string); m=length(string); c["G"]+=n-m; n=m
gsub(/[tT]+/,"",string); m=length(string); c["T"]+=n-m; n=m
gsub(/[nN]+/,"",string); m=length(string); c["N"]+=n-m; n=m
}
BEGIN{RS="\n>"; FS="\n"}
{
print ">"$1
string=substr($0,length($1)); count()
for(i in c) printf i":"c[i]" "; print ""delete c; string=""}