dssp_secondary.py

# DSSP_SECONCARY.PY
# Trajectories of Secondary Structure by DSSP
#-------------------------------------------------------------------------------
#    This source code is part of PROSTAC
#
#    Written by Naoyuki Miyashita
#    Copyright (c) 2013, Naoyuki Miyashita, Japan
#
#    This program is free software; you can redistribute it and/or
#    modify it under the terms of the GNU General Public License
#    as published by the Free Software Foundation, either 
#    version 3 of the License, or (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied 
#    warranty of MERCHANTABILITY or FITNESS FOR 
#    A PARTICULAR PURPOSE.  
#    See the GNU General Public License for more details.
#
#    You should have received a copy of the GNU General 
#    Public License along with dssp_secondary.py.  
#    If not, see <http://www.gnu.org/licenses/>.
#-------------------------------------------------------------------------------
# This script uses "numpy", "MDAnalysis" and DSSP program
# import numpy, MDAnalysis, subprocess, os, sys,shutil, random
#    The keys  [" ","B","E","G","H","I","T","S"]
#    change to   0   1   2   3   4   5   6   7 
# 
# python dssp_secondary.py [psf filename] [dcd header] [dcd extention] \ 
#                          [initial file number] [final number]
# ex)
# python dssp_secondary.py initial.psf o.replica. .crdtrj 1 1
#
# written by Naoyuki Miyashita on 09/30/2013

import subprocess
import numpy as np
from MDAnalysis import *
import numpy.linalg
import os,sys,shutil,random

def chk_secondary_structure(strx):
    list_x = [" ","B","E","G","H","I","T","S"]
    icnt=0
    for item in list_x:
        if strx == item:
            break
        icnt=icnt+1
    return(icnt)

def copy_file(orig,copy):
    shutil.copyfile(orig,copy)
    return

def initial():
    psf = sys.argv[1]
    dcd_head = sys.argv[2]
    dcd_tail = sys.argv[3]
    ini_fn = int(sys.argv[4])
    fin_fn = int(sys.argv[5])
    number_of_files = fin_fn - ini_fn +1
    return (psf,dcd_head,dcd_tail,ini_fn,number_of_files)

def argv_read(dcd):
    filex=([])
    filex=dcd.split(".")
    if filex[-1] != "dcd":
      dcdfile = dcd.replace(filex[-1], "dcd")
      shutil.move(dcd, dcdfile)
    else:
      dcdfile = dcd
    return (dcdfile)

def looping():
    (psf,dcd_head,dcd_tail,ini_fn,number_of_files)=initial()
    print psf, dcd_head, dcd_tail
    for i in range(number_of_files):
        ii       = i + ini_fn
        filename = dcd_head + str(ii) + dcd_tail
        outfile  = dcd_head + str(ii) + ".out"
        ri       = random.randint(10,1000)
        temp_file = "tmp." + str(ri) + ".pdb"
        print i,filename
        print temp_file
        dcdfile=argv_read(filename)
        print dcdfile
        print outfile
        main_file(psf,dcdfile,temp_file,outfile)
        os.remove(str(temp_file))
    return

def main_file(psf,dcd,temp_file,outfile):
    u = Universe(psf,dcd)
    protein = u.selectAtoms("protein")
    output=np.array([])
    line1=np.array([])
    sam=np.array([0.0])
    #sam = 0.0
    lsw=0
    icnt=0
    f=open(str(outfile),"w")
    for ts in u.trajectory:
        second=([])
        pdb = Writer(str(temp_file), multiframe=False)
        pdb.write(protein) 
        pdb.close()
        output = subprocess.check_output(['dssp', str(temp_file)])
        line1=output.split("\n")
        for line in line1:
           if line.find("#  RESIDUE") > 0:
               lsw=1
               continue
           if lsw == 1:
               resn=line[0:5]
               seco=line[16:17]
               nseco=chk_secondary_structure(seco)
               second.append(nseco)
               #k=0
               #print nseco
               #if nseco == 4:
               #    aa=aa+1
               #    sam[k]=sam[k]+1
               #    print aa,sam[k],k,nseco
               #k=k+1
        second.insert(0,icnt)
        for i in second:
            f.write(str(i))
            f.write(" ")
        f.write("\n")
        del second
        lsw = 0
        icnt=icnt+1
    f.close()
    #for i in arange(sam):
    #    sam[i]=float(sam[i])/float(icnt)
    #print sam
    return

if __name__ == '__main__':
     looping()