#floatcontents
* はっぴ〜☆Scripting!! [#j6377599]
#floatcontents
ということで、ここではちょっとの工夫で便利な機能、そんな簡単なスクリプトを紹介していきたいと思います。ルーチンなコマンドはスクリプトでハッピーになりましょう！~
メインに使用するスクリプト言語はシェルスクリプトのつもりですが、ちょっと複雑なものになるとPerlをつかうと思います。~
なお、単純なスクリプトの場合、エラー処理は省略しているのでその辺りは適当にヨロシク&huh;
- ''2009-05-08'' -- スクリプトの表示方法を変えてみました。表示が崩れるときはリロードしてみてください。
** ちょっと本気なスクリプト [#gc090844]
企画から使用までちょっと時間がかかる本気のスクリプト。
- [[コマンドラインでPRODRG>./コマンドラインでPRODRG]]
- [[cshをまぁまぁ使ってみよう]]
- [[pythonを思いっきり頑張って勉強してみよう]]
- [[CCP4ライブラリを使う>./CCP4ライブラリを使う]]
- [[TTCファイルを分解する>./TTCファイルを分解する]]
- [[PDB座標を操作する>./PDBファイルを操作しよう]]
** 特定のプログラム用のスクリプト [#r5f7edc3]
#ls2_1(Happy Scripting/CNS,title,new,display=hierarchy,relative=true,date)
** CROSSECでF', F"を求める -- 2009-05-07 [#yf9b4757]
*** 核種と波長から [#l9055cec]
異常分散項(F',F")はCCP4のCROSSECで簡単に求められますが、いちいち入力するのも面倒なものです。そこで簡単なシェルスクリプトを作ってみます。
#script(calcaf)
核種と波長を入力すればCROSSECの出力からF'とF"を読み取って出力します。
*** イメージから値を読み取り求める [#i1f4dae8]
calcafは核種と波長を指定しましたが、これを応用して、カレントディレクトリ以下から001.imgで終わるファイルを探し、ファイルヘッダ中のWAVELENGTHから波長を取得、異常分散項を求めてみます。WAVELENGTHパラメータのないイメージ、例えばR-AXISのファイルではヘッダがバイナリで読み取りが面倒なので対応していません。ADSCのイメージなら問題ないと思います。
#script(getafs)
- 使用方法
% 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サイトからファイルを取得する [#l729905b]
: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 [#s56805ca]
: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ファイル中の特定の値以上の原子を得る [#v6093bac]
テキスト処理言語の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の結果の確認>CNS/CNSを使ってみよう/コンポジットオミットマップを描く#t29ddd14]]でもawkを使って処理しています。
** MAPMANを使って電子密度を変換する [#mfebc4ae]
:cnsmap2ccp4|
電子密度を変換するのにいちいちMAPMANを起動するのもめんどくさい。そんなわけでシェルスクリプトにしてみました。~
もうちょっと親切にするなら引数のチェックぐらい入れたらいいかもしれませんね。
- スクリプト
#!/bin/sh
mapman <<EOS
re m1 $1 cns
wr m1 $2 ccp4
quit
EOS
- 使用方法
% cnsmap2ccp4 refine_2fofc.map 2fofc.map
** SEGIDを付加する [#d7f08940]
: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
*** ワンライナーで書いてみた [#b1541427]
せっかくなのでワンライナー(コマンドライン上でスクリプトを書くことだよ)で書いてみた。
perl -lane 'substr($_,72,1,substr($_,21,1)) if(/^(ATOM|HETATM)/); print' < hogehoge.pdb
これで上記と同じ結果が得られる。すごいなぁ。~
ワンライナーってめちゃ便利かも。SedとかAwkとか使えない私にとっては勉強の価値ありと思いました。
** MSEをMETに変換する [#i74d6a33]
: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をアミノ酸シークエンスに変換する [#rba6b71d]
#anno2(2011-12-15 : 公開)
プログラムに依存せずにPDBファイルをアミノ酸シークエンスに変換するスクリプトです。過去にざっくり作ったのでどの程度正確かは保証できませんが・・・&huh;~
インサーションなどある場合は怪しいかも知れません。また、ギャップは無視されるかも・・・~
改良提案お待ちしています(^^)。
#script(pdb2seq)
- 使用方法
% pdb2seq 2dsy.pdb
> Chain A
GMGTLTRYLEEAMARARYELIADEEPYYGEIPDLPGVWATGKSLKECEANLQAALEDWLLFLLSRGETPPPLGEVRIELPxxxxxxxxxxxxxxxxxxxxxxxxxxxx
> Chain B
MGTLTRYLEEAMARARYELIADEEPYYGEIPDLPGVWATGKSLKECEANLQAALEDWLLFLLSRGETPPPLGEVRIExxxxxxxxxxxxxxxxxxxxxxxxxxxx
同一Chain ID内のアミノ酸以外の原子(水など)はxになります。
