如何借助NumPy、Sklearn按列频率替换基因型数组字符串为0/2数值

问题背景 #

你持有尺寸为1826*5000的numpy.array矩阵，行对应样本、列对应特征，每行存储一条基因型数据，单个核苷酸以字符串形式存储，原始数据结构示例如下：

[['G' 'G' 'G' ... 'T' 'T' 'A']
 ['G' 'G' 'G' ... 'A' 'T' 'A']
 ['A' 'G' 'A' ... 'A' 'T' 'A']
 ...
 ['G' 'A' 'G' ... 'T' 'T' 'A']
 ['G' 'G' 'A' ... 'A' 'T' 'A']
 ['G' 'G' 'G' ... 'A' 'T' 'C']]

已知约束：每一列仅存在2种不同的核苷酸。
替换规则：逐列统计核苷酸出现频率，同列内出现频率更高的核苷酸替换为0，出现频率更低的核苷酸替换为2。例如第一列中'G'出现频率更高，替换为0，'A'替换为2，最终期望输出结构如下：

[['0' '0' '0' ... '2' '0' '0']
 ['0' '0' '0' ... '0' '0' '0']
 ['2' '0' '2' ... '0' '0' '0']
 ...
 ['0' '2' '0' ... '2' '0' '0']
 ['0' '0' '2' ... '0' '0' '0']
 ['0' '0' '0' ... '0' '0' '2']]

要求借助Numpy、Sklearn相关函数完成上述替换。

实现方法 #

核心逻辑分两步：

• 逐列统计频次，标记每列的高频核苷酸
• 逐元素匹配，高频位赋值0，低频位赋值2

纯Numpy实现（性能最优） #

针对5000列的规模，列级循环+列内向量化操作的运行速度极快，普通CPU上耗时不到1秒：

import numpy as np

# geno_arr 为你的原始字符串矩阵
# 初始化全0矩阵，用int8类型节省内存
res = np.zeros(geno_arr.shape, dtype=np.int8)

for col in range(geno_arr.shape[1]):
    current_col = geno_arr[:, col]
    # 统计当前列两种核苷酸的出现次数
    alleles, counts = np.unique(current_col, return_counts=True)
    # 定位高频核苷酸
    major_allele = alleles[counts.argmax()]
    # 非高频位置赋值2
    res[:, col][current_col != major_allele] = 2

# 如果需要输出字符串格式的'0'/'2'，执行类型转换即可
res_str = res.astype(str)

结合Sklearn的实现（适配流水线场景） #

如果后续需要接入Sklearn的预处理、建模流水线，可以用OrdinalEncoder实现，逻辑和规则完全一致：

import numpy as np
from sklearn.preprocessing import OrdinalEncoder

# 逐列生成类别顺序：[高频核苷酸, 低频核苷酸]
col_category_order = []
for col in range(geno_arr.shape[1]):
    current_col = geno_arr[:, col]
    alleles, counts = np.unique(current_col, return_counts=True)
    # 按频次从高到低排序，保证编码后0对应高频
    sorted_alleles = alleles[np.argsort(-counts)]
    col_category_order.append(sorted_alleles)

# 初始化编码器，指定每列的类别映射规则
encoder = OrdinalEncoder(categories=col_category_order, dtype=np.int8)
encoded_res = encoder.fit_transform(geno_arr)
# 编码器输出低频值为1，替换为2即可匹配要求
encoded_res[encoded_res == 1] = 2

# 需要字符串格式则做类型转换
res_str = encoded_res.astype(str)

说明：以上两种实现均严格匹配替换规则，因为题目明确每列仅存在2种核苷酸，不需要处理单核苷酸/3种以上核苷酸的异常场景。如果后续需要扩展支持等位基因频率计算、缺失值填充等操作，可以直接对接Sklearn的SimpleImputer等预处理组件。

内容的提问来源于stack exchange，提问作者Python NoHand