PDBファイルに変異を入れるプログラム(Ruby)。

研究室でちょっとしたプログラムを使う必要が出た。ので作った。

PDBファイルの特定の残基名を書き換えた後、AMBERパッケージを使って適切な構造に直して、分子動力学計算(MD計算)を行う。このPDBファイルの書き換えを自動でやってくれるプログラム（PDBも所詮はテキストファイルなので）。

BioRubyを使ってもいいんだが（というかその方が楽）、こんな時でもないとプログラム書かないから練習がてらに。まぁ何より、久しぶりにRubyってみたかった。...というわけで台風の土曜、家にこもって書き上げた変異体作成プログラム(不健康なので明日は出かけよう)。
100行足らずの短いプログラムだが、結構満足。次は変異入れたペプチドを組み合わせてダイマー化し、((20×20-20)/2)通りのヘテロ体を作る必要があるので、それも作ろうかな。

使い方。

引数として順に「変異を入れたいPDB」「残基番号」「変えたいアミノ酸」を指定する。

> pdbmutator.rb pdbfile target_residue_number mutation

mutationは大文字三字のアミノ酸略称を想定。
たとえば2zta.pdbというファイルの30残基目をグリシンにしたい場合は、

> pdbmutator.rb 2zta.pdb 30 GLY

一応変異と言っても非侵襲的で、元のPDBを残しながら新たなPDBファイルを作る。
変異体ファイルは、上の例で言うと2zta.pdb_30toGという表記になる。
自分の研究（数十残基のペプチドしか扱わない）用に作ったので、広く使うといろいろ問題も出そうだな。

#! /usr/bin/ruby
#
#   ---usage---
#
# > pdbmutator.rb PDBFILENAME TARGET_RSD MUTATION
#
#  <todo>  
#    - when raw pdb is loaded, strip header and footer...

Names = { 'ALA' => 'A', 'CYS' => 'C', 'ASP' => 'D', 'GLU' => 'E', 
          'PHE' => 'F', 'GLY' => 'G', 'HIS' => 'H', 'ILE' => 'I', 
          'LYS' => 'K', 'LEU' => 'L', 'MET' => 'M', 'ASN' => 'N', 
          'PRO' => 'P', 'GLN' => 'Q', 'ARG' => 'R', 'SER' => 'S', 
          'THR' => 'T', 'VAL' => 'V', 'TRP' => 'W', 'TYR' => 'Y' }

if ARGV[0] == nil || ARGV[1] == nil || ARGV[2] == nil
print("usage:"+"\n"+ 
 "> pdbmutator.rb PDBFILENAME TARGET_RSD MUTATION")
    exit    
end

require 'fileutils.rb'

filename = ARGV[0] 
workfile = filename + "_tmp"
FileUtils.cp(filename, workfile)

target_rsd = ARGV[1].to_i
mutation = ARGV[2] 

open(workfile){|a|

    residno = -1
    prev_residno = -2
    uni = File.open("#{filename}_#{target_rsd}to#{Names[mutation]}", "w")

########################    line loop start    ########################

    while line = a.gets
        column = line.split(/\s+/)

           ### a sample of read "column" ####
           ### ["ATOM", "12", "N", "MET", "A", "2", 
           ### "37.447", "16.947", "15.901", "1.00", "33.89", "N"]

        ### unless line is "TER" or "END"#
        if column[0] == "ATOM" then 
                
            ### Terminal resid names may be joined w/ the atom name.#
            ### In such a case, you need to fix it. #
            if column[2].length > 4 then
                tmpatom = column[2][0..2]
                tmpresidue = column[2][3..column[2].length]

                ### insert the "atom" in column[2] #
                ###  and replace column[3] to the real resid.#
                column[2,0] = tmpatom
                column[3] = tmpresidue
            end
        ### just a output ajustment.#
        ### inform the sequence now processing.#
        residno = column[5].to_i
        if residno != prev_residno
            print "#{residno}: #{column[3]}  "
            if (residno % 5 == 0) || residno == target_rsd
               print "\n"
            end
        end
#######################################################################
######################## when target resid No. ########################
#######################################################################

        if residno == target_rsd
             resid0 = column[3]
             line = line.sub(/#{resid0}/, "#{mutation}")
             puts "ATOM No. #{column[1]}: #{resid0} has been mutated to #{mutation}"
        end
#######################################################################

        end  # unless line is "TER" or "END"
        uni.puts <<-"EOB"
#{line.strip}
        EOB
        prev_residno = residno
    end
########################    line loop end     ########################
}