蛋白质序列-kappa参数计算算法解读

初学生物化学的人，第一次学到氨基酸的时候，首先涉及到的一般都是20种基本氨基酸侧链官能团的各种性质，

涉及到酸碱电离，一个逃不过的性质就是在不同PH环境条件下氨基酸的带电性质，到底是带正电还是带负电。

而kappa参数，就是这么一个用于描述氨基酸序列带电性质的参数。

kappa参数，首次提出是在2013年的PNAS，

华盛顿大学的pappu lab在业界内也是很有名的了

https://www.pnas.org/doi/10.1073/pnas.1304749110

这个参数其实很简单：

用于量化蛋白质序列变异中的电荷不对称性。

简答来说，kappa是一个归一化参数（范围在 [0, 1]），反映序列局部电荷涨落相对于整体电荷不对称性的偏差程度。

1，关键概念和一些符号说明

2，算法步骤解释

（1）计算整体电荷不对称性

首先是计算整段蛋白质序列的，还没有涉及到滑动窗口

分子是不是很眼熟，其实就是一段序列的NCPR（net charge per residue）；

同理，分母其实就是FCR（fraction of charged residue）

（2）划分序列为重叠 blob 段

这里其实就开始划分滑动窗口了，而且这个滑动窗口的大小size还是有意义的

（3）计算每个 blob 段的电荷不对称性

这里和（1）是对应的，（1）中是计算整段蛋白质序列的电荷不对称性，

此处（3）是计算每一小段滑动窗口的子序列的电荷不对称性

（4）计算平均平方偏差

其实就是计算方差，这个公式就是类比方差来理解，当然全局不对称性并不是局端不对称性list的均值，但是逻辑上可以类比方差来理解。

（5）计算归一化参数

疑点在于前面一步计算得到的平均平方偏差，问题在于这个值哪里来的那么多的取值可能，就是从哪一步让这个值的计算会有一个分布的输出，然后有所谓的这个值的max？

其实看定义就可以发现了，一个给定的全长序列，比如说长度为L，就可以有L！个其实就是L的全排列数个的序列可能性，也就是这里说的sequence variants，其实就是指的是我们的假设空间，或者说是所有可能的排列空间，也就是所有可能的序列世界，然后依据“极大似然”，我们只观测到了我们所有可能的序列当中的当前这个序列。

举一个简单的例子：

比如说一段全长为10的氨基酸序列，这么一条多肽，然后假设其中带正电的电荷数目是3，带负电的电荷数目也是3，不带电的电荷数目是4，

那么理论上这条多肽链，我是不是有10！种排列可能性，当然，如果我们只关注电荷排布的模式，也就是电荷分布的不对称性的话，那么正电荷3个相对如何全排列其实效果（至少对于这个参数计算而言）是等价的，同理负电荷的3个，不带电荷的4个，所以实际上的排列数是10！/（3！* 3！* 4！）种；

总之，我们有这么多种序列排列的可能性，然后在这些可能性之间，我们每一次只选择一种排列方式来计算这个偏差参数；

对于每一种排列方式，我再划分滑动窗口blob，然后再在对应的滑动窗口内计算局部blob序列的电荷不对称性，有了这个参数之后，我再计算所有blob与整体电荷不对称性的一个偏差，总之对于每一个序列排布的可能性，我就有一个偏差值。

然后因为我有很多个可能的序列排布（也就是可能性的逻辑模态），我就相当于有了一个分布，那么自然当前被观测到的（极大似然）作为分子上面的偏差参数取值，然后其中最大的作为分母上面的偏差参数取值。

（6）计算最终参数

综上，总的来说就是：

3，一些疑点、难点再强调

（1）平均偏差的取值分布来源

这一步的整体流程如下：

虽然原始文献中使用的是g=5、6，g对偏差参数的计算结果影响是次要的，主要还是序列排列对偏差参数的计算影响比较大。

逻辑上来说，就是：

为什么要遍历所有可能性的排列呢？

我们可以深入代码实现细节来看，此处选取一个计算序列氨基酸性质相关的工具，涉及到kappa的计算，

参考：https://github.com/Pappulab/localCIDER/blob/6eea6b5b570e55c2431b2a063135c30c6400467a/localcider/backend/sequence.py#L1144

可以发现关键在于delta和deltaMax这两个子函数的实现；

其中delta的实现：

我们可以发现，其实这里就已经按照最终算法中的实现，即文献中的描述，将blob=5和6的结果取了个均值，我们知道delta只是分子，delta/2再计算出deltaMax之后（delta/2）/deltaMax，其实效果和（delta/deltaMax）是一致的。

这个就是电荷不对称性，果然就是NCPR的二次方再除以FCR；

然后整体的电荷不对称性，就是累加每个blob计算下来的偏差；

看上去我们的流程以及确定到了计算一个排列下的偏差参数的步骤，那么如何从分布中获取max的偏差参数呢，这一步又是否和我们上面设想的一致？
也就是我们跟着步骤计算出了delta，验证了我们前面想法的正确性；

那么接下来的delta max呢？这个其实才是关键：

参考工具中对于这个deltaMax的实现其实有点长，其实长是情有可原的，因为我们现在计算delta Max的问题是一个组合问题，序列越长其实越容易组合爆炸，我们一般实现的时候基本上都是启发式算法、用数值解析逼近的，直接穷举组合是不现实的，尤其是在算法实现中；

具体实现细节上，该工具代码实现中通过4种计算技巧，逐步逼近deltaMax，优先级从高到低：

1. 单电荷序列优化（仅含正电或负电残基 + 中性残基）。
2. 无中性残基的最大电荷分离（正负残基完全分离）。
3. 高中性残基数（≥18）的优化（中性残基集中在特定区域）。
4. 通用情况的中性残基分布搜索（有限遍历中性残基位置）。

核心思想：直接构造理论上 delta 最大的序列模式（如正负电荷完全分离），而非穷举所有排列。

技巧1：单电荷序列（仅正电或负电）

• 适用条件：序列中只有正电 或 负电残基（另一类电荷数为0）。
• 优化策略：
  • 将带电残基聚集成一个连续块（如 ++++），中性残基分布在其余位置。
  • 通过滑动带电块的位置，计算不同排列的 delta，保留最大值。
• 示例：
  • 序列：+++000（3正电，3中性）。
  • 测试排列：+++000, 0+++00, 00+++0, 000+++，选择 delta 最大的一个。

技巧2：无中性残基的最大电荷分离

• 适用条件：序列中 无中性残基（仅正电和负电残基）。
• 优化策略：
  • 将正电和负电残基分别聚集成连续块（如 ++++----）。
  • 滑动较短的电荷块（如正电块在负电块左侧或右侧），计算 delta。
• 示例：
  • 序列：+++---（3正电，3负电）。
  • 测试排列：+++---, ---+++，选择 delta 更大的一个。

技巧3：高中性残基数（≥18）的优化

• 适用条件：中性残基数 ≥18（长中性区域）。
• 优化策略：
  • 将中性残基分配到序列的 起始、中间、末尾 三个区域（每区域最多6个中性残基）。
  • 正电和负电残基分别聚集成块，置于中性区域之间。
• 示例：
  • 序列：000++++000000----000（中性残基分布在头、中、尾）。
  • 仅测试中性残基分布的有限组合，而非全排列。

技巧4：通用情况的中性残基分布搜索

• 适用条件：上述条件均不满足时的默认方法。
• 优化策略：
  • 将中性残基分为 起始、中间、末尾 三部分，但允许更灵活的分布。
  • 遍历中性残基在中间区域的数量（midNeuts），构造序列并计算 delta。
• 示例：
  • 序列：00+++00----00（中性残基分散在正负块之间）。

与理论方法的对比

理论方法	官方代码实现
遍历所有可能的序列排列（组合爆炸）。	仅构造电荷最不均匀的少数排列（启发式）。
计算精确的 delta_max。	逼近 delta_max，可能非全局最优。
适用于任意序列长度。	对长序列（尤其是高中性残基）更高效。

一致性

• 目标一致：寻找使 delta 最大的序列排列。
• 物理意义一致：delta_max 对应电荷分布最不均匀的排列（如正负电荷完全分离）。

差异性

• 穷举 vs. 启发式：
  理论需遍历所有排列（计算不可行），代码通过构造极值排列逼近最大值。
• 精确性：
  代码可能遗漏某些高 delta 排列（如复杂电荷交替模式），但对生物物理问题足够

（2）对电荷不对称性的理解

正如我前面所说的，无论是整体的电荷不对称性，还是局部blob序列的电荷不对称性，公式本质意义上都是一样的，分子是NCPR，分母是FCR。

我们从NCPR和FCR的意义上去理解这个所谓的电荷不对称性，

首先：NCPR是一段序列的净电荷net，FCR是带电氨基酸比例（如果带电单位是±1的话，正负电荷比例差/正负电荷比例和）

总的来说，电荷不对成性是一个有量纲的物理量，从结果上看

总之来说，这个物理量（参数），实际上是量化了电荷局部涨落（分布模式）对于整体电荷不对称性的一个影响指标，

也就是局部电荷分布不对称，导致整体的kappa参数接近于1.

4，算法伪代码补充

Algorithm: Calculate_Kappa
Input: sequence: 氨基酸序列 (长度 L)g_list = [5, 6]  // 两种 blob 尺寸
Output: κ: 归一化电荷不对称性参数Begin:// 1. 计算整体电荷不对称性 σN+ = count_positive(sequence)N- = count_negative(sequence)FCR = (N+ + N-)/LNCPR = (N+ - N-)/Lσ = (NCPR²)/FCRκ_sum = 0for each g in g_list:  // 对 g=5 和 g=6 分别计算// 2. 计算当前序列的 δδ_current = Calculate_Delta(sequence, g, σ)// 3. 计算最大 δ_max (优化策略)δ_max = Calculate_Delta_Max(L, N+, N-, g)// 4. 计算归一化 κ_gκ_g = δ_current / δ_maxκ_sum += κ_g// 5. 取平均得最终 κκ = κ_sum / len(g_list)return κ
End

其中组合优化的技巧方面：

Function: Calculate_Delta(seq, g, σ)
Begin:N_blob = L - g + 1total_sq_dev = 0for i = 0 to N_blob-1:blob = seq[i:i+g]  // 取 blob 子序列n+ = count_positive(blob)n- = count_negative(blob)FCR_i = (n+ + n-)/gNCPR_i = (n+ - n-)/gσ_i = (NCPR_i²)/FCR_i  // 局部电荷不对称性total_sq_dev += (σ_i - σ)²return total_sq_dev / N_blob
EndFunction: Calculate_Delta_Max(L, N+, N-, g)
Begin:// 优化策略：生成电荷分离序列if N+ == 0 or N- == 0:  // 单电荷类型charge_block = ('+' * N+) or ('-' * N-)return max over pos: Calculate_Delta('0'*pos + charge_block + '0'*(L-len(charge_block)-pos), g, σ)else if N_neut == 0:  // 无中性残基return max(max over sep: Calculate_Delta('+'*sep + '-'*N- + '+'*(N+-sep), g, σ),max over sep: Calculate_Delta('-'*sep + '+'*N+ + '-'*(N--sep), g, σ))else if N_neut ≥ 18:  // 高中性残基return max over start,end: Calculate_Delta('0'*start + '+'*N+ + '0'*mid + '-'*N- + '0'*end, g, σ)// 约束: start,end ≤ 6, mid = N_neut - start - endelse:  // 通用情况return max over mid,start: Calculate_Delta('0'*start + '+'*N+ + '0'*mid + '-'*N- + '0'*end, g, σ)// end = N_neut - start - mid
End