【发布时间】：2023-03-11 14:04:01
【问题描述】：
首先，关于我的问题的一些背景知识。
我是一名生物信息学家，这意味着我进行信息学处理以试图回答一个生物学问题。在我的问题中，我必须操作一个名为 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>将每个核苷酸映射到一个值，但这比 boost 解决方案要慢。 - 使用算法库的
std::count，但这也太慢了。
我在互联网上搜索了解决方案，但如果序列数量很少，我找到的每个解决方案都很好，而我的情况并非如此。你有什么想法可以帮助我加快速度吗？
编辑 1 ： 我也试过这个版本，但它比 boost 慢 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 ：感谢这个线程中的每个人，与 boost 原始解决方案相比，我能够获得大约 30 倍的速度提升。这是代码：
#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 = "";
}
}
【问题讨论】：
-
我不知道 boost::count 是做什么的，但似乎你遍历同一个字符串 5 次就足够了
-
只需单独查看每个字符并增加正确的计数器。
-
迭代序列一次，增加每个核苷酸的映射值，但不要使用访问速度很慢的
map，而是使用unordered_map或使用vector的大小26，并意识到，例如'C'-'A' == 2。 （当然你会在vector中浪费 19 个未使用的空间，但我敢打赌它很快，而且几乎没有多少内存可以浪费）。 -
如果数据集非常庞大，枚举所有可能的 4 个字符组合并一次比较 4 个字节可能会更快。
-
顺便说一句，我注意到您没有使用
unordered_map，而在打开文件时使用.c_str()，这让我怀疑您使用的是不是非常旧的C++版本。你知道你用的是哪个版本吗？