はっぴ〜☆Scripting!!
ということで、ここではちょっとの工夫で便利な機能、そんな簡単なスクリプトを紹介していきたいと思います。ルーチンなコマンドはスクリプトでハッピーになりましょう！
メインに使用するスクリプト言語はシェルスクリプトのつもりですが、ちょっと複雑なものになるとPerlをつかうと思います。
なお、単純なスクリプトの場合、エラー処理は省略しているのでその辺りは適当にヨロシク ![[huh]](image/face/huh.png)
- 2009-05-08 -- スクリプトの表示方法を変えてみました。表示が崩れるときはリロードしてみてください。
ちょっと本気なスクリプト
企画から使用までちょっと時間がかかる本気のスクリプト。
特定のプログラム用のスクリプト
- CNS -- 2009-05-18 (月) 18:45:43
CROSSECでF', F"を求める -- 2009-05-07
核種と波長から
異常分散項(F',F")はCCP4のCROSSECで簡単に求められますが、いちいち入力するのも面倒なものです。そこで簡単なシェルスクリプトを作ってみます。
-
1 2 3 4 5 6 7 8 9 10 11 12
#!/bin/bash # 異常分散項(F',F")をCROSSECで求める # 2009-05-08 by Nob-rin SPECIES=$1 SPECIES=`echo $SPECIES | tr "[a-z]" "[A-Z]"` WLEN=$2 echo "SPECIES WLEN F' F\"" crossec <<EOS | grep "^$SPECIES" ATOM $SPECIES NWAV 1 $WLEN EOS - 使用方法
calcaf ATOM WAVELENGTH
- 使用例
% calcaf se 0.9898 SPECIES WLEN F' F" SE 0.9898 -4.1403 0.5089
calcaf核種と波長を入力すればCROSSECの出力からF'とF"を読み取って出力します。
イメージから値を読み取り求める
calcafは核種と波長を指定しましたが、これを応用して、カレントディレクトリ以下から001.imgで終わるファイルを探し、ファイルヘッダ中のWAVELENGTHから波長を取得、異常分散項を求めてみます。WAVELENGTHパラメータのないイメージ、例えばR-AXISのファイルではヘッダがバイナリで読み取りが面倒なので対応していません。ADSCのイメージなら問題ないと思います。
-
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29
#!/bin/bash # カレントディレクトリ以下から001.imgを検索して、 # WAVELENGTHパラメータから異常分散項を求める # 2009-05-14 by Nob-rin declare -a FACTORS function set_factors () { SPECIES=$1 SPECIES=`echo $SPECIES | tr "[a-z]" "[A-Z]"` WLEN=$2 LST=$( crossec <<EOS | grep "^$SPECIES" ATOM $SPECIES NWAV 1 $WLEN EOS ) FACTORS=($LST) } [ "$1" ] || { echo "Usage: getafs SPECIES"; exit 1; } echo FILE/SP/WLEN/F\'/F\" for fn in `find ./ -name "*001.img"` do wlen=`strings $fn | sed -n 's/^WAVELENGTH=\([0-9\.]\+\);/\1/p'` [ "$wlen" ] || { echo "$fn has no WAVELENGTH parameter. Skip."; continue; } set_factors $1 $wlen echo $fn ${FACTORS[0]} $wlen ${FACTORS[2]} ${FACTORS[3]} done
- 使用方法
% getafs pt FILE/SP/WLEN/F'/F" ./edge/10pt_1_001.img PT 1.07207 -23.0902 10.1966 ./peak/10pt_1_001.img PT 1.07180 -19.2404 10.1923 ./hrem/10pt_1_001.img PT 1.05385 -10.9851 9.9026
カレントディレクトリ以下のファイルを検索して解析後出力します。この時、核種はコマンドラインで指定しておきます。
PDBコードでPDBサイトからファイルを取得する
- getpdb
- PDBファイルが必要なとき、Webブラウザをわざわざ開いて取得する手間を省きます。wgetをあらかじめインストールしておきます。
# apt-get install wget
- スクリプト
#!/bin/sh PDBSITE=http://www.pdb.org/pdb/download/downloadFile.do wget -O $1.pdb "$PDBSITE?fileFormat=pdb&compression=NO&structureId=$1"
- 使用方法
% getpdb 1PDB
PDBファイルから特定のChainのみ取得する -- 2009-06-23
- getchain
- PDBファイル中からCRYST1と指定したChain IDを切り出します。この時、水分子は無視するようにします。
- スクリプト
#!/bin/sh CID=$2 grep -E "^(CRYST1|(ATOM |HETATOM).{11}[A-Z0-9 ]{3} $CID )" $1 | grep -v " HOH "
- 使用方法
% getchain target.pdb B
PDBファイル中の特定の値以上の原子を得る
テキスト処理言語のawkを使ってPDBファイル中で温度因子が60を超える原子を出力してみます。
% awk '/^ATOM / {if($11>60){print}}' target.pdb
こんな感じで。システムによってはgawkを使うのかもしれません(Vine 4.2ではawkはgawkへのシンボリックリンクでした)。
awkのスクリプト部分は/^ATOM / {if($11>60){print}}の部分になります。先頭部分の正規表現/^ATOM /で原子レコードのみを処理することを指定します。HETATMも入れるなら/^(ATOM |HETATM)/にすれば全ての原子を対象にすることが可能です。続きの{if($11>60){print}}でアイテム11が60を超えるなら出力、という処理をしています。アイテムの11とはスペースで区切った値の11番目のものということで、PDBの温度因子が相当します。
応用例としてCNSのComposite omit mapの結果の確認でもawkを使って処理しています。
MAPMANを使って電子密度を変換する
- cnsmap2ccp4
- 電子密度を変換するのにいちいちMAPMANを起動するのもめんどくさい。そんなわけでシェルスクリプトにしてみました。
もうちょっと親切にするなら引数のチェックぐらい入れたらいいかもしれませんね。
- スクリプト
#!/bin/sh mapman <<EOS re m1 $1 cns wr m1 $2 ccp4 quit EOS
- 使用方法
% cnsmap2ccp4 refine_2fofc.map 2fofc.map
SEGIDを付加する
- addsegid
- CNSではSEGIDというPDB ver.2.3以降廃止されたフィールドを使用しますが、CCP4ではSEGIDは基本的には使わないので付加されません。それをスクリプトで付加します。
- スクリプト
#!/usr/bin/perl use strict; open(my $fh,$ARGV[0]) || die "Can't open $ARGV[0]: $!\n"; while(<$fh>){ if(/^(ATOM|HETATM)/){ substr($_,72,1) = sprintf("%1s",substr($_,21,1) } print; } close $fh; - 使用方法
% addsegid ccp4.pdb > cns.pdb
ワンライナーで書いてみた
せっかくなのでワンライナー(コマンドライン上でスクリプトを書くことだよ)で書いてみた。
perl -lane 'substr($_,72,1,substr($_,21,1)) if(/^(ATOM|HETATM)/); print' < hogehoge.pdb
これで上記と同じ結果が得られる。すごいなぁ。
ワンライナーってめちゃ便利かも。SedとかAwkとか使えない私にとっては勉強の価値ありと思いました。
MSEをMETに変換する
- mse2met
- PDBファイル中のSe-Met残基をMetに変換します。
- スクリプト
#!/usr/bin/perl use strict; my $filename = shift(@ARGV); open(my $fh,$filename) || die "Can't open $filename [$!]\n"; while(my $line = <$fh>){ if(substr($line,17,3) eq "MSE"){ $line =~ s/^HETATM/ATOM /; substr($line,17,3) = "MET"; if(substr($line,12,4) =~ /SE /i){ substr($line,12,4) = " SD "; substr($line,76,2) = " S"; } } print $line; } close $fh; 1; - 使用方法
% mse2met protein.pdb > out.pdb
PDBをアミノ酸シークエンスに変換する
2011-12-15 : 公開
プログラムに依存せずにPDBファイルをアミノ酸シークエンスに変換するスクリプトです。過去にざっくり作ったのでどの程度正確かは保証できませんが・・・
インサーションなどある場合は怪しいかも知れません。また、ギャップは無視されるかも・・・
改良提案お待ちしています(^^)。
-
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32
#!/usr/bin/perl use strict; my %AA = ( ALA => 'A', ARG => 'R', ASN => 'N', ASP => 'D', CYS => 'C', GLN => 'Q', GLU => 'E', GLY => 'G', HIS => 'H', ILE => 'I', LEU => 'L', LYS => 'K', MET => 'M', PHE => 'F', PRO => 'P', SER => 'S', THR => 'T', TRP => 'W', TYR => 'Y', VAL => 'V', MSE => 'M', ); my $file = $ARGV[0]; open(my $fh, $file) || die "Can't open $file [$!]\n"; my ($lastrno,$lastcid) = (-1,""); my %chain; while(<$fh>){ /^(ATOM |HETATM)/ || next; my $res = substr($_,17,3); my $rno = substr($_,22,4); my $cid = substr($_,21,1); if($lastrno != $rno || $lastcid ne $cid){ exists($chain{$cid}) || ($chain{$cid} = []); push(@{$chain{$cid}},exists($AA{$res}) ? $AA{$res} : 'x'); $lastrno = $rno; $lastcid = $cid; } } close $fh; foreach my $key (sort(keys(%chain))) { print "> Chain $key\n"; print join("",@{$chain{$key}}) . "\n"; }
- 使用方法
% pdb2seq 2dsy.pdb > Chain A GMGTLTRYLEEAMARARYELIADEEPYYGEIPDLPGVWATGKSLKECEANLQAALEDWLLFLLSRGETPPPLGEVRIELPxxxxxxxxxxxxxxxxxxxxxxxxxxxx > Chain B MGTLTRYLEEAMARARYELIADEEPYYGEIPDLPGVWATGKSLKECEANLQAALEDWLLFLLSRGETPPPLGEVRIExxxxxxxxxxxxxxxxxxxxxxxxxxxx
同一Chain ID内のアミノ酸以外の原子(水など)はxになります。
