dca_frustratometer.py

import random
import numpy
import os.path
import argparse

from pathlib import Path
basedir=Path(__file__).parent.absolute()
scratchdir=f'{basedir}/scratch'

def get_pfamID(pdbID, chain):
    import pandas as pd
    #ftp://ftp.ebi.ac.uk/pub/databases/msd/sifts/flatfiles/csv/pdb_chain_pfam.csv.gz
    df = pd.read_csv('%s/pdb_chain_pfam.csv' % basedir, header=1)
    if sum((df['PDB'] == pdbID.lower()) & (df['CHAIN'] == chain.upper())) != 0:
        pfamID=df.loc[(df['PDB'] == pdbID.lower()) & (df['CHAIN'] == chain.upper())]["PFAM_ID"].values[0]
    else:
        print('cant find pfamID')
        pfamID='null'
    return(pfamID)

def get_uniprotID(pdbID, chain):
    import urllib.request
    response = urllib.request.urlopen('http://www.bioinf.org.uk/cgi-bin/pdbsws/query.pl?plain=1&qtype=pdb&id=%s&chain=%s' % (pdbID, chain))
    html = response.read()
    info=html.split('\n')
    uniprotID=info[3].split(' ')[1]
    return uniprotID

def get_pfam_map(pdbID, chain):
    import pandas as pd
    #ftp://ftp.ebi.ac.uk/pub/databases/Pfam/mappings/pdb_pfam_mapping.txt
    df = pd.read_table('%s/pdb_pfam_mapping.txt' % basedir, header=0)
    if sum((df['PDB_ID'] == pdbID.upper()) & (df['CHAIN_ID'] == chain.upper())) != 0:
        start=df.loc[(df['PDB_ID'] == pdbID.upper()) & (df['CHAIN_ID'] == chain.upper())]["PdbResNumStart"].values[0]
        end=df.loc[(df['PDB_ID'] == pdbID.upper()) & (df['CHAIN_ID'] == chain.upper())]["PdbResNumEnd"].values[0]
    else:
        print('data not found')
        pfamID='null'
    return int(start), int(end)

def download_pfam(pfamID):
    import urllib.request
    urllib.request.urlretrieve('http://pfam.xfam.org/family/%s/alignment/full' % pfamID, "%s%s.stockholm" % (directory, pfamID))


def download_pdb(pdbID):
    import urllib.request
    urllib.request.urlretrieve('http://www.rcsb.org/pdb/files/%s.pdb' % pdbID, "%s%s.pdb" % (directory, pdbID))

def stockholm2fasta(pfamID):
    from Bio import AlignIO
    # rewrite Stockholm alignment in FASTA format
    input_handle = open("%s%s.stockholm" % (directory, pfamID), "rU")
    output_handle = open("%s%s.fasta" % (directory, pfamID), "w")
    alignments = AlignIO.parse(input_handle, "stockholm")
    AlignIO.write(alignments, output_handle, "fasta")
    output_handle.close()
    input_handle.close()

def filter_fasta(gap_threshold, pfamID, pdbID, chain, seq, resnos):
    from Bio import AlignIO
    import numpy
    import subprocess
    #gap_threshold=0.25
    pfam_start, pfam_end = get_pfam_map(pdbID, chain)
    mapped_seq = seq[resnos.index(pfam_start):resnos.index(pfam_end)+1]

    #print mapped fasta file
    f = open('%s%s%s_pfam_mapped.fasta' % (directory, pdbID, chain), 'w')
    f.write('>%s:%s pdb mapped to pfam\n' % (pdbID, chain))
    f.write(mapped_seq)
    f.close()

    submit=("muscle -profile -in1 %s%s.fasta -in2 %s%s%s_pfam_mapped.fasta -out %s%s%s.fasta" % (directory,pfamID,directory,pdbID,chain,directory,pdbID,chain))

    #print(submit)

    process = subprocess.Popen(submit.split(), stdout=subprocess.PIPE)
    process.communicate()
    
    # Filter sequences based on gaps in input sequence and gap threshold
    alignment = AlignIO.read(open("%s%s%s.fasta" % (directory, pdbID, chain)), "fasta")
    targetseq=alignment[-1].seq
    targetname=alignment[-1].name
    if targetseq=='':
        print("targetseq not found")

    output_handle = open("%s%s_msa_filtered.fasta" % (directory, pdbID), "w")
    target_array = numpy.array([list(targetseq)], numpy.character)
    bools = target_array != b'-'
    sequences_passed_threshold = 0
    for i, record in enumerate(alignment):
            record_array = numpy.array([list(record.seq)], numpy.character)
            aligned_sequence = record_array[bools]
            if float(sum(aligned_sequence==b'-'))/len(aligned_sequence) < gap_threshold:
                    output_handle.write(">%s\n" % record.id + aligned_sequence.tostring().upper().decode("utf-8") +'\n')
                    sequences_passed_threshold += 1
    output_handle.close()

    fastaseq=target_array.tostring().upper().decode("utf-8")
    
    stat_output = open(stat_output_file_name, "w")
    stat_output.write("FASTA_alignments " + str(len(alignment)) + "\n")
    stat_output.write("Filtered_alignments " + str(sequences_passed_threshold) + "\n")
    stat_output.close()

    return fastaseq, sequences_passed_threshold

def calc_plm(pdbID):
    import matlab.engine
    eng = matlab.engine.start_matlab()

    #import StringIO
    #out = StringIO.StringIO()

    fastafile=("%s%s_msa_filtered.fasta" % (directory, pdbID))
    outputfile=("%soutputfile.%s" % (directory, pdbID))
    lambda_h=0.01
    lambda_J=0.01
    reweighting_threshold=0.1
    nr_of_cores=1
    outputDistribution=("%soutputDistribution.%s" % (directory, pdbID))
    outputMatrix=("%soutputMatrix.%s" % (directory, pdbID))

    eng.addpath('%s/plm' % basedir, nargout=0)
    eng.addpath('%s/plm/functions' % basedir , nargout=0)
    eng.addpath('%s/plm/3rd_party_code/minFunc' % basedir, nargout=0)
    print('plmDCA_symmetric_mod7',fastafile,outputfile,
          lambda_h,lambda_J,reweighting_threshold,nr_of_cores,
          outputDistribution,outputMatrix)
    eng.plmDCA_symmetric_mod7(fastafile,outputfile,lambda_h,lambda_J,
                              reweighting_threshold,nr_of_cores,
                              outputDistribution,outputMatrix, nargout=0)#, stdout=out )


def read_dca(pdbID):
    import scipy.io
    import numpy
    mat = scipy.io.loadmat("%soutputMatrix.%s" % (directory, pdbID))
    h=mat['h']
    J=mat['J']
    q=mat['h'].shape[1]
    N=mat['h'].shape[0]
    fields = numpy.zeros((N, q))
    couplings = numpy.swapaxes(J,1,2)
    fields=h
    return fields, couplings

def vector(p1, p2):
    return [p2[0]-p1[0], p2[1]-p1[1], p2[2]-p1[2]]

def vabs(a):
    from math import sqrt
    return sqrt(pow(a[0],2)+pow(a[1],2)+pow(a[2],2))

def calc_dis(p1, p2):
    v = vector(p1, p2)
    return vabs(v)

def three2one(prot):
    """ translate a protein sequence from 3 to 1 letter code"""

    code = {"GLY" : "G", "ALA" : "A", "LEU" : "L", "ILE" : "I",
            "ARG" : "R", "LYS" : "K", "MET" : "M", "CYS" : "C",
            "TYR" : "Y", "THR" : "T", "PRO" : "P", "SER" : "S",
            "TRP" : "W", "ASP" : "D", "GLU" : "E", "ASN" : "N",
            "GLN" : "Q", "PHE" : "F", "HIS" : "H", "VAL" : "V",
            "M3L" : "K", "MSE" : "M", "CAS" : "C" }

    newprot = ""
    for a in prot:
        newprot += code.get(a, "?")

    return newprot

def calc_distances(pdbID, chainID):
    import numpy
    from Bio.PDB.PDBParser import PDBParser
    import warnings
    from Bio.PDB.PDBExceptions import PDBConstructionWarning
    warnings.simplefilter('ignore', PDBConstructionWarning)

    p = PDBParser(PERMISSIVE=1)
    struct_id = pdbID

    filename = directory + struct_id + ".pdb"

    s = p.get_structure(struct_id, filename)
    chains = s[0].get_list()
    sequence = []
    dis = []
    all_res = []
    resnos = []

    count = 0

    for chain in chains:
        if chain.get_id() == chainID.upper():
            for resid in chain:
                is_regular_res = resid.has_id('CA') and resid.has_id('O')
                res_id = resid.get_id()[0]

                if count==0:
                    pdbstart=resid.get_id()[1]
                count=count + 1

                if (res_id==' ' or res_id=='H_MSE' or res_id=='H_M3L' or res_id=='H_CAS') and is_regular_res:
                    all_res.append(resid)
                    sequence.append(resid.get_resname())
                    resnos.append(resid.get_id()[1])

    number_of_pdb_amino_acids = len(all_res)
    native_distances = numpy.zeros((number_of_pdb_amino_acids, number_of_pdb_amino_acids))

    for i in range(0, len(all_res)):
        dis.append([]);
        ires = all_res[i]
        xyz_CAi = ires['CA'].get_coord()
        for j in range(0, len(all_res)):
            jres = all_res[j]
            xyz_CAj = jres['CA'].get_coord()
            r = calc_dis(xyz_CAi, xyz_CAj)
            dis[i].append(r);
            native_distances[i,j]=r

    sequence=three2one(sequence)
    return native_distances, sequence, resnos, pdbstart

def map_fasta_to_pdb(pdbseq, fastaseq):
    if(fastaseq in pdbseq):
        pdbstart=pdbseq.find(fastaseq)
        fastastart=0
    elif(pdbseq in fastaseq):
        fastastart=fastaseq.find(pdbseq)
        pdbstart=0
        
    else:
        import re
        from Bio import pairwise2
        alignments = pairwise2.align.globalxx(fastaseq, pdbseq)

        fastastart=re.search("[A-Z]", alignments[0][1]).start()
        pdbstart=re.search("[A-Z]", alignments[0][0]).start()

    N=min(len(fastaseq),len(pdbseq))

    pdb_indices=range(pdbstart,N+pdbstart)
    dca_indices=range(fastastart,N+fastastart)
    map_to_dca = dict(zip(pdb_indices, dca_indices))
    map_to_pdb = dict(zip(dca_indices, pdb_indices))

    return dca_indices, pdb_indices, map_to_dca, map_to_pdb, len(pdb_indices)

def evaluate_configurational_frustration(sequence, number_of_decoys, number_of_residues):
    include_fields = True
    native_totals = numpy.zeros((number_of_residues, number_of_residues))
    decoy_means = numpy.zeros((number_of_residues, number_of_residues))
    decoy_stds = numpy.zeros((number_of_residues, number_of_residues))
    configurational_frustration_indices = numpy.zeros((number_of_residues, number_of_residues))
    for i in range(number_of_residues):
        i_dca = dca_indices[i]
        i_pdb = pdb_indices[i]
        amino_acid_i = amino_acid_type_dictionary[sequence[i_dca]]
        for j in range(i+minimum_sequence_separation, number_of_residues):
            j_dca = dca_indices[j]
            j_pdb = pdb_indices[j]
            amino_acid_j = amino_acid_type_dictionary[sequence[j_dca]]
            native_distance = native_distances[i_pdb, j_pdb]
            if native_distance < max_distance_threshold:
                native_coupling = couplings[i_dca, amino_acid_i, j_dca, amino_acid_j]
                native_totals[i, j] = native_coupling
                if include_fields:
                    native_fields = fields[i_dca, amino_acid_i] + fields[j_dca, amino_acid_j]
                    native_totals[i, j] += native_fields
                # Perform decoy calculations
                decoy_fields = numpy.zeros(number_of_decoys)
                decoy_couplings = numpy.zeros(number_of_decoys)
                decoy_totals = numpy.zeros(number_of_decoys)
                for decoy_i in range(number_of_decoys):
                    i_sequence_index = i_dca
                    j_sequence_index = j_dca
                    
                    #randomize indices
                    decoy_distance = max_distance_threshold + 1.0
                    while decoy_distance > max_distance_threshold:
                        # Draw random indices based on number of pdb residues
                        i_decoy,i_sequence_index=random.choice(list(enumerate(dca_indices)))
                        j_decoy,j_sequence_index=random.choice(list(enumerate(dca_indices)))
                        # Choose the pdb indices based on those random indices
                        i_sequence_index_pdb = map_to_pdb[i_sequence_index]
                        j_sequence_index_pdb = map_to_pdb[j_sequence_index]
                        # Calculate the distance to make sure it is below the threshold
                        decoy_distance = native_distances[i_sequence_index_pdb, j_sequence_index_pdb]
                        
                    #randomize amino acid types
                    i_amino_acid_type = random.choice(dca_indices)
                    j_amino_acid_type = random.choice(dca_indices)
                    amino_acid_i = amino_acid_type_dictionary[sequence[i_amino_acid_type]]
                    amino_acid_j = amino_acid_type_dictionary[sequence[j_amino_acid_type]]
                    
                    decoy_totals[decoy_i] = couplings[i_sequence_index, amino_acid_i, j_sequence_index, amino_acid_j]
                    if include_fields:
                        decoy_totals[decoy_i] += fields[i_sequence_index, amino_acid_i] + fields[j_sequence_index, amino_acid_j]
                        
                decoy_means[i, j] = numpy.mean(decoy_totals)
                decoy_stds[i, j] = numpy.std(decoy_totals)
                configurational_frustration_indices[i, j] = (native_totals[i, j]-decoy_means[i, j])/decoy_stds[i, j]
                
    return native_totals, decoy_means, decoy_stds, configurational_frustration_indices

def evaluate_mutational_frustration(sequence, number_of_decoys, number_of_residues):
    native_totals = numpy.zeros((number_of_residues, number_of_residues))
    decoy_means = numpy.zeros((number_of_residues, number_of_residues))
    decoy_stds = numpy.zeros((number_of_residues, number_of_residues))
    mutational_frustration_indices = numpy.zeros((number_of_residues, number_of_residues))
    for i in range(number_of_residues):
        i_dca = dca_indices[i]
        i_pdb = pdb_indices[i]
        amino_acid_i = amino_acid_type_dictionary[sequence[i_dca]]
        for j in range(i+minimum_sequence_separation, number_of_residues):
            j_dca = dca_indices[j]
            j_pdb = pdb_indices[j]
            amino_acid_j = amino_acid_type_dictionary[sequence[j_dca]]
            native_distance = native_distances[i_pdb, j_pdb]
            if native_distance < max_distance_threshold:
                native_coupling = couplings[i_dca, amino_acid_i, j_dca, amino_acid_j]
                native_totals[i, j] = native_coupling
                if include_fields:
                    native_fields = fields[i_dca, amino_acid_i] + fields[j_dca, amino_acid_j]
                    native_totals[i, j] += native_fields
                # Perform decoy calculations
                decoy_fields = numpy.zeros(number_of_decoys)
                decoy_couplings = numpy.zeros(number_of_decoys)
                decoy_totals = numpy.zeros(number_of_decoys)
                
                for decoy_i in range(number_of_decoys):
                    i_sequence_index = i_dca
                    j_sequence_index = j_dca
                    
                    #get random amino acids for i and j
                    i_amino_acid_type = random.choice(dca_indices)
                    j_amino_acid_type = random.choice(dca_indices)
                    amino_acid_i = amino_acid_type_dictionary[sequence[i_amino_acid_type]]
                    amino_acid_j = amino_acid_type_dictionary[sequence[j_amino_acid_type]]
                    
                    #compute energy terms for (i,j) pair
                    decoy_totals[decoy_i] += couplings[i_sequence_index, amino_acid_i, j_sequence_index, amino_acid_j]
                    if include_fields:
                        decoy_totals[decoy_i] += fields[i_sequence_index, amino_acid_i] + fields[j_sequence_index, amino_acid_j]
                        
                    #in mutational mode, all (i,k) and (j,k) pairs also contribute to the decoy energy
                    #let's compute those pairs now
                    for k in range(number_of_residues):
                        amino_acid_k = amino_acid_type_dictionary[sequence[k]]
                        #make sure to exclude native interactions that we already counted
                        if (k==i_sequence_index or k==j_sequence_index):
                            continue
                            
                        rik = native_distances[pdb_indices[i], pdb_indices[k]]
                        if (rik < max_distance_threshold):
                            decoy_totals[decoy_i] += couplings[i_sequence_index, amino_acid_i, k, amino_acid_k]
                        rjk = native_distances[pdb_indices[j], pdb_indices[k]]
                        if (rik < max_distance_threshold):
                            decoy_totals[decoy_i] += couplings[j_sequence_index, amino_acid_j, k, amino_acid_k]
                            
                decoy_means[i, j] = numpy.mean(decoy_totals)
                decoy_stds[i, j] = numpy.std(decoy_totals)
                mutational_frustration_indices[i, j] = (native_totals[i, j]-decoy_means[i, j])/decoy_stds[i, j]
                
    return native_totals, decoy_means, decoy_stds, mutational_frustration_indices


def evaluate_single_residue_frustration(sequence, number_of_decoys, number_of_residues):
    native_totals = numpy.zeros(number_of_residues)
    decoy_means = numpy.zeros(number_of_residues)
    decoy_stds = numpy.zeros(number_of_residues)
    single_residue_frustration_indices = numpy.zeros(number_of_residues)
    for i in range(number_of_residues):
        i_dca = dca_indices[i]
        i_pdb = pdb_indices[i]
        amino_acid_i = amino_acid_type_dictionary[sequence[i_dca]]
        
        #compute native energy
        native_totals[i] = evaluate_singleres_native_energy(i_dca, i_pdb, amino_acid_i, sequence)
        
        #compute decoy energies
        decoy_totals = evaluate_singleres_decoy_energy(i_dca, i_pdb, amino_acid_i, sequence, number_of_decoys)
                            
        decoy_means[i] = numpy.mean(decoy_totals)
        decoy_stds[i] = numpy.std(decoy_totals)
        single_residue_frustration_indices[i] = (native_totals[i] - decoy_means[i])/decoy_stds[i]
                
    return native_totals, decoy_means, decoy_stds, single_residue_frustration_indices

def evaluate_singleres_decoy_energy(i_dca, ipdb, amino_acid_i, sequence, number_of_decoys):
    total_energy = 0
    decoy_totals = numpy.zeros(number_of_decoys)
    for decoy_i in range(number_of_decoys):
        #get random amino acids for i
        i_amino_acid_type = random.choice(dca_indices)
        amino_acid_i = amino_acid_type_dictionary[sequence[i_amino_acid_type]]
        
        #compute decoy_energy
        decoy_totals[decoy_i] = evaluate_singleres_native_energy(i_dca, ipdb, amino_acid_i, sequence)

    return decoy_totals

def evaluate_singleres_native_energy(i_dca, i_pdb, amino_acid_i, sequence):
    total_energy = 0
    for j in range(0, number_of_residues):
        j_dca = dca_indices[j]
        j_pdb = pdb_indices[j]
        if i_dca == j_dca:
            continue
        amino_acid_j = amino_acid_type_dictionary[sequence[j_dca]]
        native_distance = native_distances[i_pdb, j_pdb]
        if (native_distance < max_distance_threshold) and (abs(i_dca-j_dca)>minimum_sequence_separation):
            total_energy = couplings[i_dca, amino_acid_i, j_dca, amino_acid_j]
        if include_fields:
            total_energy += fields[i_dca, amino_acid_i]
            
    return total_energy

########MAIN PROGRAM##########
if __name__=='__main__':
    # Parse arguments
    parser = argparse.ArgumentParser()
    parser.add_argument("pdbID", help="id of pdb")
    parser.add_argument("chain", help="chain id of pdb")
    args = parser.parse_args()

    #Parameters for frustration calculations
    randomize_amino_acid_types=True
    randomize_indices=True
    include_fields=True
    compute_configurational_frustration=True
    compute_mutational_frustration=True
    compute_single_residue_frustration=True
    n_decoys=5000
    minimum_sequence_separation=2
    max_distance_threshold=9.5

    pdbID = args.pdbID#'5pti'
    chain = args.chain#'a'

    directory = ('%s/dca-frustratometer/automated/%s%s/' % (scratchdir, pdbID, chain))
    if not os.path.exists(directory):
        os.makedirs(directory)
        
    configurational_frustration_output_file_name=("%sconfigurational_frust_output.%s" % (directory, pdbID))
    mutational_frustration_output_file_name=("%smutatational_frust_output.%s" % (directory, pdbID))
    single_residue_frustration_output_file_name=("%ssing_frust_output.%s" % (directory, pdbID))
    stat_output_file_name = ("%s%s.info" % (directory, pdbID))
    status_file_name = ("%s%s.status" % (directory, pdbID))

    status_output = open(status_file_name, "w")

    #get the pfam and uniprot ID from the pdb info
    pfamID = get_pfamID(pdbID, chain)
    #uniprotID = get_uniprotID(pdbID, chain)

    status_output.write("pdbID: " + pdbID + "\n")
    print("pdbID: " + pdbID + "\n")
    status_output.write("pfamID: " + pfamID + "\n")
    print("pfamID: " + pfamID + "\n")
    #status_output.write("uniprotID: " + uniprotID + "\n")
    #print("uniprotID: " + uniprotID + "\n")

    #download pdb file if it doesn't exist
    if os.path.isfile(directory + pdbID + '.pdb') == False:
        download_pdb(pdbID)
        status_output.write("PDB downloaded" + "\n")
        print("PDB downloaded" + "\n")


    #download pfam msa if it doesn't exist
    if os.path.isfile(directory + pfamID + '.stockholm') == False:
        download_pfam(pfamID)
        status_output.write("Pfam MSA downloaded" + "\n")
        print("Pfam MSA downloaded" + "\n")

    if os.path.isfile(directory + pfamID + '.fasta') == False:
        stockholm2fasta(pfamID)
        status_output.write("stockholm coverted" + "\n")
        print("stockholm coverted" + "\n")

    #obtain information from PDB file (pairwise separation(native_distances) and pdb_sequences)
    native_distances, pdb_sequence, resnos, pdb_start = calc_distances(pdbID, chain)

    #filter sequences and get the target sequence and number of sequences remaining in the alignment
    #set filter threshold
    filter_threshold = 0.25
    msa_query_sequence, n_sequences = filter_fasta(filter_threshold, pfamID, pdbID, chain, pdb_sequence, resnos)
    status_output.write("FASTA filtered" + "\n")
    print("FASTA filtered" + "\n")

    #run the plm calculation which saves appropriate files
    if os.path.isfile("%soutputMatrix.%s" % (directory, pdbID)) == False:
        calc_plm(pdbID)
        status_output.write("PLM computed" + "\n")
        print("PLM computed" + "\n")

    #read in plm calculations fields and couplings
    fields, couplings = read_dca(pdbID)

    #map the fasta sequence on to the sequence of the pdb file
    dca_indices, pdb_indices, map_to_dca,map_to_pdb, number_of_residues = map_fasta_to_pdb(pdb_sequence, msa_query_sequence)
    status_output.write("parameter initialization complete" + "\n")
    print("parameter initialization complete" + "\n")


    #print some frustratometer statistics
    stat_output = open(stat_output_file_name, "a")
    stat_output.write("msa_query_sequence " + msa_query_sequence + "\n")
    stat_output.write("msa_query_sequence_length " + str(len(msa_query_sequence)) + "\n")
    stat_output.write("pdb_sequence " + pdb_sequence + "\n")
    stat_output.write("pdb_sequence_length " + str(len(pdb_sequence)) + "\n")
    stat_output.close()

    # Define parameters
    number_of_amino_acid_types = 21
    amino_acid_type_dictionary = {
    '-':0,
    'A':1,
    'C':2,
    'D':3,
    'E':4,
    'F':5,
    'G':6,
    'H':7,
    'I':8,
    'K':9,
    'L':10,
    'M':11,
    'N':12,
    'P':13,
    'Q':14,
    'R':15,
    'S':16,
    'T':17,
    'V':18,
    'W':19,
    'Y':20,
    'X':0
    }


    # Compute and write out configurational frustration indices
    target = open(("%s%s_configurational_frustration.tcl" % (directory, pdbID)), 'w')
    atomselect=0
    if compute_configurational_frustration:
        status_output.write("Computing configurational frustration...\n")
        print("Computing configurational frustration...\n")
        configurational_native_totals, configurational_decoy_means, configurational_decoy_stds, configurational_frustration_indices = evaluate_configurational_frustration(msa_query_sequence, n_decoys, number_of_residues)
        configurational_frustration_output_file = open(configurational_frustration_output_file_name, "w")
        for i in range(number_of_residues):
            for j in range(i+minimum_sequence_separation, number_of_residues):
                if native_distances[pdb_indices[i], pdb_indices[j]] < max_distance_threshold:
                    configurational_frustration_output_file.write("%4d %s %4d %s %8.2f %8.2f %8.2f %8.2f \n" % (pdb_indices[i]+pdb_start,
                                                                                                         msa_query_sequence[dca_indices[i]],
                                                                                                         pdb_indices[j]+pdb_start,
                                                                                                         msa_query_sequence[dca_indices[j]],
                                                                                                         configurational_native_totals[i, j],
                                                                                                         configurational_decoy_means[i, j],
                                                                                                         configurational_decoy_stds[i, j],
                                                                                                         configurational_frustration_indices[i, j]))
                    
                    if (configurational_frustration_indices[i,j] > 0.78 or configurational_frustration_indices[i,j] < -1):
                        pdbi=pdb_indices[i]+pdb_start-1
                        pdbj=pdb_indices[j]+pdb_start-1
                        target.write("set sel%d [atomselect top \"resid %d and name CA\"]\n" % (pdbi, pdbi+1))
                        target.write("set sel%d [atomselect top \"resid %d and name CA\"]\n" % (pdbj, pdbj+1))
                        target.write("lassign [atomselect%d get {x y z}] pos1\n" % atomselect)
                        atomselect = atomselect + 1
                        target.write("lassign [atomselect%d get {x y z}] pos2\n" % atomselect)
                        atomselect = atomselect + 1
                        if configurational_frustration_indices[i,j] > 0.78:
                            target.write("draw color green\n")
                        else:
                            target.write("draw color red\n")
                        target.write("draw line $pos1 $pos2 style solid width 1\n")
        target.write("mol modselect 0 top \"all\"\n")
        target.write("mol modstyle 0 top newcartoon\n")
        target.write("mol modcolor 0 top colorid 15\n")
        configurational_frustration_output_file.close()
    target.close()

    # Compute and write out mutational frustration indices
    target = open(("%s%s_mutational_frustration.tcl" % (directory, pdbID)), 'w')
    atomselect=0
    if compute_mutational_frustration:
        status_output.write("Computing mutational frustration...\n")
        print("Computing mutational frustration...\n")
        mutational_native_totals, mutational_decoy_means, mutational_decoy_stds, mutational_frustration_indices = evaluate_mutational_frustration(msa_query_sequence, n_decoys, number_of_residues)
        mutational_frustration_output_file = open(mutational_frustration_output_file_name, "w")
        for i in range(number_of_residues):
            for j in range(i+minimum_sequence_separation, number_of_residues):
                if native_distances[pdb_indices[i], pdb_indices[j]] < max_distance_threshold:
                    mutational_frustration_output_file.write("%4d %s %4d %s %8.2f %8.2f %8.2f %8.2f \n" % (pdb_indices[i]+pdb_start,
                                                                                                         msa_query_sequence[dca_indices[i]],
                                                                                                         pdb_indices[j]+pdb_start,
                                                                                                         msa_query_sequence[dca_indices[j]],
                                                                                                         mutational_native_totals[i, j],
                                                                                                         mutational_decoy_means[i, j],
                                                                                                         mutational_decoy_stds[i, j],
                                                                                                         mutational_frustration_indices[i, j]))
                    if (mutational_frustration_indices[i,j] > 0.78 or mutational_frustration_indices[i,j] < -1):
                        pdbi=pdb_indices[i]+pdb_start-1
                        pdbj=pdb_indices[j]+pdb_start-1
                        target.write("set sel%d [atomselect top \"resid %d and name CA\"]\n" % (pdbi, pdbi+1))
                        target.write("set sel%d [atomselect top \"resid %d and name CA\"]\n" % (pdbj, pdbj+1))
                        target.write("lassign [atomselect%d get {x y z}] pos1\n" % atomselect)
                        atomselect = atomselect + 1
                        target.write("lassign [atomselect%d get {x y z}] pos2\n" % atomselect)
                        atomselect = atomselect + 1
                        if mutational_frustration_indices[i,j] > 0.78:
                            target.write("draw color green\n")
                        else:
                            target.write("draw color red\n")
                        target.write("draw line $pos1 $pos2 style solid width 1\n")
        target.write("mol modselect 0 top \"all\"\n")
        target.write("mol modstyle 0 top newcartoon\n")
        target.write("mol modcolor 0 top colorid 15\n")
        mutational_frustration_output_file.close()
    target.close()

    # Compute and write out single residue frustration indices
    if compute_single_residue_frustration:
        status_output.write("Computing single residue frustration...\n")
        print("Computing single residue frustration...\n")
        target = open(("%s%s_singleres_frustration.tcl" % (directory, pdbID)), 'w')
        atomselect = 0
        single_residue_native_totals, single_residue_decoy_means, single_residue_decoy_stds, single_residue_frustration_indices = evaluate_single_residue_frustration(msa_query_sequence, n_decoys, number_of_residues)
        single_residue_frustration_output_file = open(single_residue_frustration_output_file_name, "w")
        for i in range(number_of_residues):
            single_residue_frustration_output_file.write("%4d %s %8.2f %8.2f %8.2f %8.2f \n" % (pdb_indices[i]+pdb_start,
                                                                                                msa_query_sequence[dca_indices[i]],
                                                                                                single_residue_native_totals[i],
                                                                                                single_residue_decoy_means[i],
                                                                                                single_residue_decoy_stds[i],
                                                                                                single_residue_frustration_indices[i]))
            atomselect = atomselect + 1
            #temp = 0.5*abs(single_residue_frustration_indices[i])
            target.write("mol addrep 0\n")
            target.write("mol modselect %d 0 resid %d\n" % (atomselect, pdb_indices[i]+pdb_start))
            target.write("mol modstyle %d 0 VDW %f 12.000000\n" % (atomselect, 0.5*abs(single_residue_frustration_indices[i])))
            target.write("mol modmaterial %d 0 Transparent\n" % atomselect)
            if single_residue_frustration_indices[i] > 0:
                target.write("mol modcolor %d 0 ColorID 7\n" % atomselect)
            else:
                target.write("mol modcolor %d 0 ColorID 1\n" % atomselect)
        target.write("mol modselect 0 top \"all\"\n")
        target.write("mol modstyle 0 top newcartoon\n")
        target.write("mol modcolor 0 top colorid 15\n")
        target.close()
        single_residue_frustration_output_file.close()

    status_output.close()