我正在努力提高用C++编写的计算(生物)模型的速度(以前的版本在我的github:原核生物上)。最耗时的功能是计算单个基因组上转录因子和结合位点之间的结合亲缘关系。
背景:在我的模型中,结合亲和力是由转录因子结合区( 20-bool阵列)与结合位点序列(也是20-bool阵列)之间的汉明距离决定的。对于单个基因组,我需要计算所有活性转录因子(通常为5-10)与所有结合位点(通常为10-50)之间的亲缘关系。我每隔一步就对人口中的一万多个细胞进行基因表达状态的更新。为了模拟我的模型中的一个时间步长,需要对20个bool阵列进行多达50万次的比较,这意味着典型的实验需要几个月(200万到10m的时间步骤)。
对于以前版本的模型(链接以上),基因组仍然相当小,所以我可以计算每一个细胞的结合亲缘关系(在出生时),并存储和重复使用这些数字在细胞的生命。然而,在最新版本中,基因组大幅扩张,多个基因组驻留在同一细胞内。因此,在细胞中存储所有转录因子结合位点对的亲缘关系变得不切实际.
在当前的实现中,我定义了一个属于Bead类的内联函数(它是转录因子类“调节器”和绑定站点类“Bsite”的基类)。它直接写在头文件Bead.hh中:
inline int Bead::BindingAffinity(const bool* sequenceA, const bool* sequenceB, int seqlen) const
{
int affinity = 0;
for (int i=0; i<seqlen; i++)
{
affinity += (int)(sequenceA[i]^sequenceB[i]);
}
return affinity;
}上面的函数接受两个指向布尔数组(sequenceA和sequenceB)的指针,以及一个指定长度(seqlen)的整数。使用简单的for-循环,然后检查数组不同的位置(sequenceA[i]^sequenceB[i]),并将其加到变量affinity中。
给定基因组上的一个结合位点(bsite),然后我们可以遍历基因组,然后对每个转录因子(reg)计算它与这个特定结合位点的亲和力,如下所示:
affinity = (double)reg->BindingAffinity(bsite->sequence, reg->sequence);所以,这就是我是如何做到的;由于我没有编程背景,我想知道是否有更好的方法来编写上面的函数或构造代码(例如,BindingAffinity应该是基Bead类的函数)?非常感谢您的建议。
发布于 2022-01-10 10:17:57
感谢“保尔麦肯齐”和“艾克”,谢谢你的建议。我根据我以前的实现测试了这两种想法。以下是调查结果。总之,这两种答案都很有效。
我以前的实现为模型的1000个时间步骤提供了5m40 +/- 7 (n=3)的平均运行时。用GPROF分析显示,函数BindingAffinity()占运行时总数的24.3%。代码见问题。
位集实现的平均运行时为5m11 +/- 6 (n=3),相当于速度增长了~9%。在BindingAffinity()中只花费了总运行时的3.5%。
//Function definition in Bead.hh
inline int Bead::BindingAffinity(std::bitset<regulator_length> sequenceA, const std::bitset<regulator_length>& sequenceB) const
{
return (int)(sequenceA ^= sequenceB).count();
}
//Function call in Genome.cc
affinity = (double)reg->BindingAffinity(bsite->sequence, reg->sequence);位集实现的主要缺点是,与布尔数组(我以前的实现)不同,我必须指定进入函数的位集的长度。我偶尔会比较不同长度的位集,因此现在我必须指定不同的函数(根据https://www.cplusplus.com/doc/oldtutorial/templates/,模板不适用于多文件项目)。
对于整型实现,我尝试了@eike建议的std::popcount(seq1^seq2)函数的两种替代方法,因为我使用的是不包含此功能的较旧版本的C++。
备选案文1:
inline int Bead::BindingAffinity(int sequenceA, int sequenceB) const
{
int i = sequenceA^sequenceB;
std::bitset<32> bi (i);
return ((std::bitset<32>)i).count();
}备选案文2:
inline int Bead::BindingAffinity(int sequenceA, int sequenceB) const
{
int i = sequenceA^sequenceB;
//SWAR algorithm, copied from https://stackoverflow.com/questions/109023/how-to-count-the-number-of-set-bits-in-a-32-bit-integer
i = i - ((i >> 1) & 0x55555555); // add pairs of bits
i = (i & 0x33333333) + ((i >> 2) & 0x33333333); // quads
i = (i + (i >> 4)) & 0x0F0F0F0F; // groups of 8
return (i * 0x01010101) >> 24; // horizontal sum of bytes
}它们分别产生了5m06 +/- 6 (n=3)和5m06 +/- 3 (n=3)的平均运行时间,与我以前的实现相比,速度增加了10%。我只分析了备选方案2,它显示只有2.2%的运行时是在BindingAffinity()中使用的。将整数用于位字符串的缺点是,每当我更改任何代码时都必须非常小心。单位突变是绝对有可能的,正如@eike所提到的,但是每件事都有点棘手。
结论:用于比较位字符串的位集和整数实现都实现了令人印象深刻的速度改进。因此,BindingAffinity()不再是我代码中的瓶颈。
https://stackoverflow.com/questions/70436663
复制相似问题