如何修改Python脚本实现按个体对汇总IBDLENGTH列值？

需求说明 #

对输入文件中第1列与第3列构成的个体对（无论重复出现多少次），汇总第9列IBDLENGTH的数值；个体对需唯一（A-B与B-A视为同一对），未重复的个体对也需保留在结果中。

输入文件格式示例 #

Sindhi_HGDP00171    0   Tunisian_39T    0   1   120437718   147097266   3.02    7.111
Sindhi_HGDP00183    1   Sindhi_HGDP00206    2   1   242708729   244766624   7.41    3.468
Sindhi_HGDP00183    1   Sindhi_HGDP00206    2   1   242708729   244766624   7.41    4.468
IBS_HG01768 2   Moroccan_MRA46  1   1   34186193    36027711    30.46   3.108
IBS_HG01710 1   Sardinian_HGDP01065 2   1   246117191   249120684   7.53    3.258
IBS_HG01768 2   Moroccan_MRA46  2   1   34186193    37320967    43.4    4.418

现有脚本问题 #

原脚本错误拆分了个体ID（将Sindhi_HGDP00171拆分为种群组Sindhi和样本ID），导致按种群组而非个体对汇总。若直接替换种群逻辑为个体逻辑却输出为空，通常是未处理A-B与B-A视为同一对的情况，或错误修改了列名引用。

修改后的完整脚本 #

import pandas as pd

# 定义列名
cols = ['ID1', 'HAP1', 'ID2', 'HAP2', 'CHR', 'STARTPOS', 'ENDPOS', 'LOD', 'IBDLENGTH']

# 加载数据（确保文件路径正确）
data = pd.read_csv("./Roma_Ref_All_sorted.txt", sep='\t', names=cols)

# 生成有序个体对：将每个个体对的两个ID按字典序排序，避免A-B和B-A被视为不同对
def get_sorted_pair(id1, id2):
    return tuple(sorted([id1, id2]))

data['sorted_pair'] = data.apply(lambda row: get_sorted_pair(row['ID1'], row['ID2']), axis=1)

# 按有序个体对分组，汇总IBDLENGTH的和
result = data.groupby('sorted_pair')['IBDLENGTH'].sum().round(3).reset_index()

# 拆分sorted_pair列为ID1和ID2列，整理输出格式
result[['ID1', 'ID2']] = pd.DataFrame(result['sorted_pair'].tolist(), index=result.index)
result = result.drop('sorted_pair', axis=1)[['ID1', 'ID2', 'IBDLENGTH']]
result = result.rename(columns={'IBDLENGTH': 'IBD_sum'})

# 保存结果
result.to_csv("./individual_pairs_sum_IBD.txt", sep='\t', index=False)

关键修改点说明 #

• 保留完整个体ID：移除原脚本中拆分ID为种群组的代码，直接使用原始的ID1和ID2列（完整个体标识）。
• 生成有序个体对：通过sorted()将每个个体对的两个ID按字典序排序，确保A-B和B-A被识别为同一对，避免重复汇总。
• 高效分组求和：使用pandas的groupby直接按有序个体对分组求和，替代原脚本的循环逻辑，更高效且不易出错。
• 结果格式化：将分组后的有序对拆分为ID1和ID2列，重命名汇总列为IBD_sum，匹配期望输出格式。

预期输出示例 #

ID1 ID2 IBD_sum
IBS_HG01710 Sardinian_HGDP01065 3.258
IBS_HG01768 Moroccan_MRA46 7.526
Sindhi_HGDP00171 Tunisian_39T 7.111
Sindhi_HGDP00183 Sindhi_HGDP00206 7.936

内容的提问来源于stack exchange，提问作者Gf.Ena