#!/usr/bin/perl -w

# damaver_series_average.pl
# -------------------------------------------------------
# This script is a wrapper for DAMAVER to compare and average 
# a series of ab inito structure reconstruction models 
# obtained by the program DAMMIN.
# The idea is to compute pairwise Normalized Spatial Discrepancy (NSD)
# values for all pairs of models. Furthermore, the models
# will be averaged and a "filtered" consensus model is generated.
# This script is design to process the output of the script
# dammin_reconstruction_series.pl
#
# Author: Jan Lipfert
# Date: 11/1/2007
#
#
# 
### --- ACKNOWLEDGEMENT AND CITATION --- ###
# If you find this script useful, please cite it as
#
# Jan Lipfert, Rhiju Das, Vincent B. Chu, Madhuri Kudaravalli, Nathan Boyd, 
# Daniel Herschlag, and Sebastian Doniach   
# "Structural Transitions and Thermodynamics of a Glycine-Dependent Riboswitch from Vibrio cholerae" 
# J. Mol. Biol. 365:1393-1406 (2007)
#
# This script calls the programs DAMAVER and SUPCOMB:
# M. B.Kozin and D. I. Svergun
# "Automated mathcing of high- and low-resolution strutural models"
# J. Appl. Cryst. 34:33-41 (2001)
#
# V. V. Volkov and D. I. Svergun
# "Uniqueness of ab inito shape determination in small-angle scattering"
# J. Appl. Cryst.36:860-864 (2003)
#
###--- BEFORE THE SCRIPT CAN BE USED ---###
# The script requires that working copies of DAMAVER and SUPCOMB 
# are installed and executable on the system. 
# DAMAVER and SUPCOMB can be downloaded from
# http://www.embl-hamburg.de/ExternalInfo/Research/Sax/damaver.html
# http://www.embl-hamburg.de/ExternalInfo/Research/Sax/supcomb.html
# 
# To let the script know which version of DAMAVER to run 
# (and possibly the path to where dammin is installed), 
# set the $DAMAVER_exec variable (in the script after this header) 
# Examples: 
# $DAMAVER_exec = 'damaver.osx';
# (Example for use with the Mac-OS and with the executable in the path)
# $DAMAVER_exec = '/home/lipfert/software/damaver/damaver.l86';
# (Example for use on a Linux operating system, with specification of the full path)
#
# In addition, you need to make sure that the SUPCOMB executable is in the path
# (see DAMAVER documentation for details).
#
###--- USAGE OF THE SCRIPT AND OUTPUT ---###
# After the $DAMAVER_exec has been set,
# the script can be run from the command line with 2 arguments:
#
# 1) <project name> is the base name for the series of structures.
# The program assumes that there are files called <project name>?-1.pdb 
# (where ? is an index, see below) in the working directory.
# These are the models generated from repeat runs of DAMMIN by running the script
# dammin_reconstruction_series.pl
#
# 2) Number of models (a number). The program assumes that there are files with indices 
# from 0 to (Number of models -1) in the working directory.
# 
# 
# The script creates a subdirectory called 
# ave_<project name> and copies the model PDB files into this subdirectory,
# naming them foo?.pdb, where ? is the running index.
# It then runs the program DAMAVER to compute pairwise NSD values for all models in the
# subdirectory, and to generated averaged and filtered models.
# For more information about these "consensus" models see the DAMAVER documentation.
# http://www.embl-hamburg.de/ExternalInfo/Research/Sax/damaver.html
# 
# The consensus "averaged" model (which is the convex hull of all models)
# is moved to a file called ave_<project name>.pdb in the ave_<project name> subdirectory.
# The filtered model (which represents the best model after averaging)
# is moved to a file called filt_<project name>.pdb in the ave_<project name> subdirectory
#
# One caveat is that the filt_* and ave_* consensus models are only meaningful if
# the pairwise NSD values of the averaged models are less than one.






###---   Set this variable to the name and path of the DAMAVER executable   ---###
$DAMAVER_exec = "damaver.exe";


###---   START OF SCRIPT    ---###
if(!( defined($ARGV[0])) || !( defined($ARGV[1])))
{
    die "Usage: ./damaver_series_average.pl <project name> <Number of structures> \n (Creates a directory called ave_<project name>, and tries to copy files <project name>\$i-1.pdb)  \n";
}

$avedir = "ave_" . $ARGV[0];
$num_files = $ARGV[1];

### Check whether the ave_ directory exists already ###
if (-e $avedir)
{ 
    print "The directory $avedir exists already \n";
    while(1)
    {
	print "Do you want to overwrite it? (Y/N):";
	$_ = <STDIN>;         
	chomp;
	if($_ =~ /Y/ || $_ =~ /y/)
	{
	    print "Overwriting directory $avedir !\n";
	    system("rm -r $avedir");
	    last;
	}
	elsif($_ =~ /N/ || $_ =~ /n/)
	{
	    die "The directory $avedir exists and you don't want to overwrite! \n ... Exiting ... \n";
	}
	else
	{
	    print "You need to answer Y/N! \n";
	}
    }
}


### Make a new directory and copy the files over
print "Making the directory $avedir and copying files \n";
system("mkdir $avedir");

for($i=0; $i<$num_files; $i++)
{
    $fname = $ARGV[0] . $i . "-1.pdb";
    system("cp $fname $avedir/foo$i.pdb");
}


### Run DAMAVER
chdir($avedir);
system("$DAMAVER_exec /a");
print "Done running DAMAVER \n";

### Rename files, etc.
$avename = "ave_" . $ARGV[0] . ".pdb";
system("mv damaver.pdb $avename");

$filtname = "filt_" . $ARGV[0] . ".pdb";
system("mv damfilt.pdb $filtname");

print "Results of superposition are as follows \n";
system("more damsel.log");