text: #!/usr/bin/env python2
  
  # The MIT License (MIT)
  #
  # Copyright (c) 2015 Noah Ollikainen, Shane O'Connor
  #
  # Permission is hereby granted, free of charge, to any person obtaining a copy
  # of this software and associated documentation files (the "Software"), to deal
  # in the Software without restriction, including without limitation the rights
  # to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  # copies of the Software, and to permit persons to whom the Software is
  # furnished to do so, subject to the following conditions:
  #
  # The above copyright notice and this permission notice shall be included in
  # all copies or substantial portions of the Software.
  #
  # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
  # THE SOFTWARE.
  
  """\
  This script acts as a wrapper for the analysis scripts written by Noah Ollikainen. It accepts a target directory
  containing a benchmark run and: i) calculates the metrics; and ii) creates the R figures.
  
  input_directory specifies where the PDB files for analysis are contained and should meet the following criteria:
    * the results for each domain should be contained within separate folders e.g. PF00013_1WVN_A;
    * input_directory should contain a JSON file, benchmarks.json, which lists the different runs generated in the folder
      along with a Python regular expression used to select those files.
  The scripts used to run the Rosetta protocols on the Sun Grid Engine cluster system meet these criteria so you will not
  need to create the JSON files using that configuration. The JSON file is a mapping from benchmark descriptions to regular
  expressions used to identify the output files. An example benchmarks.json can be found in this folder (benchmarks.json.example).
  
  
  Usage:
      analyze.py [options] (-o OUTPUT_DIRECTORY) [--benchmark=BENCHMARK ...] [--include_published_data=METHOD ...]  (JSON_FILES ...)
  
  Options:
  
      --out OUTPUT_DIRECTORY -o OUTPUT_DIRECTORY
          The folder where the analysis files will be created.
  
      --prefix PREFIX
         If this option is supplied then output filenames will use this as a prefix.
  
      --benchmark BENCHMARK -b BENCHMARK
          Benchmark runs are defined in the benchmarks.json file. These values define a list of files which are then used as input
          for the analysis scripts. For example, in benchmarks.json.example, one valid argument to --run would be "kT = 0.9".
          If this option is omitted then analysis will be run for all benchmarks list in benchmarks.json.
  
      --overwrite
          If files generated during this analysis script (e.g. .fasta and .mi files) already exist then these are not
          regenerated by default. To force regeneration, use this flag.
  
      --expectn EXPECTED_CASES
          If a certain number of sequences is expected e.g. if 500 sequences were expected to be created then this flag
          can be used to ensure that no sequences are missing.
  
      --include_published_data METHOD
          Include the published data for the specified method. Valid values are 'All', 'Fixed', '0.3', '0.6', '0.9', '1.2', '1.8', '2.4'.
  
      --published_data_backrub_method_prefix PREFIX
          If published data is included, add this string as a prefix to identify the backrub methods. This is useful if the methods are
          named differently in your benchmark scripts.
  
      --published_data_backrub_method_suffix SUFFIX
          See the published_data_backrub_method_prefix option above.
  
      --scatterplot_colors COLORS_FILE
          Specifies a JSON file to be used to color the scatterplot series. The keys are names of benchmarks i.e. the top-level keys in the benchmarks.json files.
          If you specify a file that does not yet exist, one will be created for you to edit. The colors should be in a format R can understand e.g. hex triplets
          (#ff0000) specifying red, green, and blue values or predefined colors e.g. "orange", "blue".
  
  
  Authors:
      Noah Ollikainen: Original mathematical functions used in analysis
      Shane O'Connor: This script and modifications to the analysis functions
  """
  
  import collections
  import contextlib
  import math
  import numpy
  import shutil
  import time
  import shlex
  import subprocess
  import traceback
  import copy
  import os
  import re
  import glob
  import sys
  import json
  import pprint
  
  from libraries import docopt
  from libraries import colortext
  from utils.fsio import read_file, write_file, get_file_lines
  from utils.pdb import PDB
  from utils.get_mi import create_mi_file, SequenceMatrix
  from utils.rplot import run_r_script, color_wheel
  from covariation_similarity.covariation_similarity import compute_overlap
  from profile_similarity.profile_similarity import get_covarying_pairs, get_rosetta_sequence, get_natural_sequences, background
  from published_data.published_data import get_data_frame, get_median_metric_data
  from published_data.published_data import get_method_ids as get_published_method_ids
  from published_data.published_data import published_methods as published_backrub_methods
  
  canonical_aas = set(['A', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'K', 'L', 'M', 'N', 'P', 'Q', 'R', 'S', 'T', 'V', 'W', 'Y'])
  
  class Analyzer(object):
  
  
      def __init__(self, input_jsons, benchmark_ids, output_directory, overwrite_files, expectn, published_data_methods = [], backrub_method_prefix = '', backrub_method_suffix = '', analysis_file_prefix = '', scatterplot_colors_file = None):
  
          # Set up the object
          self.output_directory = os.path.abspath(output_directory)
          self.overwrite_files = overwrite_files
          self.expectn = expectn
          self.domain_sequences = None
          self.native_sequences = None
          self.valid_positions = {} # a mapping from Pfam domain to residue positions which never contain gaps or non-canonicals in the native sequences
          self.natural_domain_entropies = {} # the entropies of valid positions in the natural domains
          self.native_sequence_matrices = {} # a mapping from Pfam domain to SequenceMatrix objects
          self.entropy_levels_by_position = {} # a mapping from pairs (domain, position) to 'low', 'medium', or 'high' indicating the relative level of entropy at that position
          self.entropy_levels_ranges = {} # the definition of the ranges of entropy for the three classes (low, medium, and high)
          self.analysis_file_prefix = analysis_file_prefix
          if 'All' in published_data_methods:
              self.published_data_methods = ['Fixed', '0.3', '0.6', '0.9', '1.2', '1.8', '2.4']
          else:
              self.published_data_methods = published_data_methods
  
          self.backrub_method_prefix = backrub_method_prefix
          self.backrub_method_suffix = backrub_method_suffix
          self.series_colors = {}
  
          if not os.path.exists(self.output_directory):
              try:
                  os.makedirs(self.output_directory)
              except: pass
              if not os.path.exists(self.output_directory):
                  raise Exception('An exception occurred creating the output directory {0}.'.format(output_directory))
  
          # Read in the domain and native sequence data
          self.get_domain_sequences(read_file('domain_sequences.txt'))
          self.determine_positions_for_entropy_calculation()
          self.get_native_sequences(os.path.join('sequence_recovery', 'native_sequences'))
  
          # Set up which benchmarks we will be analyzing
          self.benchmark_details = Analyzer.load_benchmark_json_files(input_jsons, benchmark_ids)
          self.summary_stats = {}
          for k in self.benchmark_details.keys():
              self.summary_stats[k] = {}
          if not input_jsons:
              raise Exception('No JSON files were supplied to the analysis object.')
          if not self.benchmark_details:
              if benchmark_ids:
                  raise Exception('The benchmark ids were not found in the JSON files.')
              else:
                  raise Exception('No benchmark ids were found in the JSON files.')
  
          # Read in the custom colors
          if scatterplot_colors_file:
              self.read_colors(scatterplot_colors_file)
  
  
      @staticmethod
      def load_benchmark_json_files(input_jsons, benchmark_ids):
          d = {}
          benchmark_ids = sorted(set(benchmark_ids))
          for input_json in input_jsons:
              root_path = os.path.abspath(os.path.split(input_json)[0])
              if not os.path.exists(input_json):
                  raise Exception('Could not find JSON file "{0}".'.format(input_json))
              try:
                  j = json.loads(read_file(input_json))
              except:
                  raise Exception('{0} is not a valid JSON file.'.format(input_json))
              try:
                  if benchmark_ids:
                      b = {}
                      for id in benchmark_ids:
                          id = id.split('@')
                          if len(id) == 1:
                              benchmark_name = id[0]
                              if j.get(benchmark_name):
                                  for method, details in j[benchmark_name].iteritems():
                                      benchmark_key = '{0}@{1}'.format(benchmark_name, method)
                                      if d.get(benchmark_key):
                                          raise Exception('Ambiguity: Multiple records exists for {0} using {1}'.format(benchmark_name, method))
                                      else:
                                          d[benchmark_key] = copy.deepcopy(j[benchmark_name][method])
                                          d[benchmark_key]['benchmark'] = benchmark_name
                                          d[benchmark_key]['method_id'] = method
                                          d[benchmark_key]['root_path'] = root_path
                          elif len(id) == 2:
                              benchmark_name = id[0]
                              method = id[1]
                              if j.get(benchmark_name, {}).get(method, {}):
                                  benchmark_key = '{0}@{1}'.format(benchmark_name, method)
                                  if d.get(benchmark_key):
                                      raise Exception('Ambiguity: Multiple records exists for {0} using {1}'.format(benchmark_name, method))
                                  d[benchmark_key] = copy.deepcopy(j[benchmark_name][method])
                                  d[benchmark_key]['benchmark'] = benchmark_name
                                  d[benchmark_key]['method_id'] = method
                                  d[benchmark_key]['root_path'] = root_path
                          else:
                              raise Exception('The benchmark id {0} is invalid.')
                  else:
                      for benchmark_name, methods in j.iteritems():
                          for method, details in methods.iteritems():
                              benchmark_key = '{0}@{1}'.format(benchmark_name, method)
                              if d.get(benchmark_key):
                                  raise Exception('Ambiguity: Multiple records exists for {0} using {1}'.format(benchmark_name, method))
                              else:
                                  d[benchmark_key] = copy.deepcopy(details)
                                  d[benchmark_key]['benchmark'] = benchmark_name
                                  d[benchmark_key]['root_path'] = root_path
                                  d[benchmark_key]['method_id'] = method
              except Exception, e:
                  print(traceback.format_exc())
                  raise Exception('An exception occurred processing the JSON file: {0}'.format(str(e)))
          for benchmark_id, benchmark_details in sorted(d.iteritems()):
              assert('file_filter' in benchmark_details)
              assert('method' in benchmark_details)
          return d
  
  
      @staticmethod
      def get_normalized_run_file(benchmark_id, extension):
          return re.sub("([^\w\d\-_~,;:\[\]\(\).])", "_", benchmark_id) + extension
  
  
      # Setup scripts
  
  
      def read_colors(self, scatterplot_colors_file):
          '''Read the scatterplot color scheme or create a new one.'''
          if not os.path.exists(scatterplot_colors_file):
              colortext.error('The series colors file {0} does not exist so one is being created with the default colors. Edit the color scheme and rerun this command.'.format(scatterplot_colors_file))
              series_colors = {}
              all_benchmark_ids = sorted(set([v['benchmark'] for k, v in self.benchmark_details.iteritems()]))
              rgbcolors = color_wheel(len(all_benchmark_ids), start = 15, saturation_adjustment = None)
              for x in range(len(all_benchmark_ids)):
                  series_colors[all_benchmark_ids[x]] = '#' + rgbcolors[x]
              write_file(scatterplot_colors_file, json.dumps(series_colors, sort_keys=True, indent=4))
              sys.exit(1)
          else:
              try:
                  self.series_colors = json.loads(read_file(scatterplot_colors_file))
              except:
                  colortext.error('The file {0} could not be parsed correctly as a JSON file.'.format(scatterplot_colors_file))
                  sys.exit(1)
  
  
      def determine_positions_for_entropy_calculation(self):
          '''Computes the natural sequence entropy for the natural domains.'''
  
          domain_sequences = self.domain_sequences
          natural_domain_entropies = {}
          native_sequence_matrices = {}
  
          expected_positions = []
          for l in sorted(get_file_lines(os.path.join('published_data', 'entropy_table', 'valid_positions'))):
              if l.strip():
                  ts = l.strip().split('_')
                  assert(len(ts) == 2)
                  expected_positions.append((ts[0], int(ts[1])))
  
          entropy_levels_by_position = {}
          published_natural_entropies = {}
          _natural_entropies = []
          for l in sorted(get_file_lines(os.path.join('published_data', 'sequence_entropy_backrub.txt'))[1:]):
              if l.strip() and l.startswith('PF'):
                  ts = l.split('\t')
                  natural_entropy = float(ts[2])
                  position_ = int(ts[1])
                  published_natural_entropies[(ts[0], position_)] = natural_entropy
                  _natural_entropies.append((natural_entropy, ts[0], position_))
  
          # Partition the set of all positions in all domains into three parts of roughly equal size
          _natural_entropies = sorted(_natural_entropies)
          first_cut = len(_natural_entropies) / 3
          second_cut = (2 * len(_natural_entropies)) / 3
          while _natural_entropies[first_cut] == _natural_entropies[first_cut + 1]:
              first_cut += 1
          while _natural_entropies[second_cut] == _natural_entropies[second_cut + 1]:
              second_cut += 1
          for x in _natural_entropies[:first_cut]:
              entropy_levels_by_position[(x[1], x[2])] = 'low'
          for x in _natural_entropies[first_cut:second_cut]:
              entropy_levels_by_position[(x[1], x[2])] = 'medium'
          for x in _natural_entropies[second_cut:]:
              entropy_levels_by_position[(x[1], x[2])] = 'high'
          self.entropy_levels_by_position = entropy_levels_by_position
  
          # Remember the cut-offs used for the three classes to be used in plots later
          self.entropy_levels_ranges = dict(
              low = (_natural_entropies[0][0], _natural_entropies[first_cut - 1][0]),
              medium = (_natural_entropies[first_cut][0], _natural_entropies[second_cut - 1][0]),
              high = (_natural_entropies[second_cut][0], _natural_entropies[-1][0]),
          )
  
          colortext.wgreen('\nUnit-testing sequence entropy calculator: ')
          try:
              determined_positions = []
              for f in sorted(glob.glob(os.path.join('natural_alignments', '*.align.80'))):
                  sequences = []
                  num_conserved_positions = 0
                  num_varying_positions = 0
                  lines = [l.strip() for l in get_file_lines(f) if l.strip() and not(l.startswith('>'))]
                  domain_id = os.path.split(f[:f.find('.align.80')])[1]
                  pdb_chain_id = domain_sequences[domain_id]
                  indices_file = os.path.join('indices', '%s.indices' % pdb_chain_id)
                  assert(os.path.exists(indices_file))
                  indices_lines = [l.strip() for l in get_file_lines(indices_file) if l.strip()]
                  expected_length = int(indices_lines[-1].split()[0]) + 1
                  conserved_positions = [True] * expected_length
  
                  positions = range(expected_length)
                  for l in lines:
                      assert(len(l) == expected_length)
                      for x in positions:
                          if l[x] == '-' or l[x] not in canonical_aas:
                              conserved_positions[x] = False
                  for x in positions:
                      if conserved_positions[x]:
                          determined_positions.append((domain_id, x))
  
                  for l in lines:
                      sequence = ''
                      for x in positions:
                          if conserved_positions[x]:
                              sequence += l[x]
                          else:
                              sequence += '-'
                      sequences.append(sequence)
  
                  native_sequence_matrix = SequenceMatrix(sequences)
                  entropies = native_sequence_matrix.get_sequence_entropy()
                  native_sequence_matrices[domain_id] = native_sequence_matrix
                  for position_id, entropy in sorted(entropies.iteritems()):
                      if not numpy.isnan(entropy):
                          assert(abs(published_natural_entropies[(domain_id, position_id)] - entropy) < 0.0001)
                  natural_domain_entropies[domain_id] = entropies
              assert(sorted(determined_positions) == sorted(expected_positions))
          except Exception, e:
              colortext.error('Failed')
              print(str(e))
              print(traceback.format_exc())
              print('')
              sys.exit(1)
          colortext.message('Passed.')
  
          # Check against the expected positions
          valid_positions = {}
          valid_positions_list = [[(domain_id, position_id)for position_id, entropy in sorted(entropies.iteritems()) if not numpy.isnan(entropy)] for domain_id, entropies in sorted(natural_domain_entropies.iteritems())]
          valid_positions_list = [i for sublist in valid_positions_list for i in sublist]
          num_valid_positions = len(valid_positions_list)
          for tpl in valid_positions_list:
              valid_positions[tpl[0]] = valid_positions.get(tpl[0], set())
              valid_positions[tpl[0]].add(tpl[1])
          if not(num_valid_positions == 3009):
              colortext.error('Expected 3009 positions but only read in entropies for {0} positions.'.format(num_valid_positions))
          print('')
  
          self.natural_domain_entropies = natural_domain_entropies
          self.valid_positions = valid_positions
          self.native_sequence_matrices = native_sequence_matrices
  
  
      def get_domain_sequences(self, domain_sequences_file_content):
          domain_sequences = {}
          for line in domain_sequences_file_content.split('\n'):
              line = line.strip()
              if line:
                  domain_sequences[line.split()[0]] = line.split()[1]
          self.domain_sequences = domain_sequences
  
  
      def get_native_sequences(self, native_sequences_directory):
          native_sequences = {}
          for filename in glob.glob(os.path.join(native_sequences_directory, '*.fasta.txt')):
              lines = [l.strip() for l in get_file_lines(filename) if l.strip() and l[0] != ">"]
              assert(len(lines) == 1)
              native_sequences[os.path.split(filename)[1].split('.')[0]] = lines[0]
          self.native_sequences = native_sequences
  
  
      def create_fasta_file(self, benchmark_id, input_directory, file_filter):
          expectn = self.expectn
          print('Folder: {0}'.format(os.path.split(input_directory)[1]))
          normalized_run_id = Analyzer.get_normalized_run_file(benchmark_id, '.fasta')
          output_fasta_filepath = os.path.join(input_directory, normalized_run_id)
          if os.path.exists(output_fasta_filepath) and not self.overwrite_files:
              if expectn:
                  num_fasta_headers = len([1 for l in get_file_lines(output_fasta_filepath) if l.startswith('>')])
                  if expectn != num_fasta_headers:
                      raise colortext.Exception('Expected {0} records in {1} but read {2}.'.format(expectn, output_fasta_filepath, num_fasta_headers))
              print('FASTA file exists. Skipping generation.')
              return output_fasta_filepath
  
          fasta_records = []
          c = 0
          tc = 0
          input_files = [os.path.abspath(os.path.join(input_directory, f)) for f in sorted(os.listdir(input_directory)) if re.match(file_filter, f)] # since we cut off at expectn, sorting the files is important for reliability
          input_files = input_files[:expectn]
          progress_step = (len(input_files)/(20.0))
          if input_files:
              print('[' + ('=' * 6) + 'Progress' + ('=' * 6) + ']')
              sys.stdout.write(' ')
              for pdb_file in input_files:
                  p = PDB.from_filepath(pdb_file)
                  p.pdb_id = os.path.split(pdb_file)[1]
                  if not (p.atom_sequences):
                      raise Exception('Could not extract any ATOM sequence from the PDB file.')
                  fasta_records.append(p.create_fasta(75).strip())
                  c += 1
                  tc += 1
                  if c >= progress_step:
                      c -= progress_step
                      sys.stdout.write('.'); sys.stdout.flush()
              sys.stdout.write('\n'); sys.stdout.flush()
          if fasta_records:
              if expectn and (tc != expectn):
                  raise colortext.Exception('Expected {0} records in {1} but read {2}.'.format(expectn, input_directory, tc))
              write_file(output_fasta_filepath, '\n'.join(fasta_records))
              print('{0} FASTA records written.'.format(len(fasta_records)))
              return output_fasta_filepath
          else:
              colortext.error('No FASTA files were created for directory "{0}".'.format(input_directory))
              return False
  
  
      # Metric functions
  
  
      def compute_covariation_similarity(self, domain, mutual_information_filepath):
          '''Computes the covariation similarity for a domain based on the mutual information file.'''
          struct = self.domain_sequences[domain]
          natural_covariation_file = os.path.join('natural_covariation', domain + '_80.mi')
          indices_file = os.path.join('indices', struct + ".indices")
          if not os.path.exists(natural_covariation_file):
              raise Exception('Error: The file "{0}" does not exist.'.format(natural_covariation_file))
          if not os.path.exists(indices_file):
              raise Exception('Error: The file "{0}" does not exist.'.format(indices_file))
          overlap = compute_overlap(natural_covariation_file, mutual_information_filepath, indices_file)
          return overlap
  
  
      def compute_profile_similarity(self, domain, mutual_information_filepath, fasta_file, profile_similarities_by_position):
          '''Computes the profile similarity for a domain based on the mutual information and FASTA file.'''
  
          struct = self.domain_sequences[domain]
          natural_covariation_file = os.path.join('natural_covariation', domain + '_80.mi')
          indices_file = os.path.join('indices', struct + ".indices")
          natural_alignments_file = os.path.join('natural_alignments', domain + ".align.80")
          if not os.path.exists(natural_covariation_file):
              raise Exception('Error: The file "{0}" does not exist.'.format(natural_covariation_file))
          if not os.path.exists(indices_file):
              raise Exception('Error: The file "{0}" does not exist.'.format(indices_file))
          if not os.path.exists(natural_alignments_file):
              raise Exception('Error: The file "{0}" does not exist.'.format(natural_alignments_file))
  
          natural_pairs, positions = get_covarying_pairs(natural_covariation_file, mutual_information_filepath)
          designed_sequences, designed_indices = get_rosetta_sequence(fasta_file, indices_file)
          natural_sequences = get_natural_sequences(natural_alignments_file)
  
          similarity_score_mean = []
          sorted_positions = sorted(list(positions))
  
          P = {}
          Q = {}
  
          aa = 'ACDEFGHIKLMNPQRSTVWY'
          single = {}
          for a in aa:
              single[a] = 0.0
  
          for i in sorted_positions:
              P[i] = single.copy()
              Q[i] = single.copy()
  
          count_Q = 0.0
          for index in range(0, len(designed_sequences)):
              count_Q += 1.0
              for i in sorted_positions:
                  Q[i][designed_sequences[index][designed_indices[i]]] += 1.0
  
          count_P = 0.0
          for index in range(0, len(natural_sequences)):
              count_P += 1.0
              for i in sorted_positions:
                  P[i][natural_sequences[index][i]] += 1.0
  
          overall_sum = 0
          overall_count = 0.0
          for i in P:
              sum_PR = 0.0
              sum_QR = 0.0
              sum_bg_RR = 0.0
              sum_R_RR = 0.0
              for ii in P[i]:
                  P_i = float(P[i][ii]) / float(count_P)
                  Q_i = float(Q[i][ii]) / float(count_Q)
                  R_i = P_i * 0.5 + Q_i * 0.5
  
                  bg_i = background[ii]
                  RR_i = bg_i * 0.5 + R_i * 0.5
  
                  if R_i != 0 and Q_i != 0:
                      sum_PR += math.log( (Q_i / R_i), 2 ) * Q_i
  
                  if R_i != 0 and P_i != 0:
                      sum_QR += math.log( (P_i / R_i), 2) * P_i
  
                  if RR_i != 0 and bg_i != 0:
                      sum_bg_RR += math.log( (bg_i / RR_i), 2) * bg_i
  
                  if RR_i != 0 and R_i != 0:
                      sum_R_RR += math.log( (R_i / RR_i), 2) * R_i
  
              divergence = float(sum_PR + sum_QR) * 0.5
              significance = float(sum_bg_RR + sum_R_RR) * 0.5
              position_profile_similarity = ((1 - divergence) * (1 + significance)) * 0.5
              assert(profile_similarities_by_position.get((domain, i)) == None)
              profile_similarities_by_position[(domain, i)] = position_profile_similarity
              overall_sum += position_profile_similarity
              overall_count += 1.0
  
          return float(overall_sum) / float(overall_count)
  
  
      def compute_sequence_recovery(self, domain, fasta_file):
  
          struct = self.domain_sequences[domain]
          lines = get_file_lines(fasta_file)
  
          predicted_sequences = []
          seq = ''
          for l in lines:
              if l.strip():
                  if l[0] == '>':
                      if seq:
                          predicted_sequences.append(seq)
                          seq = ''
                  else:
                      seq += l.strip()
          predicted_sequences.append(seq)
          #pprint.pprint(predicted_sequences)
  
          #predicted_sequences = [l.strip() for l in lines if l.strip() and l[0] != '>']
          #print(len(predicted_sequences))
          assert(len(set(map(len, predicted_sequences))) == 1)
          if self.expectn != None:
              assert(len(predicted_sequences) == self.expectn)
  
          native_sequence = self.native_sequences[struct]
          # todo: These are hacks for the general Rosetta behavior for the Rosetta protocols in this capture where both
          #       the N-terminal residue of 1TEN_A and the C-terminal residue of 1UNQ are discarded due to missing ATOM
          #       records. This removal can be avoided using the '-ignore_zero_occupancy false' Rosetta option.
          #       If the length of the native sequence does not match the length of the predicted sequences, an error will
          #       be printed out below.
          if struct == '1TEN_A':
              native_sequence = native_sequence[1:]
          elif struct == '1UNQ_A':
              native_sequence = native_sequence[:-1]
  
          #print(native_sequence)
          if not len(native_sequence) == len(predicted_sequences[0]):
              colortext.error('The natural sequence length ({0}) does not match the predicted sequence lengths ({1}).'.format(len(native_sequence), len(predicted_sequences[0])))
              colortext.error('This may indicate an error in your method or could be due to truncations of the natural sequence in the analysis code. Please check the code near this error message in {0}.'.format(os.path.basename(__file__)))
              print('Native sequence    : {0}'.format(native_sequence))
              for x in range(min(10, len(predicted_sequences))):
                  print('Predicted sequence : {0}'.format(predicted_sequences[x]))
              native_sequence = native_sequence[:-1]
          sequence_length = min(len(native_sequence), len(predicted_sequences[0]))
          sum = 0
          for predicted_sequence in predicted_sequences:
              match = 0
              total = 0
              for i in range(0, sequence_length):
                  if predicted_sequence[i] == native_sequence[i]:
                      match += 1
                  total += 1
              sum += float(match) / float(total)
          #colortext.warning('Sequence recovery: {0}'.format(float(sum) / float(len(predicted_sequences))))
          return float(sum) / float(len(predicted_sequences))
  
  
      # Main function
  
  
      def read_published_profile_similarities_by_position(self):
          '''Returns a mapping from benchmark name -> position pairs (domain, position) -> profile similarity at that position
             for the data used in the original publication.'''
  
          # We only have the raw data for two methods at present
          profile_similarities_by_position = {}
          if self.published_data_methods:
              colortext.wgreen('Calculating profile similarity from the publication data: ')
              if '0.9' in self.published_data_methods:
                  profile_similarities_by_position['Published (score12)@backrub_0.9'] = {}
                  for mutual_information_filepath in sorted(glob.glob(os.path.join('..', 'output', 'backrub', '*_0.9.mi'))):
                      fasta_file = os.path.splitext(mutual_information_filepath)[0] + '.fasta'
                      assert(os.path.exists(fasta_file))
                      domain = os.path.split(mutual_information_filepath)[1].split('_')[0]
                      self.compute_profile_similarity(domain, mutual_information_filepath, fasta_file, profile_similarities_by_position['Published (score12)@backrub_0.9'])
              if 'Fixed' in self.published_data_methods:
                  profile_similarities_by_position['Published (score12)@fixbb'] = {}
                  for mutual_information_filepath in sorted(glob.glob(os.path.join('..', 'output', 'fixed_backbone', '*_0.0.mi'))):
                      fasta_file = os.path.splitext(mutual_information_filepath)[0] + '.fasta'
                      assert(os.path.exists(fasta_file))
                      domain = os.path.split(mutual_information_filepath)[1].split('_')[0]
                      self.compute_profile_similarity(domain, mutual_information_filepath, fasta_file, profile_similarities_by_position['Published (score12)@fixbb'])
              colortext.message('Done.')
          return profile_similarities_by_position
  
  
      def run(self):
          covariation_similarities = {}
          profile_similarities = {}
          sequence_recoveries = {}
          average_entropies = {}
          residue_entropies = {}
          profile_similarities_by_position = self.read_published_profile_similarities_by_position()
  
          for benchmark_id, benchmark_details in sorted(self.benchmark_details.iteritems()):
              covariation_similarities[benchmark_id] = {}
              profile_similarities[benchmark_id] = {}
              sequence_recoveries[benchmark_id] = {}
              average_entropies[benchmark_id] = {}
              residue_entropies[benchmark_id] = {}
              profile_similarities_by_position[benchmark_id] = {}
              input_directory = benchmark_details['root_path']
              assert(os.path.exists(input_directory))
              file_filter = benchmark_details['file_filter']
              # Create the input files (.mi and .fasta)
  
              colortext.message('\n== Creating input files for benchmark: {0} ==\n'.format(benchmark_id))
              colortext.message('Creating FASTA files and mutual information files')
              for f in sorted(os.listdir(input_directory)):
                  sub_dir = os.path.join(input_directory, f)
                  if os.path.isdir(sub_dir):
                      fasta_file = self.create_fasta_file(benchmark_id, sub_dir, file_filter)
                      if fasta_file:
                          domain = os.path.split(sub_dir)[1].split('_')[0]
                          indices_directory = os.path.abspath('indices')
                          mutual_information_filepath = os.path.join(sub_dir, Analyzer.get_normalized_run_file(benchmark_id, '.mi'))
                          residue_entropies_filepath = os.path.join(sub_dir, Analyzer.get_normalized_run_file(benchmark_id, '.residue.entropies'))
  
                          if self.overwrite_files or not(os.path.exists(mutual_information_filepath)) or not(os.path.exists(residue_entropies_filepath)):
                              mi_file, entropies, natural_indexed_residue_entropies = create_mi_file(domain, read_file(fasta_file), self.domain_sequences, indices_directory)
                              write_file(mutual_information_filepath, mi_file)
                              write_file(residue_entropies_filepath, json.dumps(natural_indexed_residue_entropies))
  
                          # Compute the benchmark metrics for each domain
                          # todo: this is where most of the time is taken in the loop - we could cache these values
                          covariation_similarities[benchmark_id][domain] = self.compute_covariation_similarity(domain, mutual_information_filepath)
                          profile_similarities[benchmark_id][domain] = self.compute_profile_similarity(domain, mutual_information_filepath, fasta_file, profile_similarities_by_position[benchmark_id])
                          sequence_recoveries[benchmark_id][domain] = self.compute_sequence_recovery(domain, fasta_file)
                          residue_entropies[benchmark_id][domain] = json.loads(read_file(residue_entropies_filepath))
                          average_entropies[benchmark_id][domain] = numpy.mean(residue_entropies[benchmark_id][domain].values())
                      else:
                          print('Mutual information (.mi) file exists. Skipping generation.' + input_directory)
                      print('')
  
          common_domains = set(covariation_similarities[covariation_similarities.keys()[0]].keys())
          for benchmark_id, covariation_similarity_values in covariation_similarities.iteritems():
              common_domains = common_domains.intersection(set(covariation_similarity_values.keys()))
          if not common_domains:
              raise Exception('There were no domains common to the specified benchmark runs.')
          common_domains = sorted(common_domains)
  
          colortext.message('Creating analysis plots in {0}.'.format(self.output_directory))
          self.cross_analyze_methods_and_benchmarks_by_position(profile_similarities_by_position, common_domains)
          self.cross_analyze_methods_and_benchmarks(covariation_similarities, profile_similarities, sequence_recoveries, common_domains)
          #self.analyze_residue_sequence_entropy(residue_entropies, common_domains)
          self.analyze_covariation_similarity(covariation_similarities, common_domains)
          self.analyze_profile_similarity(profile_similarities, common_domains)
          self.analyze_sequence_recovery(sequence_recoveries, common_domains)
          self.analyze_mean_sequence_entropy(average_entropies, common_domains)
  
          # File clean-up: move the R files and txt files into a subdirectory - it is likely that they will not be used
          self.cleanup_output_directory()
  
          # Compute the best and worst values for each metric. In all of the metrics, it happens that larger values indicate
          # an improvement.
          min_max = {}
          for k, v in self.summary_stats.iteritems():
              for metric, value in sorted(v.iteritems()):
                  if not min_max.get(metric):
                      min_max[metric] = (value, value)
                  else:
                      min_max[metric] = (min(min_max[metric][0], value), max(min_max[metric][1], value))
  
          # Print the summary (mean) value for each metric for all of the benchmarks
          # We color-code the values to highlight the best (green) and worst (red) values
          colortext.message('\nMetric values for the benchmarks.')
          for k, v in sorted(self.summary_stats.iteritems()):
              colortext.warning(k)
              metric_name_length = max(map(len, v.keys()))
              for metric, value in sorted(v.iteritems()):
                  value_str = str(value)
                  if len(self.summary_stats) != 1:
                      if min_max[metric][0] == value:
                          value_str = colortext.mred(value_str)
                      elif min_max[metric][1] == value:
                          value_str = colortext.mgreen(value_str)
                  print('  {0} : {1}'.format(metric.ljust(metric_name_length), value_str))
          print('')
  
  
      def cleanup_output_directory(self):
          '''Move the R files and txt files into a subdirectory of the output directory.'''
          rawdata_dir = os.path.join(self.output_directory, 'rawdata')
          try:
              os.mkdir(rawdata_dir)
          except: pass
          for f in glob.glob(os.path.join(self.output_directory, '*.R')):
              try:
                  print(f)
                  print(os.path.join(rawdata_dir, os.path.split(f)[1]))
                  shutil.move(f, os.path.join(rawdata_dir, os.path.split(f)[1]))
              except: pass
          for f in glob.glob(os.path.join(self.output_directory, '*.txt')):
              try:
                  print(f)
                  print(os.path.join(rawdata_dir, os.path.split(f)[1]))
                  shutil.move(f, os.path.join(rawdata_dir, os.path.split(f)[1]))
              except: pass
  
  
      def get_scatterplot_title_(self, benchmark_id):
          details = self.benchmark_details[benchmark_id]
          if details.get('kT'):
              method_details = 'kT={0}'.format(str(details.get('kT', '')))
          else:
              method_details = str(details.get('kT', ''))
          title = details['benchmark']
          return title, (details['method'], method_details)
  
  
      def cross_analyze_methods_and_benchmarks_by_position(self, profile_similarities_by_position, common_domains):
          '''Create scatterplots comparing the same metrics across different benchmark runs.'''
  
          benchmark_metrics = [
              # Metric name, dict name, name in published_data.py
              ['Profile similarity', profile_similarities_by_position, 'profile_similarity', 0.01],
          ]
          filename_prefixes, similarity_ranges = {}, {}
          for benchmark_metric in benchmark_metrics:
              filename_prefixes[benchmark_metric[0]] = benchmark_metric[2]
              similarity_ranges[benchmark_metric[0]] = benchmark_metric[3]
  
          # Read in the published and benchmark run data
          data_by_method = {}
          for benchmark_id, positional_data in profile_similarities_by_position.iteritems():
              if benchmark_id.startswith('Published (score12)@'):
                  title = benchmark_id.split('@')[0]
                  if benchmark_id.endswith('Fixed'):
                      method = ('Fixed', '')
                  else:
                      method = ('backrub', 'kT={0}'.format(benchmark_id.split('@')[1].replace('backrub_', '')))
              else:
                  title, method = self.get_scatterplot_title_(benchmark_id)
              data_by_method[method] = data_by_method.get(method, {})
              for benchmark_metric in benchmark_metrics:
                  data_by_method[method][benchmark_metric[0]] = data_by_method[method].get(benchmark_metric[0], {})
                  if title in data_by_method[method][benchmark_metric[0]]:
                      colortext.error('Ambiguity: The selected benchmarks do not have unique titles ("{0}"). The scatterplots for comparing benchmark runs over the same methods cannot be created.'.format(title))
                      return
                  data_by_method[method][benchmark_metric[0]][title] = benchmark_metric[1][benchmark_id]
  
          # For all comparable sets (same method, same metric), create a scatterplot of the values
          all_benchmark_ids = set()
          for method, method_data in data_by_method.iteritems():
              for metric_name, benchmark_values in method_data.iteritems():
                  all_benchmark_ids = all_benchmark_ids.union(set(benchmark_values.keys()))
          all_benchmark_ids.remove('Published (score12)')
          all_benchmark_ids = sorted(all_benchmark_ids)
  
          # Use a fixed color for the published data (this can be overridden with colors.json)
          series_colors = {
              'Published (score12)' : '#008B8B',
          }
          # Use colors spaced around the HSV wheel for the remaining benchmarks
          rgbcolors = color_wheel(len(all_benchmark_ids), start = 15, saturation_adjustment = None)
          for x in range(len(all_benchmark_ids)):
              series_colors[all_benchmark_ids[x]] = '#' + rgbcolors[x]
          # Override any colors with those specified in the command line, if any
          series_colors.update(self.series_colors)
  
          # Create a subdirectory to hold the individual plots
          if not os.path.exists(os.path.join(self.output_directory, 'positional_scatterplots')):
              os.mkdir(os.path.join(self.output_directory, 'positional_scatterplots'))
  
          # For all comparable sets (same method, same metric), create a scatterplot of the values
          for method, method_data in data_by_method.iteritems():
              # method e.g. = ('Fixed', '')
              for metric_name, benchmark_values in method_data.iteritems():
                  # metric_name e.g. = ('Covariation similarity', '') benchmark_values: benchmark_id -> domain -> value
                  similarity_range = similarity_ranges[metric_name]
                  benchmark_ids = sorted(benchmark_values.keys())
                  for x in range(len(benchmark_ids)):
                      for y in range(x + 1, len(benchmark_ids)):
                          benchmark_id_1 = benchmark_ids[x]
                          benchmark_id_2 = benchmark_ids[y]
  
                          x_color = series_colors[benchmark_id_1]
                          y_color = series_colors[benchmark_id_2]
  
                          # Get the list of positions common to both runs
                          xy_keys = sorted(set(benchmark_values[benchmark_id_1].keys()).intersection(set(benchmark_values[benchmark_id_2].keys())))
  
                          plot_title = '{0}, {1}'.format(metric_name, '@'.join([a for a in method if a]))
                          x_axis_label = benchmark_id_1
                          y_axis_label = benchmark_id_2
                          data_table_headers = ('Position', benchmark_id_1, benchmark_id_2, 'Entropy', 'Classification')
                          data_table = []
                          for domain_position in xy_keys:
                              entropy_classification = self.entropy_levels_by_position[domain_position]
  
                              xydiff = benchmark_values[benchmark_id_1][domain_position] - benchmark_values[benchmark_id_2][domain_position]
                              if -similarity_range <= xydiff <= similarity_range:
                                  similarity_classification = 'Similar'
                              elif xydiff > 0:
                                  similarity_classification = 'X'
                              elif xydiff < 0:
                                  similarity_classification = 'Y'
  
                              data_table.append(('{0}_{1}'.format(domain_position[0], domain_position[1]), benchmark_values[benchmark_id_1][domain_position], benchmark_values[benchmark_id_2][domain_position], entropy_classification, similarity_classification))
                          file_prefix = '_'.join([a.lower().replace(' ', '_').replace('=', '_') for a in [self.analysis_file_prefix, filename_prefixes[metric_name], method[0], method[1], x_axis_label, 'vs', y_axis_label] if a])
                          self.create_positional_scatterplot(os.path.join(self.output_directory, 'positional_scatterplots', file_prefix), data_table_headers, data_table, 1, 2, x_color, y_color, plot_title, x_axis_label, y_axis_label)
  
          # Combine the plots into a PDF file
          colortext.message('Creating a PDF containing all of the positional scatterplots.')
          try:
              if os.path.exists(os.path.join(self.output_directory, 'benchmark_vs_benchmark_positional_scatterplots.pdf')):
                  os.remove(os.path.join(self.output_directory, 'benchmark_vs_benchmark_positional_scatterplots.pdf'))
              p = subprocess.Popen(shlex.split('convert {0} benchmark_vs_benchmark_positional_scatterplots.pdf'.format(os.path.join('positional_scatterplots', '*.png'))), cwd = self.output_directory)
              stdoutdata, stderrdata = p.communicate()
              if p.returncode != 0: raise Exception('')
          except:
              colortext.error('An error occurred while combining the positional scatterplots using the convert application (ImageMagick).')
  
  
      def cross_analyze_methods_and_benchmarks(self, covariation_similarities, profile_similarities, sequence_recoveries, common_domains):
          '''Create scatterplots comparing the same metrics across different benchmark runs.'''
  
          benchmark_metrics = [
              # Metric name, dict name, name in published_data.py
              ['Covariation similarity', covariation_similarities, 'covariation_similarity', 0.01],
              ['Sequence recovery', sequence_recoveries, 'sequence_recovery', 0.01],
              ['Profile similarity', profile_similarities, 'profile_similarity', 0.01],
          ]
          filename_prefixes, similarity_ranges = {}, {}
          for benchmark_metric in benchmark_metrics:
              filename_prefixes[benchmark_metric[0]] = benchmark_metric[2]
              similarity_ranges[benchmark_metric[0]] = benchmark_metric[3]
  
  
          # Read in the published data
          #published_data = {}
          data_by_method = {}
          if self.published_data_methods:
              for published_data_method in self.published_data_methods:
                  for benchmark_metric in benchmark_metrics:
                      #published_data[benchmark_metric[0]] = published_data.get(benchmark_metric[0], {})
                      median_scores_by_domain = get_median_metric_data(published_data_method, benchmark_metric[2])
                      if published_data_method == 'Fixed':
                          method = ('Fixed', '')
                          #published_data[benchmark_metric[0]][('Fixed', '')] = median_scores_by_domain
                      else:
                          method = ('backrub', 'kT={0}'.format(str(published_data_method)))
                      data_by_method[method] = data_by_method.get(method, {})
                      assert(benchmark_metric[0] not in data_by_method[method])
                      data_by_method[method][benchmark_metric[0]] = {'Published (score12)' : median_scores_by_domain}
  
          # Read in the benchmark run data
          for benchmark_id, benchmark_details in self.benchmark_details.iteritems():
              title, method = self.get_scatterplot_title_(benchmark_id)
              data_by_method[method] = data_by_method.get(method, {})
              for benchmark_metric in benchmark_metrics:
                  data_by_method[method][benchmark_metric[0]] = data_by_method[method].get(benchmark_metric[0], {})
                  if title in data_by_method[method][benchmark_metric[0]]:
                      colortext.error('Ambiguity: The selected benchmarks do not have unique titles ("{0}"). The scatterplots for comparing benchmark runs over the same methods cannot be created.'.format(title))
                      return
                  data_by_method[method][benchmark_metric[0]][title] = benchmark_metric[1][benchmark_id]
  
          # For all comparable sets (same method, same metric), create a scatterplot of the values
          all_benchmark_ids = set()
          for method, method_data in data_by_method.iteritems():
              for metric_name, benchmark_values in method_data.iteritems():
                  all_benchmark_ids = all_benchmark_ids.union(set(benchmark_values.keys()))
          all_benchmark_ids.remove('Published (score12)')
          all_benchmark_ids = sorted(all_benchmark_ids)
  
          # Use a fixed color for the published data (this can be overridden with colors.json)
          series_colors = {
              'Published (score12)' : '#008B8B',
          }
          # Use colors spaced around the HSV wheel for the remaining benchmarks
          rgbcolors = color_wheel(len(all_benchmark_ids), start = 15, saturation_adjustment = None)
          for x in range(len(all_benchmark_ids)):
              series_colors[all_benchmark_ids[x]] = '#' + rgbcolors[x]
          # Override any colors with those specified in the command line, if any
          series_colors.update(self.series_colors)
  
          # Create a subdirectory to hold the individual plots
          if not os.path.exists(os.path.join(self.output_directory, 'scatterplots')):
              os.mkdir(os.path.join(self.output_directory, 'scatterplots'))
  
          # For all comparable sets (same method, same metric), create a scatterplot of the values
          for method, method_data in data_by_method.iteritems():
              # method e.g. = ('Fixed', '')
              for metric_name, benchmark_values in method_data.iteritems():
                  # metric_name e.g. = ('Covariation similarity', '') benchmark_values: benchmark_id -> domain -> value
                  similarity_range = similarity_ranges[metric_name]
                  benchmark_ids = sorted(benchmark_values.keys())
                  for x in range(len(benchmark_ids)):
                      for y in range(x + 1, len(benchmark_ids)):
                          benchmark_id_1 = benchmark_ids[x]
                          benchmark_id_2 = benchmark_ids[y]
                          x_color = series_colors[benchmark_id_1]
                          y_color = series_colors[benchmark_id_2]
  
                          # Get the list of domains common to both runs
                          xy_keys = sorted(set(benchmark_values[benchmark_id_1].keys()).intersection(set(benchmark_values[benchmark_id_2].keys())))
  
                          plot_title = '{0}, {1}'.format(metric_name, '@'.join([a for a in method if a]))
                          x_axis_label = benchmark_id_1
                          y_axis_label = benchmark_id_2
                          data_table_headers = ('Domain', benchmark_id_1, benchmark_id_2, 'Classification')
                          data_table = []
                          for domain in xy_keys:
                              xydiff = benchmark_values[benchmark_id_1][domain] - benchmark_values[benchmark_id_2][domain]
                              if -similarity_range <= xydiff <= similarity_range:
                                  classification = 'Similar'
                              elif xydiff > 0:
                                  classification = 'X'
                              elif xydiff < 0:
                                  classification = 'Y'
                              data_table.append((domain, benchmark_values[benchmark_id_1][domain], benchmark_values[benchmark_id_2][domain], classification))
                          file_prefix = '_'.join([a.lower().replace(' ', '_').replace('=', '_') for a in [self.analysis_file_prefix, filename_prefixes[metric_name], method[0], method[1], x_axis_label, 'vs', y_axis_label] if a])
                          self.create_scatterplot(os.path.join(self.output_directory, 'scatterplots', file_prefix), data_table_headers, data_table, 1, 2, x_color, y_color, plot_title, x_axis_label, y_axis_label)
  
          # Combine the plots into a PDF file
          colortext.message('Creating a PDF containing all of the scatterplots.')
          try:
              if os.path.exists(os.path.join(self.output_directory, 'benchmark_vs_benchmark_scatterplots.pdf')):
                  os.remove(os.path.join(self.output_directory, 'benchmark_vs_benchmark_scatterplots.pdf'))
              p = subprocess.Popen(shlex.split('convert {0} benchmark_vs_benchmark_scatterplots.pdf'.format(os.path.join('scatterplots', '*.png'))), cwd = self.output_directory)
              stdoutdata, stderrdata = p.communicate()
              if p.returncode != 0: raise Exception('')
          except:
              colortext.error('An error occurred while combining the scatterplots using the convert application (ImageMagick).')
  
  
      def create_positional_scatterplot(self, file_prefix, data_table_headers, data_table, x_series_index, y_series_index, x_color, y_color, plot_title = '', x_axis_label = '', y_axis_label = ''):
          # Create the R script
          boxplot_r_script = '''
  library(ggplot2)
  library(gridExtra)
  library(scales)
  library(qualV)
  
  # PNG generation
  png('%(file_prefix)s.png', width=2560, height=2048, bg="white", res=600)
  txtalpha <- 0.5
  redtxtalpha <- 0.8
  
  %(png_plot_commands)s
      '''
  
          xy_table_file = file_prefix + '.txt'
          write_file(file_prefix + '.txt', '\n'.join(','.join(map(str, line)) for line in [data_table_headers] + data_table))
  
          high_range = '{0:.2f} - {1:.2f}'.format(self.entropy_levels_ranges['high'][0], self.entropy_levels_ranges['high'][1])
          medium_range = '{0:.2f} - {1:.2f}'.format(self.entropy_levels_ranges['medium'][0], self.entropy_levels_ranges['medium'][1])
          low_range = '{0:.2f} - {1:.2f}'.format(self.entropy_levels_ranges['low'][0], self.entropy_levels_ranges['low'][1])
  
          single_plot_commands = '''
  # Set the margins
  par(mar=c(5, 5, 1, 1))
  
  xy_data <- read.csv('%(xy_table_file)s', header=T)
  
  names(xy_data)[%(x_series_index)d + 1] <- "xvalues"
  names(xy_data)[%(y_series_index)d + 1] <- "yvalues"
  
  # coefs contains two values: (Intercept) and yvalues
  coefs <- coef(lm(xvalues~yvalues, data = xy_data))
  fitcoefs = coef(lm(xvalues~0 + yvalues, data = xy_data))
  fitlmv_yvalues <- as.numeric(fitcoefs[1])
  lmv_intercept <- as.numeric(coefs[1])
  lmv_yvalues <- as.numeric(coefs[2])
  lm(xy_data$yvalues~xy_data$xvalues)
  
  xlabel <- "%(x_axis_label)s"
  ylabel <- "%(y_axis_label)s"
  plot_title <- "%(plot_title)s"
  rvalue <- cor(xy_data$yvalues, xy_data$xvalues)
  
                              #pprint.pprint(self.entropy_levels_ranges)
                              #{'high': (0.770027052352, 0.949615557346),
                              #'low': (0.0, 0.535423817387),
                              #'medium': (0.535768318026, 0.769974734857)}
  
  
  # Alphabetically, "Similar" < "X" < "Y" so the logic below works
  count_high <- "High\n(%(high_range)s)"
  count_medium <- "Medium\n(%(medium_range)s)"
  count_low <- "Low\n(%(low_range)s)"
  count_high <- "High"
  count_medium <- "Medium"
  count_low <- "Low"
  
  # Set graph limits and the position for the correlation value
  minx <- 0.0
  maxx <- min(1.0, max(xy_data$xvalues) + 0.1)
  miny <- 0.0
  maxy <- min(1.0, max(xy_data$yvalues) + 0.1)
  maxx <- max(maxx, maxy)
  maxy <- maxx
  xpos <- maxx / 10.0
  ypos <- maxy - (maxy / 10.0)
  
  plot_scale <- scale_color_manual(
      "Entropy",
      values = c( "high" = '#770000', "medium" = '#007777', "low" ='#0000aa'),
      labels = c(  "high" = count_high, "medium" = count_medium, "low" = count_low) )
  
  p <- qplot(main="", xvalues, yvalues, data=xy_data, xlab=xlabel, ylab=ylabel, size=I(1), shape = I(19), alpha = I(txtalpha), col=factor(Entropy)) +
          plot_scale +
          labs(title = "%(plot_title)s") +
          theme(plot.title = element_text(color = "#555555", size=rel(0.75))) +
          guides(col = guide_legend()) +
          geom_abline(size = 0.125, color="orange", intercept = lmv_intercept, slope = lmv_yvalues, alpha=0.2) +
          geom_abline(slope=1, intercept=0, linetype=3, size=0.25, alpha=0.4) + # add a diagonal (dotted)
          coord_cartesian(xlim = c(0.0, maxx), ylim = c(0.0, maxy)) + # set the graph limits
          geom_text(hjust=0, size=4, colour="black", aes(xpos, ypos, fontface="plain", family = "sans", label=sprintf("R = %%0.2f", round(rvalue, digits = 4)))) # add correlation text; hjust=0 sets left-alignment
  
  # Plot graph
  p
  
  dev.off()
          '''
  
          # Create the R script
          plot_type = 'png'
          png_plot_commands = single_plot_commands % locals()
          boxplot_r_script = boxplot_r_script % locals()
          r_script_filepath = os.path.join(self.output_directory, '{0}.R'.format(file_prefix))
          write_file(r_script_filepath, boxplot_r_script)
  
          # Run the R script
          run_r_script(r_script_filepath, cwd = os.path.join(self.output_directory, 'positional_scatterplots'))
  
  
      def create_scatterplot(self, file_prefix, data_table_headers, data_table, x_series_index, y_series_index, x_color, y_color, plot_title = '', x_axis_label = '', y_axis_label = ''):
  
          # Create the R script
          boxplot_r_script = '''
  library(ggplot2)
  library(gridExtra)
  library(scales)
  library(qualV)
  
  # PNG generation
  png('%(file_prefix)s.png', width=2560, height=2048, bg="white", res=600)
  txtalpha <- 0.8
  redtxtalpha <- 0.8
  
  %(png_plot_commands)s
      '''
  
          xy_table_file = file_prefix + '.txt'
          write_file(file_prefix + '.txt', '\n'.join(','.join(map(str, line)) for line in [data_table_headers] + data_table))
  
          single_plot_commands = '''
  # Set the margins
  par(mar=c(5, 5, 1, 1))
  
  xy_data <- read.csv('%(xy_table_file)s', header=T)
  
  names(xy_data)[%(x_series_index)d + 1] <- "xvalues"
  names(xy_data)[%(y_series_index)d + 1] <- "yvalues"
  
  # coefs contains two values: (Intercept) and yvalues
  coefs <- coef(lm(xvalues~yvalues, data = xy_data))
  fitcoefs = coef(lm(xvalues~0 + yvalues, data = xy_data))
  fitlmv_yvalues <- as.numeric(fitcoefs[1])
  lmv_intercept <- as.numeric(coefs[1])
  lmv_yvalues <- as.numeric(coefs[2])
  lm(xy_data$yvalues~xy_data$xvalues)
  
  xlabel <- "%(x_axis_label)s"
  ylabel <- "%(y_axis_label)s"
  plot_title <- "%(plot_title)s"
  rvalue <- cor(xy_data$yvalues, xy_data$xvalues)
  
  # Alphabetically, "Similar" < "X" < "Y" so the logic below works
  countsim <- paste("X ~ Y =", dim(subset(xy_data, Classification=="Similar"))[1])
  countX <- paste("X > Y =", dim(subset(xy_data, Classification=="X"))[1])
  countY <- paste("Y > X =", dim(subset(xy_data, Classification=="Y"))[1])
  
  # Set graph limits and the position for the correlation value
  minx <- 0.0
  maxx <- min(1.0, max(xy_data$xvalues) + 0.1)
  miny <- 0.0
  maxy <- min(1.0, max(xy_data$yvalues) + 0.1)
  maxx <- max(maxx, maxy)
  maxy <- maxx
  xpos <- maxx / 10.0
  ypos <- maxy - (maxy / 10.0)
  
  plot_scale <- scale_color_manual(
      "Counts",
      values = c( "Similar" = '#444444', "X" = '%(x_color)s', "Y" ='%(y_color)s'),
      labels = c(  "Similar" = countsim, "X" = countX, "Y" = countY) )
  
  p <- qplot(main="", xvalues, yvalues, data=xy_data, xlab=xlabel, ylab=ylabel, shape = I(18), alpha = I(txtalpha), col=factor(Classification)) +
          plot_scale +
          labs(title = "%(plot_title)s") +
          theme(plot.title = element_text(color = "#555555", size=rel(0.75))) +
          guides(col = guide_legend()) +
          geom_abline(size = 0.125, color="orange", intercept = lmv_intercept, slope = lmv_yvalues, alpha=0.2) +
          geom_abline(slope=1, intercept=0, linetype=3, size=0.25, alpha=0.4) + # add a diagonal (dotted)
          coord_cartesian(xlim = c(0.0, maxx), ylim = c(0.0, maxy)) + # set the graph limits
          geom_text(size=1.5, color="#000000", alpha=0.4, data=subset(xy_data, yvalues - xvalues > 0.15), aes(xvalues, yvalues+0.015, label=Domain)) + # label outliers
          geom_text(hjust=0, size=4, colour="black", aes(xpos, ypos, fontface="plain", family = "sans", label=sprintf("R = %%0.2f", round(rvalue, digits = 4)))) # add correlation text; hjust=0 sets left-alignment
  
  # Plot graph
  p
  
  dev.off()
          '''
  
          # Create the R script
          plot_type = 'png'
          png_plot_commands = single_plot_commands % locals()
          boxplot_r_script = boxplot_r_script % locals()
          r_script_filepath = os.path.join(self.output_directory, '{0}.R'.format(file_prefix))
          write_file(r_script_filepath, boxplot_r_script)
  
          # Run the R script
          run_r_script(r_script_filepath, cwd = os.path.join(self.output_directory, 'scatterplots'))
  
  
      def analyze_covariation_similarity(self, covariation_similarities, common_domains):
          '''This function creates a text file with the covariation similarities by benchmark and domain and runs R to
             generate boxplots of the covariation similarities.'''
          self.generic_analysis_function(covariation_similarities, common_domains, 'covariation_similarity', 'Covariation similarity')
  
  
      def analyze_profile_similarity(self, profile_similarities, common_domains):
          self.generic_analysis_function(profile_similarities, common_domains, 'profile_similarity', 'Profile similarity')
  
  
      def analyze_sequence_recovery(self, sequence_recoveries, common_domains):
          self.generic_analysis_function(sequence_recoveries, common_domains, 'sequence_recovery', 'Sequence recovery')
  
  
      def analyze_mean_sequence_entropy(self, average_entropies, common_domains):
          self.generic_analysis_function(average_entropies, common_domains, 'sequence_entropy', 'Sequence entropy')
  
  
      def analyze_residue_sequence_entropy(self, residue_entropies, common_domains):
          '''This function creates a text file with the metric data over the benchmarks and domains and runs R to
             generate boxplots of the metric data. The mean of the metrics over each benchmark is computed and added to
             the summary_stats dictionary.
          '''
          # todo: merge this function with generic_analysis_function
  
          metric_data = residue_entropies
          metric_column_name = 'sequence_entropy'
          metric_column_title = 'Sequence entropy'
  
          # Create the R data frame
          data_frame_values = []
          benchmark_ids = set()
          method_ids = set()
          domain_indices = set()
          for benchmark_id, domain_values in sorted(metric_data.iteritems()):
              benchmark_ids.add(benchmark_id)
              for domain, metric_values in sorted(domain_values.iteritems()):
                  for natural_sequence_index, metric_value in sorted(metric_values.iteritems()):
                      benchmark_details = self.benchmark_details[benchmark_id]
                      if benchmark_details.get('title'):
                          # If the user specifies a title, that should be used
                          method = benchmark_details['title']
                      else:
                          method = benchmark_details['method']
                          if benchmark_details.get('kT') != None:
                              method += ', kT {0}'.format(benchmark_details.get('kT'))
                      method_ids.add(method)
                      domain_indices.add(domain + '_' + natural_sequence_index)
                      data_frame_values.append((benchmark_details['benchmark'], method, domain, int(natural_sequence_index), metric_value))
  
          benchmark_column = "','".join(d[0] for d in data_frame_values)
          method_column = "','".join(d[1] for d in data_frame_values)
          domain_column = "','".join(d[2] for d in data_frame_values)
          natural_sequence_index_column = "','".join(d[3] for d in data_frame_values)
          metric_data_column = ",".join(str(d[4]) for d in data_frame_values)
  
          # Create the R script
          boxplot_generation_commands = '''
  p <- ggplot(data = benchmark_data, aes(factor(method), {0}, fill=method))
  p <- p + stat_boxplot(geom ='errorbar', linetype="solid", width = 0.25, position = "dodge")
  p <- p + stat_boxplot(geom = "boxplot", linetype="solid", position = "dodge", width = 0.60, na.rm = TRUE) + xlab("Method") + ylab("{1}") + facet_grid(.~benchmark) + theme(legend.position="none")
  p
  '''.format(metric_column_name, metric_column_title)
  
          analysis_file_prefix = self.analysis_file_prefix
          boxplot_r_script = '''
  library(ggplot2)
  library(grid)
  
  '''
          published_data_dataframe = None
          if self.published_data_methods:
              published_data_dataframe = get_data_frame(self.published_data_methods, metric_column_name, backrub_method_prefix = self.backrub_method_prefix, backrub_method_suffix = self.backrub_method_suffix)
              if published_data_dataframe:
                  boxplot_r_script += published_data_dataframe
              method_ids = method_ids.union(set(get_published_method_ids(self.published_data_methods, backrub_method_prefix = self.backrub_method_prefix, backrub_method_suffix = self.backrub_method_suffix)))
              benchmark_ids.add('_plos_one_published_data')
  
  
          boxplot_r_script += '''
  benchmark_data <- data.frame(benchmark=c('%(benchmark_column)s'),
                               method=c('%(method_column)s'),
                               domain=c('%(domain_column)s'),
                               %(metric_column_name)s=c(%(metric_data_column)s))
  ''' % locals()
          if published_data_dataframe:
              boxplot_r_script += '''benchmark_data <- rbind(plos_benchmark_data, benchmark_data)\n'''
  
          boxplot_r_script += '''
  pdf('%(analysis_file_prefix)s%(metric_column_name)s.pdf')
  %(boxplot_generation_commands)s
  dev.off()
  ''' % locals()
  
          r_script_filepath = os.path.join(self.output_directory, '{0}{1}.R'.format(self.analysis_file_prefix, metric_column_name))
          write_file(r_script_filepath, boxplot_r_script)
  
          # Run the R script
          run_r_script(r_script_filepath, cwd = self.output_directory)
  
          # Create a readable tabular file with the raw data
          metric_data_xy_file_lines = []
          benchmark_ids = sorted(metric_data.keys())
          metric_data_xy_file_lines.append('\t'.join(['Domain'] + benchmark_ids))
          for domain in common_domains:
              columns = [domain]
              for benchmark_id in benchmark_ids:
                  columns.append(metric_data[benchmark_id][domain])
              metric_data_xy_file_lines.append('\t'.join(map(str, columns)))
          write_file(os.path.join(self.output_directory, '{0}{1}_xy.txt'.format(self.analysis_file_prefix, metric_column_name)), '\n'.join(metric_data_xy_file_lines))
  
          # Compute the mean of the values per benchmark
          for benchmark_id, metric_values in sorted(metric_data.iteritems()):
              self.summary_stats[benchmark_id][metric_column_title] = numpy.median(metric_values.values())
  
  
      def generic_analysis_function(self, metric_data, common_domains, metric_column_name, metric_column_title):
          '''This function creates a text file with the metric data over the benchmarks and domains and runs R to
             generate boxplots of the metric data. The mean of the metrics over each benchmark is computed and added to
             the summary_stats dictionary.
          '''
  
          # Create the R data frame
          data_frame_values = []
          benchmark_ids = set()
          method_ids = set()
          for benchmark_id, domain_values in sorted(metric_data.iteritems()):
              benchmark_ids.add(benchmark_id)
              for domain, metric_value in sorted(domain_values.iteritems()):
                  benchmark_details = self.benchmark_details[benchmark_id]
                  if benchmark_details.get('title'):
                      # If the user specifies a title, that should be used
                      method = benchmark_details['title']
                  else:
                      method = benchmark_details['method']
                      if benchmark_details.get('kT') != None:
                          method += ', kT {0}'.format(benchmark_details.get('kT'))
                  method_ids.add(method)
                  data_frame_values.append((benchmark_details['benchmark'], method, domain, metric_value))
  
          benchmark_column = "','".join(d[0] for d in data_frame_values)
          method_column = "','".join(d[1] for d in data_frame_values)
          domain_column = "','".join(d[2] for d in data_frame_values)
          metric_data_column = ",".join(str(d[3]) for d in data_frame_values)
  
          # Create the R script
          boxplot_generation_commands = '''
  p <- ggplot(data = benchmark_data, aes(factor(method), {0}, fill=method))
  p <- p + stat_boxplot(geom ='errorbar', linetype="solid", width = 0.25, position = "dodge")
  p <- p + stat_boxplot(geom = "boxplot", linetype="solid", position = "dodge", width = 0.60, na.rm = TRUE) + xlab("Method") + ylab("{1}") + facet_grid(.~benchmark) + theme(legend.position="none")
  p
  '''.format(metric_column_name, metric_column_title)
  
          analysis_file_prefix = self.analysis_file_prefix
          boxplot_r_script = '''
  library(ggplot2)
  library(grid)
  
  '''
          published_data_dataframe = None
          if self.published_data_methods and metric_column_name != 'sequence_entropy':
              published_data_dataframe = get_data_frame(self.published_data_methods, metric_column_name, backrub_method_prefix = self.backrub_method_prefix, backrub_method_suffix = self.backrub_method_suffix)
              if published_data_dataframe:
                  boxplot_r_script += published_data_dataframe
              method_ids = method_ids.union(set(get_published_method_ids(self.published_data_methods, backrub_method_prefix = self.backrub_method_prefix, backrub_method_suffix = self.backrub_method_suffix)))
              benchmark_ids.add('_plos_one_published_data')
  
  
          boxplot_r_script += '''
  benchmark_data <- data.frame(benchmark=c('%(benchmark_column)s'),
                               method=c('%(method_column)s'),
                               domain=c('%(domain_column)s'),
                               %(metric_column_name)s=c(%(metric_data_column)s))
  ''' % locals()
          if published_data_dataframe:
              boxplot_r_script += '''benchmark_data <- rbind(plos_benchmark_data, benchmark_data)\n'''
  
          boxplot_r_script += '''
  pdf('%(analysis_file_prefix)s%(metric_column_name)s_boxplots.pdf')
  %(boxplot_generation_commands)s
  dev.off()
  ''' % locals()
  
          r_script_filepath = os.path.join(self.output_directory, '{0}{1}.R'.format(self.analysis_file_prefix, metric_column_name))
          write_file(r_script_filepath, boxplot_r_script)
  
          # Run the R script
          run_r_script(r_script_filepath, cwd = self.output_directory)
  
          # Create a readable tabular file with the raw data
          metric_data_xy_file_lines = []
          benchmark_ids = sorted(metric_data.keys())
          metric_data_xy_file_lines.append('\t'.join(['Domain'] + benchmark_ids))
          for domain in common_domains:
              columns = [domain]
              for benchmark_id in benchmark_ids:
                  columns.append(metric_data[benchmark_id][domain])
              metric_data_xy_file_lines.append('\t'.join(map(str, columns)))
          write_file(os.path.join(self.output_directory, '{0}{1}_xy.txt'.format(self.analysis_file_prefix, metric_column_name)), '\n'.join(metric_data_xy_file_lines))
  
          # Compute the mean of the values per benchmark
          for benchmark_id, metric_values in sorted(metric_data.iteritems()):
              self.summary_stats[benchmark_id][metric_column_title] = numpy.median(metric_values.values())
  
  
  if __name__ == '__main__':
  
      try:
          arguments = docopt.docopt(__doc__.format(**locals()))
          input_jsons = sorted(set(arguments['JSON_FILES']))
          output_directory = arguments['--out']
          overwrite_files = arguments['--overwrite']
          expectn = arguments['--expectn'] or None
          if expectn:
              try:
                  assert(expectn.isdigit())
                  expectn = int(expectn)
              except:
                  raise Exception('The --expectn argument should be an integer.')
          benchmark_ids = []
          if arguments['--benchmark']:
              benchmark_ids = arguments['--benchmark']
          analysis_file_prefix = arguments['--prefix'] or ''
  
          if 'all' in [s.lower() for s in arguments['--include_published_data']]:
              arguments['--include_published_data'] = published_backrub_methods
  
          Analyzer(input_jsons, benchmark_ids, output_directory, overwrite_files, expectn,
                   analysis_file_prefix = analysis_file_prefix,
                   published_data_methods = arguments['--include_published_data'],
                   backrub_method_prefix = arguments.get('--published_data_backrub_method_prefix') or '',
                   backrub_method_suffix = arguments.get('--published_data_backrub_method_suffix') or '',
                   scatterplot_colors_file = arguments.get('--scatterplot_colors')).run()
      except KeyboardInterrupt:
          print
  