如何用Bash/Awk脚本修改PDB文件TER行后的第二列数值？

给定以下PDB文件片段：

ATOM     52  OXT GLU     5     -23.455 -21.595  12.691  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

需要实现：将第一个TER行的前一行第二列数值加1，以此作为TER之后所有行第二列的起始值，后续每行依次递增1，目标效果如下：

ATOM 52 OXT GLU 5 -23.455 -21.595 12.691 0.00 0.00
TER
ATOM 53 N MET A 1 -9.976 22.279 65.378 1.00 37.35 N
ATOM 54 H MET A 1 -9.180 21.915 65.882 1.00 37.35 H

尝试了以下Bash脚本但无法正常运行：

#!/bin/bash

# Get the value of the second column in the last line before "TER"
last_col2=$(grep -v "TER" model.pdb | awk '{if($2>max){max=$2}}END{print max}')
new_col2=$((last_col2+1))

# Loop through the lines of the input file
while read line; do
  if [[ "$line" == "TER" ]]; then
    # Skip "TER" line, but print it to the output file
    echo "$line" >> model_new.pdb
  else
    # Increment the second and fifth columns for lines after "TER"
    # and print the updated line to the output file
    col2=$(echo "$line" | cut -c 7-11 | tr -d '[:space:]')
    if [[ "$col2" -eq 1 ]]; then
      new_line=$(echo "$line" | sed "s/$col2/$new_col2/")
      echo "$new_line" >> model_new.pdb
    else
      echo "$line" >> model_new.pdb
    fi
  fi
done < model.pdb

用Awk可以高效完成需求，避免Bash循环的低效和逻辑漏洞，脚本如下：

#!/usr/bin/awk -f

BEGIN {
    in_ter_section = 0
    next_id = 0
}

/^TER/ {
    print $0
    next_id = prev_id + 1
    in_ter_section = 1
    next
}

!in_ter_section {
    print $0
    prev_id = $2
    next
}

in_ter_section {
    $2 = next_id++
    # 若需严格保持PDB固定列宽格式，替换为下面的printf语句：
    # printf "%-6s%5d %-3s %-3s %1s%4d    %8.3f%8.3f%8.3f %6.2f%6.2f          %s\n", $1, $2, $3, $4, $5, $6, $7, $8, $9, $10, $11, $12
    print $0
}

脚本说明 #

1. 初始化：in_ter_section标记是否进入TER后的段落，next_id存储TER后续的起始编号。
2. TER行处理：打印TER行，计算后续起始编号为TER前一行的ID+1，标记进入后续段落。
3. TER前行处理：直接输出行内容，同时记录当前行的第二列值到prev_id。
4. TER后行处理：将第二列替换为next_id，随后next_id自增，输出修改后的行。

若需要严格匹配PDB文件的固定列宽格式，取消注释脚本中的printf语句即可，该语句对应示例片段的字段格式，可根据实际PDB结构调整。

使用方法 #

1. 将脚本保存为fix_pdb_ids.awk，赋予执行权限：

chmod +x fix_pdb_ids.awk

2. 运行脚本处理目标PDB文件：

./fix_pdb_ids.awk model.pdb > model_new.pdb

原脚本问题分析 #

你之前的Bash脚本存在几个关键问题：

• 用grep -v "TER" | awk找TER前的最大ID，逻辑错误——若TER前行的ID并非递增，会取全局最大值而非TER前一行的ID。
• 用cut -c 7-11提取第二列依赖固定列位置，格式变化时会出错；且sed替换可能误匹配其他位置的数字（比如第五列的1）。
• 仅处理了第二列为1的行，未实现后续行的连续递增编号。

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