[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[ccp4bb]: Re:



Satish Nair wrote:
> 
> Hello,
> 
> We have high resolution data (2.0 ?A) and a good molecular replacement model
> (Rfree starting is 35%) which we'd like to throw into wARP for rebuilding/
> refinement.  The problem is that the crystal are of a protein-DNA complex
> and protin is not very friendly to DNA.  Is there any easy way to avoid
> having to generate a complete DNA dictionary for protin?
> 
> Thank you,
> Satish Nair
> (Burley lab)
Hello,

use the newest version of refmac: you don't need to use protin anymore !
download it from : 
http://www.ysbl.york.ac.uk/~garib/refmac/latest_refmac.html

attached is an example script for use of refmac5 with arp (it worked in
my case, no garanty for your !)

good luck
-- 


-------------------------------------------------------------------
Laurent Maveyraud               
Groupe de Cristallographie Biologique  -  IPBS-CNRS

205, route de Narbonne       	 Tel (+33) 05 61 17 54 35
31077 TOULOUSE CEDEX        	 Fax (+33) 05 61 17 54 48
France                       	 mailto:laurent.maveyraud@ipbs.fr
-------------------------------------------------------------------
#!/bin/ksh -f
#
# Refmac refinement with the possibility to run arp_warp for
# water molecules detection and refinement.
#
# ccp4 version 4 
# refmac verion 5
# arp_warp verion 5
#
# Laurent Maveyraud, 6 Jul. 2000
#############################################################

#########
# files #
#########

# starting coordinate file
pdb="ref1-modified.pdb"

# structure factor file (MTZ format, uniqueified)
mtz="../6bHOPr-125mM-4min.mtz"

# column labels for FP, SIGFP and FREE
fp="F"
sigfp="SIGF"
free="FreeR_flag"

# root name for output filenames
outfile="ref1"

# binaries
#myrefmac=refmac5
myarp=arp_warp

#########
# misc. #
#########

# maximum number of cycles to run 
totalcycles="10"

# do you want to stop if Rfree stops decreasing?
#   not implemented yet !
rfreestop="no"

######################
# solvent correction #
######################

# bulk solvent correction strategy < NONE | MASK | BABINET >
solvcorr="MASK"

##########
# refmac #
##########

# refmac binary
myrefmac="/usr/mnt/frames1/maveyrau/refmac.5.0.25/refmac5"

# how many refinement cycles within refmac ?
refmac_cycle="3"

# refinement method < CGMA | CGRA | CDIR >
refi_meth="CGMAT"

# refinement type < REST | UNRE | IDEA | RIGI | TLSR >
refi_type="REST"

# refinement residual function < LSQF | MLKF >
refi_resi="MLKF"

# B-factor refinement < OVER | ISOT | ANIS | MIXE >
refi_bref="ISOT"

# apply weighting to < MATR | GRAD >
weight_scheme="MATR"

# weight to get reasonable geometry:
weight_value="0.10"

# use of experiment sigma in weighting < EXPE | NOEX >:
weight_expe="EXPE"

# use of experiment sigma in scaling < EXPE | "blank" >:
scal_expe="EXPE"

# scaling using FREE set or WORK set ?
scal_set="FREE"

#######
# arp # 
#######

# should we run arp < yes | no >
doarp="yes"

# how many refmac refinement cycles before starting arp ?
delayarp="2"

# grid for arp
grid="128 160 200"

# asu for arp
xyzlim="0. 1. 0. 0.5 0. 1."

# no of atoms to remove:
natom_rem="60"

# no of atoms to add (total nb of atoms*0.08/(resolmax)^3):
natom_add="60"

# should arp work for water or for all atoms < WATERS | ALLATOMS >
arpmode="water"

# merge distance (0.6 for allatoms, 2.2 for waters)
merge="2.2"

# Chain-Id for water:
chain="W"

# remove all waters below $remove sigma in 2Fo-Fc
remove="1.0"
  
# find waters with density above $find in Fo-Fc (could be "AUTO")
find="AUTO"

#####################
# end of user input #
#####################

# setup CCP4 if not already done
if [[ "${CCP4}" = "" ]]; then
  . /usr/mnt/programs/xtal_4.0.1/ccp4/include/ccp4.setup
fi

# check that require files are present

if [ ! -f $pdb ]; then
 echo "ERROR: file not found = $pdb"
 exit
fi

if [ ! -f $mtz ]; then
  echo "ERROR: file not found = $mtz"
  exit
fi

# read some information in FOBS

mtzdump hklin ${mtz} <<EOF > mtzdump_$$.log
NREF 0
GO
EOF
  if [ $? -ne 0 ]; then
    echo "ERROR in mtzdump !!! Maybe ${mtz} is not an MTZ file !!!"
    rm -f ${CCP4_SCR}/$$*
    exit
  else
    spacegroup=`grep "Space group"  mtzdump_$$.log |  cut -d" " -f6 `
    nref=`grep "Number of Reflections"  mtzdump_$$.log |  cut -d" " -f8 `
    a=`awk '/Cell Dimensions/ {getline;getline;print $1}'  mtzdump_$$.log`
    b=`awk '/Cell Dimensions/ {getline;getline;print $2}'  mtzdump_$$.log`
    c=`awk '/Cell Dimensions/ {getline;getline;print $3}'  mtzdump_$$.log`
    alpha=`awk '/Cell Dimensions/  {getline;getline;print $4}'  mtzdump_$$.log`
    beta=`awk '/Cell Dimensions/  {getline;getline;print $5}'  mtzdump_$$.log`
    gamma=`awk '/Cell Dimensions/ {getline;getline;print $6}'  mtzdump_$$.log`
    low_res=`awk '/OVERALL FILE STATISTICS for resolution range /  {getline;getline;getline;getline;getline;getline;getline;getline;
getline;getline;print $9}'  mtzdump_$$.log`
    high_res=`awk '/OVERALL FILE STATISTICS for resolution range /  {getline;getline;getline;getline;getline;getline;getline;getline
;getline;getline;print $10}'  mtzdump_$$.log`

    rm -f  mtzdump_$$.log
fi

########################################################################
#
echo " "
echo " =============== REFMAC / ARP  REFINEMENT ==================="
echo " "
echo " user ............................... `whoami`"
echo " date ............................... `date`"
echo " directory .......................... `pwd`"
echo " "
echo " coordinates ........................ ${pdb}"
echo " Fobs ............................... ${mtz}"
echo " spacegroup ......................... ${spacegroup}"
echo " cell ............................... ${a} ${b} ${c} ${alpha} ${beta} ${gamma}"
echo " resolution ......................... ${low_res} ${high_res}"
echo " "
echo " # of refinement cycles ............. ${totalcycles}"
echo " refmac mode ........................ ${refi_type}"
echo " run of arp_warp .................... ${doarp}"
if [ ${doarp} = "yes" ]; then
  echo " arp mode ........................... ${arpmode}"
fi
echo " ============================================================"
echo " "
########################################################################

# set solvent correction parameters

if [ $solvcorr = "NONE" ]; then
  scal_type="SIMP"
  solv_meth="NO" 
  echo " no solvent correction requested "
fi

if [ $solvcorr = "BABINET" ]; then
  scal_type="BULK"
  solv_meth="NO" 
  echo " solvent correction according to Babinet's principle "
fi

if [ $solvcorr = "MASK" ]; then
  scal_type="SIMP"
  solv_meth="YES" 
  echo " solvent correction with constant value outside protein"
fi

echo " "

# set defaults and starting values
exit="false"
previous="0"
currentcycle="1"
nextcycle="2"

ln -s ${pdb} refmac_1.pdb
ln -s ${mtz} refmac_1.mtz

while [[ ${exit} = 'false' ]] ;do

if [[ ${doarp} = 'yes' && ${nextcycle} > ${delayarp} ]]; then 
# computes maps : 2fofcwt and fofcwt, and put them in the appropriate ASU

fft hklin refmac_${currentcycle}.mtz \
    mapout 2fofcwt_${currentcycle}.map \
<<EOF >fft_${currentcycle}.log
labi F1=FWT PHI=PHWT
titl 2FoFc map, cycle ${currentcycle}
grid ${grid}
reso ${low_res} ${high_res}
binm
end
EOF
  if [ $? -ne 0 ]; then
    echo "ERROR in fft (1), cycle ${currentcycle} !!!"
    exit
  fi

mapmask mapin 2fofcwt_${currentcycle}.map \
        mapout 2fofcwt_${currentcycle}.ext \
<<EOF >>fft_${currentcycle}.log
xyzlim ${xyzlim}
EOF
  if [ $? -ne 0 ]; then
    echo "ERROR in mapmask (1), cycle ${currentcycle} !!!"
    exit
  fi

rm 2fofcwt_${currentcycle}.map

fft hklin refmac_${currentcycle}.mtz \
    mapout fofcwt_${currentcycle}.map \
<<EOF >>fft_${currentcycle}.log
labi F1=DELFWT PHI=PHDELWT
titl FoFc map, cycle ${currentcycle}
grid ${grid}
reso ${low_res} ${high_res}
binm
end
EOF
  if [ $? -ne 0 ]; then
    echo "ERROR in fft (2), cycle ${currentcycle} !!!"
    exit
  fi

mapmask mapin fofcwt_${currentcycle}.map \
        mapout fofcwt_${currentcycle}.ext \
<<EOF >>fft_${currentcycle}.log
xyzlim ${xyzlim}
EOF
  if [ $? -ne 0 ]; then
    echo "ERROR in mapmask (2), cycle ${currentcycle} !!!"
    exit
  fi

rm fofcwt_${currentcycle}.map 
rm fft_${previous}.log

# run arp to locate water molecules

echo " arp_warp run, cycle #" ${currentcycle}

arp_warp xyzin refmac_${currentcycle}.pdb \
      mapin1 2fofcwt_${currentcycle}.ext \
      mapin2 fofcwt_${currentcycle}.ext \
      xyzout arp_${currentcycle}.pdb \
<<EOF >arp_${currentcycle}.log
mode update $arpmode
cell $a $b $c $alpha $beta $gamma
symm $spacegroup
remo atoms $natom_rem cutsigma $remove merge $merge
find atoms $natom_add chain $chain cutsigma $find
reso $low_res $high_res
refi waters
end
EOF
  if [ $? -ne 0 ]; then
    echo "ERROR in arp_warp, cycle ${currentcycle} !!!"
    exit
  fi

rm arp_${previous}.log
rm fofcwt_${currentcycle}.ext
rm 2fofcwt_${currentcycle}.ext

# find stats in arp_warp logfile and print them

 grep "removed during merging"  arp_${currentcycle}.log
 grep "will be removed:" arp_${currentcycle}.log
 grep "Number of new atoms:" arp_${currentcycle}.log
 grep "Density threshold:" arp_${currentcycle}.log
 echo " "

else
  mv refmac_${currentcycle}.pdb arp_${currentcycle}.pdb
fi
# run refmac to refine


echo " ============================================================"
echo "                 grand refinement cycle #" $currentcycle
echo " "

$myrefmac xyzin arp_${currentcycle}.pdb \
       xyzout refmac_${nextcycle}.pdb \
       hklin $mtz \
       hklout refmac_${nextcycle}.mtz \
<<EOF > refmac_${currentcycle}.log
LABI FP=${fp} SIGFP=${sigfp} FREE=${free}
LABO FWT=FWT PHWT=PHWT DELFWT=DELFWT PHDELWT=PHDELWT

! general keywords
NCYC ${refmac_cycle}
MONI medium
BINS 20

!refinement parameters
REFI TYPE ${refi_type}
REFI RESI ${refi_resi}
REFI BREF ${refi_bref}
REFI METH ${refi_meth}
REFI RESO $low_res $high_res

!scaling/bulksolvent correction
SCAL TYPE ${scal_type}
SCAL LSSC ANIS 
SCAL LSSC ${scal_expe} 
SCAL LSSC ${scal_set}
SOLV ${solv_meth} 

!weighting
WEIG ${weight_expe} ${weight_scheme} ${weight_value} 

EOF
  if [ $? -ne 0 ]; then
    echo "ERROR in refmac, cycle ${currentcycle} !!!"
    exit
  fi

echo " refmac refinement:"
# get some stats from refmac logfile and print them
 rworki=`grep "Overall R factor " refmac_${currentcycle}.log | head -1 | cut -c 44-49`
 rworkf=`grep "Overall R factor " refmac_${currentcycle}.log | tail -1 | cut -c 44-49`
 rfreei=`grep "Free R factor " refmac_${currentcycle}.log | head -1 | cut -c 44-49`
 rfreef=`grep "Free R factor " refmac_${currentcycle}.log | tail -1 | cut -c 44-49`
 ccworki=`grep "Overall correlation" refmac_${currentcycle}.log | head -1 | cut -c 44-49`
 ccworkf=`grep "Overall correlation" refmac_${currentcycle}.log | tail -1 | cut -c 44-49`
 ccfreei=`grep "Free correlation" refmac_${currentcycle}.log | head -1 | cut -c 44-49`
 ccfreef=`grep "Free correlation" refmac_${currentcycle}.log | tail -1 | cut -c 44-49`
 fomi=`grep "Overall figure " refmac_${currentcycle}.log | head -1 | cut -c 44-49`
 fomf=`grep "Overall figure " refmac_${currentcycle}.log | tail -1 | cut -c 44-49`
 rmsbondi=`grep "distances: ref" refmac_${currentcycle}.log | head -1 | cut -c57-61`
 rmsbondf=`grep "distances: ref" refmac_${currentcycle}.log | tail -1 | cut -c57-61`
 rmsanglei=`grep "angles  : ref" refmac_${currentcycle}.log | head -1 | cut -c57-61`
 rmsanglef=`grep "angles  : ref" refmac_${currentcycle}.log | tail -1 | cut -c57-61`
 rmschirali=`grep "Chiral centres:" refmac_${currentcycle}.log | head -1 | cut -c57-61`
 rmschiralf=`grep "Chiral centres:" refmac_${currentcycle}.log | tail -1 | cut -c57-61`
# echo "      | 0.7867 0.7867 | 0.7867 0.7867 | 0.7867 | 0.312 0.312 0.312 |" 
 echo "      |    Rfactor    |  Correlation  |        |  rms deviations   |"
 echo "      | work     free | work     free |  fom   | bonds angl. chir. |"
 echo "--------------------------------------------------------------------"
 echo "init. | $rworki $rfreei | $ccworki $ccfreei | $fomi | $rmsbondi $rmsanglei $rmschirali |"
 echo "final | $rworkf $rfreef | $ccworkf $ccfreef | $fomf | $rmsbondf $rmsanglef $rmschiralf |"
 echo "--------------------------------------------------------------------"
 echo " "
# should we stop cycling ?

  rm -f refmac_${previous}.log
  rm -f refmac_${previous}.pdb
  rm -f arp_${previous}.pdb
  rm -f refmac_${previous}.mtz


  if [[ ${currentcycle} = ${totalcycles} ]]; then
    echo " "
    echo " maximum number of cycle reached..."
    exit=true
  else
    export currentcycle=$((${currentcycle}+1))
    export nextcycle=$((${nextcycle}+1))
    export previous=$((${previous}+1))
  fi
done

mv refmac_${nextcycle}.mtz ${outfile}.mtz
mv refmac_${nextcycle}.pdb ${outfile}.pdb


fft hklin ${outfile}.mtz \
    mapout ${outfile}_2fofcwt.map \
<<EOF > fft_${nextcycle}.log
labi F1=FWT PHI=PHWT
titl 2FoFc map, cycle ${nextcycle}, R=$rworkf Rfree=$rfreef
end
EOF
  if [ $? -ne 0 ]; then
    echo "ERROR in fft (last) !!!"
    exit
  fi

mapmask mapin ${outfile}_2fofcwt.map \
        mapout ${outfile}_2fofcwt.ext \
        xyzin ${outfile}.pdb \
<<EOF >> fft_${nextcycle}.log
border 5
EOF
  if [ $? -ne 0 ]; then
    echo "ERROR in mapamsk (last) !!!"
    exit
  fi

rm ${outfile}_2fofcwt.map

fft hklin ${outfile}.mtz \
    mapout ${outfile}_fofcwt.map \
<<EOF >> fft_${nextcycle}.log
labi F1=DELFWT PHI=PHDELWT
titl FoFc map, cycle ${nextcycle}, R=$rworkf Rfree=$rfreef
end
EOF
  if [ $? -ne 0 ]; then
    echo "ERROR in fft (last) !!!"
    exit
  fi

mapmask mapin ${outfile}_fofcwt.map \
        mapout ${outfile}_fofcwt.ext \
        xyzin ${outfile}.pdb \
<<EOF >> fft_${nextcycle}.log
border 5        
EOF
  if [ $? -ne 0 ]; then
    echo "ERROR in mapamsk (last) !!!"
    exit
  fi

rm ${outfile}_fofcwt.map


exit