蛋白质ATOM与配体HETATM距离计算及Perl脚本问题排查

让我帮你一步步排查原脚本的问题，并且实现你需要的功能：

一、原脚本无输出的核心原因 #

你的脚本里有个致命逻辑错误：当遇到HETATM行时，你执行了$part++; next;，这直接跳过了后面的坐标提取代码，导致所有HETATM的坐标都没有被存入@points数组。最后遍历@points的时候数组是空的，自然没有任何输出。

另外还有个小问题：提取的坐标字符串带有空格，直接计算的话可能会因为字符串转数字不彻底出现错误，需要先清理空格。

二、修复并增强后的脚本 #

下面是修改后的脚本，实现了：

• 正确读取ATOM和HETATM的坐标
• 计算每一对ATOM与HETATM的距离
• 过滤掉所有与其他类型原子距离都>5的原子行（保留至少有一个配对距离≤5的行）

#!/usr/local/bin/perl
use strict;
use warnings;

# 检查输入参数
die "Usage: $0 <pdb_file>\n" unless @ARGV == 1;

open(my $in_fh, '<', $ARGV[0]) or die "Cannot open file: $!";

my (@atoms, @hetatms);
my (%keep_atom, %keep_hetatm); # 标记需要保留的原子行

# 第一步：读取所有ATOM和HETATM的行与坐标
while (my $line = <$in_fh>) {
    chomp $line;
    if ($line =~ /^ATOM/) {
        # 提取坐标并清理空格转成数字
        my $x = substr($line, 30, 8);
        my $y = substr($line, 38, 8);
        my $z = substr($line, 46, 8);
        # 转成数字自动去除前后空格
        push @atoms, { line => $line, x => +$x, y => +$y, z => +$z };
    } elsif ($line =~ /^HETATM/) {
        my $x = substr($line, 30, 8);
        my $y = substr($line, 38, 8);
        my $z = substr($line, 46, 8);
        push @hetatms, { line => $line, x => +$x, y => +$y, z => +$z };
    }
}
close $in_fh;

# 第二步：计算所有ATOM与HETATM的距离，标记需要保留的行
foreach my $atom (@atoms) {
    my $has_close_hetatm = 0;
    foreach my $hetatm (@hetatms) {
        my $dist = sqrt(
            ($atom->{x} - $hetatm->{x})**2 +
            ($atom->{y} - $hetatm->{y})**2 +
            ($atom->{z} - $hetatm->{z})**2
        );
        # 可选：打印距离信息（如果需要调试）
        # printf "ATOM line: %s\nHETATM line: %s\nDistance: %.2f\n\n", $atom->{line}, $hetatm->{line}, $dist;
        
        if ($dist <= 5) {
            $has_close_hetatm = 1;
            $keep_hetatm{$hetatm->{line}} = 1; # 标记该HETATM需要保留
        }
    }
    if ($has_close_hetatm) {
        $keep_atom{$atom->{line}} = 1; # 标记该ATOM需要保留
    }
}

# 第三步：输出保留的行（先ATOM，后HETATM，顺序和原文件一致）
print "保留的ATOM行:\n";
foreach my $atom (@atoms) {
    print "$atom->{line}\n" if exists $keep_atom{$atom->{line}};
}

print "\n保留的HETATM行:\n";
foreach my $hetatm (@hetatms) {
    print "$hetatm->{line}\n" if exists $keep_hetatm{$hetatm->{line}};
}

三、脚本说明 #

1. 数据存储：用哈希数组存储每个原子的原行信息和坐标，方便后续关联原行和计算距离。
2. 距离计算：遍历所有ATOM-HETATM对，计算欧氏距离。
3. 过滤逻辑：
   • 如果某个ATOM和至少一个HETATM的距离≤5，就保留这个ATOM行；
   • 如果某个HETATM和至少一个ATOM的距离≤5，就保留这个HETATM行；
   • 这样就避免了删除那些有近距离配对的原子行，同时去掉了所有没有近距离交互的原子。
4. 输出：分别输出保留的ATOM和HETATM行，保持原文件中的顺序。

四、测试脚本 #

用你提供的示例输入测试的话，脚本会计算每个ATOM和4个HETATM的距离，然后保留那些距离≤5的原子对对应的行。

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