#! /bin/tcsh -f
#
#	Compute the correlation coefficient between a bunch of pdb files
#	and a map file						-James Holton  9-8-11
#
#	usage:
#	local_corr.com rotamer1.pdb rotamer2.pdb observed.map
#
#	will compute the correlation coefficient of each pdb with observed.map
#	for all points in space occupied by any atom in all the pdbs
#

set mapfile   = observed.map
set lap_B     = 0

set logfile  = local_corr_details.log
set tempfile = local_corr_temp

set MAPMAN = /programs/o/rave/lx_mapman


if($#argv == 0) goto Help
goto Setup
Help:
cat << EOF
usage: 
local_corr.com observed.map rotamer1.pdb rotamer2.pdb ...

where: 
observed.map	- is the map file you want to compare to the pdbs
rotamer?.pdb	- are the pdb files you want to compare to the map

eg:
$0 2fofc.map rotamer1.pdb rotamer2.pdb
will print out the correlation coefficient of 2fofc.map with the calculated
map produced by rotamer1.pdb and then rotamer2.pdb at points in space occupied
by either pdb file.

EOF
exit 9
Return_from_Setup:

if($?debug) then
    echo "mapfile:  $mapfile"
    echo "pdbfiles: $pdbfiles"
endif

# quit if no input pdb is specified
if(($pdbfiles == 0)||(! -e "$mapfile")) goto Help
echo "" >! $logfile


# retrieve header from the input map
echo "GO" | mapdump MAPIN $mapfile >! ${tempfile}.mapdump

# determine space group and grid spacing from the map header
set SG     = `awk '/Space-group/{print $NF; exit}' ${tempfile}.mapdump`
set CELL   = `awk '/Cell dimensions/{print $4, $5, $6, $7, $8, $9; exit}' ${tempfile}.mapdump`
set GRID   = `awk '/Grid sampling/{print $(NF-2), $(NF-1), $NF; exit}' ${tempfile}.mapdump`
set AXIS   = `awk '/Fast, medium, slow/{print $(NF-2), $(NF-1), $NF}' ${tempfile}.mapdump | awk '! /[^ XYZ]/'`
set LIMITS = `awk '/Fast, medium, slow/{o[$(NF-2)]=1;o[$(NF-1)]=2;o[$NF]=3; print b[o["X"]],e[o["X"]], b[o["Y"]],e[o["Y"]], b[o["Z"]],e[o["Z"]]; exit} /Start and stop points/{b[1]=$(NF-5); e[1]=$(NF-4); b[2]=$(NF-3); e[2]=$(NF-2); b[3]=$(NF-1); e[3]=$NF}' ${tempfile}.mapdump`
if($?debug) echo "SG = $SG"
set sfallaxis = ""

#################################################################
# create a "probe" pdb for scoring the correlation space
cat ${tempfile}.mapdump |\
awk '/Cell dimensions/&& NF>6{printf "CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f\n",\
        $(NF-5),$(NF-4),$(NF-3),$(NF-2),$(NF-1),$NF}' |\
cat >! ${tempfile}probe.pdb

# pile all the atoms from all the structures into here
cat ${tempfile}pdbfiles |\
awk 'NF==1{file=$1; while(getline < file) print; close(file)}' |\
awk '/^ATOM/ || /^HETATM/{\
printf "ATOM      1  C   ALA     1    %s\n", substr($0, 31) }' |\
cat >> ${tempfile}probe.pdb
# this will be the "probe file" from now on

rm -f ${tempfile}.mapdump >& /dev/null


#################################################################
# calculate the labeled-space "probe map" from the probe pdb
sfall XYZIN ${tempfile}probe.pdb MAPOUT ${tempfile}flag.map << EOF | tee ${tempfile}sfall.log >> $logfile
MODE ATMMAP RESMOD
SYMM $SG
GRID $GRID
CELL $CELL
VDWR 3
EOF
rm -f ${tempfile}probe.pdb   >& /dev/null

if(! $?NO_RESIZE_MAPS) then
    # now get the size of the standard SFALL asymmetric unit
    set LIMITS = `awk '/Map dimensions are/{print $6, $8, $11, $13, $16, $18}' ${tempfile}sfall.log`
    set AXIS   = `awk '/Fast, medium, slow/{print $(NF-2), $(NF-1), $NF}' ${tempfile}sfall.log | awk '! /[^ XYZ]/'`

    echo "AXIS $AXIS" |\
    mapmask mapin $mapfile \
	    mapout ${tempfile}temp.map >> $logfile
    echo "XYZLIM $LIMITS" |\
    mapmask mapin ${tempfile}temp.map \
	    mapout ${tempfile}observed.map >> $logfile
    rm -f ${tempfile}temp.map
else
    ln -sf $mapfile ${tempfile}observed.map
endif
rm -f ${tempfile}sfall.log


# calculate the Laplacian of the observed map
if(-x $MAPMAN) then
    setenv MAPSIZE `ls -l ${tempfile}observed.map | awk '{printf "%d", $5/3.5}'`
    $MAPMAN << EOF >>& $logfile
read map1 ${tempfile}observed.map CCP4
filter laplace map1
write map1 ${tempfile}del_observed_mapman.map CCP4
Quit
EOF
endif




# calculate the Laplacian of the observed map
sfall mapin ${tempfile}observed.map hklout ${tempfile}sfall.mtz << EOF | tee ${tempfile}.log >> $logfile
MODE SFCALC MAPIN
EOF
if($status) then
    echo "fixing that..."
    # check axis order
    set temp = `cat ${tempfile}.log | awk '/Check Iuvw/{printf " %c %c %c", 87+$6, 87+$7, 87+$8}' | awk '! /[^ XYZ]/'`
    if($#temp == 3) set sfallaxis = "AXIS $temp"
    mapmask mapin ${tempfile}observed.map mapout ${tempfile}sfallme.map << EOF >> $logfile
    $sfallaxis
    xyzlim cell
EOF
    sfall mapin ${tempfile}sfallme.map hklout ${tempfile}sfall.mtz << EOF >> $logfile
MODE SFCALC MAPIN
EOF
endif


# apply the laplacian filter in reciprocal space
rm -f ${tempfile}lap.mtz
sftools << EOF >> $logfile
read ${tempfile}sfall.mtz
calc ( col Flap PHIlap ) = ( col FC PHIC ) ( STOLSQ 3.14159265358979 ) *
write ${tempfile}lap.mtz col Flap PHIlap
exit
y 
EOF
fft hklin ${tempfile}lap.mtz mapout ${tempfile}fft.map << EOF >> $logfile
LABIN F1=Flap PHI=PHIlap
scale F1 1 $lap_B
GRID $GRID
EOF
echo "AXIS $AXIS" |\
mapmask mapin ${tempfile}fft.map \
	mapout ${tempfile}temp.map >> $logfile
echo "XYZLIM $LIMITS" |\
mapmask mapin ${tempfile}temp.map \
	mapout ${tempfile}del_observed.map >> $logfile
rm -f ${tempfile}sfallme.map
rm -f ${tempfile}sfall.mtz
rm -f ${tempfile}lap.mtz
rm -f ${tempfile}fft.map
rm -f ${tempfile}temp.map



set pdb = 0
while ( $pdb < $pdbfiles )
    @ pdb = ( $pdb + 1 )
    
    # get the file name from the stored list
    set pdbfile = `awk "NR==$pdb" ${tempfile}pdbfiles`
    echo -n "$pdb $pdbfile "
    
    pdbset xyzin $pdbfile xyzout ${tempfile}calc.pdb << EOF >> $logfile
    CELL $CELL
EOF
    
    #####################################################################
    # calculate a full map from the input pdb (on same grid as $mapfile)
    sfall XYZIN ${tempfile}calc.pdb MAPOUT ${tempfile}calc.map << EOF >> $logfile
    MODE ATMMAP 
    SYMM $SG
    CELL $CELL
    GRID $GRID
EOF

    # correlate the indicated regions
    echo "correlate residue" |\
    overlapmap mapin1 ${tempfile}observed.map \
	       mapin2 ${tempfile}calc.map \
               mapin3 ${tempfile}flag.map mapout /dev/null |\
    tee -a $logfile | awk '/Total:/{printf "%s ", $2}'


  if(-x $MAPMAN) then
	setenv MAPSIZE `ls -l ${tempfile}observed.map | awk '{printf "%d", $5/3.5}'`
	$MAPMAN << EOF >>& $logfile
read map1 ${tempfile}calc.map CCP4
filter laplace map1
write map1 ${tempfile}del_calc.map CCP4
quit
EOF


    echo "correlate residue" |\
    overlapmap mapin1 ${tempfile}del_observed_mapman.map \
	       mapin2 ${tempfile}del_calc.map \
               mapin3 ${tempfile}flag.map mapout /dev/null |\
    tee -a $logfile | awk '/Total:/{printf "%s ", $2}'
  endif


    # calculate the Laplacian of the calculated map
    mapmask mapin ${tempfile}calc.map mapout ${tempfile}sfallme.map << EOF >> $logfile
    $sfallaxis
    xyzlim cell
EOF
    sfall mapin ${tempfile}sfallme.map hklout ${tempfile}sfall.mtz << EOF >> $logfile
    MODE SFCALC MAPIN
EOF
    # apply the laplacian filter in reciprocal space
    rm -f ${tempfile}lap.mtz
    sftools << EOF >> $logfile
    read ${tempfile}sfall.mtz
    calc ( col Flap PHIlap ) = ( col FC PHIC ) ( STOLSQ 3.14159265358979 ) *
    write ${tempfile}lap.mtz col Flap PHIlap
    exit
    y 
EOF
    fft hklin ${tempfile}lap.mtz mapout ${tempfile}fft.map << EOF >> $logfile
    LABIN F1=Flap PHI=PHIlap
    scale F1 1 $lap_B
    GRID $GRID
EOF
    echo "AXIS $AXIS" |\
    mapmask mapin ${tempfile}fft.map \
           mapout ${tempfile}temp.map >> $logfile
    echo "XYZLIM $LIMITS" |\
    mapmask mapin ${tempfile}temp.map \
           mapout ${tempfile}del_calc.map >> $logfile
    rm -f ${tempfile}sfallme.map
    rm -f ${tempfile}sfall.mtz
    rm -f ${tempfile}lap.mtz
    rm -f ${tempfile}fft.map
    rm -f ${tempfile}temp.map
    rm -f ${tempfile}.log


    echo "correlate residue" |\
    overlapmap mapin1 ${tempfile}del_observed.map \
	       mapin2 ${tempfile}del_calc.map \
               mapin3 ${tempfile}flag.map mapout /dev/null |\
    tee -a $logfile | awk '/Total:/{printf "%s ", $2}'

    echo ""
end

if($?debug) exit
# clean up 
rm -f ${tempfile}calc.pdb   >& /dev/null
rm -f ${tempfile}calc.map   >& /dev/null
rm -f ${tempfile}flag.map   >& /dev/null
rm -f ${tempfile}observed.map   >& /dev/null
rm -f ${tempfile}del_observed.map   >& /dev/null
rm -f ${tempfile}del_observed_mapman.map   >& /dev/null
rm -f ${tempfile}del_calc.map   >& /dev/null
rm -f ${tempfile}pdbfiles   >& /dev/null


exit

Setup:

set i = 0
set pdbfiles = 0
echo -n "" >! ${tempfile}pdbfiles
while ( $i < $#argv )
    @ i = ( $i + 1 )
    set arg = "$argv[$i]"

    if("$arg" =~ *.map) then
	set mapfile = "$argv[$i]"
	continue
    endif
    if("$arg" =~ B=*) then
	set lap_B = `echo $arg | awk -F "=" '{print $2}'`
	continue
    endif
    if(-d "$arg") then
	echo "looking for PDB files in $arg"
	ls -1 $arg |\
	awk -v dir=$arg '/.pdb$/{print dir"/"$1}' |\
	cat >> ${tempfile}pdbfiles
    endif
    if(("$arg" =~ *.pdb)||("$arg" =~ *.brk)) then
	if(-e "$arg") then
	    echo "$arg" >> ${tempfile}pdbfiles
	else
	    echo "WARNING: $arg does not exitst!"
	endif
	continue
    endif
end

set pdbfiles = `cat ${tempfile}pdbfiles | wc -l`

goto Return_from_Setup