Быстрый подсчет типов нуклеотидов в большом количестве последовательностей - PullRequest
0 голосов
/ 05 ноября 2018

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

>Header 1  
ATGACTGATCGNTGACTGACTGTAGCTAGC  
>Header 2  
ATGCATGCTAGCTGACTGATCGTAGCTAGC  
ATCGATCGTAGCT

Таким образом, файл FASTA - это просто заголовок, перед которым стоит символ «>», затем последовательность в одной или нескольких строках, состоящая из нуклеотидов. Нуклеотиды - это символы, которые могут принимать 5 возможных значений: A, T, C, G или N.

Я хотел бы подсчитать, сколько раз встречается каждый тип нуклеотидов, поэтому, если мы рассмотрим этот фиктивный файл FASTA:

>Header 1  
ATTCGN

В результате я должен был получить:
A:1 T:2 C:1 G:1 N:1

Вот что я получил до сих пор:

ifstream sequence_file(input_file.c_str());
string line;
string sequence = "";
map<char, double> nucleotide_counts;

while(getline(sequence_file, line)) {
    if(line[0] != '>') {
        sequence += line;
    }
    else {
        nucleotide_counts['A'] = boost::count(sequence, 'A');
        nucleotide_counts['T'] = boost::count(sequence, 'T');
        nucleotide_counts['C'] = boost::count(sequence, 'C');
        nucleotide_counts['G'] = boost::count(sequence, 'G');
        nucleotide_counts['N'] = boost::count(sequence, 'N');
        sequence = "";
    }
}

Таким образом, он читает файл построчно, если он встречает '>' в качестве первого символа строки, он знает, что последовательность завершена, и начинает считать. Теперь проблема, с которой я сталкиваюсь, состоит в том, что у меня есть миллионы последовательностей с несколькими миллиардами нуклеотидов в общей сложности. Я вижу, что мой метод не оптимизирован, потому что я вызываю boost::count пять раз в одной и той же последовательности.

Другие вещи, которые я пробовал:

  • Анализ последовательности для увеличения счетчика для каждого типа нуклеотидов. Я попытался использовать map<char, double>, чтобы отобразить каждый нуклеотид на значение, но это было медленнее, чем раствор для повышения.
  • Использование std::count библиотеки алгоритмов, но это тоже было слишком медленно.

Я искал решения в Интернете, но каждое найденное мной решение было бы хорошим, если число последовательностей было низким, что не в моем случае. Не могли бы вы придумать, как можно ускорить процесс?

РЕДАКТИРОВАТЬ 1 : Я также пробовал эту версию, но она была в 2 раза медленнее, чем буст:

ifstream sequence_file(input_file.c_str());
string line;
string sequence = "";
map<char, double> nucleotide_counts;

while(getline(sequence_file, line)) {
    if(line[0] != '>') {
        sequence += line;
    }
    else {
        for(int i = 0; i < sequence.size(); i++) {
           nucleotide_counts[sequence[i]]++;
        }
        sequence = "";
    }
}

РЕДАКТИРОВАТЬ 2 : Спасибо всем в этой ветке, я смог добиться ускорения примерно в 30 раз по сравнению с оригинальным решением Boost. Вот код:

#include <map> // std::array
#include <fstream> // std::ifstream
#include <string> // std::string  

void count_nucleotides(std::array<double, 26> &nucleotide_counts, std::string sequence) {
    for(unsigned int i = 0; i < sequence.size(); i++) {
        ++nucleotide_counts[sequence[i] - 'A'];
    }
}  

std::ifstream sequence_file(input_file.c_str());
std::string line;
std::string sequence = "";
std::array<double, 26> nucleotide_counts;

while(getline(sequence_file, line)) {
    if(line[0] != '>') {
        sequence += line;
    }
    else {
        count_nucleotides(nucleotide_counts, sequence);
        sequence = "";
    }
}

Ответы [ 5 ]

0 голосов
/ 05 ноября 2018

В порядке значимости:

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

  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 как минимум в этой версии .

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

0 голосов
/ 05 ноября 2018

Не используйте карту, если вы хотите скорость и можете использовать массив. Кроме того, std::getline может использовать пользовательский разделитель (вместо \n).

ifstream sequence_file(input_file.c_str());
string sequence = "";
std::array<int, 26> nucleotide_counts;

// For one sequence
getline(sequence_file, sequence, '>');
for(auto&& c : sequence) {
    ++nucleotide_counts[c-'A'];
}

// nucleotide_counts['X'-'A'] contains the count of nucleotide X in the sequence

Демо

0 голосов
/ 05 ноября 2018

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

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;
    }
}

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

0 голосов
/ 05 ноября 2018

Если это основная задача, которую вам нужно выполнить, вас может заинтересовать решение 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:

Скрипт 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=""
}
0 голосов
/ 05 ноября 2018

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

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

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

switch(char)
{
case 'A':
    ++a;
    break;
case 'G':
    ++g;
    break;
}
...
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...