基于SNV列表的基因组重组Python实现性能优化问询

问题背景 #

编写种群模拟程序时，个体基因组的两条染色体以有序单核苷酸变异（SNV）列表（实际用numpy浮点数组存储）表示，重组规则如下：

示例：SNV列表[0, 3, 6, 9]和[1, 2, 5, 11]，交叉点[4,10]。第一次交叉后得到[0,3,5,11]和[1,2,6,9]；第二次交叉后得到[0,3,5]和[1,2,6,9,11]。

现有两种Python实现，但随着世代推进SNV数量增多，函数性能持续下降成为瓶颈，需优化实现。

现有实现的瓶颈分析 #

无论初始的列表拼接还是优化后的切片交换，核心问题在于每次交叉操作都涉及大量元素的内存复制/移动：

• bisect操作是O(log n)复杂度，但切片/拼接操作是O(k)（k为切片长度）
• 多次交叉循环累积后，时间复杂度会趋近于O(m*n)（m为交叉点数量，n为SNV总数），当n增大时性能急剧下降

高效优化方案 #

1. Numpy向量化操作（优先推荐） #

利用numpy的C底层实现替代Python原生列表操作，内存连续性更好，向量化操作效率远超Python循环。

优化思路

• 预排序交叉点，避免重复排序开销
• 用np.searchsorted一次性计算所有交叉点在SNV数组中的分割索引
• 一次性合并所有分段，减少多次内存复制的开销

代码实现

import numpy as np

def recombine_numpy(snvs_0, snvs_1, crossovers):
    # 预排序交叉点，确保处理顺序正确
    crossovers = np.sort(crossovers)
    m = len(crossovers)
    
    # 获取所有分割索引（包含首尾边界）
    idx0 = np.concatenate([[0], np.searchsorted(snvs_0, crossovers), [len(snvs_0)]])
    idx1 = np.concatenate([[0], np.searchsorted(snvs_1, crossovers), [len(snvs_1)]])
    
    res0_segments = []
    res1_segments = []
    
    # 交替取对应亲本的分段，构建结果
    for i in range(m + 1):
        if i % 2 == 0:
            # 偶数段：取原亲本的对应区间
            res0_segments.append(snvs_0[idx0[i]:idx0[i+1]])
            res1_segments.append(snvs_1[idx1[i]:idx1[i+1]])
        else:
            # 奇数段：交换取另一个亲本的对应区间
            res0_segments.append(snvs_1[idx1[i]:idx1[i+1]])
            res1_segments.append(snvs_0[idx0[i]:idx0[i+1]])
    
    # 合并所有分段，生成最终结果
    return np.concatenate(res0_segments), np.concatenate(res1_segments)

2. 使用SortedList数据结构 #

若需频繁进行插入、删除和二分查找操作，sortedcontainers库的SortedList比普通列表更高效——它采用分段存储，切片和拼接时的内存复制开销远低于原生列表。

代码实现

from sortedcontainers import SortedList

def recombine_sortedlist(snvs_0, snvs_1, crossovers):
    sl0 = SortedList(snvs_0)
    sl1 = SortedList(snvs_1)
    
    for crossover_point in crossovers:
        idx0 = sl0.bisect_left(crossover_point)
        idx1 = sl1.bisect_left(crossover_point)
        
        # 交换后半段，SortedList的分段存储减少内存复制
        part0 = sl0[idx0:]
        part1 = sl1[idx1:]
        
        del sl0[idx0:]
        del sl1[idx1:]
        
        sl0.update(part1)
        sl1.update(part0)
    
    return list(sl0), list(sl1)

3. Cython/C扩展（极致性能） #

对于性能要求极高的场景，可将核心逻辑（二分查找、数组操作）用Cython或纯C实现，彻底规避Python循环的开销。通过直接操作内存数组，能将性能提升数倍至数十倍。

方案选择建议 #

• 若已使用numpy存储SNV，优先选择numpy向量化方案，无需额外依赖，性能提升显著
• 若需频繁动态修改SNV列表，选择SortedList方案，平衡查找和修改效率
• 若模拟规模极大（百万级SNV），考虑Cython/C扩展方案，压榨极致性能

内容的提问来源于stack exchange，提问作者Noah Jornet Calomarde