Python/Bash实现蛋白质残基互作统计与展示及报错修复

问题概述 #

需要处理记录蛋白质残基间相互作用的文本文件，每行包含互作类型（如sb、pc、vdw，hb子类如hbbb需统一归为hb）和两个残基，最终按指定格式统计每组残基对的各类互作次数。用户原有Python脚本因键生成逻辑错误触发ValueError。

原始数据示例：

sb ASP-11 LYS-15
sb GLU-309 HIS-46
pc HIS-46 TYR-45
vdw ASP-270 LEU-273

期望输出格式：

sb    pc    vdw    hb    residue1    residue2
2     0     0      0     ASP-11      LYS-15

报错原因 #

原有脚本用"-".join(sorted([residue1, residue2]))生成残基对的键，但残基本身包含-（如ASP-11），后续用split("-")拆分时会得到4个元素，无法赋值给2个变量，触发ValueError: too many values to unpack (expected 2)。

修复后的Python脚本 #

改用元组作为字典的键（元组可哈希且无需分隔符），同时处理hb子类型的归类逻辑：

# 定义存储统计结果的字典，键为排序后的残基对元组，值为各互作类型的计数
interaction_counts = {}

# 定义需要归为hb的子类型，可根据实际需求补充
HB_SUBTYPES = {"hbbb", "hbbs", "hbss", "hbsb"}

with open("file1.txt", "r") as file:
    for line in file:
        line = line.strip()
        if not line:
            continue  # 跳过空行
        parts = line.split()
        if len(parts) != 3:
            continue  # 跳过格式错误的行
        interaction_type, res1, res2 = parts
        
        # 处理hb子类型，统一归为hb
        if interaction_type in HB_SUBTYPES:
            interaction_type = "hb"
        
        # 生成排序后的残基对元组作为唯一键，确保无序残基对被视为同一组
        res_pair = tuple(sorted([res1, res2]))
        
        # 更新计数：初始化残基对的计数模板，避免后续get操作
        if res_pair not in interaction_counts:
            interaction_counts[res_pair] = {"sb":0, "pc":0, "vdw":0, "hb":0}
        interaction_counts[res_pair][interaction_type] += 1

# 打印对齐的表头
print(f"{'sb':<6}{'pc':<6}{'vdw':<6}{'hb':<6}{'residue1':<12}{'residue2':<12}")

# 遍历输出统计结果
for (res1, res2), counts in interaction_counts.items():
    print(f"{counts['sb']:<6}{counts['pc']:<6}{counts['vdw']:<6}{counts['hb']:<6}{res1:<12}{res2:<12}")

脚本说明 #

1. 用元组存储排序后的残基对，彻底避免分隔符冲突，同时保证无序残基对（如ASP-11 LYS-15和LYS-15 ASP-11）被视为同一组；
2. 新增HB_SUBTYPES集合，统一处理hb子类型的归类；
3. 初始化每个残基对的计数为0，简化后续计数更新逻辑；
4. 用格式化字符串保证输出列对齐，更贴合示例格式。

可选Bash脚本实现 #

用awk实现相同功能，适合处理超大型文本文件：

#!/bin/bash

awk '
BEGIN {
    # 打印对齐的表头
    printf "%-6s%-6s%-6s%-6s%-12s%-12s\n", "sb", "pc", "vdw", "hb", "residue1", "residue2"
    # 定义hb子类型映射
    hb_subtypes["hbbb"]=1; hb_subtypes["hbbs"]=1; hb_subtypes["hbss"]=1; hb_subtypes["hbsb"]=1
}

# 处理有效数据行
NF == 3 {
    type=$1; res1=$2; res2=$3
    
    # 统一hb子类型
    if (type in hb_subtypes) type="hb"
    
    # 生成排序后的残基对键，确保无序对归一
    key = (res1 < res2) ? res1 "\t" res2 : res2 "\t" res1
    
    # 更新对应类型的计数
    counts[key][type]++
}

# 输出最终统计结果
END {
    for (key in counts) {
        split(key, res, "\t")
        # 为空的计数补0
        sb = counts[key]["sb"] ? counts[key]["sb"] : 0
        pc = counts[key]["pc"] ? counts[key]["pc"] : 0
        vdw = counts[key]["vdw"] ? counts[key]["vdw"] : 0
        hb = counts[key]["hb"] ? counts[key]["hb"] : 0
        printf "%-6d%-6d%-6d%-6d%-12s%-12s\n", sb, pc, vdw, hb, res[1], res[2]
    }
}
' file1.txt

内容的提问来源于stack exchange，提问作者Rohan Nath