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

按重要性排序：

1. 此任务的好代码将100%被I/O绑定。您的处理器计算字符的速度比磁盘将字符泵送到CPU的速度快得多。因此，我面临的第一个问题是：您的存储介质的吞吐量是多少？您理想的RAM和缓存吞吐量是多少？这些都是上限。如果你已经碰到了它们，那么进一步研究你的代码就没有多大意义了。你的助推解决方案可能已经存在了。

2. CCD_ 5查找相对昂贵。是的，它是O(log(N))，但你的N=5很小，而且是恒定的，所以这个告诉你什么。对于5个值，映射将不得不为每次查找寻找大约三个指针(更不用说这对于分支预测器来说是多么不可能)。您的count解决方案有5个映射查找和每个字符串的5次遍历，而您的手动解决方案对每个核苷酸都有一个映射查找(但字符串只有一次遍历)。

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

3. 考虑std::string_view(这将需要不同的读取文件的方法)以避免创建数据的副本。你可以将整个数据从磁盘加载到内存中，然后对每个序列进行复制。这并不是真正必要的，而且(取决于你的编译器有多聪明)可能会让你陷入困境。特别是因为你一直在字符串后面追加，直到下一个标题(这是不必要的复制——你可以在每一行之后计数)。

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

顺便说一句，至少在这个版本中，boost::count是std::count的薄包装。

不过，我认为你在这里做了正确的事情：编写好的可读代码，然后将其识别为性能瓶颈，并检查是否可以让它运行得更快(可能会让它稍微难看一点)。

如果这是您必须执行的主要任务，那么您可能对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($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的基准：

• 测试文件：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=""
}

如果你想要速度并且可以使用数组，就不要使用映射。此外，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

演示

它之所以如此缓慢，是因为您一直都有间接访问或对同一字符串进行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;
}
}

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