蛋白质序列-omega参数计算算法解读
参考我的上一篇博客:https://blog.csdn.net/weixin_62528784/article/details/149234912?spm=1001.2014.3001.5501
同样omega参数,也是一个衡量或描述氨基酸序列带电性质的参数,类似于kappa。
首次提出是在2016年的Journal of the American Chemical Society上,
首先是原文中对于该参数的描述,其实在算法上我们可以看到和kappa参数是非常相似的:
从物理意义上看,其实omega参数和kappa参数一样,都是描述氨基酸序列分布混合模式的参数,都是看混合得均匀不均匀、对称不对称,
只不过kappa参数是描述带正电的氨基酸和带负电的氨基酸混合得是否均匀,
然后omega参数是描述带电氨基酸/Pro与其他氨基酸的序列混合模式。
一,从计算流程也就是算法上看
1,都是先计算分布不对称性的指标,类似于kappa中的电荷分布不对称性质;
首先全局的可以计算一个分布不对称性,
然后局部的,都是选定1个序列窗口,类似于kappa中的blob,也就是滑动窗口,然后在滑动窗口里计算
带电氨基酸+pro的比例,以及其他氨基酸的比例,然后同样的可以计算一个分布不对称性;
2,然后就是局部不对称性与全局不对称性的一个偏差,然后计算一个类似于方差的项,用blob滑动窗口数做一个平均
3,然后就是组合问题,对于目前观测的序列,我们可以计算一个偏差值,对于所有排列情况下的序列,我们可以取到一个max的序列偏差,然后就是用这个max做一个类似于归一化的操作
4,最后和kappa参数计算类似,就是blob滑动窗口的设计有g=5或者是6的两种情况,所以再分别对归一化之后的计算取一个平均
在代码实现的细节上,其实我们可以看到:
实现过程中首先将一段序列降维表示成只有2字母的序列,其中RKDEP使用X表示,然后其余的序列只用O表示,这样计算比例参数的时候,以及后面组合优化的时候,其实就简洁很多了,
后续的组合优化以及算法实现,其实都可以直接照搬kappa参数实现的算法。
二,然后同样的算法伪代码实现部分
Algorithm: Calculate_Omega
算法:计算蛋白序列的Omega参数
Input: sequence: 氨基酸序列 (长度 L)g_list = [5, 6] // 两种blob尺寸
Output: Ω: 脯氨酸/电荷分布参数 [0,1]Begin:// 1. 残基分类class_A = {D,E,R,K,P} // 带电残基+脯氨酸class_B = 其他15种氨基酸// 2. 计算全局参数N_A = count(sequence, class_A)σ_G = (N_A/L - (L-N_A)/L)^2 / (N_A/L + (L-N_A)/L)Ω_sum = 0for g in g_list:// 3. 计算δδ = 0for i=0 to L-g:blob = sequence[i:i+g]n_A = count(blob, class_A)σ_i = (n_A/g - (g-n_A)/g)^2 / (n_A/g + (g-n_A)/g)δ += (σ_i - σ_G)^2δ /= (L-g+1)// 4. 计算δ_maxδ_max = max_delta(L, N_A, g)// 5. 计算Ω_gΩ_sum += δ/δ_maxreturn Ω_sum/2
End
综合算法的blob=5和6的考虑之后,获得的该参数omega的数学表达式为: