求助：PDB文件TER行后原子与残基编号批量修改问题

需求说明 #

要修改肽-蛋白复合物的PDB文件，第一个TER行之后的所有行必须满足：

• 第6列（残基编号）：原起始值1改为从6开始递增
• 第2列（原子编号）：原起始值1改为从53开始递增

自己写的bash脚本改完后，原子编号没更新，而且Pymol里显示不正常，求修复。

输入样本 #

ATOM     51  O   ARG     4      18.189  21.505 -30.356  0.00  0.00
ATOM     52  OXT ARG     5      19.822  21.322 -27.773  0.00  0.00
TER
ATOM      1  N   MET A   1      -9.976  22.279  65.378  1.00 37.35           N
ATOM      2  H   MET A   1      -9.180  21.915  65.882  1.00 37.35           H
ATOM      3  N   LYS A   2     -11.970  21.837  62.804  1.00 40.65           N
ATOM      4  H   LYS A   2     -11.194  21.438  62.295  1.00 40.65           H

期望输出 #

ATOM 51 O ARG 4 18.189 21.505 -30.356 0.00 0.00
ATOM 52 OXT ARG 5 19.822 21.322 -27.773 0.00 0.00
TER
ATOM 53 N MET A 6 -9.976 22.279 65.378 1.00 37.35 N
ATOM 54 H MET A 6 -9.180 21.915 65.882 1.00 37.35 H
ATOM 55 N LYS A 7 -11.970 21.837 62.804 1.00 40.65 N
ATOM 56 H LYS A 7 -11.194 21.438 62.295 1.00 40.65 H

现有脚本 #

#!/bin/bash

# read input file and output file names from command line arguments
input_file=complex.pdb
output_file=renum.pdb

# initialize residue counter and flag for tracking first "TER" occurrence
residue_num=1
ter_found=false

# loop through the lines of the input file
while read line
do
    # check if the line contains "TER"
    if [[ "$line" == "TER" ]]
    then
        # if it does, reset the residue counter to 5 and set the flag to true
        residue_num=5
        ter_found=true
    else
        # if it doesn't, extract the residue name and chain ID from the line
        residue_name=$(echo $line | awk '{print $4}')
        chain_id=$(echo $line | awk '{print $5}')

        # if the residue name or chain ID has changed, increment the residue counter
        if [[ "$residue_name" != "$prev_residue_name" || "$chain_id" != "$prev_chain_id" ]]
        then
            residue_num=$((residue_num+1))
        fi

        # if the first "TER" has been found, replace the 6th column with the new residue number
        if [[ "$ter_found" == true ]]
        then
            line=$(echo $line | awk -v num="$residue_num" '{$6=num; print}')
        fi

        # save the current residue name and chain ID for comparison in the next iteration
        prev_residue_name=$residue_name
        prev_chain_id=$chain_id
    fi

    # write the modified line to the output file
    echo $line >> $output_file
done < $input_file

实际输出 #

ATOM 51 O ARG 4 18.189 21.505 -30.356 0.00 0.00
ATOM 52 OXT ARG 5 19.822 21.322 -27.773 0.00 0.00
TER
ATOM 1 N MET A 6 -9.976 22.279 65.882 1.00 37.35 N
ATOM 2 H MET A 6 -9.180 21.915 65.882 1.00 37.35 H
ATOM 3 N LYS A 7 -11.970 21.837 62.804 1.00 40.65 N
ATOM 4 H LYS A 7 -11.194 21.438 62.295 1.00 40.65 H

原脚本为啥不行 #

1. 完全没处理原子编号：这是最明显的问题，脚本里根本没碰第2列的原子编号
2. 残基编号逻辑绕且不稳定：依赖空格分割列提取字段，PDB文件字段间空格数量不固定，容易出错
3. 破坏PDB格式：用echo $line | awk修改列会把固定空格变成任意空格，PDB是列位置敏感的格式，格式乱了Pymol肯定解析出错

修复后的脚本 #

#!/bin/bash

input_file="complex.pdb"
output_file="renum.pdb"

# 初始化计数器：TER前最后一个原子编号是52，残基编号初始为5（下一个残基从6开始）
atom_num=52
residue_num=5
ter_found=false
prev_residue_id=""
prev_chain=""

# 先清空输出文件，避免重复追加
> "$output_file"

while IFS= read -r line; do
    # 处理TER行
    if [[ "$line" == "TER"* ]]; then
        echo "$line" >> "$output_file"
        ter_found=true
        # 重置残基计数器，确保下一个残基从6开始
        residue_num=5
        continue
    fi

    # 只处理ATOM行，其他行直接输出
    if [[ "$line" == "ATOM"* ]]; then
        if [[ "$ter_found" == false ]]; then
            # TER之前的ATOM行直接输出，同时记录最后一个原子编号
            echo "$line" >> "$output_file"
            atom_num=$(echo "$line" | awk '{print $2}')
            continue
        fi

        # TER之后的ATOM行，开始修改编号
        # 用PDB固定列位置提取字段：链ID在第22位，残基编号在23-26位
        current_chain=${line:21:1}
        current_residue_id=${line:22:4}
        current_residue_id=$(echo "$current_residue_id" | xargs) # 去掉前后空格

        # 残基或链变化时，递增残基编号
        if [[ "$current_residue_id" != "$prev_residue_id" || "$current_chain" != "$prev_chain" ]]; then
            residue_num=$((residue_num + 1))
            prev_residue_id="$current_residue_id"
            prev_chain="$current_chain"
        fi

        # 原子编号每次都递增
        atom_num=$((atom_num + 1))

        # 按PDB标准格式重写行，保证列位置完全正确
        printf "%-6s%5d %-4s %3s %1s%4d    %8.3f%8.3f%8.3f%6.2f%6.2f          %s\n" \
            "ATOM" \
            "$atom_num" \
            "${line:12:4}" \
            "${line:17:3}" \
            "$current_chain" \
            "$residue_num" \
            "${line:30:8}" \
            "${line:38:8}" \
            "${line:46:8}" \
            "${line:54:6}" \
            "${line:60:6}" \
            "${line:76:2}" >> "$output_file"
    else
        # 非ATOM/TER行直接输出
        echo "$line" >> "$output_file"
    fi
done < "$input_file"

修复点说明 #

1. 原子编号处理：TER前先记录最后一个原子编号（52），TER后每一行原子编号自动+1，从53开始递增
2. 残基编号处理：TER后初始设为5，遇到新残基（残基ID或链ID变化）就+1，确保第一个残基是6
3. 格式兼容：用printf严格按照PDB的固定列宽输出，彻底解决Pymol解析异常的问题
4. 字段提取更可靠：直接通过字符串切片提取PDB固定位置的字段，比按空格分割准得多，不会因为字段间空格数量出问题

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