Perl のハッシュを利用する例題

PDBをFASTAに変換する

PDB のファイルはタンパク質の構造情報やそのアミノ酸配列が書かれている。次の例は、PDB のファイルからアミノ酸配列の情報だけを抽出し、FASTA 形式で出力する処理である。

アミノ酸情報は PDB ファイルの中で、SEQRES レコードと ATOM レコードに保存されている。SEQRES レコードにはタンパク質を構成するすべてのアミノ酸が記載されている。これに対し、ATOM レコードには、立体構造(空間座標)のわかっているアミノ酸だけが記載されている。従って、両者を比べると、しばしば ATOM レコードに記載されるアミノ酸の数が少ない。

次のスクリプトは、PDB ファイルの ATOM レコードからアミノ酸を取得して、FASTA ファイルに関する Perl スクリプトである。タンパク質を構成する全アミノ酸を取得できないが、立体情報を持つアミノ酸のみを取得できるため、ホモロジーモデリングソフトウェア MODELLER の入力ファイルなどを作成するときに役立つ。

use strict;

 
my $pdb = '1ALK.pdb';
my $fa = {};
 
my $amino_code = {
     'ALA' => 'A', 'ARG' => 'R', 'ASN' => 'N', 'ASP' => 'D',
     'CYS' => 'C', 'GLN' => 'Q', 'GLU' => 'E', 'GLY' => 'G',
     'ILE' => 'I', 'LEU' => 'L', 'LYS' => 'K', 'MET' => 'M',
     'PHE' => 'F', 'PRO' => 'P', 'SER' => 'S', 'THR' => 'T',
     'TRP' => 'W', 'TYR' => 'Y', 'VAL' => 'V', 'HIS' => 'H',
     'ASX' => 'B', 'GLX' => 'Z', 'UNK' => 'K'
   };
 

open (my $fh, '<', $pdb) or die;
 
while (my $buff = <$fh>) {
  if ($buff !~ /^ATOM.{9}CA/) {
    next;
  }
 
  my $chain = substr($buff, 21, 1);
  my $num = substr($buff, 22, 4); 
  my $amino = substr($buff, 17, 3);
 
  if (!exists $fa->{$chain}) {
    $fa->{$chain} = [];
  }
 
  if (exists $amino_code->{$amino}) {
    $fa->{$chain}->[$num] = $amino_code->{$amino};
  } else {
    $fa->{$chain}->[$num] = 'X';
  }
}
 
close($fh);
 
$pdb =~ /(.+)\.[^\.]+$/;
my $pdb_id = $1;
 
 
open(my $fh, '>', $pdb_id . '.fa') or die;
foreach my $key (sort keys %{$fa}) {
  print $fh ">$pdb_id;$key\n";
    for (my $i = 1; $i < @{$fa->{$key}}; $i++) {
      if (!defined $fa->{$key}->[$i]) {
        print STDERR "$key:  $i\n";
      }
      print $fh $fa->{$key}->[$i];
    }
  print $fh "\n";
}

上のコードを pdb2fa.pl として保存して、ターミナル次のように実行できる。なお、PDB ファイルは PDB 1ALK からダウンロードできる。

perl pdb2fa.pl