"""Performs comparisons between site-specific amino-acid preferences of influenza NP measured in two sequence contexts in biological replicates.
Written by Mike Doud, 2014.

This script is run with a single argument at the command line which specifies a configuration file that 
lists some paths to scripts and preference files which are used by this script, eg:

python compare_preferences.py compare_prefs_config.txt

The configuration file specifies paths to various directories containing scripts and preference files.
Each line of the configuration file should have a pathname and a path, separated by a space, for the following pathnames:

preferencefiles_dir: path to the preference files made by ``dms_tools``. likely ends in '/dmstools_output/'
calculate_rmsd_script_path: path to ``calculate_rmsd.py``, which is called by this script.
pdb_file_path: path to ``2IQH_monomerC.pdb`` for calculating distances between sites
output_directory: path to a directory for output from this script. Will be created if it doesn't already exist.

example config file:

preferencefiles_dir  /home/user/NP_analysis/dmstools_output/
calculate_rmsd_script_path /home/user/NP_analysis/scripts/calculate_rmsd.py
pdb_file_path /home/user/NP_analysis/2IQH_monomerC.pdb
output_directory /home/user/NP_analysis/compare_prefs_output/

Description of the analysis
---------------------------

Each inference of amino-acid preferences from an individual deep mutation scanning experiment contains noise, so differences in 
site-specific amino-acid preferences *between* two groups of replicate experiments is adjusted for the average differences *within*
each set of experimental replicates. 

For more information about how the RMSD_corrected value measures distance in amino-acid preferences between groups of replicate experiments, 
see the documentation for the `calculate_RMSD.py` script, which is used by this master script.

The options for running this analysis are all hardcoded at the beginning of the main function. 
These options are currently set to make comparisons using the amino-acid preference files that are generated by `run_dmstools.py`.

First the site-specific RMSD_within, RMSD_between, and RMSD_corrected are calculated for the comparision between two groups of experiments.

Then, null distributions of site-specific RMSD_corrected are generated by two methods:

    1) Simulated data is drawn from a Dirichlet distribution centered on the mean preferences of the two groups, with variance tuned using
    a scaling parameter that is empirically determined to give between-replicate correlations equal to the between-replicate correlations
    observed in the actual data. 

    2) Randomization is used to shuffle which individual experiment is assigned to which group of experiments, and site-specific RMSD_corrected values
    are calculated for all possible shuffles. If there are no true differences in the amino-acid preferences between
    the two groups of experiments, then the RMSD_corrected distribution should not change during shuffling.


The RMSD_corrected distances calculated for the two groups of replicate experiments can then be compared to the set of
RMSD_corrected distances generated under these two null models which assume that there are no true differences in amino-acid preferences, 
and p-values are assigned to sites based on either null distribution after applying a Bonferonni correction.

Figures 4 and 5B from the paper are generated by this script.
"""


import os
import sys
import time
import string
import rmsdtools
import mapmuts.io
import numpy as np
import mapmuts.sequtils
import dms_tools.file_io
import itertools
import sitegroups
import scipy
import operator


def RunScript(rundir, run_name, script_name, commands, use_sbatch, sbatch_cpus, walltime=None):
    """Runs a `python` script with a single argument being an infile composed of key/value pairs
    specified in the dictionary *commands* like so:

    python script_name infile

    *rundir* is the directory in which we run the job. Created if it does
    not exist.

    *run_name* is the name of the run, which should be a string without
    spaces. The input file has this prefix followed by ``_infile.txt``.

    *script_name* is the name of the script that we run. This should contain
    the full path to the script, which should be hardcoded once in the main function below.

    *commands* contains the commands written to the input file. Itis a list 
    of 2-tuples giving the key / value pairs.
    Both keys and values should be strings.

    *use_sbatch* is a Boolean switch specifying whether we use ``sbatch``
    to run the script. If *False*, the script is just run with the command
    line instruction. If *True*, then ``sbatch`` is used, and the command file
    has the prefix *run_name* followed by the suffix ``.sbatch``.

    *sbatch_cpus* is an option that is only meaningful if *use_sbatch* is 
    *True*. It gives the integer number of CPUs that are claimed via
    ``sbatch`` using the option ``sbatch -c``. 

    *waltime* is an option that is only meaningful if *use_sbatch* is
    *True*. If so, it should be an integer giving the number of hours 
    to allocate for the job. If *walltime* has its default value of 
    *None*, no wall time for the job is specified.

    Returns *runfailed*: *True* if run failed, and *False* otherwise.
    """
    print "Running %s for %s in directory %s..." % (script_name, run_name, rundir)
    currdir = os.getcwd()
    if not os.path.isdir(rundir):
        os.mkdir(rundir)
    os.chdir(rundir)
    if (not run_name) or not all([x not in string.whitespace for x in run_name]):
        raise ValueError("Invalid run_name of %s" % run_name)
    infile = '%s_infile.txt' % run_name
    open(infile, 'w').write('# input file for running script %s for %s\n%s' % (script_name, run_name, '\n'.join(['%s %s' % (key, value) for (key, value) in commands])))
    if use_sbatch:
        sbatchfile = '%s.sbatch' % run_name # sbatch command file
        jobidfile = 'sbatch_%s_jobid' % run_name # holds sbatch job id
        jobstatusfile = 'sbatch_%s_jobstatus' % run_name # holds sbatch job status
        joberrorsfile = 'sbatch_%s_errors' % run_name # holds sbatch job errors
        sbatch_f = open(sbatchfile, 'w')
        sbatch_f.write('#!/bin/sh\n#SBATCH\n')
        if walltime:
            sbatch_f.write('#PBS -l walltime=%d:00:00\n' % walltime)
        sbatch_f.write('python %s %s' % (script_name, infile))
        sbatch_f.close()
        os.system('sbatch -c %d -e %s %s > %s' % (sbatch_cpus, joberrorsfile, sbatchfile, jobidfile))
        time.sleep(60) # short 1 minute delay
        jobid = int(open(jobidfile).read().split()[-1])
        nslurmfails = 0
        while True:
            time.sleep(60) # delay 1 minute
            returncode = os.system('squeue -j %d > %s' % (jobid, jobstatusfile))
            if returncode != 0:
                nslurmfails += 1
                if nslurmfails > 180: # error is squeue fails at least 180 consecutive times
                    raise ValueError("squeue is continually failing, which means that slurm is not working on your system. Note that although this script has crashed, many of the jobs submitted via slurm may still be running. You'll want to monitor (squeue) or kill them (scancel) -- unfortunately you can't do that until slurm starts working again.")
                continue # we got an error while trying to run squeue
            nslurmfails = 0
            lines = open(jobstatusfile).readlines()
            if len(lines) < 2:
                break # no longer in slurm queue
        errors = open(joberrorsfile).read().strip()
    else:
        errors = os.system('python %s %s' % (script_name, infile))
    os.chdir(currdir)
    if errors:
        print "ERROR running %s for %s in directory %s." % (script_name, run_name, rundir)
        return True
    else:
        print "Successfully completed running %s for %s in directory %s." % (script_name, run_name, rundir)
        return False


def RunCommandLineOptionsScript(rundir, run_name, script_name, commands, use_sbatch, sbatch_cpus, walltime=None, use_full_partition=True):
    """Runs an installed script with specified command-line options.

    Written by Jesse Bloom and modified by Mike Doud to run scripts with command-line
    arguments (eg dms_tools scripts) instead of infiles (eg mapmuts scripts).

    *rundir* is the directory in which we run the job. Created if it does
    not exist.
    *run_name* is the name of the run, which should be a string without
    spaces.
    *script_name* is the name of the script that we run.
    *commands* is a list of commands for the script (command-line arguments).
    *use_sbatch* is a Boolean switch specifying whether we use ``sbatch``
    to run the script. If *False*, the script is just run with the command
    line instruction. If *True*, then ``sbatch`` is used, and the command file
    has the prefix *run_name* followed by the suffix ``.sbatch``.
    *sbatch_cpus* is an option that is only meaningful if *use_sbatch* is 
    *True*. It gives the integer number of CPUs that are claimed via
    ``sbatch`` using the option ``sbatch -c``. 
    *waltime* is an option that is only meaningful if *use_sbatch* is
    *True*. If so, it should be an integer giving the number of hours 
    to allocate for the job. If *walltime* has its default value of 
    *None*, no wall time for the job is specified.
    It is assumed that the script can be run at the command line using::
        script_name (command-line options)
    Returns *runfailed*: *True* if run failed, and *False* otherwise.
    """
    print "Running %s for %s in directory %s..." % (script_name, run_name, rundir)
    currdir = os.getcwd()
    if not os.path.isdir(rundir):
        os.mkdir(rundir)
    os.chdir(rundir)
    if (not run_name) or not all([x not in string.whitespace for x in run_name]):
        raise ValueError("Invalid run_name of %s" % run_name)
    if use_sbatch:
        sbatchfile = '%s.sbatch' % run_name # sbatch command file
        jobidfile = 'sbatch_%s_jobid' % run_name # holds sbatch job id
        jobstatusfile = 'sbatch_%s_jobstatus' % run_name # holds sbatch job status
        joberrorsfile = 'sbatch_%s_errors' % run_name # holds sbatch job errors
        sbatch_f = open(sbatchfile, 'w')
        sbatch_f.write('#!/bin/sh\n#SBATCH\n')
        if walltime:
            sbatch_f.write('#PBS -l walltime=%d:00:00\n' % walltime)
        sbatch_f.write('%s ' % (script_name))
        for command in commands:
            sbatch_f.write('%s ' % command)
        sbatch_f.close()

        if use_full_partition:
            os.system('sbatch -c %d -p full -e %s %s > %s' % (sbatch_cpus, joberrorsfile, sbatchfile, jobidfile))
        else:
            os.system('sbatch -c %d -e %s %s > %s' % (sbatch_cpus, joberrorsfile, sbatchfile, jobidfile))
        time.sleep(60) # short 1 minute delay
        jobid = int(open(jobidfile).read().split()[-1])
        nslurmfails = 0
        while True:
            time.sleep(60) # delay 1 minute
            returncode = os.system('squeue -j %d > %s' % (jobid, jobstatusfile))
            if returncode != 0:
                nslurmfails += 1
                if nslurmfails > 180: # error is squeue fails at least 180 consecutive times
                    raise ValueError("squeue is continually failing, which means that slurm is not working on your system. Note that although this script has crashed, many of the jobs submitted via slurm may still be running. You'll want to monitor (squeue) or kill them (scancel) -- unfortunately you can't do that until slurm starts working again.")
                continue # we got an error while trying to run squeue
            nslurmfails = 0
            lines = open(jobstatusfile).readlines()
            if len(lines) < 2:
                break # no longer in slurm queue
        errors = open(joberrorsfile).read().strip()
    else: 
        full_command = '%s ' % (script_name) + ' '.join(commands)
        errors = os.system(full_command)
    os.chdir(currdir)
    if errors:
        print "ERROR running %s for %s in directory %s." % (script_name, run_name, rundir)
        return True
    else:
        print "Successfully completed running %s for %s in directory %s." % (script_name, run_name, rundir)
        return False



def main(): 

    # Parse config file:
    configfilename = sys.argv[1]
    if not os.path.isfile(configfilename):
        raise IOError("Failed to find configuration file %s" % configfilename)
    d = mapmuts.io.ParseInfile(open(configfilename))
    preferencefiles_dir = mapmuts.io.ParseStringValue(d, 'preferencefiles_dir')
    if preferencefiles_dir[-1] != '/':
        preferencefiles_dir = preferencefiles_dir + '/'
    print "Using preferencefiles_dir %s" % preferencefiles_dir
    calculate_rmsd_script_path = mapmuts.io.ParseStringValue(d, 'calculate_rmsd_script_path')
    print "Using calculate_rmsd_script_path %s" % calculate_rmsd_script_path
    output_directory = mapmuts.io.ParseStringValue(d, 'output_directory')
    print "Using output_directory %s" % output_directory
    pdb_file_path = mapmuts.io.ParseStringValue(d, 'pdb_file_path')
    print "Using pdb_file_path %s" % pdb_file_path

    # seed the random number generator used to simulate preferences
    # by drawing from the Dirichlet distribution
    np.random.seed(1)

    # make output directory if doesn't exist.
    if not os.path.isdir(output_directory):
        os.mkdir(output_directory)
    # perform all analysis from this script within this directory:
    os.chdir(output_directory)

    # The script can use existing simulation and randomization data if 
    # it has already been generated in the appropriate output directory
    # by setting do_simulations and/or do_randomizations to False.
    do_simulations = False
    number_of_simulated_experiments = 1000
    do_randomizations = False

    logfile = "%s/compare_preferences_log.txt" % output_directory
    logout = open(logfile, 'w')

    distance_function = "Jensen-Shannon"     
    sites = range(2,499)
    aas = mapmuts.sequtils.AminoAcids(includestop=False)
    all_sites = sites
    variable_sites = sitegroups.VariableSites()
    rna_sites = sitegroups.RNABindingGrooveSites()
    conserved_contact_sites, all_contact_sites = sitegroups.SitesContactingVariableSites()
    site_groups = [all_sites, rna_sites, variable_sites, conserved_contact_sites]
    site_group_labels = ["All sites", "RNA-binding", "Variable", "Contacting variable"]

    # Specify the names of the preference files. 
    # All these files should be in the preferencefiles directory after running the `run_dmstools.py` script.
    PR8_preference_files = ['PR8_replicate_1_prefs.txt',
                            'PR8_replicate_2_prefs.txt',
                            'PR8_replicate_3_prefs.txt']
    Aichi68A_preference_files = [   'Aichi68A_replicate_1_prefs.txt',
                                    'Aichi68A_replicate_2_prefs.txt',
                                    'Aichi68A_replicate_3_prefs.txt',
                                    'Aichi68A_replicate_4_prefs.txt'
                                    ]
    Aichi68B_preference_files = [   'Aichi68B_replicate_1_prefs.txt',
                                    'Aichi68B_replicate_2_prefs.txt',
                                    'Aichi68B_replicate_3_prefs.txt',
                                    'Aichi68B_replicate_4_prefs.txt'
                                    ]
    Aichi68C_preference_files = [   'Aichi68C_replicate_1_prefs.txt',
                                    'Aichi68C_replicate_2_prefs.txt']
    WSN_HA_files = [ 'WSN_HA_rep1_prefs.txt', 'WSN_HA_rep2_prefs.txt', 'WSN_HA_rep3_prefs.txt' ]


    # Because the WSN-HA preference files don't contain a preference for site 1, we need to
    # make NP files that are stripped of site 1 preferences to allow taking the means of the two
    # while simulating data for the null model comparing NP to WSN-HA. We also need to strip sites beyond 498
    # from the WSN-HA files.
    PR8_preference_files_stripped = []
    for f in PR8_preference_files:
        these_prefs = rmsdtools.ReadDMSToolsFormattedPrefsToMMFormattedDict(preferencefiles_dir+f)
        rmsdtools.WriteEquilFreqs(these_prefs, aas, sites, "%s/%s_stripped"%(preferencefiles_dir,f))
        PR8_preference_files_stripped.append("%s_stripped"%f)

    Aichi68_preference_files_stripped = []
    for f in Aichi68A_preference_files + Aichi68B_preference_files + Aichi68C_preference_files:
        these_prefs = rmsdtools.ReadDMSToolsFormattedPrefsToMMFormattedDict(preferencefiles_dir+f)
        rmsdtools.WriteEquilFreqs(these_prefs, aas, sites, "%s/%s_stripped"%(preferencefiles_dir,f))
        Aichi68_preference_files_stripped.append("%s_stripped"%f)

    WSN_HA_preference_files_stripped = []
    for f in WSN_HA_files:
        these_prefs = rmsdtools.ReadDMSToolsFormattedPrefsToMMFormattedDict(preferencefiles_dir+f)
        rmsdtools.WriteEquilFreqs(these_prefs, aas, sites, "%s/%s_stripped"%(preferencefiles_dir,f))
        WSN_HA_preference_files_stripped.append("%s_stripped"%f)

    # Define the groups of preference files to do RMSD calculations on.
    # Specify tuples of (groupname, list of preference files for this group)
    PR8_group = ('PR/1934', PR8_preference_files)
    Aichi68_previousstudy_group = ('Aichi/1968 previous study', Aichi68A_preference_files + Aichi68B_preference_files)
    Aichi68_currentstudy_group = ('Aichi/1968 current study', Aichi68C_preference_files)
    Aichi68_AllMeasurements_group = ('Aichi/1968', Aichi68A_preference_files + Aichi68B_preference_files + Aichi68C_preference_files)
    WSN_HA_group = ('HA', WSN_HA_preference_files_stripped)
    PR8_group_stripped = ('PR/1934', PR8_preference_files_stripped)
    Aichi68_group_stripped = ('Aichi/1968', Aichi68_preference_files_stripped)

    # Now define the specific comparisons to make between groups. 
    # These comparisons are defined as a list of tuples in the form of (group1, group2):
    # The first is a comparison between previously published amino-acid preferences for Aichi68 and a new set of data from an additional set of two replicates.
    # We do not expect differences in the amino-acid preferences that have been inferred from these separate experiments on the same Aichi68 NP.
    # The second is a comparison between PR8 NP preferences and all Aichi68 NP preferences.
    # The third is a comparison between PR8 NP and the non-homologous WSN/1933(H1N1) HA.
    comparisons = [ (Aichi68_previousstudy_group, Aichi68_currentstudy_group),
                    (Aichi68_AllMeasurements_group, PR8_group),
                    (Aichi68_group_stripped, WSN_HA_group) ]

    logout.write("Running compare_preferences with the following options:\n")
    logout.write("Using distance function %s\n" % distance_function)
    logout.write("Using sites %s\n" % sites)
    logout.write("Running in output directory %s\n" % output_directory)
    logout.write("Performing comparisons between the following pairs of groups:\n\n")
    for i in range(len(comparisons)):
        logout.write("Comparison #%d:\n" % (i+1))
        logout.write("Comparing %s vs %s\n" % (comparisons[i][0][0], comparisons[i][1][0]) )
        logout.write("Preference files for %s: %s\n" % (comparisons[i][0][0], comparisons[i][0][1]) )
        logout.write("Preference files for %s: %s\n" % (comparisons[i][1][0], comparisons[i][1][1]) )
        logout.write("----------------\n\n")

    # keep RMSD calculations for each of the 3 comparisons performed, for later analysis:
    RMSD_dictionary_list = [] # entries in this list will be 2-tuples of (label, RMSD_dictionary)

    # Now perform each comparison between two groups of replicates in a new subdirectory:
    for comparison in comparisons:
        (group1, group2) = comparison
        group_1_name = group1[0]
        group_1_files = list(group1[1])
        group_2_name = group2[0]
        group_2_files = list(group2[1])
        label = "%s vs %s" % (group_1_name,group_2_name)
        pathlabel = label.replace("/", "").replace(" ", "_")
        logout.write("Beginning comparision: %s\n" % label)

        # Add the base directory path to all filenames:
        for i,f in enumerate(group_1_files):
            group_1_files[i] = preferencefiles_dir + '/' + f
        for i,f in enumerate(group_2_files):
            group_2_files[i] = preferencefiles_dir + '/' + f

        # make a subdirectory for this comparison and go to it.
        this_comparison_subdir = output_directory + ('/%s' % pathlabel)
        if not os.path.isdir(this_comparison_subdir):
            os.mkdir(this_comparison_subdir)
        os.chdir(this_comparison_subdir)

        # Specify the commands for `calculate_rmsd.py`:
        command_d = {   'number_of_groups' : '2',
                        'comparisons' : '1,2',
                        'outfileprefix' : 'RMSD',
                        'distance_function' : distance_function,
                        'sites_type' : 'range',
                        'sites' : '2 498',
                        'group_1_files' : ' '.join( group_1_files ),
                        'group_2_files' : ' '.join( group_2_files ),
                        'group_1_name' : group_1_name,
                        'group_2_name' : group_2_name
                        }

        # Run the `calculate_rmsd.py` script:
        RunScript(this_comparison_subdir, pathlabel + '_calculate_RMSD', calculate_rmsd_script_path, 
                    list(command_d.items()), False, 1)

        # read in the RMSD analysis that was saved by ``calculate_rmsd.py``:
        filename = this_comparison_subdir + '/RMSD_' + pathlabel + '_RMSD_calcs.txt'
        (rmsds_between, rmsds_within, rmsds_corrected) = rmsdtools.ReadRMSDFile(filename)

        # convert to a dictionary of rmsd dictionaries, keyed by rmsd type.
        this_comparison_data = {    'rmsds_between' : dict(rmsds_between),
                                    'rmsds_within' : dict(rmsds_within),
                                    'rmsds_corrected' : dict(rmsds_corrected),
                                }
        RMSD_dictionary_list.append((label,this_comparison_data))

        # Plot scatter of RMSD_between vs. RMSD_within for all sites:
        rmsdtools.PlotDistancesBetweenWithinWithSubgroup(rmsds_between, rmsds_within, [], label, [], "Scatter_nosubgroup", enforce_limits=True)
        
        # Plot scatter of RMSD_between vs. RMSD_within for all sites and highlight each sitegroup:
        for site_group, site_group_label in zip(site_groups, site_group_labels):
            rmsdtools.PlotDistancesBetweenWithinWithSubgroup(rmsds_between, rmsds_within, [site_group], label, [site_group_label], "Scatter_%s" % site_group_label, subgroup_colors = ['Blue'], enforce_limits=True)

        # Plot scatter of RMSD_between vs. RMSD_within for all sites and highlight RNA-binding and Variable sites, specifically:
        rmsdtools.PlotDistancesBetweenWithinWithSubgroup(rmsds_between, rmsds_within,
                                                        [rna_sites, variable_sites], label, ["RNA-binding", "Variable"],
                                                        "Scatter_RNA_Variable", subgroup_colors=['Magenta','Lime'], enforce_limits=True)

        
        if do_simulations:
            logout.write("Simulating data")
            # make a null distribution of RMSD scores by performing many simulated experiments:
            #
            # The null hypothesis here is that there are no true differences in the amino-acid preferences between the
            # two groups we are comparing and that the individual observations within each group are sampling 
            # amino-acid preferences from one true set of amino-acid preferences. Each group has some level of sampling variance
            # between replicates, and we quantify this noise by calculating the average Pearson's R
            # for all pairwise comparisons of amino-acid preferences within that group.
            #
            # To simulate replicate observations of amino-acid preferences sampled from one true set of amino-acid preferences
            # with a defined level of noise, we draw preference vectors from a Dirichlet distribution with mean equal to the mean
            # of the amino-acid preferences for the two groups we are comparing, and with variance adjusted for each simulated group
            # to give an average simulated Pearson's R equal to the average Pearson's R for the actual data of that group.
            #
            # This simulation_null distribution will be generated in a subdirectory:
            this_simulation_null_subdir = this_comparison_subdir + '/simulation_null'
            if not os.path.isdir(this_simulation_null_subdir):
                os.mkdir(this_simulation_null_subdir)
            os.chdir(this_simulation_null_subdir)

            # Average the preferences across both groups we are comparing. The mean preferences here are
            # assumed to be the true set of preferences for both groups under the null model, and we calculate these mean
            # preferences in order to draw simulated amino-acid preference observations from a Dirichlet distribution with
            # mean equal to this mean.
            # Find the average for each group:
            group_1_mean_outfile = pathlabel + '_group1meanprefs.txt'
            group_2_mean_outfile = pathlabel + '_group2meanprefs.txt'

            commands = ['--excludestop', group_1_mean_outfile, 'average'] + group_1_files
            RunCommandLineOptionsScript(this_simulation_null_subdir, 'dmsmerge_group1', 'dms_merge', commands, False, 1)

            commands = ['--excludestop', group_2_mean_outfile, 'average'] + group_2_files
            RunCommandLineOptionsScript(this_simulation_null_subdir, 'dmsmerge_group2', 'dms_merge', commands, False, 1)
            
            # Average evenly across the two groups:
            means_outfile = pathlabel + '_overallmeanprefs.txt'
            means_infiles = [group_1_mean_outfile, group_2_mean_outfile]
            commands = ['--excludestop', means_outfile, 'average'] + means_infiles
            RunCommandLineOptionsScript(this_simulation_null_subdir, 'dmsmerge_overallMean', 'dms_merge', commands, False, 1)

            # Read the mean preferences into memory for use in the simulation below:
            mean_prefs = rmsdtools.ReadDMSToolsFormattedPrefsToMMFormattedDict(means_outfile)

            # perform many simulated experiments drawing amino-acid preferences from the mean, with noise similar
            # to the noise in each respective group as determined by average pairwise Pearson's R. We will save the
            # simulated amino-acid preference files in sub-directories for each simulated experiment, and we will run
            # calculate_rmsd.py for each simulated experiment in that subdirectory too.
            logout.write("Performing %d simulated experiments\n" % number_of_simulated_experiments)
            # First determine the average Pearson's R for pairwise comparisons within each group.
            group_1_average_r = rmsdtools.AveragePairwisePearsonsR(group_1_files, sites)
            group_2_average_r = rmsdtools.AveragePairwisePearsonsR(group_2_files, sites)
            logout.write("Average pairwise Pearson's R in group %s = %f\n" % (group_1_name, group_1_average_r))
            logout.write("Average pairwise Pearson's R in group %s = %f\n" % (group_2_name, group_2_average_r))
            # define the range of scaling parameters to vary across in order to tune the sampling variance to match the actual data:
            scaling_parameters = np.arange(1,20,0.02)
            # define how close we want the simulated data's average pearson's R to match the actual data's:
            r_diff = 0.01
            logout.write("Searching for a scaling parameter for the Dirichlet that results in data with an average pairwise Pearson's R within %f of the actual data\n" % r_diff)
            group_1_scaling_parameter = 0
            group_2_scaling_parameter = 0

            # find the scaling parameter to match group_1_average_r:
            for scale in scaling_parameters:
                simulated_prefs = [] # a list of simulated observations of amino-acid preferences
                for i in range(len(group_1_files)): # simulate the same number of observations as in the actual data for this group
                    simulated_prefs.append( rmsdtools.DrawDirichletEquilFreqs(mean_prefs, aas, scale_conc = scale) )
                logout.write("Found r_diff of %f for scaling parameter of %f.\n" % ( abs( rmsdtools.AveragePairwisePearsonsR(simulated_prefs, sites, read_from_file = False) - group_1_average_r ), scale ))
                if abs( rmsdtools.AveragePairwisePearsonsR(simulated_prefs, sites, read_from_file = False) - group_1_average_r ) <= r_diff:
                    logout.write("Found scaling parameter of %f to match average R of %f for group %s\n" % (scale, group_1_average_r, group_1_name))
                    group_1_scaling_parameter = scale
                    break

             # find the scaling parameter to match group_2_average_r:
            for scale in scaling_parameters:
                simulated_prefs = [] # a list of simulated observations of amino-acid preferences
                for i in range(len(group_2_files)): # simulate the same number of observations as in the actual data for this group
                    simulated_prefs.append( rmsdtools.DrawDirichletEquilFreqs(mean_prefs, aas, scale_conc = scale) )
                logout.write("Found r_diff of %f for scaling parameter of %f.\n" % ( abs( rmsdtools.AveragePairwisePearsonsR(simulated_prefs, sites, read_from_file = False) - group_2_average_r ), scale ))
                if abs( rmsdtools.AveragePairwisePearsonsR(simulated_prefs, sites, read_from_file = False) - group_2_average_r ) <= r_diff:
                    logout.write("Found scaling parameter of %f to match average R of %f for group %s\n" % (scale, group_2_average_r, group_2_name))
                    group_2_scaling_parameter = scale
                    break

            assert group_1_scaling_parameter > 0
            assert group_2_scaling_parameter > 0

            # Now that scaling parameters are defined for each group, start simulating experiments:
            # the simulated data will be stored in "sims"
            sims = [] # each item in sims will be be a single simulation's analysis represented as a nested dictionary keyed first by the rmsd score type, and then by the site.
            # for instance, the RMSD_corrected score for site 45 in simulation 1 can be accessed by sims[0]['rmsds_corrected']['45']

            for simulation_i in range(1, number_of_simulated_experiments + 1):
                # for each simulated set of experiments, make a new subdirectory:
                this_simulation_i_subdir = this_simulation_null_subdir + '/sim_%d_data' % simulation_i
                if not os.path.isdir(this_simulation_i_subdir):
                    os.mkdir(this_simulation_i_subdir)
                os.chdir(this_simulation_i_subdir)

                # keep track of simulated preference files as they are generated, for each group:
                group_1_simulated_pref_files = []
                group_2_simulated_pref_files = []

                # simulate group 1 preference observations:
                for i in range(len(group_1_files)):
                    outfile = "sim_%s_scalingparameter_%.2f_replicate_%d_prefs.txt" % (group_1_name.replace("/", "").replace(" ", "_"), group_1_scaling_parameter, i+1)
                    rmsdtools.WriteEquilFreqs(rmsdtools.DrawDirichletEquilFreqs(mean_prefs, aas, scale_conc = group_1_scaling_parameter), aas, sites, outfile)
                    group_1_simulated_pref_files.append(outfile)

                # simulate group 2 preference observations:
                for i in range(len(group_2_files)):
                    outfile = "sim_%s_scalingparamter_%.2f_replicate_%d_prefs.txt" % (group_2_name.replace("/", "").replace(" ", "_"), group_2_scaling_parameter, i+1)
                    rmsdtools.WriteEquilFreqs(rmsdtools.DrawDirichletEquilFreqs(mean_prefs, aas, scale_conc = group_2_scaling_parameter), aas, sites, outfile)
                    group_2_simulated_pref_files.append(outfile)

                # run `calculate_rmsd.py` on this simulated experiment:
                command_d = {   'number_of_groups' : '2',
                                'comparisons' : '1,2',
                                'outfileprefix' : 'sim_%d_RMSD' % simulation_i,
                                'distance_function' : distance_function,
                                'sites_type' : 'range',
                                'sites' : '2 498',
                                'group_1_files' : ' '.join( group_1_simulated_pref_files ),
                                'group_2_files' : ' '.join( group_2_simulated_pref_files ),
                                'group_1_name' : 'sim_%s' % group_1_name.replace("/","").replace(" ","_"),
                                'group_2_name' : 'sim_%s' % group_2_name.replace("/","").replace(" ","_")
                                }
                RunScript(this_simulation_i_subdir, 'calculate_rmsd_simulation_%d' % simulation_i, calculate_rmsd_script_path, list(command_d.items()), False, 1)

                # Read in the RMSD analysis from the file that was just generated by calculate_rmsd.py:
                filename = 'sim_%d_RMSD_sim_%s_vs_sim_%s_RMSD_calcs.txt' % (simulation_i, group_1_name.replace("/","").replace(" ","_"), group_2_name.replace("/","").replace(" ","_"))
                (rmsds_between, rmsds_within, rmsds_corrected) = rmsdtools.ReadRMSDFile(filename)

                # Append all these calculations in `sims` so that we can go through them later
                this_sim_dict = {   'rmsds_between' : dict(rmsds_between),
                                    'rmsds_within' : dict(rmsds_within),
                                    'rmsds_corrected' : dict(rmsds_corrected),
                                    }
                
                sims.append(this_sim_dict)

            # done running simulations
    
        else:
            logout.write("Skipping data simulation and loading from previously-existing data")
            this_simulation_null_subdir = this_comparison_subdir + '/simulation_null'
            # load all the simulated data into `sims`
            sims = []
            for simulation_i in range(1, number_of_simulated_experiments+1):
                this_simulation_i_subdir = this_simulation_null_subdir + '/sim_%d_data' % simulation_i
                if not os.path.isdir(this_simulation_i_subdir):
                    assert ValueError("Attempting to load previous simulation data that does not exist!")
                os.chdir(this_simulation_i_subdir)
                filename = 'sim_%d_RMSD_sim_%s_vs_sim_%s_RMSD_calcs.txt' % (simulation_i, group_1_name.replace("/","").replace(" ","_"), group_2_name.replace("/","").replace(" ","_"))
                (rmsds_between, rmsds_within, rmsds_corrected) = rmsdtools.ReadRMSDFile(filename)
                # Append all these calculations in `sims` so that we can go through them later
                this_sim_dict = {   'rmsds_between' : dict(rmsds_between),
                                    'rmsds_within' : dict(rmsds_within),
                                    'rmsds_corrected' : dict(rmsds_corrected),
                                    }
                sims.append(this_sim_dict)

        # Now aggregate and analyze the simulated data in sims[]:
        # We have `number_of_simulated_experiments` independent simulations of the experiment under the null hypothesis that there is one true 
        # set of amino-acid preferences and that the differences between groups is due to group-specific noise.
        # Each of these simulations has an RMSD analysis associated with it.
        # Append all the RMSD_corrected scores from all the simulations together to have one null distribution of scores. 
        simulated_null_distributions = {} # keyed by [score_type] and contains ALL simulated scores for that rmsd score type.
        for score_type in ['rmsds_corrected']:
            simulated_null_distributions[score_type] = []
            for r in sites:
                r=str(r)
                for simulation in sims:
                    simulated_null_distributions[score_type].append( simulation[score_type][r] )

        # Now output a file in this_simulation_null_subdir that lists the probability of obtaining an
        # RMSD_corrected score under the null simulation model equal or greater than than the score that 
        # was observed in the real data for each site:
        os.chdir(this_simulation_null_subdir)
        outfile_name = '%s_score_probabilities.txt' % (pathlabel)
        fileout = open(outfile_name, 'w') 
        simulation_p_vals = {} # keyed by site
        fileout.write("SITE\tRMSD_probability\n")
        for site in sites:
            site = str(site)
            rmsd_corrected_score = this_comparison_data['rmsds_corrected'][site]
            rmsd_corrected_prob = sum(x >= rmsd_corrected_score for x in simulated_null_distributions['rmsds_corrected']) / float(len(simulated_null_distributions['rmsds_corrected']))
            fileout.write("%s\t%.9f\n" % (site, rmsd_corrected_prob))
            simulation_p_vals[site] = rmsd_corrected_prob
        fileout.close()

        # Plot the actual and simulated RMSD_corrected distributions.
        # Unpack the RMSD_corrected scores for all sites from the real data into lists
        real_rmsd_corrected_distribution = []
        for site in sites:
            site = str(site)
            real_rmsd_corrected_distribution.append( this_comparison_data['rmsds_corrected'][site] )

        rmsdtools.PlotHistogram(    [simulated_null_distributions['rmsds_corrected'], real_rmsd_corrected_distribution],
                                    ['Simulated data', 'Experimental data'], 
                                    "%s_simulated_RMSD_distribution.pdf" % pathlabel, 0.008, r"$RMSD_{corrected}$", "Frequency",
                                    rwidth = 0.8, legend_loc = 7, normalize = True, title = label.replace("_",""), show_legend = True, 
                                    hist_type='bar', xlim = [-.15,.7] )

        if do_randomizations:
            logout.write("Performing exact randomization tests")
            # Under the null hypothesis that there is no difference in the amino-acid preferences between
            # these two groups, we expect the distribution of RMSD_corrected scores to look the same
            # when we shuffle the replicate datasets between the two group labels and re-calculate RMSD_corrected.
            this_randomization_null_subdir = this_comparison_subdir + '/randomization_null'
            if not os.path.isdir(this_randomization_null_subdir):
                os.mkdir(this_randomization_null_subdir)
            os.chdir(this_randomization_null_subdir)

            # we will use itertools.combinations() to enumerate all possible randomizations
            # of the individual preference files from both groups into two randomized groups 
            # with sizes equal to the actual group sizes.
            size_group_1 = len(group_1_files)
            size_group_2 = len(group_2_files)
            allfiles = group_1_files + group_2_files
            # all possible sets of preference files for group_1:
            group_1_combinations = list(itertools.combinations(allfiles, size_group_1))
            group_2_combinations = [] # assemble the corresponding group_2_combinations for each group_1_combination:
            for g1combo in group_1_combinations:
                g2combo = []
                for f in allfiles:
                    if f not in g1combo:
                        g2combo.append(f)
                group_2_combinations.append(g2combo)
                assert len(g1combo)+len(g2combo) == len(allfiles)

            # Perform RMSD analysis on each randomized sample, and store the results in randomizations[]
            randomizations = [] # each item in randomizations will be be a single randomized experiment's analysis represented as a 
            # nested dictionary keyed first by the rmsd score type, and then by the site.
            # for instance, the rmsd_corrected score for site 45 in randomized experiment 1 can be accessed by randomizations[0]['rmsds_corrected']['45']
            randomization_i = 0
            for shuffled_group_1, shuffled_group_2 in zip(group_1_combinations,group_2_combinations):
                randomization_i += 1
                command_d = {   'number_of_groups' : '2',
                                'comparisons' : '1,2',
                                'outfileprefix' : 'randomization_%d_RMSD' % randomization_i,
                                'distance_function' : distance_function,
                                'sites_type' : 'range',
                                'sites' : '2 498',
                                'group_1_files' : ' '.join( shuffled_group_1 ),
                                'group_2_files' : ' '.join( shuffled_group_2 ),
                                'group_1_name' : 'shuffled_%s_group_1' % pathlabel,
                                'group_2_name' : 'shuffled_%s_group_2' % pathlabel
                                }
                RunScript(this_randomization_null_subdir, 'calculate_RMSD_randomization_%d' % randomization_i, calculate_rmsd_script_path, list(command_d.items()), False, 1)
                # Read in the RMSD analysis from the file that was just generated by `calculate_rmsd.py`:
                filename = 'randomization_%d_RMSD_shuffled_%s_group_1_vs_shuffled_%s_group_2_RMSD_calcs.txt' % (randomization_i, pathlabel, pathlabel)
                (rmsds_between, rmsds_within, rmsds_corrected) = rmsdtools.ReadRMSDFile(filename)
                # Append all these calculations in `randomizations` so that we can go through them later
                this_rand_dict = {  'rmsds_between' : dict(rmsds_between),
                                    'rmsds_within' : dict(rmsds_within),
                                    'rmsds_corrected' : dict(rmsds_corrected),
                                    }
                randomizations.append(this_rand_dict)

        else:
            logout.write("Skipping randomization and loading from previously-existing data")
            this_randomization_null_subdir = this_comparison_subdir + '/randomization_null'
            if not os.path.isdir(this_randomization_null_subdir):
                assert ValueError("Attempting to load previous randomization data that does not exist!")
            os.chdir(this_randomization_null_subdir)
            # determine how many randomizations should have been done:
            allfiles = group_1_files + group_2_files
            num_of_combinations = len( list(itertools.combinations(allfiles, len(group_1_files))) )
            randomizations = []
            for randomization_i in range(1,num_of_combinations+1):
                filename = 'randomization_%d_RMSD_shuffled_%s_group_1_vs_shuffled_%s_group_2_RMSD_calcs.txt' % (randomization_i, pathlabel, pathlabel)
                assert os.path.isfile(filename)
                (rmsds_between, rmsds_within, rmsds_corrected) = rmsdtools.ReadRMSDFile(filename)
                this_rand_dict = {  'rmsds_between' : dict(rmsds_between),
                                    'rmsds_within' : dict(rmsds_within),
                                    'rmsds_corrected' : dict(rmsds_corrected),
                                    }
                randomizations.append(this_rand_dict)

        # Append all the scores from the randomizations
        randomized_null_distributions = {}
        for score_type in ['rmsds_corrected']:
            randomized_null_distributions[score_type] = []
            for r in sites:
                r = str(r)
                for randomization in randomizations:
                    randomized_null_distributions[score_type].append( randomization[score_type][r] )

        # Calculate probabilities of observing scores from randomization equal or greater than the actual score.
        outfile_name = 'randomized_%s_score_probabilities.txt' % (pathlabel)
        fileout = open(outfile_name, 'w')
        randomization_p_vals = {} # keyed by site
        fileout.write("SITE\tRMSD_probability\n")
        for site in sites:
            site = str(site)
            rmsd_corrected_score = this_comparison_data['rmsds_corrected'][site]
            rmsd_corrected_prob = sum(x >= rmsd_corrected_score for x in randomized_null_distributions['rmsds_corrected']) / float(len(randomized_null_distributions['rmsds_corrected']))
            fileout.write("%s\t%.8f\n" % (site, rmsd_corrected_prob))
            randomization_p_vals[site] = rmsd_corrected_prob
        fileout.close()

        # Plot RMSD_corrected distribution along with the corresponding randomized distribution:
        rmsdtools.PlotHistogram(    [randomized_null_distributions['rmsds_corrected'], real_rmsd_corrected_distribution ],
                                    ['Randomized data', 'Experimental data'], "%s_randomization_RMSD_distribution.pdf" % (pathlabel), 0.008, r"$RMSD_{corrected}$", "Frequency",
                                    rwidth = 0.8, legend_loc = 7, normalize = True, title = label.replace("_",""), show_legend = True, hist_type='bar', xlim = [-.15,.7] )

        # move up directory
        os.chdir(this_comparison_subdir)



        # Using the Benjamini-Hochberg (BH) procedure to control the False Discovery Rate (FDR),
        # enumerate the significant sites using a FDR = 5% using each null distribution (simulation or randomization)
        rand_sig_sites = rmsdtools.BenjaminiHochbergCorrection(randomization_p_vals.items(), 0.05)
        sim_sig_sites = rmsdtools.BenjaminiHochbergCorrection(simulation_p_vals.items(), 0.05)

        sig_both = list( set(sim_sig_sites) & set(rand_sig_sites) ) # significant under both null models
        sig_both.sort()
        sig_either = list ( set(sim_sig_sites) | set(rand_sig_sites) ) # significant under either null model
        sig_either.sort()
        sig_only_sim = [s for s in sim_sig_sites if s not in rand_sig_sites] # significant only under simulation and not under randomization
        sig_only_rand = [s for s in rand_sig_sites if s not in sim_sig_sites] # significant only under randomization and not under simulation
        variable_sig = list( set(variable_sites) & set(sig_either) ) # which variable sites are significant 
        variable_sig.sort()
        variable_nonsig = [s for s in variable_sites if s not in sig_either] # which variable sites are non-significant
        sig_nonvariable = [s for s in sig_either if s not in variable_sites] # which significant sites are conserved

        # test clustering between significant sites found with either randomization or simulation nulls:
        cluster_test_rand_sig_sites = rmsdtools.TestSiteClustering(rand_sig_sites, pdbfile=pdb_file_path)
        cluster_test_sim_sig_sites = rmsdtools.TestSiteClustering(sim_sig_sites, pdbfile=pdb_file_path)

        # test distances to variable sites for significant sites found with either randomization or simulation nulls:
        proximity_to_variable_test_rand_sig_sites = rmsdtools.TestSiteProximityToVariableSites(rand_sig_sites, pdbfile=pdb_file_path)
        proximity_to_variable_test_sim_sig_sites = rmsdtools.TestSiteProximityToVariableSites(sim_sig_sites, pdbfile=pdb_file_path)

        # output summary of significant sites under both null distributions and results of statistical tests for
        # clustering of those sites in the NP crystal structure and distances of those sites which are evolutionarily conserved 
        # to evolutionarily variable sites
        outfile_name = 'summary_significantsites_%s.txt' % (pathlabel)
        fileout = open(outfile_name, 'w')
        fileout.write("Sites significant using both null distributions:\n")
        for s in sig_both:
            fileout.write("%d, " % int(s))
        fileout.write("\nSites significant using either null distribution:\n")
        for s in sig_either:
            fileout.write("%d, " % int(s))
        fileout.write("\nSites significant only under simulation null (and not under randomization):\n")
        for s in sig_only_sim:
            fileout.write("%d, " % int(s))
        fileout.write("\nSites significant only under randomization null (and not under simulation):\n")
        for s in sig_only_rand:
            fileout.write("%d, " % int(s))
        fileout.write("\n\n\nAll sites significant under simulation null:\n")
        for s in sim_sig_sites:
            fileout.write("%d, " % int(s))
        fileout.write("\n\n\nAll sites significant under randomization null:\n")
        for s in rand_sig_sites:
            fileout.write("%d, " % int(s))
        fileout.write("\n\n\nVariable sites with significant (using either null) RMSD_corrected:\n")
        for s in variable_sig:
            fileout.write("%d, " % int(s))
        fileout.write("\n\n\nVariable sites with non-significant RMSD_corrected:\n")
        for s in variable_nonsig:
            fileout.write("%d, " % int(s))
        fileout.write("\n\n\nSignificant (using either null) RMSD sites that are conserved between strains:\n")
        for s in sig_nonvariable:
            fileout.write("%d, " % int(s))
        fileout.write("\n\n")

        fileout.write("Tests for clustering of significant sites:\n")
        fileout.write("Of the top %d sites with largest preference shifts (null hypotheses rejected using exact randomization),\n" % len(rand_sig_sites))
        fileout.write("%d of these sites are resolved in the crystal structure and tested for spatial clustering here.\n" % cluster_test_rand_sig_sites[0])
        fileout.write("the median distance to nearest neighbor among these sites is %f,\n" % cluster_test_rand_sig_sites[1])
        fileout.write("the median distance to nearest neighbor among 1000 random draws of %d sites is %f,\n" % (cluster_test_rand_sig_sites[0], cluster_test_rand_sig_sites[2]))
        fileout.write("the one-sided Mann-Whitney U Test P-value comparing these distributions is %f\n\n" % cluster_test_rand_sig_sites[3])

        fileout.write("Of the top %d sites with largest preference shifts (null hypotheses rejected using simulation),\n" % len(sim_sig_sites))
        fileout.write("%d of these sites are resolved in the crystal structure and tested for spatial clustering here.\n" % cluster_test_sim_sig_sites[0])
        fileout.write("the median distance to nearest neighbor among these sites is %f,\n" % cluster_test_sim_sig_sites[1])
        fileout.write("the median distance to nearest neighbor among 1000 random draws of %d sites is %f,\n" % (cluster_test_sim_sig_sites[0], cluster_test_sim_sig_sites[2]))
        fileout.write("the one-sided Mann-Whitney U Test P-value comparing these distributions is %f\n\n\n" % cluster_test_sim_sig_sites[3])

        fileout.write("Test for distance of significant sites to variable sites:\n")
        fileout.write("Of the top %d sites with largest preference shifts (null hypotheses rejected using exact randomization),\n" % len(rand_sig_sites))
        fileout.write("%d of these sites are evolutionarily conserved between PR/1934 and Aichi/1968 and are resolved in the crystal structure and tested here.\n" % proximity_to_variable_test_rand_sig_sites[0])
        fileout.write("the median distance to nearest variable site among these sites is %f,\n" % proximity_to_variable_test_rand_sig_sites[1])
        fileout.write("the median distance to nearest variable site among 1000 random draws of %d sites is %f,\n" % (proximity_to_variable_test_rand_sig_sites[0],proximity_to_variable_test_rand_sig_sites[2]))
        fileout.write("the one-sided Mann-Whitney U Test P-value comparing these distribution is %f\n\n" % proximity_to_variable_test_rand_sig_sites[3])

        fileout.write("Of the top %d sites with largest preference shifts (null hypotheses rejected using simulation),\n" % len(sim_sig_sites))
        fileout.write("%d of these sites are evolutionarily conserved between PR/1934 and Aichi/1968 and are resolved in the crystal structure and tested here.\n" % proximity_to_variable_test_sim_sig_sites[0])
        fileout.write("the median distance to nearest variable site among these sites is %f,\n" % proximity_to_variable_test_sim_sig_sites[1])
        fileout.write("the median distance to nearest variable site among 1000 random draws of %d sites is %f,\n" % (proximity_to_variable_test_sim_sig_sites[0],proximity_to_variable_test_sim_sig_sites[2]))
        fileout.write("the one-sided Mann-Whitney U Test P-value comparing these distribution is %f\n\n" % proximity_to_variable_test_sim_sig_sites[3])

        fileout.close()


        # Make a csv file containing RMSD calculations and top preferences for each site in NP:
        if pathlabel == "Aichi1968_vs_PR1934":
            mean_PR8_prefs_dict = rmsdtools.ReadDMSToolsFormattedPrefsToMMFormattedDict(preferencefiles_dir+'/PR8_mean_prefs.txt')
            mean_Aichi68_prefs_dict = rmsdtools.ReadDMSToolsFormattedPrefsToMMFormattedDict(preferencefiles_dir+'/mean_Aichi68_both_studies_prefs.txt')
            threshold = 0.65 # minimal total fraction of the preferences to list in each table entry
            max_aa_per_entry = 3
            # A table to show all significant RMSD sites:
            outfile_name = 'RMSD_and_Preferences_Summary_%s.csv' % (pathlabel)
            fileout = open(outfile_name, 'w')

            # Header:
            fileout.write("Site, Amino Acid, RMSD_between, RMSD_within, RMSD_corrected, Significant by Randomization, Significant by Simulation, PR/1934 Preferences, Aichi/1968 Preferences\n")
            for s in sites:
                # site:
                fileout.write("%s, " % (s) )
                # AA identity/identities:
                if int(s) in variable_sites:
                    fileout.write("variable (%s / %s), " % (sitegroups.pr8_identity(s), sitegroups.aichi68_identity(s)))
                else:
                    fileout.write("conserved (%s), " % (sitegroups.pr8_identity(s)))
                #RMSD_between:
                fileout.write("%.3f, " % (this_comparison_data['rmsds_between'][str(s)]))
                #RMSD_within:
                fileout.write("%.3f, " % (this_comparison_data['rmsds_within'][str(s)]))
                #RMSD_corrected:
                fileout.write("%.3f, " % (this_comparison_data['rmsds_corrected'][str(s)]))

                if s in rand_sig_sites:
                    fileout.write("*, ")
                else:
                    fileout.write(", ")

                if s in sim_sig_sites:
                    fileout.write("*, ")
                else:
                    fileout.write(", ")

                # PR8 prefs:
                sorted_prefs = sorted([p for p in mean_PR8_prefs_dict[s].items() if 'PI_' in p[0]], key=operator.itemgetter(1), reverse=True)
                preference_sum = 0
                num_aa_this_entry = 0
                table_entry = ''
                for pref in sorted_prefs:
                    table_entry = table_entry + '%s (%.2f); ' % (pref[0][-1], pref[1])
                    preference_sum += pref[1]
                    num_aa_this_entry += 1
                    if preference_sum > threshold:
                        break
                    elif num_aa_this_entry == max_aa_per_entry:
                        break
                table_entry = table_entry[0:-2]
                fileout.write("%s, " % table_entry)
                # Aichi prefs:
                sorted_prefs = sorted([p for p in mean_Aichi68_prefs_dict[s].items() if 'PI_' in p[0]], key=operator.itemgetter(1), reverse=True)
                preference_sum = 0
                num_aa_this_entry = 0
                table_entry = ''
                for pref in sorted_prefs:
                    table_entry = table_entry + '%s (%.2f); ' % (pref[0][-1], pref[1])
                    preference_sum += pref[1]
                    num_aa_this_entry += 1
                    if preference_sum > threshold:
                        break
                    elif num_aa_this_entry == max_aa_per_entry:
                        break
                table_entry = table_entry[0:-2]
                fileout.write("%s" % table_entry)
                # End the row:
                fileout.write("\n")
            fileout.close()

        # End of this comparison loop

    # Make a grouped box plot for the distribution of RMSD_corrected scores at each sitegroup of interest, for the first two comparisions (Aichi previous vs. Aichi current, PR8 vs Aichi)
    os.chdir(output_directory)
    comp_1_bars = []
    comp_2_bars = []
    for site_group, site_group_label in zip(site_groups, site_group_labels):
        comp_1_scores = [] #eg, aichi-aichi
        comp_2_scores = [] #eg, aichi-pr8
        for r in site_group:
            r=str(r)
            comp_1_scores.append( RMSD_dictionary_list[0][1]['rmsds_corrected'][r] )
            comp_2_scores.append( RMSD_dictionary_list[1][1]['rmsds_corrected'][r] )
        comp_1_bars.append(comp_1_scores)
        comp_2_bars.append(comp_2_scores)
        (u, p) = scipy.stats.mannwhitneyu(comp_1_scores, comp_2_scores)
        logout.write("Mann-Whitney U-test between Aichi/Aichi and PR8/Aichi comparison RMSD for site group %s: p=%f\n" % (site_group_label, p))
    list_of_groups = [comp_1_bars, comp_2_bars]
    ylabel = r"$RMSD_{corrected}$"
    rmsdtools.PlotGroupedBoxPlot(list_of_groups, site_group_labels, ylabel, "RMSD_boxplot", title = False, y_lim = [-0.1,0.6])
    

main()
