快速计数大量序列中的核苷酸类型

Question

首先，介绍一下我的问题的背景。
我是一名生物信息学家，这意味着我进行信息学处理以尝试回答生物学问题。在我的问题中，我必须操作一个名为 FASTA 文件的文件，如下所示：

>Header 1  
ATGACTGATCGNTGACTGACTGTAGCTAGC  
>Header 2  
ATGCATGCTAGCTGACTGATCGTAGCTAGC  
ATCGATCGTAGCT

所以 FASTA 文件基本上只是一个 header，前面有一个 '>' 字符，然后是一行或多行由核苷酸组成的序列。核苷酸是可以取 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：我也试过这个版本，但它比 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 = "";
    }
}

Answer 1

之所以这么慢，是因为您一直在间接访问或对同一字符串进行 5 次扫描。

不需要映射，使用5个整数，分别递增。那么它应该比 boost::count 版本快，因为你不会遍历字符串 5 次，它会比 map 或 unordered_map 增量更快，因为你不会n 间接访问。

所以使用类似的东西：

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

Answer 2

像评论里的人建议的那样试试

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

并为每一行内容调用此函数。（未测试代码。）

Answer 3

如果这是您必须执行的主要任务，您可能会对 awk 解决方案感兴趣。使用 awk 可以很容易地解决 FASTA 文件的各种问题：

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([=10=]);++i) a[substr([=10=],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 的基准测试：

测试文件：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=[=12=]
    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 [=13=] }
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 
  string=substr([=14=],length()); 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 ">"
  string=substr([=15=],length()); count()
  for(i in c) printf i":"c[i]" "; print ""
  delete c; string=""
}

Answer 4

如果您想要速度并且可以使用数组，请不要使用地图。此外，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

Demo

Answer 5

按重要性排序：

此任务的好代码将 100% I/O-bound。您的处理器计算字符的速度比您的磁盘将字符泵送到 CPU 的速度快得多。因此，我的第一个问题是：您的存储介质的吞吐量是多少？您理想的 RAM 和缓存吞吐量是多少？这些是上限。如果您已经找到了它们，那么进一步查看您的代码就没有多大意义了。您的增强解决方案可能已经存在。
std::map 查找相对昂贵。是的，它是 O(log(N))，但是你的 N=5 很小而且不变，所以这告诉你什么都没有。对于 5 个值，地图每次查找都必须追踪大约三个指针（更不用说这对于分支预测器来说是多么不可能）。您的 count 解决方案对每个字符串进行 5 次映射查找和 5 次遍历，而您的手动解决方案对每个核苷酸 进行映射查找（但仅对字符串进行一次遍历）。

严重建议：为每个计数器使用一个局部变量。这些几乎肯定会放在 CPU 寄存器中，因此基本上是免费的。与 map、unordered_map、vector 等
不同，您永远不会用计数器污染缓存用像这样的重复来代替抽象通常不是一个好主意，但在这种情况下，很难想象你会需要更多的计数器，所以可伸缩性不是问题。

考虑std::string_view（这需要不同的文件读取方法）以避免创建副本 的数据。您将整个数据从磁盘加载到内存中，然后为每个序列复制它。这并不是真正必要的，并且（取决于你的编译器有多聪明）会让你陷入困境。特别是因为你一直附加到字符串直到下一个 header （这是更多不必要的复制 - 你可以在每一行之后计算）。

如果由于某种原因您没有达到理论吞吐量，请考虑多线程 and/or 矢量化。但我无法想象这是必要的。

顺便说一句，boost::count 是 std::count 的薄包装，至少 in this version。

不过我认为您在这里做的是对的：编写良好且可读的代码，然后将其确定为性能瓶颈并检查是否可以使它运行更快（可能通过使其稍微更难看）。

快速计数大量序列中的核苷酸类型

Fast counting of nucleotide types in a large number of sequences

c++

character

count

fasta