autooperator

#!/usr/bin/env phenix.python
import sys
from iotbx import mtz, pdb
from cctbx import crystal

def main():
  if len(sys.argv) < 3:
    print("Usage: auto_operator.py  ")
    sys.exit(1)
  mtz_file = sys.argv[1]
  pdb_file = sys.argv[2]

  # MTZファイルから結晶対称性情報を取得
  mtz_obj = mtz.object(mtz_file)
  mtz_sym = mtz_obj.crystal_symmetry()

  # PDBファイルから結晶対称性情報を取得
  pdb_inp = pdb.input(file_name=pdb_file)
  pdb_sym = pdb_inp.crystal_symmetry()

  # もし両者が一致していれば、変換は不要（identity operator）
  if mtz_sym.unit_cell() == pdb_sym.unit_cell() and \
     mtz_sym.space_group().type().lookup_symbol() == pdb_sym.space_group().type().lookup_symbol():
    operator_str = "a,b,c"
  else:
    try:
      # PDBからMTZへの変換演算子を自動算出
      op = pdb_sym.change_of_basis_op(mtz_sym)
      operator_str = op.as_string()
    except Exception as e:
      print("自動算出に失敗しました: ", e)
      operator_str = "a,b,c"

  # 逆変換演算子を算出
  op = crystal.change_of_basis_op(operator_str)
  inv_op = op.inverse()
  inverse_operator_str = inv_op.as_string()

  print("Determined change-of-basis operator: {}".format(operator_str))
  print("Inverse operator: {}".format(inverse_operator_str))

if __name__=="__main__":
  main()

#!/usr/bin/env phenix.python
import sys, math
from iotbx import ccp4_map
from cctbx import maptbx
from scitbx.array_family import flex

def main():
    if len(sys.argv) < 3:
        print("Usage: auto_rebox.py input_map.ccp4 output_cubic_map.ccp4")
        sys.exit(1)
    input_map = sys.argv[1]
    output_map = sys.argv[2]
    
    # 元のマップ読み込み
    map_reader = ccp4_map.map_reader(input_map)
    map_data = map_reader.map_data.as_double()  # 3D flex.double array
    header = map_reader.header
    # 単位セルパラメータ (a,b,c,alpha,beta,gamma)
    a, b, c, alpha, beta, gamma = header.unit_cell_parameters
    nx, ny, nz = map_data.all()
    
    # 各軸の格子間隔
    spacing_x = a / nx
    spacing_y = b / ny
    spacing_z = c / nz
    # 解像度として最大の格子間隔を採用（保守的に）
    spacing = max(spacing_x, spacing_y, spacing_z)
    
    # 立方体とするため，新たなボックスサイズをL = max(a,b,c)に設定
    L = max(a, b, c)
    new_n = int(math.ceil(L / spacing))
    new_grid = (new_n, new_n, new_n)
    
    # 元のマップの物理中心を計算
    # 元のオフセット（header.origin）は、一般に物理的な原点位置（Å）として与えられている
    orig_origin = header.origin  # タプル (ox, oy, oz)
    center_phys = (orig_origin[0] + (nx*spacing_x)/2.0,
                   orig_origin[1] + (ny*spacing_y)/2.0,
                   orig_origin[2] + (nz*spacing_z)/2.0)
    
    # 新たな立方体マップのグリッド間隔は元と同じとし，物理中心を保つよう新オリジンを決定
    new_spacing = spacing
    new_origin = (center_phys[0] - new_n/2.0*new_spacing,
                  center_phys[1] - new_n/2.0*new_spacing,
                  center_phys[2] - new_n/2.0*new_spacing)
    
    # 新たな立方体マップ用の空配列を作成
    new_map = flex.double(flex.grid(new_grid))
    
    # 新たなグリッド各点の物理座標を計算し、元のマップから補間して新マップに格納（trilinear補間）
    for i in range(new_n):
        for j in range(new_n):
            for k in range(new_n):
                # 新グリッド点の物理座標
                x = new_origin[0] + i * new_spacing
                y = new_origin[1] + j * new_spacing
                z = new_origin[2] + k * new_spacing
                # 対応する元のグリッド内の位置（連続値）
                fx = (x - orig_origin[0]) / spacing_x
                fy = (y - orig_origin[1]) / spacing_y
                fz = (z - orig_origin[2]) / spacing_z
                # maptbxの補間関数を利用（trilinear補間）
                val = maptbx.interpolate_map_values(map_data, (fx, fy, fz))
                new_map[i,j,k] = val

    # 新たなマップヘッダー：立方体ユニットセル（L, L, L, 90,90,90）と新グリッド，新オリジンを設定
    new_unit_cell = (L, L, L, 90, 90, 90)
    new_header = header.customized_copy(unit_cell_parameters=new_unit_cell, origin=new_origin, grid=new_grid)
    
    # 新たな立方体マップを出力
    from iotbx import ccp4_map_writer
    ccp4_map_writer.write_ccp4_map_file(file_name=output_map, map_data=new_map, header=new_header)
    
    # 適用されたシフトは new_origin - orig_origin (物理座標での平行移動)
    shift_vector = (new_origin[0] - orig_origin[0],
                    new_origin[1] - orig_origin[1],
                    new_origin[2] - orig_origin[2])
    print("Applied shift vector (Å): {} {} {}".format(*shift_vector))

if __name__=="__main__":
    main()

#!/bin/bash
# 使用例:
#   ./auto_transform.sh input.mtz input.pdb refined_model.pdb
# ※ refined_model.pdb は、shifted_model.pdbを用いて精密化済みのモデルファイルとします。

if [ "$#" -ne 3 ]; then
    echo "Usage: $0   "
    exit 1
fi

MTZ_FILE=$1
PDB_FILE=$2
REFINED_MODEL=$3

echo "=== Step 1: MTZファイルから電子密度マップ (map.ccp4) を作成 ==="
phenix.mtz2map ${MTZ_FILE} d_min=2.0 keep_origin=True map_file=map.ccp4

echo "=== Step 2: マップを立方体 (cubic_map.ccp4) に再グリッド ==="
phenix.python auto_rebox.py map.ccp4 cubic_map.ccp4 > rebox_info.txt
# 出力例の行:
# Applied shift vector (Å): 5.0  -3.2  2.7
SHIFT_LINE=$(grep "Applied shift vector" rebox_info.txt)
Sx=$(echo $SHIFT_LINE | awk '{print $4}')
Sy=$(echo $SHIFT_LINE | awk '{print $5}')
Sz=$(echo $SHIFT_LINE | awk '{print $6}')
echo "適用されたシフトベクトル: ${Sx} ${Sy} ${Sz} (Å)"

echo "=== Step 3: PDBファイルの座標を変換（再配列演算子＋平行移動） ==="
# 自動で再配列演算子を算出
OP_INFO=$(phenix.python auto_operator.py ${MTZ_FILE} ${PDB_FILE})
OPERATOR=$(echo "$OP_INFO" | grep "Determined change-of-basis operator:" | awk '{print $5}')
INVERSE_OPERATOR=$(echo "$OP_INFO" | grep "Inverse operator:" | awk '{print $3}')

echo "算出された変換演算子: ${OPERATOR}"
echo "算出された逆変換演算子: ${INVERSE_OPERATOR}"

# ※ここでは、MTZ→map変換時にkeep_origin=Trueなら座標系は一致していると仮定し、
# さらにauto_rebox.pyで立方体化した際に適用された平行移動 (shift_vector) を反映します。
# そのため、PDB座標に対しては「change_of_basis」オプションで再配列演算子を（必要なら）適用し、
# さらに sites.translate で -shift_vector を適用します。
phenix.pdbtools ${PDB_FILE} change_of_basis="${OPERATOR}" sites.translate="-${Sx},-${Sy},-${Sz}" out=shifted_model.pdb

echo "shifted_model.pdb を作成しました。"
echo "=== Step 4: （ここで shifted_model.pdb を用いてモデル精密化を実施してください） ==="
echo "    ※精密化後のモデルは、引数で指定された ${REFINED_MODEL} とします。"

echo "=== Step 5: 精密化後のモデルを元の座標系に戻す ==="
# 逆変換：まず sites.translate で +shift_vector を適用し、その後逆の再配列演算子を適用
phenix.pdbtools ${REFINED_MODEL} sites.translate="${Sx},${Sy},${Sz}" change_of_basis="${INVERSE_OPERATOR}" out=final_model.pdb

echo "最終モデル (元のMTZ座標系に戻したもの) は final_model.pdb に出力されました。"