How to use opposite method in Sure

Best Python code snippet using sure_python

twist_discrete.py

Source:twist_discrete.py

1# -*- coding: utf-8 -*-2"""3Created on Fri Oct 14 23:01:08 201645@author: tunaz6"""789import sys1011import re12import numpy13import math1415from scipy import interpolate1617from mpl_toolkits.mplot3d import Axes3D18import matplotlib.pyplot as plt192021class Tee(object):22    def __init__(self, *files):23        self.files = files24    def write(self, obj):25        for f in self.files:26            f.write(obj)27            f.flush() # If you want the output to be visible immediately28    def flush(self) :29        for f in self.files:30            f.flush()31            32            33if len(sys.argv) == 1 or len(sys.argv) == 2 or len(sys.argv) > 3  :34    print "Usage: python twist_region.py pdbfilename.pdb direction"35    sys.exit(0)363738lines_File = [lines1.rstrip('\n') for lines1 in open(sys.argv[1], 'r')]   394041try:42    direction = int(sys.argv[2]) # give direction as input43    '''44    direction= 0 means forward direction, 45    direction = 1 means backword direction 46    '''47    if direction != 0 and direction != 1 :48        print "please input either 0 or 1. 0 = forward, 1 = backword"49        sys.exit(0)50except ValueError:  51        print "please input either 0 or 1. 0 = forward, 1 = backword"52        sys.exit(0)          5354#sys.stdout = open("output_all.pdb", "w")55out_filename = 'twist_angles.txt'56outf = open(out_filename, 'w')5758sys.stdout = open('output_pdb.pdb','w')59#outf2 = open(out_filename2, 'w')60#lines_File = [lines1.rstrip('\n') for lines1 in open('1aky.pdb')]    61#lines_File = [lines1.rstrip('\n') for lines1 in open('1aop.pdb')]  62#lines_File = [lines1.rstrip('\n') for lines1 in open('2p8y.pdb')] #index problem63#lines_File = [lines1.rstrip('\n') for lines1 in open('1s04.pdb')] #complex structure (multple same chain duplicate index)64#lines_File = [lines1.rstrip('\n') for lines1 in open('1atz.pdb')] 65#lines_File = [lines1.rstrip('\n') for lines1 in open('2qtr.pdb')] 66#lines_File = [lines1.rstrip('\n') for lines1 in open('1d5t.pdb')] 67#lines_File = [lines1.rstrip('\n') for lines1 in open('1elu.pdb')]68#lines_File = [lines1.rstrip('\n') for lines1 in open('1a12.pdb')]69#lines_File = [lines1.rstrip('\n') for lines1 in open('1b3a.pdb')]7071#direction = 0 # direction= 0 means forward direction, direction = 1 means backword direction so change here        72            73def draw_interpolation(input_array):74    #print "input  ", input_array75    #print "length  ", len(input_array)76    if  len(input_array) > 2 :77        data = numpy.array(input_array)78        data = data.astype(float)79        80       # print data #ok81           82        data = data[0:100].transpose()83        #print "hi", data #ok84        #now we get all the knots and info about the interpolated spline85		#spacing = int(len(input_array)/.3)86        tck, u= interpolate.splprep(data,k=2)   87        #new = interpolate.splev(numpy.linspace(0,1,50), tck)88        #new = interpolate.splev(numpy.linspace(0,1,int(len(input_array)/.8)), tck) #unit length given 0.189        #new = interpolate.splev(numpy.linspace(0,1,int(len(input_array)/.5)), tck) #we used it for ACM_BCB90        new = interpolate.splev(numpy.linspace(0, 1, int(len(input_array) /.05)), tck)9192        #now lets plot it!93        94        #fig = plt.figure()95        #ax = Axes3D(fig)96        #ax.plot(data[0], data[1], data[2], label='originalpoints', lw =2, c='Dodgerblue')97        #ax.plot(new[0], new[1], new[2], label='fit', lw =2, c='red')98        #ax.legend()99        #plt.savefig('junk.png')100        #plt.show() 101        return new[0], new[1], new[2]             102103104#print (lines_File1)105106direction_opposite = 1 #hydrogen bond 107108count_beta_sheet = 0109count_beta_strand = 0110counter = 0111counter_index = []112problem = 0113a=[]114beta_sheet = []115atom = []116splitted_line = []117count_atom = 0118for sheet_line in lines_File: #find each line in all lines119     if sheet_line.startswith('SHEET'):     #find those lines start with "SHEET"120         #count_atom = count_atom+1121         #print (sheet_line)       #print the lines started with "SHEET"        122         #beta_sheet.append ( re.sub("\s\s+" , " ",sheet_line).split()) #ok for all except 1123         beta_sheet.append([sheet_line[0:6], sheet_line[7:10], sheet_line[11:14], sheet_line[14:16],sheet_line[17:20], sheet_line[21], sheet_line[22:27], sheet_line[28:31], sheet_line[32], sheet_line[33:37], sheet_line[37], sheet_line[38:40],sheet_line[41:45],sheet_line[45:48], sheet_line[49], sheet_line[50:55], sheet_line[56:60], sheet_line[60:63], sheet_line[64], sheet_line[65:69], sheet_line[69]])124     if sheet_line.startswith('ATOM'):  #find those lines start with "ATOM"125         count_atom = count_atom+1126         #print (sheet_line)       #print the lines started with "ATOM"127         atom.append( [sheet_line[:6], sheet_line[6:11], sheet_line[12:16], sheet_line[17:20], sheet_line[21], sheet_line[22:26], sheet_line[30:38], sheet_line[38:46], sheet_line[46:54]])128         #atom.append( [sheet_line[0:6], sheet_line[6:11], sheet_line[12:16], sheet_line[16:17],sheet_line[17:20], sheet_line[21:22], sheet_line[22:26], sheet_line[26:27], sheet_line[30:38], sheet_line[38:46], sheet_line[46:54], sheet_line[46:54], sheet_line[54:60], sheet_line[60:66], sheet_line[76:78], sheet_line[78:80]])129         #atom.append ( re.sub("\s\s+" , " ",sheet_line).split())130#print (int(atom[1836][5].strip()))131#print (atom[1836])132133strand_vector=[]134center_vector = []135CA_vector_between_strand = []136amino_acid = []137138139#common hydrogen bond140strand_vector_opposite=[]141center_vector_opposite = []142CA_vector_between_strand_opposite = []143amino_acid_opposite = []144145146poi = 0147track = 1148149poi_opposite = 0 #hydrogen bond150151152for sheet_index in range (0,len(beta_sheet)): # handle if there is less than 3 atoms in a strand153    if int(beta_sheet[sheet_index][1]) == int(beta_sheet[sheet_index][3]):154        counter = counter + 1155156#print ("number of beta sheet ",counter) #ok157out_char = "number of beta sheet " + str(counter) + '\n\n'158outf.write(out_char)            159for sheet_index in range (0,len(beta_sheet)): # handle if there is less than 3 atoms in a strand160    if abs((int(beta_sheet[sheet_index][9]))  - int(beta_sheet[sheet_index][6])) >= 2 :161        if int(beta_sheet[sheet_index][1]) == int(beta_sheet[sheet_index][3]):162            #counter = counter + 1163            if track == 0 :164                counter_index.append(int(beta_sheet[sheet_index][1]) - problem)165                track = 1166            else:   167                counter_index.append(int(beta_sheet[sheet_index][1]))168                track = 1169    else:170        #print "hi " #ok171        if  (int(beta_sheet[sheet_index][3]) > 2):           172             problem = problem + 1173             #print problem  #ok174             track = 0175            176            177#for i in range (0,len(counter_index)):178    #print ("counter_index ",counter_index[i]) #ok  179     180for sheet_index in range (0,len(beta_sheet)): # 5 should be number of beta stand in one beta sheet181    #print("love ",beta_sheet[sheet_index][10])182   183    atom_start_index=(int(beta_sheet[sheet_index][6]))184    #print("start index ",atom_start_index) #ok185    atom_end_index=(int(beta_sheet[sheet_index][9])) 186    #print("end index ",atom_end_index) #ok187    sheet_chain = beta_sheet[sheet_index][5]188    #print sheet_chain189190    mid_point = []191    array_strand_index =[]192    amino_acid_array = []193    aa_index = []194    #array_strand_index = [[0 for i in range(2)] for i in range(2)] #2 should be change according to length of one beta sheet195    array_x = []196    array_y = []197    array_z = []198    CA_x = []199    CA_y = []200    CA_z = []201    202    #print (count_atom)203    204    for i in range (0,count_atom) : 205        206        if (atom_start_index <= int(atom[i][5].strip()) <=  atom_end_index)  :207            208            if  atom[i][2].strip() == 'N' or atom[i][2].strip() == 'C': 209                #print "duplicate ", int(atom[i][5].strip())210                211                if sheet_chain == atom[i][4].strip(): #this is needed to differentiate same startindex & endindex but different chains212                    array_strand_index.append(atom[i][5].strip())                  213                    array_x.append(atom[i][6].strip())214                    array_y.append(atom[i][7].strip())215                    array_z.append(atom[i][8].strip())216                    length=len(array_strand_index)  217            if atom[i][2].strip() == 'CA':218                if sheet_chain == atom[i][4].strip():219                    aa_index.append( i )             220                    amino_acid_array.append(atom[i][3].strip())        221                    CA_x.append(atom[i][6].strip())222                    CA_y.append(atom[i][7].strip())223                    CA_z.append(atom[i][8].strip())224                    length_CA = len (CA_x)225    #print ("strand index  ",array_strand_index) #ok226   # print (aa_index)227    #print ("Amino acid  ",amino_acid_array)  #ok        228   229   230   231    232    if direction == 0: #forward direction #ok233        midpoint_x_bef = []234        midpoint_y_bef = []235        midpoint_z_bef = []  236        237      238        239        midpoint_x = []240        midpoint_y = []241        midpoint_z = []242        243        new_x = []244        new_y = []245        new_z = []  246        247       248        #coef = [1,1,1,1,1,1,1,1]249        #coef = []250        251        c_index = 1252        n_index = 2253        while c_index < length-1 and n_index <  length-1 :  #mid point calculation254            #print ("c_index ",c_index, "n_index ", n_index)255            midpoint_x_bef.append((float(array_x[c_index]) + float(array_x[n_index]))/float(2.0))    #making midpoint       256            midpoint_y_bef.append((float(array_y[c_index]) + float(array_y[n_index]))/float(2.0))   #making midpoint 257            midpoint_z_bef.append((float(array_z[c_index]) + float(array_z[n_index]))/float(2.0))  #making midpoint 258           259            c_index = c_index + 2260            n_index = n_index +2  261        262        #print ("midpoint_x_bef  ", midpoint_x_bef ) 263        #print ("midpoint_y_bef  ", midpoint_y_bef  )264        #print ("midpoint_z_bef  ", midpoint_z_bef  )265        midpoint_x=midpoint_x_bef266        midpoint_y=midpoint_y_bef267        midpoint_z=midpoint_z_bef268        '''269        midpoint_draw = []270        draw_index = 0271        while draw_index < len (midpoint_x_bef) :272            midpoint_point = []273            mid_x = midpoint_x_bef [draw_index]274            mid_y = midpoint_y_bef [draw_index]275            mid_z = midpoint_z_bef [draw_index]           276            midpoint_point.append (mid_x)277            midpoint_point.append (mid_y)278            midpoint_point.append (mid_z)279            midpoint_draw.append (midpoint_point) 280            draw_index = draw_index + 1281        #print   (midpoint_draw)  #ok282        pdb_format = '6s5s1s4s1s3s1s1s4s1s3s8s8s6s10s2s3s'283        #shape = 'ATOM 1 N LEU A 81 25.865 39.612 19.376 1.00 53.28 N '284        shape = ['ATOM', '1','N', 'LEU', 'A', '81', '25.865', '39.612', '19.376']285        if len(midpoint_draw) <= 2:286             midpoint_x=midpoint_x_bef287             midpoint_y=midpoint_y_bef288             midpoint_z=midpoint_z_bef289             290             for i in range (len( midpoint_x)):291                poi = poi +1292                                           293                #print "%-6s%5s %4s %3s %s%4s    %8s%8s%8s\n"%tuple(shape) ok294                print "%-6s%5s %4s %3s %s%4s    %8s%8s%8s"%tuple(['ATOM',poi,'S','LEU', 'A', '81',round(midpoint_x[i],3),round(midpoint_y[i],3),round(midpoint_z[i],3)]) #print new atom position in PDB file295                296             297        else:    298            value = draw_interpolation (midpoint_draw)   #making 3D cubic spline interpolation for midpoints299            #print value[0] #ok300            midpoint_x=value[0]  #new midpoint after cubic spline301            midpoint_y=value[1]  #new midpoint after cubic spline302            midpoint_z=value[2]   #new midpoint after cubic spline303        304            305           306            pdb_length = len( midpoint_x)307          308            for i in range (pdb_length):309                poi = poi +1          310                print "%-6s%5s %4s %3s %s%4s    %8s%8s%8s"%tuple(['ATOM',poi,'S','LEU', 'A', '81',round(midpoint_x[i],3),round(midpoint_y[i],3),round(midpoint_z[i],3)]) #print new atom position in PDB file311           '''     312       313       314        315        316       317        #print ("midpoint_x ", midpoint_x ) 318         319        320     321     322        length_midpoint_array = len(midpoint_x)323        vector_array =[]324        cen_array = []325        midpoint_index = 1326        while midpoint_index < length_midpoint_array:327            point = []328            cen_point = []329            temp_mid_index = midpoint_index - 1330            vec_x = midpoint_x [temp_mid_index] -  midpoint_x [temp_mid_index+1]  #making midpoint vector331            vec_y = midpoint_y [temp_mid_index] -  midpoint_y [temp_mid_index+1]  #making midpoint vector332            vec_z = midpoint_z [temp_mid_index] -  midpoint_z [temp_mid_index+1]  #making midpoint vector333            334            cen_x = (float(midpoint_x [temp_mid_index]) + float( midpoint_x [temp_mid_index+1])) / float(2.0)  #making center of midpoint vector335            cen_y =(float( midpoint_y [temp_mid_index]) + float(midpoint_y [temp_mid_index+1]))/float(2.0)   #making center of midpoint vector336            cen_z = (float(midpoint_z [temp_mid_index]) + float( midpoint_z [temp_mid_index+1]))/float(2.0)   #making center of midpoint vector337            338            #print ("vector x ",vec_x) 339           # print ("vector y ",vec_y) 340           # print ("vector z ",vec_z) 341            #make unit vector342            magnitude = math.sqrt(pow(vec_x, 2) + pow(vec_y, 2) + pow(vec_z, 2))343            344           345            #print("vector length ",magnitude)346            347            point.append(vec_x/magnitude) # unit vector348            point.append(vec_y/magnitude)  # unit vector349            point.append(vec_z/magnitude)  # unit vector350            351      352            #point.append(vec_x) #ok353            #point.append(vec_y)354            #point.append(vec_z)355      356         357            cen_point.append(cen_x) #ok358            cen_point.append(cen_y)359            cen_point.append(cen_z)360            361            362            363            364            vector_array.append(point)365            366            cen_array.append(cen_point)367            368            midpoint_index = midpoint_index +1369        #print "vector array ",vector_array370        371        372        aa = []373        aa_index = 1374        while aa_index < len(amino_acid_array) - 1:   # to get amino acid value375             aa.append(amino_acid_array[aa_index]) 376             aa_index = aa_index +1377        amino_acid.append(aa)378        379        #print ("x cordinate of CA in each strand" CA_x )380    381        #print("length ", length)382        #print ("CA length ",length_CA)383        CA_x_new = []384        CA_y_new = []385        CA_z_new = []386        CA_vector = []387        CA_index = 1388        while CA_index < length_CA - 1:389            CA_point = []390            CA_x_new=(float(CA_x[CA_index]))391            CA_y_new=(float(CA_y[CA_index]))392            CA_z_new=(float(CA_z[CA_index]))393            CA_point.append (CA_x_new)394            CA_point.append (CA_y_new)395            CA_point.append (CA_z_new)396            CA_vector.append (CA_point)   397            CA_index = CA_index + 1398        #draw_interpolation (CA_vector) #draw spline according to carbon alpha 399400        CA_vector_between_strand.append(CA_vector)401        #draw_interpolation (CA_vector_between_strand)402403         404    if direction_opposite == 1:   #backword direction405    406        midpoint_x_bef_opposite = []407        midpoint_y_bef_opposite = []408        midpoint_z_bef_opposite = []  409410        midpoint_x_opposite = []411        midpoint_y_opposite = []412        midpoint_z_opposite = [] 413        414       415        416        c_index_opposite = length-3417        n_index_opposite = length-2418        while c_index_opposite > 0  and n_index_opposite > 0 : 419            #print ("c_index ",c_index, "n_index ", n_index)420            midpoint_x_bef_opposite.append((float(array_x[c_index_opposite]) + float(array_x[n_index_opposite]))/float(2.0))421            midpoint_y_bef_opposite.append((float(array_y[c_index_opposite]) + float(array_y[n_index_opposite]))/float(2.0))422            midpoint_z_bef_opposite.append((float(array_z[c_index_opposite]) + float(array_z[n_index_opposite]))/float(2.0)) 423            c_index_opposite = c_index_opposite - 2424            n_index_opposite = n_index_opposite - 2   425       426       427        midpoint_x_opposite = midpoint_x_bef_opposite428        midpoint_y_opposite = midpoint_y_bef_opposite429        midpoint_z_opposite = midpoint_z_bef_opposite430       431        '''432        midpoint_draw_opposite = []433        draw_index_opposite = len (midpoint_x_bef_opposite) -1434        while draw_index_opposite >= 0 :435            midpoint_point_opposite = []436            mid_x_opposite = midpoint_x_bef_opposite [draw_index_opposite]437            mid_y_opposite = midpoint_y_bef_opposite [draw_index_opposite]438            mid_z_opposite = midpoint_z_bef_opposite [draw_index_opposite]           439            midpoint_point_opposite.append (mid_x_opposite)440            midpoint_point_opposite.append (mid_y_opposite)441            midpoint_point_opposite.append (mid_z_opposite)442            midpoint_draw_opposite.append (midpoint_point_opposite) 443            draw_index_opposite = draw_index_opposite - 1444            445        #draw_interpolation (midpoint_draw)446       447        if len(midpoint_draw_opposite) <= 2:448             midpoint_x_opposite = midpoint_x_bef_opposite 449             midpoint_y_opposite = midpoint_y_bef_opposite450             midpoint_z_opposite = midpoint_z_bef_opposite451             452             for i in range (len( midpoint_x_opposite)):453                poi_opposite = poi_opposite +1454                                           455                #print "%-6s%5s %4s %3s %s%4s    %8s%8s%8s\n"%tuple(shape) ok456               # print "%-6s%5s %4s %3s %s%4s    %8s%8s%8s"%tuple(['ATOM',poi_opposite,'S','LEU', 'A', '81',round(midpoint_x_opposite[i],3),round(midpoint_y_opposite[i],3),round(midpoint_z_opposite[i],3)]) #print new atom position in PDB file457                458             459        else:    460            value_opposite = draw_interpolation (midpoint_draw_opposite)   #making 3D cubic spline interpolation for midpoints461            #print value_opposite[0] #ok462            midpoint_x_opposite = value_opposite[0]  #new midpoint after cubic spline463            midpoint_y_opposite = value_opposite[1]  #new midpoint after cubic spline464            midpoint_z_opposite = value_opposite[2]   #new midpoint after cubic spline465        466            467           468            pdb_length_opposite = len( midpoint_x_opposite)469          470            for i in range (pdb_length_opposite):471                poi_opposite = poi_opposite +1          472                #print "%-6s%5s %4s %3s %s%4s    %8s%8s%8s"%tuple(['ATOM',poi_opposite,'S','LEU', 'A', '81',round(midpoint_x_opposite[i],3),round(midpoint_y_opposite[i],3),round(midpoint_z_opposite[i],3)]) #print new atom position in PDB file473                474       475       476        midpoint_x_opposite = midpoint_x_bef_opposite477        midpoint_y_opposite = midpoint_y_bef_opposite478        midpoint_z_opposite = midpoint_z_bef_opposite479        '''480    481        #print ("midpoint_x ", midpoint_x )      482        length_midpoint_array_opposite = len (midpoint_x_opposite)   483        vector_array_opposite =[]484        cen_array_opposite = []485        486        #midpoint_index_opposite = length_midpoint_array487        #while midpoint_index > 1:488        489        midpoint_index_opposite = 1 #big change in this loop490        while midpoint_index_opposite < length_midpoint_array_opposite:491            point_opposite = []492            cen_point_opposite = []493            494            temp_mid_index_opposite = midpoint_index_opposite - 1495            #vec_x = midpoint_x [temp_mid_index] -  midpoint_x [temp_mid_index-1]496            #vec_y = midpoint_y [temp_mid_index] -  midpoint_y [temp_mid_index-1]497            #vec_z = midpoint_z [temp_mid_index] -  midpoint_z [temp_mid_index-1]498            499            vec_x_opposite = midpoint_x_opposite [temp_mid_index_opposite] -  midpoint_x_opposite [temp_mid_index_opposite+1]500            vec_y_opposite = midpoint_y_opposite [temp_mid_index_opposite] -  midpoint_y_opposite [temp_mid_index_opposite+1]501            vec_z_opposite = midpoint_z_opposite [temp_mid_index_opposite] -  midpoint_z_opposite [temp_mid_index_opposite+1]502            503            #print ("vector ",vec_x,vec_y,vec_z) 504            505            #cen_x = (float(midpoint_x [temp_mid_index]) + float( midpoint_x [temp_mid_index-1])) / float(2.0)  #making center of midpoint vector506            #cen_y =(float( midpoint_y [temp_mid_index]) + float(midpoint_y [temp_mid_index-1]))/float(2.0)   #making center of midpoint vector507            #cen_z = (float(midpoint_z [temp_mid_index]) + float( midpoint_z [temp_mid_index-1]))/float(2.0)   #making center of midpoint vector508            509            cen_x_opposite = (float(midpoint_x_opposite [temp_mid_index_opposite]) + float( midpoint_x_opposite [temp_mid_index_opposite+1])) / float(2.0)  #making center of midpoint vector510            cen_y_opposite =(float( midpoint_y_opposite [temp_mid_index_opposite]) + float(midpoint_y_opposite [temp_mid_index_opposite+1]))/float(2.0)   #making center of midpoint vector511            cen_z_opposite = (float(midpoint_z_opposite [temp_mid_index_opposite]) + float( midpoint_z_opposite [temp_mid_index_opposite+1]))/float(2.0)   #making center of midpoint vector512            513              #make unit vector514            magnitude_opposite = math.sqrt(pow(vec_x_opposite, 2) + pow(vec_y_opposite, 2) + pow(vec_z_opposite, 2))515            #print("vector length opposite ",magnitude_opposite)516            #point.append(vec_x/magnitude) #not ok517            #point.append(vec_y/magnitude)518            #point.append(vec_z/magnitude)519            520            '''521            point_opposite.append(vec_x_opposite) #ok522            point_opposite.append(vec_y_opposite)523            point_opposite.append(vec_z_opposite)524            '''525            point_opposite.append(vec_x_opposite/magnitude_opposite) #ok526            point_opposite.append(vec_y_opposite/magnitude_opposite)527            point_opposite.append(vec_z_opposite/magnitude_opposite)528            529            530            531            cen_point_opposite.append(cen_x_opposite) #ok532            cen_point_opposite.append(cen_y_opposite)533            cen_point_opposite.append(cen_z_opposite)534     535            vector_array_opposite.append(point_opposite)536            537            cen_array_opposite.append(cen_point_opposite)538            539            midpoint_index_opposite = midpoint_index_opposite + 1540            541            #midpoint_index = midpoint_index - 1 542        543        aa_opposite = []544        aa_index_opposite = len(amino_acid_array) - 2545        while aa_index_opposite > 0:   # to get amino acid value546             aa_opposite.append(amino_acid_array[aa_index_opposite]) 547             aa_index_opposite = aa_index_opposite - 1 548        amino_acid_opposite.append(aa_opposite)549   550   551        #print ("x cordinate of CA in each strand", CA_x )552    553        #print("length ", length)554        #print ("CA length ",length_CA)555        CA_x_new_opposite = []556        CA_y_new_opposite = []557        CA_z_new_opposite = []558        CA_vector_opposite = []559        CA_index_opposite = length_CA - 2560        while CA_index_opposite > 0 :561            CA_point_opposite = []562            CA_x_new_opposite=(float(CA_x[CA_index_opposite]))563            CA_y_new_opposite=(float(CA_y[CA_index_opposite]))564            CA_z_new_opposite=(float(CA_z[CA_index_opposite]))565            CA_point_opposite.append (CA_x_new_opposite)566            CA_point_opposite.append (CA_y_new_opposite)567            CA_point_opposite.append (CA_z_new_opposite)568            CA_vector_opposite.append (CA_point_opposite)   569            CA_index_opposite = CA_index_opposite - 1570        CA_vector_between_strand_opposite.append(CA_vector_opposite)571        572        573    strand_vector.append(vector_array)  574    575    center_vector.append(cen_array)576    577    #for hydrogen bond calculation578    strand_vector_opposite.append(vector_array_opposite)  579    580    center_vector_opposite.append(cen_array_opposite)581    582   #print(beta_sheet[sheet_index][1])  583    584    585586 587   588    589    '''590    print ("vector ",vector_array)591    print (atom_start_index)592    print (array_strand_index) 593    print(length)594    595    print (array_x) 596    print (array_y) 597    print (array_z) 598    print (midpoint_x)599    print (midpoint_y)600    print (midpoint_z)601    print (length_midpoint_array)602    '''603print "\n"604length_strand_vector = len(strand_vector)605#print (strand_vector)606#print(amino_acid) # ok607#print ("vector strand length ",length_strand_vector)608#print ("CA_vector_between_strand ",len(CA_vector_between_strand)) # ok.  vector_strand length and CA_vector same609#print ("CA_vector_between_strand ",CA_vector_between_strand) 610#draw_interpolation (CA_vector_between_strand)611612def distance(x,y):   613    return numpy.sqrt(numpy.sum((x-y)**2))614615616def find_nearest_distance(a,B, pair_aa1, pair_aa2):617    618    619    #print ("a ", a)620    min_point = B[0]621    #min_dist = abs(numpy.linalg.norm(a- numpy.array(B[0])))622    min_dist = math.sqrt ( pow((a[0]-B[0][0]),2) + pow((a[1]-B[0][1]),2) + pow((a[2]-B[0][2]),2)  )623    #min_dist = 12000624   # print ("mindistance ", min_dist)625    #print("len of B",len(B)) #ok626    aa_min_index = 0627    for i in range (0, len(B)): 628        #dist = []629                       630        b = numpy.array(B[i])        631        #print ("b ", b) #ok632        #dist=(numpy.linalg.norm(a-b))633        try_dist = math.sqrt ( pow((a[0]-b[0]),2) + pow((a[1]-b[1]),2) + pow((a[2]-b[2]),2)  )634        #print ("distance ", try_dist)635        #dist=abs(distance(a,b))636        if try_dist < min_dist:637            min_dist = try_dist638            min_point = b639            aa_min_index = i640 641    642    return (min_point, aa_min_index)643644645def calculate_twist_angle(a,minimum_dist_point,j,aa_min_index2):646        647    #print ("mindist point ", minimum_dist_point)648    dot_product = numpy.dot(a,minimum_dist_point) # A.B649    #print (dot_product)650       651    a_sum = pow(a[0],2)+pow(a[1],2)+pow(a[2],2)652    a_root=math.sqrt(a_sum)653    b_sum = pow(minimum_dist_point[0],2)+pow(minimum_dist_point[1],2)+pow(minimum_dist_point[2],2)654    b_root=math.sqrt(b_sum)655    a_b = a_root*b_root656    657    angle_dot = math.acos(dot_product/a_b)658    #print (angle_dot)659    twist_angle = math.degrees(angle_dot)660    661    '''662    ### cross product ###663    cross_product = numpy.cross(a,minimum_dist_point) # A*B664    #print (cross_product)665    sum_square_cross_product = pow(cross_product[0],2) +  pow(cross_product[1],2)+ pow(cross_product[2],2)666    #print (sum_square_cross_product)667    root_sum_square_cross_product = math.sqrt(sum_square_cross_product)668    #print (root_sum_square_cross_product)669   670   671   672    angle_cross = math.asin(root_sum_square_cross_product/a_b)673    #print (angle_cross)674    print(math.degrees(angle_cross))675    ######## end of cross product ###########676    '''677    if twist_angle > 90.:678        twist_angle = 180.0 - twist_angle679   # sys.stdout = original680   # print(" Amino Acid",pair_aa1[j], " and ", pair_aa2[aa_min_index2],":", "Twist angle::" , twist_angle) #necessary681    #print("Twist angle::" , round(twist_angle,3)) #necessary682    out_char = "Twist angle:: " + str(round(twist_angle,3)) + '\n'683    outf.write(out_char)684    return twist_angle685       # print ( "amino acid ", pair_aa1[j], " ", pair_aa2[aa_min_index2])686minimum_dist_point = []  687minimum_dist_point_dup = []688flag = 0689all_avg_twist_angle = []  690avg_count = 0    691avg_track = 1692pair_array = []693len_array = []694i = 0695while i < length_strand_vector : # pairwise calculation 696    #print ("beta sheet index", beta_sheet[i][1])697    #print ("beta strand", beta_sheet[i])698    699    #print(int(beta_sheet[i][1]), " : ", int(beta_sheet[i][3]))700    #print ("i ", i)701    #print ("flag ",flag)702    avg_twist_count = []703            704    if(int(beta_sheet[i][1])>=int(beta_sheet[i][3])):705        flag = int(beta_sheet[i][1])706        avg_count = 0707       708        #print ("avg count ",avg_count)709    else:710        flag = flag + 1711        avg_count = avg_count + 1712        A_t = []713        B_t= []714        A = []715        B = []716        717        #for hydrogen bond718        A_t_opposite = []719        B_t_opposite  = []720        A_opposite  = []721        B_opposite  = []722        723        alpha_A_t = []724        alpha_B_t = []725        alpha_A = []726        alpha_B = []727        728          #for hydrogen bond729        alpha_A_t_opposite = []730        alpha_B_t_opposite = []731        alpha_A_opposite = []732        alpha_B_opposite = []733        734        pair_aa1_t = []735        pair_aa2_t = []736        pair_aa1 = []737        pair_aa2 = []738        739         740        #for hydrogen bond741        pair_aa1_t_opposite = []742        pair_aa2_t_opposite = []743        pair_aa1_opposite = []744        pair_aa2_opposite = []745        746        747        index_vec = i748        A_t = strand_vector[i] #copy strand to A array749        B_t = strand_vector[index_vec+1] #copy strand to B array 750        751        #for hydrogen bond752        A_t_opposite = strand_vector_opposite[i] #copy strand to A array753        B_t_opposite = strand_vector_opposite[index_vec+1] #copy strand to B array 754        755        #alpha carbon756        alpha_A_t = CA_vector_between_strand[i] #copy strand to A array757        alpha_B_t = CA_vector_between_strand[index_vec+1] #copy strand to B array 758        759         #for hydrogen bond760        alpha_A_t_opposite = CA_vector_between_strand_opposite[i] #copy strand to A array (alpha-carbon) backword direction761        alpha_B_t_opposite = CA_vector_between_strand_opposite[index_vec+1] #copy strand to B array (alpha-carbon) backword direction762        763        '''764        #center vector765        alpha_A_t = center_vector[i] #copy strand to A array766        alpha_B_t = center_vector[index_vec+1] #copy strand to B array 767        768         #for hydrogen bond769        alpha_A_t_opposite = center_vector_opposite[i] #copy strand to A array (alpha-carbon) backword direction770        alpha_B_t_opposite = center_vector_opposite[index_vec+1] #copy strand to B array (alpha-carbon) backword direction771        '''772        773        pair_aa1_t = amino_acid[i] #copy strand to A array774        pair_aa2_t = amino_acid[index_vec+1] #copy strand to B array 775        776        #hydrogen bond777        pair_aa1_t_opposite = amino_acid_opposite[i] #copy strand to A array778        pair_aa2_t_opposite = amino_acid_opposite[index_vec+1] #copy strand to B array 779        780        781      #  print ("A_t ",A_t)782       783      #  print ("B_t ",B_t)784       # sys.stdout = original785        #print ("beta strand 1", beta_sheet[i]) #necessary786        #print ("beta strand 2", beta_sheet[i+1]) #necessary787        out_char = "beta strand 1" + str(beta_sheet[i]) + '\n'788        outf.write(out_char)789        out_char = "beta strand 2" + str(beta_sheet[i+1]) + '\n'790        outf.write(out_char)791        #print ("A_t ",len(A_t))792       793        #print ("B_t ",len(B_t))794        if len(A_t) > len(B_t):795            A = B_t796            B = A_t797            798            #hydrogen bond799            A_opposite = B_t_opposite800            B_opposite = A_t_opposite801            802            alpha_A = alpha_B_t803            alpha_B = alpha_A_t804            805             #hydrogen bond806            alpha_A_opposite = alpha_B_t_opposite807            alpha_B_opposite = alpha_A_t_opposite808            809            pair_aa1 = pair_aa2_t810            pair_aa2 = pair_aa1_t811            812            #hydrogen bond813            pair_aa1_opposite = pair_aa2_t_opposite814            pair_aa2_opposite = pair_aa1_t_opposite815        else:816            A = A_t817            B = B_t818            819             #hydrogen bond820            A_opposite = A_t_opposite821            B_opposite = B_t_opposite822            823            alpha_A = alpha_A_t824            alpha_B = alpha_B_t 825            826            #hydrogen bond827            alpha_A_opposite = alpha_A_t_opposite828            alpha_B_opposite = alpha_B_t_opposite829            830            pair_aa1 = pair_aa1_t831            pair_aa2 = pair_aa2_t832            833             #hydrogen bond834            pair_aa1_opposite = pair_aa1_t_opposite835            pair_aa2_opposite = pair_aa2_t_opposite836        837        #aa_min_index2 = 0838        k = 0839        twist_angle_return = 0.0840        add_1 = 0 841        add_n = 0842        add_hp = 0 #hydrogen bond843        pick_index = 0844        845        '''for calculation of avg twist angle between two beta strands'''846        twist_count = 0847        sum_twist_angle = 0.0848        avg_twist_angle = []849        850        #taking minimum of 4 consecutive  angles' sum851        cons_twist_angle = []852       853        854        for j in range (0, len(A_opposite)):855            a = numpy.array(A[j])  856            857            #hydrogen bond858            a_opposite = numpy.array(A_opposite[j])859            860            alpha_a = numpy.array(alpha_A[j])  #this is needed for calculating center of vector in 3D space861            862            #hydrogen bond863            alpha_a_opposite = numpy.array(alpha_A_opposite[j])864           865               866            #minimum_dist_point, aa_min_index2 = find_nearest_distance(a, B, pair_aa1[j], pair_aa2)867            #print ("k ", k)868            if k == 0  : # determining first pair869              870                minimum_dist_point, aa_min_index2 = find_nearest_distance(alpha_a, alpha_B, pair_aa1[j], pair_aa2)#nearest neighbour checking it will work for determining first pair871                #minimum_dist_point, aa_min_index2 = find_nearest_distance(a, B, pair_aa1[j], pair_aa2)872               #hydrogen bond873                minimum_dist_point_opposite, aa_min_index2_opposite = find_nearest_distance(alpha_a_opposite, alpha_B_opposite, pair_aa1_opposite[j], pair_aa2_opposite)#nearest neighbour checking it will work for determining first pair874                #minimum_dist_point_opposite, aa_min_index2_opposite = find_nearest_distance(a_opposite, B_opposite, pair_aa1_opposite[j], pair_aa2_opposite)875               876               # sys.stdout = original877                print(" MIN_INDEX_Forward", aa_min_index2)   878                print(" MIN_INDEX_Backword", aa_min_index2_opposite)    879                #print(" MIN_POINT ", minimum_dist_point)880                #print(" MIN_B ", B[aa_min_index2])881                882                '''883                twist_angle_return = calculate_twist_angle(a,B[aa_min_index2],j,aa_min_index2)884                sum_twist_angle = sum_twist_angle + twist_angle_return885                twist_count = twist_count + 1             886                pick_index = aa_min_index2887                '''888                889                if aa_min_index2 == 0 and len(A) <= aa_min_index2_opposite: #hydrogen bond parallel and anti-parallel case equal length strand or small strand is inside of big strand (complete overlap)890                    add_1 = 1891                    892                    twist_angle_return = calculate_twist_angle(a,B[aa_min_index2],j,aa_min_index2)893                    cons_twist_angle.append(twist_angle_return)   894                    sum_twist_angle = sum_twist_angle + twist_angle_return895                    twist_count = twist_count + 1             896                    pick_index = aa_min_index2897                    898                if aa_min_index2 == 0  and len(A) > aa_min_index2_opposite: #hydrogen bond parallel and anti-parallel case   small strand is outside (partial overlap) of big strand899                    add_hp = 1900                    901                    twist_angle_return = calculate_twist_angle(a_opposite,B_opposite[aa_min_index2_opposite],j,aa_min_index2)902                    cons_twist_angle.append(twist_angle_return)   903                    sum_twist_angle = sum_twist_angle + twist_angle_return904                    twist_count = twist_count + 1             905                    pick_index = aa_min_index2_opposite906                    907                if aa_min_index2 != 0:908                    add_n = 1909                    twist_angle_return = calculate_twist_angle(a,B[aa_min_index2],j,aa_min_index2)910                    cons_twist_angle.append(twist_angle_return)   911                    sum_twist_angle = sum_twist_angle + twist_angle_return912                    twist_count = twist_count + 1             913                    pick_index = aa_min_index2914                915            else: #next pairs916                if add_1 == 1 and k < len(B):917                    minimum_dist_point_dup = B[k]918                    pick_index = k919                    twist_angle_return = calculate_twist_angle(a,minimum_dist_point_dup,j,pick_index)920                    cons_twist_angle.append(twist_angle_return)   921                    cons_twist_angle.append(twist_angle_return)   922                    sum_twist_angle = sum_twist_angle + twist_angle_return923                    twist_count = twist_count + 1924                    925                if add_hp == 1 and  int(beta_sheet[i+1][11])!=-1  and  pick_index < len(B_opposite) - 1:		#parallel  small strand is outside (partial overlap) of big strand                         926                   # j = pick_index - 1927                    pick_index = pick_index + 1928                  929                    print "k ",k930                    print "pick_index",pick_index931                    minimum_dist_point_dup = B_opposite[pick_index]932                    #print "minimum_dist_point_dup ",minimum_dist_point_dup #ok933                   # print "a_opposite ",a_opposite934            935                    twist_angle_return = calculate_twist_angle(a_opposite,minimum_dist_point_dup,j,pick_index)936                    cons_twist_angle.append(twist_angle_return)   937                    sum_twist_angle = sum_twist_angle + twist_angle_return938                    twist_count = twist_count + 1939                    940                if add_hp == 1 and  int(beta_sheet[i+1][11])==-1  and  pick_index != 0: # anti-parallel  small strand is outside (partial overlap) of big strand                         941                   # j = pick_index - 1942                    pick_index = pick_index - 1 #may be need to change later 943                  944                    print "k ",k945                    print "pick_index",pick_index946                    minimum_dist_point_dup = B_opposite[pick_index]947                    #print "minimum_dist_point_dup ",minimum_dist_point_dup #ok948                   # print "a_opposite ",a_opposite949                    950                   951                    twist_angle_return = calculate_twist_angle(a_opposite,minimum_dist_point_dup,j,pick_index)952                    cons_twist_angle.append(twist_angle_return)   953                    sum_twist_angle = sum_twist_angle + twist_angle_return954                    twist_count = twist_count + 1955                956                957                if add_n == 1 and  int(beta_sheet[i+1][11])!=-1 and k < len(B) - aa_min_index2: # Parallel beta strands 	958                    pick_index = pick_index + 1959                    minimum_dist_point_dup = B[pick_index]960                    twist_angle_return = calculate_twist_angle(a,minimum_dist_point_dup,j,pick_index)961                    cons_twist_angle.append(twist_angle_return)   962                    sum_twist_angle = sum_twist_angle + twist_angle_return963                    twist_count = twist_count + 1964                    965                if add_n == 1 and   int(beta_sheet[i+1][11])==-1  and pick_index!=0: # anti-Parallel beta strands 		966                    pick_index = pick_index - 1           967                    minimum_dist_point_dup = B[pick_index]968                    twist_angle_return = calculate_twist_angle(a,minimum_dist_point_dup,j,pick_index)969                    cons_twist_angle.append(twist_angle_return)   970                    sum_twist_angle = sum_twist_angle + twist_angle_return971                    twist_count = twist_count + 1972            973                                                                   974            k = k + 1975                 976        if  twist_count != 0:  977            avg_twist_angle.append (float(sum_twist_angle) / float(twist_count))  978           979            cons_sum_array = []980            if len(cons_twist_angle) > 4:981                for c in range (0,len(cons_twist_angle) - 3):982                    cons_sum = 0.0983                    for cc in range (c,c+4):       984                         cons_sum = cons_twist_angle[cc] + cons_sum985                    cons_sum_array.append(cons_sum)986                pair_array.append((round(min(cons_sum_array)  / 4.0,3)))987                len_array.append(len(cons_sum_array))988                #print "sum of consecutive 4 angles ", cons_sum_array989                #print "min of consecutive 4 angles ",min(cons_sum_array)990                out_char = "min of sum of consecutive 4 angles :: " + str(round(min(cons_sum_array),3)) + '\n'991                outf.write(out_char)992                out_char = "Twist per angle :: " + str(round(min(cons_sum_array)  / 4.0 ,3)) + '\n'993                outf.write(out_char)994            995            elif len(cons_twist_angle) == 4:996                for c in range (0,len(cons_twist_angle) - 2):997                    cons_sum = 0.0998                    for cc in range (c,c+3):       999                         cons_sum = cons_twist_angle[cc] + cons_sum1000                    cons_sum_array.append(cons_sum)1001                pair_array.append((round(min(cons_sum_array)  / 3.0,3)))1002                len_array.append(len(cons_sum_array))1003                #print "sum of consecutive 4 angles ", cons_sum_array1004                #print "min of consecutive 4 angles ",min(cons_sum_array)1005                out_char = "min of sum of consecutive 3 angles :: " + str(round(min(cons_sum_array),3)) + '\n'1006                outf.write(out_char)  1007                out_char = "Twist per angle :: " + str(round(min(cons_sum_array) / 3.0 ,3) ) + '\n'1008                outf.write(out_char)1009                1010            elif len(cons_twist_angle) == 2 or len(cons_twist_angle) == 3 or len(cons_twist_angle) == 1:1011               #print "smallest twist ",min(cons_twist_angle)1012                pair_array.append((round(min(cons_twist_angle) , 3 )))1013                len_array.append(len(cons_sum_array))1014                out_char = "smallest angles :: " + str(round(min(cons_twist_angle),3)) + '\n'1015                outf.write(out_char)1016            avg_track = 11017        else:1018            avg_track = 01019        1020        for t in range (len(avg_twist_angle)):   1021            #print ("avg_twist_angle between strands ", round(avg_twist_angle[t],3)) #necessary1022            #print ( round(avg_twist_angle[t],3)) #necessary1023            out_char = "avg_twist_angle between strands  " + str( round(avg_twist_angle[t],3) ) + '\n'1024            outf.write(out_char)1025            all_avg_twist_angle.append(avg_twist_angle[t])            1026     1027        #print ("...........................")#necessary1028        #print cons_twist_angle#ok1029      1030        outf.write("................................................"+"\n")1031   1032    i = i +11033#print ("avg_twist_angle between strands ", all_avg_twist_angle) #ok     103410351036'''for calculation of avg twist angle in a beta sheet'''1037outf.write("\n\n")10381039twist_angle_beta_sheet = []1040avg_avg_twist_angle_beta_sheet = []1041jc = 0   1042for i in range (len(counter_index)):  1043    sum_avg_twist_angle = 0.0  1044    j = jc1045    #print jc, " ", counter_index[i]-1, " ",j #ok1046    while j < (jc + counter_index[i]-1):1047        #print "j ",j1048        sum_avg_twist_angle = sum_avg_twist_angle + all_avg_twist_angle[j]  1049        sum_avg_twist_angle = round(sum_avg_twist_angle,3)       1050        j = j + 1    1051    jc = jc + counter_index[i]-11052    #avg_avg_twist_angle_beta_sheet.append(sum_avg_twist_angle)1053    twist_angle_beta_sheet.append(float(sum_avg_twist_angle)/float(counter_index[i]-1))1054    #print ("avg_twist angle in all beta sheets ", round(twist_angle_beta_sheet[i],3)) #necessary1055  1056    #print ( round(twist_angle_beta_sheet[i],3)) 1057    out_char = "avg_twist angle in all beta sheets " + str( round(twist_angle_beta_sheet[i],3)) + '\n'1058    outf.write(out_char)1059   1060#print ("avg_twist angle in all beta sheets ", float(avg_avg_twist_angle_beta_sheet[0])/float(counter_index[0]-1))1061#print pair_array #ok10621063outf.write("................................................"+"\n")10641065#print pair_array # ok1066#print len_array #ok1067'''1068#### According to minimum value pair #############10691070kc = 01071for i in range (len(counter_index)):  1072    j = kc1073    min_arr = []1074    while j < (kc + counter_index[i]-1):1075        min_arr.append (pair_array[j])1076        1077        j = j + 1  1078    kc = kc + counter_index[i]-11079    arr = numpy.array (min_arr) 1080    ind1 , ind2 = arr.argsort()[:2] 1081    min_val1 =   arr[ind1]1082    min_val2 =   arr[ind2]1083    #print min_val1, min_val21084    out_char = "Min pair val1, val2 :: " + str(min_val1) + ",  "+ str(min_val2)+'\n'1085    outf.write(out_char)1086    1087    min_avg_twist = float(min_val1 + min_val2)/2.01088    #print min_avg_twist1089    out_char = "Min Twist :: " + str(min_avg_twist) + '\n'1090    outf.write(out_char)10911092outf.write("................................................"+"\n")10931094#### According to longest pair #############1095kl = 01096for il in range (len(counter_index)):  1097    jl = kl1098    min_arr_l = []1099    while jl < (kl + counter_index[il]-1):1100        min_arr_l.append (len_array[jl])1101        1102        jl = jl + 1  1103    kl = kl + counter_index[il]-11104    arr_len = numpy.array (len_array) 1105    ind1 , ind2 = arr_len.argsort()[::-1][:2]1106    max_val1 =   pair_array[ind1]1107    max_val2 =   pair_array[ind2]1108   # print ind1, ind2 #ok1109    out_char = "longest pair val1, val2 :: " + str(max_val1) + ",  "+ str(max_val2)+'\n'1110    outf.write(out_char)1111    1112    min_avg_twist_long_pair = float(max_val1 + max_val2)/2.01113    #print min_avg_twist_long_pair1114    out_char = "Avg MinTwist according to longest pair :: " + str(min_avg_twist_long_pair) + '\n'1115    outf.write(out_char)11161117'''111811191120#print (beta_sheet)1121#print (beta_sheet[0][6])1122#print (beta_sheet[0][9])1123#print (atom[0])1124#print (atom[0][5])1125#print (a)1126#f.close()1127
...

twist_region_unitlength_spline.py

Source:twist_region_unitlength_spline.py

1# -*- coding: utf-8 -*-2"""3Created on Fri Oct 14 23:01:08 201645@author: tunaz6"""789import sys1011import re12import numpy13import math1415from scipy import interpolate1617from mpl_toolkits.mplot3d import Axes3D18import matplotlib.pyplot as plt192021class Tee(object):22    def __init__(self, *files):23        self.files = files24    def write(self, obj):25        for f in self.files:26            f.write(obj)27            f.flush() # If you want the output to be visible immediately28    def flush(self) :29        for f in self.files:30            f.flush()31            32            33if len(sys.argv) == 1 or len(sys.argv) == 2 or len(sys.argv) > 3  :34    print "Usage: python twist_region.py pdbfilename.pdb direction"35    sys.exit(0)363738lines_File = [lines1.rstrip('\n') for lines1 in open(sys.argv[1], 'r')]   394041try:42    direction = int(sys.argv[2]) # give direction as input43    '''44    direction= 0 means forward direction, 45    direction = 1 means backword direction 46    '''47    if direction != 0 and direction != 1 :48        print "please input either 0 or 1. 0 = forward, 1 = backword"49        sys.exit(0)50except ValueError:  51        print "please input either 0 or 1. 0 = forward, 1 = backword"52        sys.exit(0)          5354#sys.stdout = open("output_all.pdb", "w")55out_filename = 'twist_angles.txt'56outf = open(out_filename, 'w')5758sys.stdout = open('output_pdb.pdb','w')59#outf2 = open(out_filename2, 'w')60#lines_File = [lines1.rstrip('\n') for lines1 in open('1aky.pdb')]    61#lines_File = [lines1.rstrip('\n') for lines1 in open('1aop.pdb')]  62#lines_File = [lines1.rstrip('\n') for lines1 in open('2p8y.pdb')] #index problem63#lines_File = [lines1.rstrip('\n') for lines1 in open('1s04.pdb')] #complex structure (multple same chain duplicate index)64#lines_File = [lines1.rstrip('\n') for lines1 in open('1atz.pdb')] 65#lines_File = [lines1.rstrip('\n') for lines1 in open('2qtr.pdb')] 66#lines_File = [lines1.rstrip('\n') for lines1 in open('1d5t.pdb')] 67#lines_File = [lines1.rstrip('\n') for lines1 in open('1elu.pdb')]68#lines_File = [lines1.rstrip('\n') for lines1 in open('1a12.pdb')]69#lines_File = [lines1.rstrip('\n') for lines1 in open('1b3a.pdb')]7071#direction = 0 # direction= 0 means forward direction, direction = 1 means backword direction so change here        72            73def draw_interpolation(input_array):74    #print "input  ", input_array75    #print "length  ", len(input_array)76    if  len(input_array) > 2 :77        data = numpy.array(input_array)78        data = data.astype(float)79        80       # print data #ok81           82        data = data[0:100].transpose()83        #print "hi", data #ok84        #now we get all the knots and info about the interpolated spline85		#spacing = int(len(input_array)/.3)86        tck, u= interpolate.splprep(data,k=2)   87        #new = interpolate.splev(numpy.linspace(0,1,50), tck)88        #new = interpolate.splev(numpy.linspace(0,1,int(len(input_array)/.8)), tck) #unit length given 0.189        #new = interpolate.splev(numpy.linspace(0,1,int(len(input_array)/.5)), tck) #we used it for ACM_BCB90        new = interpolate.splev(numpy.linspace(0, 1, int(len(input_array) /.4)), tck) #vector length 1.19192        #now lets plot it!93        94        #fig = plt.figure()95        #ax = Axes3D(fig)96        #ax.plot(data[0], data[1], data[2], label='originalpoints', lw =2, c='Dodgerblue')97        #ax.plot(new[0], new[1], new[2], label='fit', lw =2, c='red')98        #ax.legend()99        #plt.savefig('junk.png')100        #plt.show() 101        return new[0], new[1], new[2]             102103104#print (lines_File1)105106direction_opposite = 1 #hydrogen bond 107108count_beta_sheet = 0109count_beta_strand = 0110counter = 0111counter_index = []112problem = 0113a=[]114beta_sheet = []115atom = []116splitted_line = []117count_atom = 0118for sheet_line in lines_File: #find each line in all lines119     if sheet_line.startswith('SHEET'):     #find those lines start with "SHEET"120         #count_atom = count_atom+1121         #print (sheet_line)       #print the lines started with "SHEET"        122         #beta_sheet.append ( re.sub("\s\s+" , " ",sheet_line).split()) #ok for all except 1123         beta_sheet.append([sheet_line[0:6], sheet_line[7:10], sheet_line[11:14], sheet_line[14:16],sheet_line[17:20], sheet_line[21], sheet_line[22:27], sheet_line[28:31], sheet_line[32], sheet_line[33:37], sheet_line[37], sheet_line[38:40],sheet_line[41:45],sheet_line[45:48], sheet_line[49], sheet_line[50:55], sheet_line[56:60], sheet_line[60:63], sheet_line[64], sheet_line[65:69], sheet_line[69]])124     if sheet_line.startswith('ATOM'):  #find those lines start with "ATOM"125         count_atom = count_atom+1126         #print (sheet_line)       #print the lines started with "ATOM"127         atom.append( [sheet_line[:6], sheet_line[6:11], sheet_line[12:16], sheet_line[17:20], sheet_line[21], sheet_line[22:26], sheet_line[30:38], sheet_line[38:46], sheet_line[46:54]])128         #atom.append( [sheet_line[0:6], sheet_line[6:11], sheet_line[12:16], sheet_line[16:17],sheet_line[17:20], sheet_line[21:22], sheet_line[22:26], sheet_line[26:27], sheet_line[30:38], sheet_line[38:46], sheet_line[46:54], sheet_line[46:54], sheet_line[54:60], sheet_line[60:66], sheet_line[76:78], sheet_line[78:80]])129         #atom.append ( re.sub("\s\s+" , " ",sheet_line).split())130#print (int(atom[1836][5].strip()))131#print (atom[1836])132133strand_vector=[]134center_vector = []135CA_vector_between_strand = []136amino_acid = []137138139#common hydrogen bond140strand_vector_opposite=[]141center_vector_opposite = []142CA_vector_between_strand_opposite = []143amino_acid_opposite = []144145146poi = 0147track = 1148149poi_opposite = 0 #hydrogen bond150151152for sheet_index in range (0,len(beta_sheet)): # handle if there is less than 3 atoms in a strand153    if int(beta_sheet[sheet_index][1]) == int(beta_sheet[sheet_index][3]):154        counter = counter + 1155156#print ("number of beta sheet ",counter) #ok157out_char = "number of beta sheet " + str(counter) + '\n\n'158outf.write(out_char)            159for sheet_index in range (0,len(beta_sheet)): # handle if there is less than 3 atoms in a strand160    if abs((int(beta_sheet[sheet_index][9]))  - int(beta_sheet[sheet_index][6])) >= 2 :161        if int(beta_sheet[sheet_index][1]) == int(beta_sheet[sheet_index][3]):162            #counter = counter + 1163            if track == 0 :164                counter_index.append(int(beta_sheet[sheet_index][1]) - problem)165                track = 1166            else:   167                counter_index.append(int(beta_sheet[sheet_index][1]))168                track = 1169    else:170        #print "hi " #ok171        if  (int(beta_sheet[sheet_index][3]) > 2):           172             problem = problem + 1173             #print problem  #ok174             track = 0175            176            177#for i in range (0,len(counter_index)):178    #print ("counter_index ",counter_index[i]) #ok  179     180for sheet_index in range (0,len(beta_sheet)): # 5 should be number of beta stand in one beta sheet181    #print("love ",beta_sheet[sheet_index][10])182   183    atom_start_index=(int(beta_sheet[sheet_index][6]))184    #print("start index ",atom_start_index) #ok185    atom_end_index=(int(beta_sheet[sheet_index][9])) 186    #print("end index ",atom_end_index) #ok187    sheet_chain = beta_sheet[sheet_index][5]188    #print sheet_chain189190    mid_point = []191    array_strand_index =[]192    amino_acid_array = []193    aa_index = []194    #array_strand_index = [[0 for i in range(2)] for i in range(2)] #2 should be change according to length of one beta sheet195    array_x = []196    array_y = []197    array_z = []198    CA_x = []199    CA_y = []200    CA_z = []201    202    #print (count_atom)203    204    for i in range (0,count_atom) : 205        206        if (atom_start_index <= int(atom[i][5].strip()) <=  atom_end_index)  :207            208            if  atom[i][2].strip() == 'N' or atom[i][2].strip() == 'C': 209                #print "duplicate ", int(atom[i][5].strip())210                211                if sheet_chain == atom[i][4].strip(): #this is needed to differentiate same startindex & endindex but different chains212                    array_strand_index.append(atom[i][5].strip())                  213                    array_x.append(atom[i][6].strip())214                    array_y.append(atom[i][7].strip())215                    array_z.append(atom[i][8].strip())216                    length=len(array_strand_index)  217            if atom[i][2].strip() == 'CA':218                if sheet_chain == atom[i][4].strip():219                    aa_index.append( i )             220                    amino_acid_array.append(atom[i][3].strip())        221                    CA_x.append(atom[i][6].strip())222                    CA_y.append(atom[i][7].strip())223                    CA_z.append(atom[i][8].strip())224                    length_CA = len (CA_x)225    #print ("strand index  ",array_strand_index) #ok226   # print (aa_index)227    #print ("Amino acid  ",amino_acid_array)  #ok        228   229   230   231    232    if direction == 0: #forward direction #ok233        midpoint_x_bef = []234        midpoint_y_bef = []235        midpoint_z_bef = []  236        237      238        239        midpoint_x = []240        midpoint_y = []241        midpoint_z = []242        243        new_x = []244        new_y = []245        new_z = []  246        247       248        #coef = [1,1,1,1,1,1,1,1]249        #coef = []250        251        c_index = 1252        n_index = 2253        while c_index < length-1 and n_index <  length-1 :  #mid point calculation254            #print ("c_index ",c_index, "n_index ", n_index)255            midpoint_x_bef.append((float(array_x[c_index]) + float(array_x[n_index]))/float(2.0))    #making midpoint       256            midpoint_y_bef.append((float(array_y[c_index]) + float(array_y[n_index]))/float(2.0))   #making midpoint 257            midpoint_z_bef.append((float(array_z[c_index]) + float(array_z[n_index]))/float(2.0))  #making midpoint 258           259            c_index = c_index + 2260            n_index = n_index +2  261        262        #print ("midpoint_x_bef  ", midpoint_x_bef ) 263        #print ("midpoint_y_bef  ", midpoint_y_bef  )264        #print ("midpoint_z_bef  ", midpoint_z_bef  )265            266        midpoint_draw = []267        draw_index = 0268        while draw_index < len (midpoint_x_bef) :269            midpoint_point = []270            mid_x = midpoint_x_bef [draw_index]271            mid_y = midpoint_y_bef [draw_index]272            mid_z = midpoint_z_bef [draw_index]           273            midpoint_point.append (mid_x)274            midpoint_point.append (mid_y)275            midpoint_point.append (mid_z)276            midpoint_draw.append (midpoint_point) 277            draw_index = draw_index + 1278        #print   (midpoint_draw)  #ok279        pdb_format = '6s5s1s4s1s3s1s1s4s1s3s8s8s6s10s2s3s'280        #shape = 'ATOM 1 N LEU A 81 25.865 39.612 19.376 1.00 53.28 N '281        shape = ['ATOM', '1','N', 'LEU', 'A', '81', '25.865', '39.612', '19.376']282        if len(midpoint_draw) <= 2:283             midpoint_x=midpoint_x_bef284             midpoint_y=midpoint_y_bef285             midpoint_z=midpoint_z_bef286             287             for i in range (len( midpoint_x)):288                poi = poi +1289                                           290                #print "%-6s%5s %4s %3s %s%4s    %8s%8s%8s\n"%tuple(shape) ok291                print "%-6s%5s %4s %3s %s%4s    %8s%8s%8s"%tuple(['ATOM',poi,'S','LEU', 'A', '81',round(midpoint_x[i],3),round(midpoint_y[i],3),round(midpoint_z[i],3)]) #print new atom position in PDB file292                293             294        else:    295            value = draw_interpolation (midpoint_draw)   #making 3D cubic spline interpolation for midpoints296            #print value[0] #ok297            midpoint_x=value[0]  #new midpoint after cubic spline298            midpoint_y=value[1]  #new midpoint after cubic spline299            midpoint_z=value[2]   #new midpoint after cubic spline300        301            302           303            pdb_length = len( midpoint_x)304          305            for i in range (pdb_length):306                poi = poi +1          307                print "%-6s%5s %4s %3s %s%4s    %8s%8s%8s"%tuple(['ATOM',poi,'S','LEU', 'A', '81',round(midpoint_x[i],3),round(midpoint_y[i],3),round(midpoint_z[i],3)]) #print new atom position in PDB file308                309       310       311        312        313       314        #print ("midpoint_x ", midpoint_x ) 315         316        317     318     319        length_midpoint_array = len(midpoint_x)320        vector_array =[]321        cen_array = []322        midpoint_index = 1323        while midpoint_index < length_midpoint_array:324            point = []325            cen_point = []326            temp_mid_index = midpoint_index - 1327            vec_x = midpoint_x [temp_mid_index] -  midpoint_x [temp_mid_index+1]  #making midpoint vector328            vec_y = midpoint_y [temp_mid_index] -  midpoint_y [temp_mid_index+1]  #making midpoint vector329            vec_z = midpoint_z [temp_mid_index] -  midpoint_z [temp_mid_index+1]  #making midpoint vector330            331            cen_x = (float(midpoint_x [temp_mid_index]) + float( midpoint_x [temp_mid_index+1])) / float(2.0)  #making center of midpoint vector332            cen_y =(float( midpoint_y [temp_mid_index]) + float(midpoint_y [temp_mid_index+1]))/float(2.0)   #making center of midpoint vector333            cen_z = (float(midpoint_z [temp_mid_index]) + float( midpoint_z [temp_mid_index+1]))/float(2.0)   #making center of midpoint vector334            335            #print ("vector x ",vec_x) 336           # print ("vector y ",vec_y) 337           # print ("vector z ",vec_z) 338            #make unit vector339            magnitude = math.sqrt(pow(vec_x, 2) + pow(vec_y, 2) + pow(vec_z, 2))340            341           342            print("vector length ",magnitude)343            344            point.append(vec_x/magnitude) # unit vector345            point.append(vec_y/magnitude)  # unit vector346            point.append(vec_z/magnitude)  # unit vector347            348      349            #point.append(vec_x) #ok350            #point.append(vec_y)351            #point.append(vec_z)352      353         354            cen_point.append(cen_x) #ok355            cen_point.append(cen_y)356            cen_point.append(cen_z)357            358            359            360            361            vector_array.append(point)362            363            cen_array.append(cen_point)364            365            midpoint_index = midpoint_index +1366        #print "vector array ",vector_array367        368        369        aa = []370        aa_index = 1371        while aa_index < len(amino_acid_array) - 1:   # to get amino acid value372             aa.append(amino_acid_array[aa_index]) 373             aa_index = aa_index +1374        amino_acid.append(aa)375        376        #print ("x cordinate of CA in each strand" CA_x )377    378        #print("length ", length)379        #print ("CA length ",length_CA)380        CA_x_new = []381        CA_y_new = []382        CA_z_new = []383        CA_vector = []384        CA_index = 1385        while CA_index < length_CA - 1:386            CA_point = []387            CA_x_new=(float(CA_x[CA_index]))388            CA_y_new=(float(CA_y[CA_index]))389            CA_z_new=(float(CA_z[CA_index]))390            CA_point.append (CA_x_new)391            CA_point.append (CA_y_new)392            CA_point.append (CA_z_new)393            CA_vector.append (CA_point)   394            CA_index = CA_index + 1395        #draw_interpolation (CA_vector) #draw spline according to carbon alpha 396397        CA_vector_between_strand.append(CA_vector)398        #draw_interpolation (CA_vector_between_strand)399400         401    if direction_opposite == 1:   #backword direction402    403        midpoint_x_bef_opposite = []404        midpoint_y_bef_opposite = []405        midpoint_z_bef_opposite = []  406407        midpoint_x_opposite = []408        midpoint_y_opposite = []409        midpoint_z_opposite = [] 410        411       412        413        c_index_opposite = length-3414        n_index_opposite = length-2415        while c_index_opposite > 0  and n_index_opposite > 0 : 416            #print ("c_index ",c_index, "n_index ", n_index)417            midpoint_x_bef_opposite.append((float(array_x[c_index_opposite]) + float(array_x[n_index_opposite]))/float(2.0))418            midpoint_y_bef_opposite.append((float(array_y[c_index_opposite]) + float(array_y[n_index_opposite]))/float(2.0))419            midpoint_z_bef_opposite.append((float(array_z[c_index_opposite]) + float(array_z[n_index_opposite]))/float(2.0)) 420            c_index_opposite = c_index_opposite - 2421            n_index_opposite = n_index_opposite - 2   422       423       424       425        midpoint_draw_opposite = []426        draw_index_opposite = len (midpoint_x_bef_opposite) -1427        while draw_index_opposite >= 0 :428            midpoint_point_opposite = []429            mid_x_opposite = midpoint_x_bef_opposite [draw_index_opposite]430            mid_y_opposite = midpoint_y_bef_opposite [draw_index_opposite]431            mid_z_opposite = midpoint_z_bef_opposite [draw_index_opposite]           432            midpoint_point_opposite.append (mid_x_opposite)433            midpoint_point_opposite.append (mid_y_opposite)434            midpoint_point_opposite.append (mid_z_opposite)435            midpoint_draw_opposite.append (midpoint_point_opposite) 436            draw_index_opposite = draw_index_opposite - 1437            438        #draw_interpolation (midpoint_draw)439       440        if len(midpoint_draw_opposite) <= 2:441             midpoint_x_opposite = midpoint_x_bef_opposite 442             midpoint_y_opposite = midpoint_y_bef_opposite443             midpoint_z_opposite = midpoint_z_bef_opposite444             445             for i in range (len( midpoint_x_opposite)):446                poi_opposite = poi_opposite +1447                                           448                #print "%-6s%5s %4s %3s %s%4s    %8s%8s%8s\n"%tuple(shape) ok449               # print "%-6s%5s %4s %3s %s%4s    %8s%8s%8s"%tuple(['ATOM',poi_opposite,'S','LEU', 'A', '81',round(midpoint_x_opposite[i],3),round(midpoint_y_opposite[i],3),round(midpoint_z_opposite[i],3)]) #print new atom position in PDB file450                451             452        else:    453            value_opposite = draw_interpolation (midpoint_draw_opposite)   #making 3D cubic spline interpolation for midpoints454            #print value_opposite[0] #ok455            midpoint_x_opposite = value_opposite[0]  #new midpoint after cubic spline456            midpoint_y_opposite = value_opposite[1]  #new midpoint after cubic spline457            midpoint_z_opposite = value_opposite[2]   #new midpoint after cubic spline458        459            460           461            pdb_length_opposite = len( midpoint_x_opposite)462          463            for i in range (pdb_length_opposite):464                poi_opposite = poi_opposite +1          465                #print "%-6s%5s %4s %3s %s%4s    %8s%8s%8s"%tuple(['ATOM',poi_opposite,'S','LEU', 'A', '81',round(midpoint_x_opposite[i],3),round(midpoint_y_opposite[i],3),round(midpoint_z_opposite[i],3)]) #print new atom position in PDB file466                467       468       469        midpoint_x_opposite = midpoint_x_bef_opposite470        midpoint_y_opposite = midpoint_y_bef_opposite471        midpoint_z_opposite = midpoint_z_bef_opposite472        473    474        #print ("midpoint_x ", midpoint_x )      475        length_midpoint_array_opposite = len (midpoint_x_opposite)   476        vector_array_opposite =[]477        cen_array_opposite = []478        479        #midpoint_index_opposite = length_midpoint_array480        #while midpoint_index > 1:481        482        midpoint_index_opposite = 1 #big change in this loop483        while midpoint_index_opposite < length_midpoint_array_opposite:484            point_opposite = []485            cen_point_opposite = []486            487            temp_mid_index_opposite = midpoint_index_opposite - 1488            #vec_x = midpoint_x [temp_mid_index] -  midpoint_x [temp_mid_index-1]489            #vec_y = midpoint_y [temp_mid_index] -  midpoint_y [temp_mid_index-1]490            #vec_z = midpoint_z [temp_mid_index] -  midpoint_z [temp_mid_index-1]491            492            vec_x_opposite = midpoint_x_opposite [temp_mid_index_opposite] -  midpoint_x_opposite [temp_mid_index_opposite+1]493            vec_y_opposite = midpoint_y_opposite [temp_mid_index_opposite] -  midpoint_y_opposite [temp_mid_index_opposite+1]494            vec_z_opposite = midpoint_z_opposite [temp_mid_index_opposite] -  midpoint_z_opposite [temp_mid_index_opposite+1]495            496            #print ("vector ",vec_x,vec_y,vec_z) 497            498            #cen_x = (float(midpoint_x [temp_mid_index]) + float( midpoint_x [temp_mid_index-1])) / float(2.0)  #making center of midpoint vector499            #cen_y =(float( midpoint_y [temp_mid_index]) + float(midpoint_y [temp_mid_index-1]))/float(2.0)   #making center of midpoint vector500            #cen_z = (float(midpoint_z [temp_mid_index]) + float( midpoint_z [temp_mid_index-1]))/float(2.0)   #making center of midpoint vector501            502            cen_x_opposite = (float(midpoint_x_opposite [temp_mid_index_opposite]) + float( midpoint_x_opposite [temp_mid_index_opposite+1])) / float(2.0)  #making center of midpoint vector503            cen_y_opposite =(float( midpoint_y_opposite [temp_mid_index_opposite]) + float(midpoint_y_opposite [temp_mid_index_opposite+1]))/float(2.0)   #making center of midpoint vector504            cen_z_opposite = (float(midpoint_z_opposite [temp_mid_index_opposite]) + float( midpoint_z_opposite [temp_mid_index_opposite+1]))/float(2.0)   #making center of midpoint vector505            506              #make unit vector507            magnitude_opposite = math.sqrt(pow(vec_x_opposite, 2) + pow(vec_y_opposite, 2) + pow(vec_z_opposite, 2))508            #print("vector length opposite ",magnitude_opposite)509            #point.append(vec_x/magnitude) #not ok510            #point.append(vec_y/magnitude)511            #point.append(vec_z/magnitude)512            513            '''514            point_opposite.append(vec_x_opposite) #ok515            point_opposite.append(vec_y_opposite)516            point_opposite.append(vec_z_opposite)517            '''518            point_opposite.append(vec_x_opposite/magnitude_opposite) #ok519            point_opposite.append(vec_y_opposite/magnitude_opposite)520            point_opposite.append(vec_z_opposite/magnitude_opposite)521            522            523            524            cen_point_opposite.append(cen_x_opposite) #ok525            cen_point_opposite.append(cen_y_opposite)526            cen_point_opposite.append(cen_z_opposite)527     528            vector_array_opposite.append(point_opposite)529            530            cen_array_opposite.append(cen_point_opposite)531            532            midpoint_index_opposite = midpoint_index_opposite + 1533            534            #midpoint_index = midpoint_index - 1 535        536        aa_opposite = []537        aa_index_opposite = len(amino_acid_array) - 2538        while aa_index_opposite > 0:   # to get amino acid value539             aa_opposite.append(amino_acid_array[aa_index_opposite]) 540             aa_index_opposite = aa_index_opposite - 1 541        amino_acid_opposite.append(aa_opposite)542   543   544        #print ("x cordinate of CA in each strand", CA_x )545    546        #print("length ", length)547        #print ("CA length ",length_CA)548        CA_x_new_opposite = []549        CA_y_new_opposite = []550        CA_z_new_opposite = []551        CA_vector_opposite = []552        CA_index_opposite = length_CA - 2553        while CA_index_opposite > 0 :554            CA_point_opposite = []555            CA_x_new_opposite=(float(CA_x[CA_index_opposite]))556            CA_y_new_opposite=(float(CA_y[CA_index_opposite]))557            CA_z_new_opposite=(float(CA_z[CA_index_opposite]))558            CA_point_opposite.append (CA_x_new_opposite)559            CA_point_opposite.append (CA_y_new_opposite)560            CA_point_opposite.append (CA_z_new_opposite)561            CA_vector_opposite.append (CA_point_opposite)   562            CA_index_opposite = CA_index_opposite - 1563        CA_vector_between_strand_opposite.append(CA_vector_opposite)564        565        566    strand_vector.append(vector_array)  567    568    center_vector.append(cen_array)569    570    #for hydrogen bond calculation571    strand_vector_opposite.append(vector_array_opposite)  572    573    center_vector_opposite.append(cen_array_opposite)574    575   #print(beta_sheet[sheet_index][1])  576    577    578579 580   581    582    '''583    print ("vector ",vector_array)584    print (atom_start_index)585    print (array_strand_index) 586    print(length)587    588    print (array_x) 589    print (array_y) 590    print (array_z) 591    print (midpoint_x)592    print (midpoint_y)593    print (midpoint_z)594    print (length_midpoint_array)595    '''596print "\n"597length_strand_vector = len(strand_vector)598#print (strand_vector)599#print(amino_acid) # ok600#print ("vector strand length ",length_strand_vector)601#print ("CA_vector_between_strand ",len(CA_vector_between_strand)) # ok.  vector_strand length and CA_vector same602#print ("CA_vector_between_strand ",CA_vector_between_strand) 603#draw_interpolation (CA_vector_between_strand)604605def distance(x,y):   606    return numpy.sqrt(numpy.sum((x-y)**2))607608609def find_nearest_distance(a,B, pair_aa1, pair_aa2):610    611    612    #print ("a ", a)613    min_point = B[0]614    #min_dist = abs(numpy.linalg.norm(a- numpy.array(B[0])))615    min_dist = math.sqrt ( pow((a[0]-B[0][0]),2) + pow((a[1]-B[0][1]),2) + pow((a[2]-B[0][2]),2)  )616    #min_dist = 12000617   # print ("mindistance ", min_dist)618    #print("len of B",len(B)) #ok619    aa_min_index = 0620    for i in range (0, len(B)): 621        #dist = []622                       623        b = numpy.array(B[i])        624        #print ("b ", b) #ok625        #dist=(numpy.linalg.norm(a-b))626        try_dist = math.sqrt ( pow((a[0]-b[0]),2) + pow((a[1]-b[1]),2) + pow((a[2]-b[2]),2)  )627        #print ("distance ", try_dist)628        #dist=abs(distance(a,b))629        if try_dist < min_dist:630            min_dist = try_dist631            min_point = b632            aa_min_index = i633 634    635    return (min_point, aa_min_index)636637638def calculate_twist_angle(a,minimum_dist_point,j,aa_min_index2):639        640    #print ("mindist point ", minimum_dist_point)641    dot_product = numpy.dot(a,minimum_dist_point) # A.B642    #print (dot_product)643       644    a_sum = pow(a[0],2)+pow(a[1],2)+pow(a[2],2)645    a_root=math.sqrt(a_sum)646    b_sum = pow(minimum_dist_point[0],2)+pow(minimum_dist_point[1],2)+pow(minimum_dist_point[2],2)647    b_root=math.sqrt(b_sum)648    a_b = a_root*b_root649    650    angle_dot = math.acos(dot_product/a_b)651    #print (angle_dot)652    twist_angle = math.degrees(angle_dot)653    654    '''655    ### cross product ###656    cross_product = numpy.cross(a,minimum_dist_point) # A*B657    #print (cross_product)658    sum_square_cross_product = pow(cross_product[0],2) +  pow(cross_product[1],2)+ pow(cross_product[2],2)659    #print (sum_square_cross_product)660    root_sum_square_cross_product = math.sqrt(sum_square_cross_product)661    #print (root_sum_square_cross_product)662   663   664   665    angle_cross = math.asin(root_sum_square_cross_product/a_b)666    #print (angle_cross)667    print(math.degrees(angle_cross))668    ######## end of cross product ###########669    '''670    if twist_angle > 90.:671        twist_angle = 180.0 - twist_angle672   # sys.stdout = original673   # print(" Amino Acid",pair_aa1[j], " and ", pair_aa2[aa_min_index2],":", "Twist angle::" , twist_angle) #necessary674    #print("Twist angle::" , round(twist_angle,3)) #necessary675    out_char = "Twist angle:: " + str(round(twist_angle,3)) + '\n'676    outf.write(out_char)677    return twist_angle678       # print ( "amino acid ", pair_aa1[j], " ", pair_aa2[aa_min_index2])679minimum_dist_point = []  680minimum_dist_point_dup = []681flag = 0682all_avg_twist_angle = []  683avg_count = 0    684avg_track = 1685pair_array = []686len_array = []687i = 0688while i < length_strand_vector : # pairwise calculation 689    #print ("beta sheet index", beta_sheet[i][1])690    #print ("beta strand", beta_sheet[i])691    692    #print(int(beta_sheet[i][1]), " : ", int(beta_sheet[i][3]))693    #print ("i ", i)694    #print ("flag ",flag)695    avg_twist_count = []696            697    if(int(beta_sheet[i][1])>=int(beta_sheet[i][3])):698        flag = int(beta_sheet[i][1])699        avg_count = 0700       701        #print ("avg count ",avg_count)702    else:703        flag = flag + 1704        avg_count = avg_count + 1705        A_t = []706        B_t= []707        A = []708        B = []709        710        #for hydrogen bond711        A_t_opposite = []712        B_t_opposite  = []713        A_opposite  = []714        B_opposite  = []715        716        alpha_A_t = []717        alpha_B_t = []718        alpha_A = []719        alpha_B = []720        721          #for hydrogen bond722        alpha_A_t_opposite = []723        alpha_B_t_opposite = []724        alpha_A_opposite = []725        alpha_B_opposite = []726        727        pair_aa1_t = []728        pair_aa2_t = []729        pair_aa1 = []730        pair_aa2 = []731        732         733        #for hydrogen bond734        pair_aa1_t_opposite = []735        pair_aa2_t_opposite = []736        pair_aa1_opposite = []737        pair_aa2_opposite = []738        739        740        index_vec = i741        A_t = strand_vector[i] #copy strand to A array742        B_t = strand_vector[index_vec+1] #copy strand to B array 743        744        #for hydrogen bond745        A_t_opposite = strand_vector_opposite[i] #copy strand to A array746        B_t_opposite = strand_vector_opposite[index_vec+1] #copy strand to B array 747        748        749        #alpha_A_t = CA_vector_between_strand[i] #copy strand to A array750        #alpha_B_t = CA_vector_between_strand[index_vec+1] #copy strand to B array 751        752         #for hydrogen bond753        #alpha_A_t_opposite = CA_vector_between_strand_opposite[i] #copy strand to A array (alpha-carbon) backword direction754        #alpha_B_t_opposite = CA_vector_between_strand_opposite[index_vec+1] #copy strand to B array (alpha-carbon) backword direction755        756        757        alpha_A_t = center_vector[i] #copy strand to A array758        alpha_B_t = center_vector[index_vec+1] #copy strand to B array 759        760         #for hydrogen bond761        alpha_A_t_opposite = center_vector_opposite[i] #copy strand to A array (alpha-carbon) backword direction762        alpha_B_t_opposite = center_vector_opposite[index_vec+1] #copy strand to B array (alpha-carbon) backword direction763        764        765        pair_aa1_t = amino_acid[i] #copy strand to A array766        pair_aa2_t = amino_acid[index_vec+1] #copy strand to B array 767        768        #hydrogen bond769        pair_aa1_t_opposite = amino_acid_opposite[i] #copy strand to A array770        pair_aa2_t_opposite = amino_acid_opposite[index_vec+1] #copy strand to B array 771        772        773      #  print ("A_t ",A_t)774       775      #  print ("B_t ",B_t)776       # sys.stdout = original777        #print ("beta strand 1", beta_sheet[i]) #necessary778        #print ("beta strand 2", beta_sheet[i+1]) #necessary779        out_char = "beta strand 1" + str(beta_sheet[i]) + '\n'780        outf.write(out_char)781        out_char = "beta strand 2" + str(beta_sheet[i+1]) + '\n'782        outf.write(out_char)783        #print ("A_t ",len(A_t))784       785        #print ("B_t ",len(B_t))786        if len(A_t) > len(B_t):787            A = B_t788            B = A_t789            790            #hydrogen bond791            A_opposite = B_t_opposite792            B_opposite = A_t_opposite793            794            alpha_A = alpha_B_t795            alpha_B = alpha_A_t796            797             #hydrogen bond798            alpha_A_opposite = alpha_B_t_opposite799            alpha_B_opposite = alpha_A_t_opposite800            801            pair_aa1 = pair_aa2_t802            pair_aa2 = pair_aa1_t803            804            #hydrogen bond805            pair_aa1_opposite = pair_aa2_t_opposite806            pair_aa2_opposite = pair_aa1_t_opposite807        else:808            A = A_t809            B = B_t810            811             #hydrogen bond812            A_opposite = A_t_opposite813            B_opposite = B_t_opposite814            815            alpha_A = alpha_A_t816            alpha_B = alpha_B_t 817            818            #hydrogen bond819            alpha_A_opposite = alpha_A_t_opposite820            alpha_B_opposite = alpha_B_t_opposite821            822            pair_aa1 = pair_aa1_t823            pair_aa2 = pair_aa2_t824            825             #hydrogen bond826            pair_aa1_opposite = pair_aa1_t_opposite827            pair_aa2_opposite = pair_aa2_t_opposite828        829        #aa_min_index2 = 0830        k = 0831        twist_angle_return = 0.0832        add_1 = 0 833        add_n = 0834        add_hp = 0 #hydrogen bond835        pick_index = 0836        837        '''for calculation of avg twist angle between two beta strands'''838        twist_count = 0839        sum_twist_angle = 0.0840        avg_twist_angle = []841        842        #taking minimum of 4 consecutive  angles' sum843        cons_twist_angle = []844       845        846        for j in range (0, len(A_opposite)):847            a = numpy.array(A[j])  848            849            #hydrogen bond850            a_opposite = numpy.array(A_opposite[j])851            852            alpha_a = numpy.array(alpha_A[j])  #this is needed for calculating center of vector in 3D space853            854            #hydrogen bond855            alpha_a_opposite = numpy.array(alpha_A_opposite[j])856           857               858            #minimum_dist_point, aa_min_index2 = find_nearest_distance(a, B, pair_aa1[j], pair_aa2)859            #print ("k ", k)860            if k == 0  : # determining first pair861              862                minimum_dist_point, aa_min_index2 = find_nearest_distance(alpha_a, alpha_B, pair_aa1[j], pair_aa2)#nearest neighbour checking it will work for determining first pair863                #minimum_dist_point, aa_min_index2 = find_nearest_distance(a, B, pair_aa1[j], pair_aa2)864               #hydrogen bond865                minimum_dist_point_opposite, aa_min_index2_opposite = find_nearest_distance(alpha_a_opposite, alpha_B_opposite, pair_aa1_opposite[j], pair_aa2_opposite)#nearest neighbour checking it will work for determining first pair866                #minimum_dist_point_opposite, aa_min_index2_opposite = find_nearest_distance(a_opposite, B_opposite, pair_aa1_opposite[j], pair_aa2_opposite)867               868               # sys.stdout = original869                print(" MIN_INDEX_Forward", aa_min_index2)   870                print(" MIN_INDEX_Backword", aa_min_index2_opposite)    871                #print(" MIN_POINT ", minimum_dist_point)872                #print(" MIN_B ", B[aa_min_index2])873                874                '''875                twist_angle_return = calculate_twist_angle(a,B[aa_min_index2],j,aa_min_index2)876                sum_twist_angle = sum_twist_angle + twist_angle_return877                twist_count = twist_count + 1             878                pick_index = aa_min_index2879                '''880                881                if aa_min_index2 == 0 and len(A) <= aa_min_index2_opposite: #hydrogen bond parallel and anti-parallel case equal length strand or small strand is inside of big strand (complete overlap)882                    add_1 = 1883                    884                    twist_angle_return = calculate_twist_angle(a,B[aa_min_index2],j,aa_min_index2)885                    cons_twist_angle.append(twist_angle_return)   886                    sum_twist_angle = sum_twist_angle + twist_angle_return887                    twist_count = twist_count + 1             888                    pick_index = aa_min_index2889                    890                if aa_min_index2 == 0  and len(A) > aa_min_index2_opposite: #hydrogen bond parallel and anti-parallel case   small strand is outside (partial overlap) of big strand891                    add_hp = 1892                    893                    twist_angle_return = calculate_twist_angle(a_opposite,B_opposite[aa_min_index2_opposite],j,aa_min_index2)894                    cons_twist_angle.append(twist_angle_return)   895                    sum_twist_angle = sum_twist_angle + twist_angle_return896                    twist_count = twist_count + 1             897                    pick_index = aa_min_index2_opposite898                    899                if aa_min_index2 != 0:900                    add_n = 1901                    twist_angle_return = calculate_twist_angle(a,B[aa_min_index2],j,aa_min_index2)902                    cons_twist_angle.append(twist_angle_return)   903                    sum_twist_angle = sum_twist_angle + twist_angle_return904                    twist_count = twist_count + 1             905                    pick_index = aa_min_index2906                907            else: #next pairs908                print aa_min_index2," ",len(B), " ",aa_min_index2_opposite909                if add_1 == 1 and k < len(B):910                    minimum_dist_point_dup = B[k]911                    pick_index = k912                    twist_angle_return = calculate_twist_angle(a,minimum_dist_point_dup,j,pick_index)913                    cons_twist_angle.append(twist_angle_return)   914                    cons_twist_angle.append(twist_angle_return)   915                    sum_twist_angle = sum_twist_angle + twist_angle_return916                    twist_count = twist_count + 1917                    918                if add_hp == 1 and  int(beta_sheet[i+1][11])!=-1  and  pick_index < len(B_opposite) - 1:		#parallel  small strand is outside (partial overlap) of big strand                         919                   # j = pick_index - 1920                    pick_index = pick_index + 1921                  922                    print "k ",k923                    print "pick_index",pick_index924                    minimum_dist_point_dup = B_opposite[pick_index]925                    #print "minimum_dist_point_dup ",minimum_dist_point_dup #ok926                   # print "a_opposite ",a_opposite927            928                    twist_angle_return = calculate_twist_angle(a_opposite,minimum_dist_point_dup,j,pick_index)929                    cons_twist_angle.append(twist_angle_return)   930                    sum_twist_angle = sum_twist_angle + twist_angle_return931                    twist_count = twist_count + 1932                    933                if add_hp == 1 and  int(beta_sheet[i+1][11])==-1  and  pick_index != 0: # anti-parallel  small strand is outside (partial overlap) of big strand                         934                   # j = pick_index - 1935                    pick_index = pick_index - 1 #may be need to change later 936                  937                    print "k ",k938                    print "pick_index",pick_index939                    minimum_dist_point_dup = B_opposite[pick_index]940                    #print "minimum_dist_point_dup ",minimum_dist_point_dup #ok941                   # print "a_opposite ",a_opposite942                    943                   944                    twist_angle_return = calculate_twist_angle(a_opposite,minimum_dist_point_dup,j,pick_index)945                    cons_twist_angle.append(twist_angle_return)   946                    sum_twist_angle = sum_twist_angle + twist_angle_return947                    twist_count = twist_count + 1948                949                950                if add_n == 1 and  int(beta_sheet[i+1][11])!=-1 and k < len(B) - aa_min_index2: # Parallel beta strands 	951                    pick_index = pick_index + 1952                    minimum_dist_point_dup = B[pick_index]953                    twist_angle_return = calculate_twist_angle(a,minimum_dist_point_dup,j,pick_index)954                    cons_twist_angle.append(twist_angle_return)   955                    sum_twist_angle = sum_twist_angle + twist_angle_return956                    twist_count = twist_count + 1957                    958                if add_n == 1 and   int(beta_sheet[i+1][11])==-1  and pick_index!=0: # anti-Parallel beta strands 		959                    pick_index = pick_index - 1           960                    minimum_dist_point_dup = B[pick_index]961                    twist_angle_return = calculate_twist_angle(a,minimum_dist_point_dup,j,pick_index)962                    cons_twist_angle.append(twist_angle_return)   963                    sum_twist_angle = sum_twist_angle + twist_angle_return964                    twist_count = twist_count + 1965            966                                                                   967            k = k + 1968                 969        if  twist_count != 0:  970            avg_twist_angle.append (float(sum_twist_angle) / float(twist_count))  971           972            cons_sum_array = []973            if len(cons_twist_angle) > 4:974                for c in range (0,len(cons_twist_angle) - 3):975                    cons_sum = 0.0976                    for cc in range (c,c+4):       977                         cons_sum = cons_twist_angle[cc] + cons_sum978                    cons_sum_array.append(cons_sum)979                pair_array.append((round(min(cons_sum_array)  / 4.0,3)))980                len_array.append(len(cons_sum_array))981                #print "sum of consecutive 4 angles ", cons_sum_array982                #print "min of consecutive 4 angles ",min(cons_sum_array)983                out_char = "min of sum of consecutive 4 angles :: " + str(round(min(cons_sum_array),3)) + '\n'984                outf.write(out_char)985                out_char = "Twist per angle :: " + str(round(min(cons_sum_array)  / 4.0 ,3)) + '\n'986                outf.write(out_char)987            988            elif len(cons_twist_angle) == 4:989                for c in range (0,len(cons_twist_angle) - 2):990                    cons_sum = 0.0991                    for cc in range (c,c+3):       992                         cons_sum = cons_twist_angle[cc] + cons_sum993                    cons_sum_array.append(cons_sum)994                pair_array.append((round(min(cons_sum_array)  / 3.0,3)))995                len_array.append(len(cons_sum_array))996                #print "sum of consecutive 4 angles ", cons_sum_array997                #print "min of consecutive 4 angles ",min(cons_sum_array)998                out_char = "min of sum of consecutive 3 angles :: " + str(round(min(cons_sum_array),3)) + '\n'999                outf.write(out_char)  1000                out_char = "Twist per angle :: " + str(round(min(cons_sum_array) / 3.0 ,3) ) + '\n'1001                outf.write(out_char)1002                1003            elif len(cons_twist_angle) == 2 or len(cons_twist_angle) == 3 or len(cons_twist_angle) == 1:1004               #print "smallest twist ",min(cons_twist_angle)1005                pair_array.append((round(min(cons_twist_angle) , 3 )))1006                len_array.append(len(cons_sum_array))1007                out_char = "smallest angles :: " + str(round(min(cons_twist_angle),3)) + '\n'1008                outf.write(out_char)1009            avg_track = 11010        else:1011            avg_track = 01012        1013        for t in range (len(avg_twist_angle)):   1014            #print ("avg_twist_angle between strands ", round(avg_twist_angle[t],3)) #necessary1015            #print ( round(avg_twist_angle[t],3)) #necessary1016            out_char = "avg_twist_angle between strands  " + str( round(avg_twist_angle[t],3) ) + '\n'1017            outf.write(out_char)1018            all_avg_twist_angle.append(avg_twist_angle[t])            1019     1020        #print ("...........................")#necessary1021        #print cons_twist_angle#ok1022      1023        outf.write("................................................"+"\n")1024   1025    i = i +11026#print ("avg_twist_angle between strands ", all_avg_twist_angle) #ok     102710281029'''for calculation of avg twist angle in a beta sheet'''1030outf.write("\n\n")10311032twist_angle_beta_sheet = []1033avg_avg_twist_angle_beta_sheet = []1034jc = 0   1035for i in range (len(counter_index)):  1036    sum_avg_twist_angle = 0.0  1037    j = jc1038    #print jc, " ", counter_index[i]-1, " ",j #ok1039    while j < (jc + counter_index[i]-1):1040        #print "j ",j1041        sum_avg_twist_angle = sum_avg_twist_angle + all_avg_twist_angle[j]  1042        sum_avg_twist_angle = round(sum_avg_twist_angle,3)       1043        j = j + 1    1044    jc = jc + counter_index[i]-11045    #avg_avg_twist_angle_beta_sheet.append(sum_avg_twist_angle)1046    twist_angle_beta_sheet.append(float(sum_avg_twist_angle)/float(counter_index[i]-1))1047    #print ("avg_twist angle in all beta sheets ", round(twist_angle_beta_sheet[i],3)) #necessary1048  1049    #print ( round(twist_angle_beta_sheet[i],3)) 1050    out_char = "avg_twist angle in all beta sheets " + str( round(twist_angle_beta_sheet[i],3)) + '\n'1051    outf.write(out_char)1052   1053#print ("avg_twist angle in all beta sheets ", float(avg_avg_twist_angle_beta_sheet[0])/float(counter_index[0]-1))1054#print pair_array #ok10551056outf.write("................................................"+"\n")10571058#print pair_array # ok1059#print len_array #ok1060'''1061#### According to minimum value pair #############10621063kc = 01064for i in range (len(counter_index)):  1065    j = kc1066    min_arr = []1067    while j < (kc + counter_index[i]-1):1068        min_arr.append (pair_array[j])1069        1070        j = j + 1  1071    kc = kc + counter_index[i]-11072    arr = numpy.array (min_arr) 1073    ind1 , ind2 = arr.argsort()[:2] 1074    min_val1 =   arr[ind1]1075    min_val2 =   arr[ind2]1076    #print min_val1, min_val21077    out_char = "Min pair val1, val2 :: " + str(min_val1) + ",  "+ str(min_val2)+'\n'1078    outf.write(out_char)1079    1080    min_avg_twist = float(min_val1 + min_val2)/2.01081    #print min_avg_twist1082    out_char = "Min Twist :: " + str(min_avg_twist) + '\n'1083    outf.write(out_char)10841085outf.write("................................................"+"\n")10861087#### According to longest pair #############1088kl = 01089for il in range (len(counter_index)):  1090    jl = kl1091    min_arr_l = []1092    while jl < (kl + counter_index[il]-1):1093        min_arr_l.append (len_array[jl])1094        1095        jl = jl + 1  1096    kl = kl + counter_index[il]-11097    arr_len = numpy.array (len_array) 1098    ind1 , ind2 = arr_len.argsort()[::-1][:2]1099    max_val1 =   pair_array[ind1]1100    max_val2 =   pair_array[ind2]1101   # print ind1, ind2 #ok1102    out_char = "longest pair val1, val2 :: " + str(max_val1) + ",  "+ str(max_val2)+'\n'1103    outf.write(out_char)1104    1105    min_avg_twist_long_pair = float(max_val1 + max_val2)/2.01106    #print min_avg_twist_long_pair1107    out_char = "Avg MinTwist according to longest pair :: " + str(min_avg_twist_long_pair) + '\n'1108    outf.write(out_char)11091110'''111111121113#print (beta_sheet)1114#print (beta_sheet[0][6])1115#print (beta_sheet[0][9])1116#print (atom[0])1117#print (atom[0][5])1118#print (a)1119#f.close()1120
...

merge_sample.py

Source:merge_sample.py

1from extractor.models import *2def compare_stretch_with_confusions(first,second, all_indexes):3    inside_if=04    '''check if my is inside of opposite5    check if opposite is inside of my6    check of my and opposite contain each other for at least 50%7    '''8    # SELECT AVG(percentage_based_on_chars) from extractor_compatibility;9    for my in first:10        overlap = []11        for opposite in second:12            state = 0 # 1 = first if is true, 2 = second if is true13            if my[1]>=opposite[1] and my[2]<=opposite[2]:14                state = 115            if my[2]>=opposite[1] and my[1]<=opposite[1] and (my[2]-opposite[1])>=((my[2]-my[1])/2):16                state = 217            if state==1 or state==2:18                if my[3]==opposite[3] \19                        or (my[3]==1 and opposite[3]==4) or (my[3]==4 and opposite[3]==1) \20                        or (my[3]==1 and opposite[3]==2) or (my[3]==2 and opposite[3]==1)\21                        or (my[3]==1 and opposite[3]==12) or (my[3]==12 and opposite[3]==1)\22                        or (my[3]==1 and opposite[3]==10) or (my[3]==10 and opposite[3]==1)\23                        or (my[3]==2 and opposite[3]==7) or (my[3]==7 and opposite[3]==2)\24                        or (my[3]==1 and opposite[3]==8) or (my[3]==8 and opposite[3]==1)\25                        or (my[3]==8 and opposite[3]==9) or (my[3]==9 and opposite[3]==8)\26                        or (my[3]==1 and opposite[3]==6) or (my[3]==6 and opposite[3]==1)\27                        or (my[3]==2 and opposite[3]==4) or (my[3]==4 and opposite[3]==2)\28                        or (my[3]==7 and opposite[3]==8) or (my[3]==8 and opposite[3]==7)\29                        or (my[3]==8 and opposite[3]==4) or (my[3]==4 and opposite[3]==8)\30                        or (my[3]==5 and opposite[3]==7) or (my[3]==7 and opposite[3]==5)\31                        or (my[3]==2 and opposite[3]==5) or (my[3]==5 and opposite[3]==2):32                    inside_if+=133                    all_indexes.update({my[0]:True})34                    all_indexes.update({opposite[0]:True})35                    break36        #     if state==1:37        #         overlap.append((my[3],opposite[3],my[0],opposite[0],my[2]-my[3]))38        #     if state==2:39        #         overlap.append((my[3],opposite[3],my[0],opposite[0],my[2]-opposite[1]))40        # if len(overlap)>0:41        #     to_save = max(overlap,key=itemgetter(4))42        #     save_confusion(to_save)43    return inside_if44def save_confusion(as_list):45    atype = KnowledgeType.objects.get(id=as_list[0])46    btype = KnowledgeType.objects.get(id=as_list[1])47    aid = MarkedUnit.objects.get(id=(max(as_list[2],as_list[3])))48    bid = MarkedUnit.objects.get(id=(min(as_list[2],as_list[3])))49    try:50        change_entry = Confusions.objects.get(idofa=aid,idofb=bid)51        change_entry.atype = atype52        change_entry.btype = btype53        change_entry.length = as_list[4]54    except Confusions.DoesNotExist:55        Confusions.objects.create(atype=atype,56                                  btype=btype,57                                  idofa=aid,58                                  idofb=bid,59                                  length=as_list[4])...

Automation Testing Tutorials

Learn to execute automation testing from scratch with LambdaTest Learning Hub. Right from setting up the prerequisites to run your first automation test, to following best practices and diving deeper into advanced test scenarios. LambdaTest Learning Hubs compile a list of step-by-step guides to help you be proficient with different test automation frameworks i.e. Selenium, Cypress, TestNG etc.