PCA and FEL through gromacs

 

Principal component analysis and Free energy landscapes of protein systems in gromacs.

>gmx covar -f md_t14cnopbc.xtc -s md_t14c.tpr -o eigenval.xvg -v eigenvect.trr -xpm covara.xpm

[least squares fit:c-alpha; covariance: c alpha]

>gmx xpm2ps -f covara.xpm -o covara.eps -do covar.m2p

>gmx anaeig -v eigenvect.trr -f md_t14cnopbc.xtc -s md_t14c.tpr -first 1 -last 2 -proj proj_eig.xvg -2d 2d_proj.xvg

>gmx sham -f 2d_proj.xvg -ls gibbs.xpm -notime #for free energy, "sham" module uses PC1 and PC2 file

gmx xpm2ps -f gibbs.xpm  -o gibbs.eps -rainbow red

python xpm2txt.py -f gibbs.xpm -o gibbs.txt 

#xpm2txt.py is mentioned here. Just copy and paste this script. Run in command line as mentioned above

#!/usr/bin/env python

import sys

"""

Utility tool to convert xpm files generated by GROMACS to a 3-column text file.

"""

USAGE = "USAGE: xpm2txt.py -f <input xpm file> -o <output txt file> [-s]\n"

USAGE+= "Options:\n"

USAGE+= "\t-s\t(int)\tSorts the output by a given column"

USAGE+= "\n" # always keep this line

# Parse arguments

read_input, read_output, sort = False, False, False

xpm_file, out_file, column_sort = None, None, None

for arg in sys.argv[1:]:

    if read_input:

        read_input = False

        xpm_file = arg

    elif read_output:

        read_output = False

        out_file = arg

    elif sort:

        sort = False

        column_sort = int(arg)

    if arg[0] == "-":

        if arg == "-f":

            read_input = True

            continue

        elif arg == "-o":

            read_output = True

            continue

        elif arg == "-s":

            sort = True

        else:

            print USAGE

            sys.stderr.write('ERROR: Option not recognized: %s\n' %arg)

            sys.exit(1)

if not xpm_file:

    print USAGE

    sys.stderr.write('ERROR: You forgot to provide an input file.\n')

    sys.exit(1)

if not out_file:

    out_file = "out.txt"

# Parse XPM file

xpm_handle = open(xpm_file)

xpm_data = []

x_axis, y_axis = [], []

letter_to_value = {}

for line in xpm_handle:

    if line.startswith("/* x-axis"):

        x_axis = map(float, line.split()[2:-2]) # We trim the last value

    if line.startswith("/* y-axis"):

        y_axis = map(float, line.split()[2:-2]) # We trim the last value

    if line.startswith('"') and x_axis and y_axis: # Read data

        xpm_data.insert(0, line.strip().strip(',')[1:-1])

    if line.startswith('"') and len(line.split()) > 4:

        letter = line.split()[0][1:]

        value = float(line.split()[-2][1:-1])

        letter_to_value[letter] = value

xpm_handle.close()

# Match x/y/data

txt_values = []

for y_index, data_value in enumerate(xpm_data):

    y_value = y_axis[y_index]

    for x_index, x_value in enumerate(x_axis):

        txt_values.append([x_value, y_value, letter_to_value[data_value[x_index]]])

# Apply sorting if requested

if column_sort:

    try:

        txt_values.sort(key=lambda x: x[column_sort-1])

    except IndexError:

        print USAGE

        sys.stderr.write('ERROR: Column not found (%s)\n' %(column_sort))

        sys.exit(1)

# Print to file

out_handle = open(out_file, 'w')

for x, y, z in txt_values:

    out_handle.write("%3.5f\t%3.5f\t%3.5f\n" %(x,y,z))

out_handle.close()

    

  For any queries regarding the script, kindly drop your message in the comments section below.