DIALS でたくさんのデータセットを一気に処理する

構造

small wedge で測定した大量のデータセットを DIALS で一気に処理してみた。まず、wedge を検出して、wedge ごとに作業フォルダを作る。

TOPDIR=/path/to/topdir/

find $TOPDIR -name "*_000001.img" -a ! -name "scan*.img" > targets.lst

i=0
while read fn; do
   echo "START $fn"
   mkdir -p data-$i
   cd data-$i
   echo dials.import `find ../${fn%000001.img}*.img ! -name "*_000000.img"` > run.sh
   i=$((i + 1))
   echo "DONE $fn"
   cd ..
done < targets.lst

これで、data-0, data-1, ... というフォルダができるので、各フォルダの中で run.sh を実行し、dials.import を行う。

for d in data-*; cd $d; do bash run.sh; cd ..; done

面倒なことに、000000 番はデータセットではないので除外する必要がある。それでも、dials.import は ? の桁数を揃える必要がないので、三文字.三文字生成器より使いやすい。

なぜか、上のスクリプトの while ループの中で dials.import を行うと、ループが 1 回しか回らないという不思議な現象が起こった。これは dials.import が標準入力を最後まで読むという仕様になっていて、read に渡るべき targets.lst の中身を全部消費してしまうためである。

   find ${fn%000001.img}*.img ! -name "*_000000.img" | dials.import

とすれば、期待した動作をする。

あとは、いつもどおり、find_spot, index, integrate, export_mtz と順番にやればよい。ここは時間がかかるので、parallel や qsub を使ってフォルダごとに並列化するのもよいだろう。

for d in data-*; do 
  cd $d
  dials.find_spots datablock.json global_threshold=200
  dials.index datablock.json strong.pickle unit_cell=100,110,120,90,90,90 space_group=P212121 indexing.method=fft1d index_assignment.method=local
  dials.integrate experiments.json indexed.pickle
  dials.export_mtz integrated_experiments.json integrated.pickle
  cd ..
done

うまくいったフォルダには hklout.mtz ができるので、これを BLEND などでクラスタリングしつつマージする。

cd blend
ls ../data-*/hklout.mtz > hkls.lst
blend -a hkls.lst <<EOF
GO
EOF

blend -s 5 <<EOF
RESOLUTION HIGH 2.5
TOLERANCE 100
GO
EOF

DIALS は精密化のときにカメラ長が変化しやすいようなので、small wedge では固定したほうがいいかもしれない。さもないと、格子定数のズレが結晶の不均一性を本当に反映しているのか分からなくなる。なお、上では TOLERANCE 100 として、格子のズレで AIMLESS が止まるのを防いでいる。

dials.export_mtz はデフォルトで partial reflection を捨ててしまうが、small wedge の場合は completeness への影響が大きい。include_partials=True をいれると、partial reflection を partiality で補正するようになるが、keep_partials=True もいれないと、partiality が 99% 未満のを捨ててしまうので、ほとんど効果がない。dials.export_mtz を改造して、partiality 0.5 以上の反射だけ出力するといった動作も可能にしたいところ。

なお、dials.export_mtz の出力する MTZ ファイル内で MPART カラム(PEAK に相当) は常に 1.00 になる。partiality 補正済みの値を出力しているから、意味的には正しいが、下流行程で partiality の情報を参照できないのは不便な気がする。