import sys import csv import math file = sys.argv[1] TotModelN = sys.argv[2] tmp=sys.argv[3] BPNUM=int(tmp) TotModelNum=int(TotModelN) print TotModelNum FileList=[] #initialize output arrays BPNUM_A=[] SHIFT_A=[] SLIDE_A=[] RISE_A=[] TILT_A = [] ROLL_A = [] TWIST_A = [] SHEAR_A = [] STRETCH_A = [] STAGGER_A = [] BUCKLE_A = [] PROPELLER_A = [] OPENING_A = [] MINPP_A = [] MAJPP_A = [] XDISP_A=[] YDISP_A=[] HRISE_A=[] INCL_A=[] TIP_A=[] HTWIST_A=[] BP_STR_A=[] #BPfile=open(file+'_numbp.dat','r') #for line in BPfile: # BPNUM_A.append(int(line)) #BPfile.close() outputInfoDict = { 'NumBPs': '''File name: Date and time: Number of base-pairs: Number of atoms: ''', 'baseRmsds': '''RMSD of the bases (----- for WC bp, + for isolated bp, x for helix change)''', 'hbondInfo': '''Detailed H-bond information: atom-name pair and length [ O N]''', 'polygonOverlapArea': ''' Overlap area in Angstrom^2 between polygons defined by atoms on successive bases. Polygons projected in the mean plane of the designed base-pair step. Values in parentheses measure the overlap of base ring atoms only. Those outside parentheses include exocyclic atoms on the ring. Intra- and inter-strand overlap is designated according to the following diagram: i2 3' 5' j2 /|\ | | | Strand I | | II | | | | | \|/ i1 5' 3' j ''', 'originAndMNVector': ''' Origin (Ox, Oy, Oz) and mean normal vector (Nx, Ny, Nz) of each base-pair in the coordinate system of the given structure ''', 'localBPPars': '''Local base-pair parameters''', 'localBPStepPars': '''Local base-pair step parameters''', 'localBPHelicalPars': '''Local base-pair helical parameters''', 'structureClass': '''Structure classification:''', 'lambda': ''' lambda: virtual angle between C1'-YN1 or C1'-RN9 glycosidic bonds and the base-pair C1'-C1' line C1'-C1': distance between C1' atoms for each base-pair RN9-YN1: distance between RN9-YN1 atoms for each base-pair RC8-YC6: distance between RC8-YC6 atoms for each base-pair ''', 'diNucStepClassRHelix': ''' Classification of each dinucleotide step in a right-handed nucleic acid structure: A-like; B-like; TA-like; intermediate of A and B, or other cases ''', 'grooveWidthsPPDist': ''' Minor and major groove widths: direct P-P distances and refined P-P distances which take into account the directions of the sugar-phosphate backbones (Subtract 5.8 Angstrom from the values to take account of the vdw radii of the phosphate groups, and for comparison with FreeHelix and Curves.) Ref: M. A. El Hassan and C. R. Calladine (1998). ``Two Distinct Modes of Protein-induced Bending in DNA.'' J. Mol. Biol., v282, pp331-343. ''', 'globalHelixAxis': ''' Global linear helical axis defined by equivalent C1' and RN9/YN1 atom pairs ''', 'mainChainAndChiAngles': ''' Main chain and chi torsion angles: Note: alpha: O3'(i-1)-P-O5'-C5' beta: P-O5'-C5'-C4' gamma: O5'-C5'-C4'-C3' delta: C5'-C4'-C3'-O3' epsilon: C4'-C3'-O3'-P(i+1) zeta: C3'-O3'-P(i+1)-O5'(i+1) chi for pyrimidines(Y): O4'-C1'-N1-C2 chi for purines(R): O4'-C1'-N9-C4 ''', 'sugarConfParameters': ''' Sugar conformational parameters: Note: v0: C4'-O4'-C1'-C2' v1: O4'-C1'-C2'-C3' v2: C1'-C2'-C3'-C4' v3: C2'-C3'-C4'-O4' v4: C3'-C4'-O4'-C1' tm: amplitude of pseudorotation of the sugar ring P: phase angle of pseudorotation of the sugar ring ''', 'intraStrandDist': '''Same strand P--P and C1'--C1' virtual bond distances''', 'helixRadius': ''' Helix radius (radial displacement of P, O4', and C1' atoms in local helix frame of each dimer) ''', 'diNucStepPosition': ''' Position (Px, Py, Pz) and local helical axis vector (Hx, Hy, Hz) for each dinucleotide step ''' } def chop(lst, n): """Chop lst into sublists of n items.""" return [lst[i:i+n] for i in range(0, len(lst), n)] def columnize(lst, n): # """Arrange the elements of lst into n columns.""" # assemble long & short columns q, r = divmod(len(lst), n) longcount = r # number of long columns longlen = q + 1 # length of long columns shortlen = q # length of short columns longcols = chop(lst[:longcount*longlen], longlen) shortcols = chop(lst[longcount*longlen:], shortlen) # pad short columns shortcols = [c + [None] for c in shortcols] # transpose rows = map(list, zip(*(longcols + shortcols))) # unpad short columns #del rows[-1][longcols:] return rows def identifyParameterBlock(block): # Identify the block by the first line # print '#',block.split('\n')[1] found = False for parameterBlockId, infoText in outputInfoDict.iteritems(): # print parameterBlockId, infoText try: if block.split('\n')[1].strip() in infoText: found = True break except IndexError: pass if found: return parameterBlockId, infoText else: return None, None def step2bp(step): ''' converts "GA/TC" to "G-C" ''' return step[0] + '-' + step[-1] #parameterBlockId="" for i in range(TotModelNum): i2 = str(i+1) FileList.append(file+"."+i2+".out") print FileList[i] #Traj_i = X3dna(FileList[i],i) x3dnaOutput = open(FileList[i], 'r').read() #print x3dnaOutput parameterBlocks = x3dnaOutput.split('****************************************************************************') for parameterBlock in parameterBlocks: #print parameterBlock parameterBlockId, infoText = None, None for paraBId, infoT in outputInfoDict.iteritems(): found = False try: if parameterBlock.split('\n')[1].strip() in infoT: found = True #print paraBId,infoT #print paraBId,'#',parameterBlock.split('\n')[1].strip() except IndexError: pass if found: parameterBlockId, infoText = paraBId, infoT #print parameterBlockId splitLines = parameterBlock.split('\n') #print parameterBlockId, "OUT OF LOOP" #Number of Base-pairs parseLine = False # if parameterBlockId == 'NumBPs': # print parameterBlockId,"CCCCCC" #print splitLines # for line in splitLines: # wordList = line.split() # print line," HEYOHEYO" # if line.strip() == '': # print "NEEXT1" # continue # elif line.startswith('File',0): # print "NEEXT" # continue # elif line.startswith('Date',0): # print "NEEXT2" # continue # elif wordList[2] == ('base-pairs:'): # print wordList[0:2] # parseLine = True # BP_Num = int(wordList[3]) # BPNUM_A.append(BP_Num) # print BPNUM,int(wordList[3]),"HERE IT ISSSS:" # elif line.strip().startswith('Number of atoms'): # print "NEEEXT3" # continue #Local Base-pair parameters if parameterBlockId == 'localBPPars': # bp Shear Stretch Stagger Buckle Propeller Opening # 1 G-C -0.44 -0.22 -0.62 -2.36 3.15 -3.95 # 2 A-T -0.09 -0.22 -0.19 12.78 -7.95 -0.93 # ...... # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # ave. 0.02 -0.24 3.33 0.20 4.00 33.41 # s.d. 0.71 0.90 0.26 3.68 8.17 6.97 #print parameterBlockId,"AAAAaaaaaa" #print splitLines for line in splitLines: wordList = line.split() #print wordList if line.strip() == '': continue elif wordList[0] == 'bp': parseLine = True continue elif wordList[0] == '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~': break elif parseLine == True: bp = int(wordList[0]) SHEAR_A.append(float(wordList[2])) STRETCH_A.append(float(wordList[3])) STAGGER_A.append(float(wordList[4])) BUCKLE_A.append(float(wordList[5])) PROPELLER_A.append(float(wordList[6])) OPENING_A.append(float(wordList[7])) if parameterBlockId == 'localBPStepPars': # step Shift Slide Rise Tilt Roll Twist # 1 GA/TC -0.86 0.54 2.99 -4.13 2.08 31.35 # 2 AA/TT -0.88 -0.27 3.34 -4.59 3.92 39.72 # ...... # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # ave. 0.02 -0.24 3.33 0.20 4.00 33.41 # s.d. 0.71 0.90 0.26 3.68 8.17 6.97 for line in splitLines: wordList = line.split() if line.strip() == '': continue #elif '----' in line: # line = line.replace('----', '0.00') # print line elif wordList[0] == 'step': parseLine = True continue elif wordList[0] == '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~': break elif parseLine == True: step = int(wordList[0]) bp = step step_str = wordList[1] bp_str = step2bp(step_str) BP_STR_A.append(bp_str) if wordList[2]=='----': SHIFT_A.append(" ") elif wordList[3]=='----': SLIDE_A.append(" ") elif wordList[4]=='----': RISE_A.append(" ") elif wordList[5]=='----': TILT_A.append(" ") elif wordList[6]=='----': ROLL_A.append(" ") elif wordList[7]=='----': TWIST_A.append(" ") else: SHIFT_A.append(float(wordList[2])) SLIDE_A.append(float(wordList[3])) RISE_A.append(float(wordList[4])) TILT_A.append(float(wordList[5])) ROLL_A.append(float(wordList[6])) TWIST_A.append(float(wordList[7])) if parameterBlockId == 'grooveWidthsPPDist': # Minor Groove Major Groove # P-P Refined P-P Refined # 1 GA/TC --- --- --- --- # 2 AA/TT --- --- --- --- # 3 AT/AT 11.2 --- 15.9 --- # 4 TT/AA 11.1 11.0 16.1 15.3 # 5 TG/CA 13.2 13.1 18.5 17.5 for line in splitLines: wordList = line.split() if line.strip() == '': continue elif wordList[0] == 'P-P': parseLine = True continue elif parseLine == True: line = line.replace('---', ' 0 ') # watch out these are not zero but NAN wordList = line.split() # NTdebug('wordList %r' % wordList) step = int(wordList[0]) #bp = step #step_str = wordList[1] #bp_str = step2bp(step_str) #step = step #step_str = step_str #bp_str = bp_str MINPP_A.append(float(wordList[2])) #minPP_ref = self.parseX3dnaFloat(wordList[3]), MAJPP_A.append(float(wordList[4])) #majPP_ref = self.parseX3dnaFloat(wordList[5]) if parameterBlockId == 'localBPHelicalPars': #step X-disp Y-disp h-Rise Incl. Tip h-Twist #1 CC/GG -2.72 -1.28 3.35 4.26 8.22 27.43 #2 CT/AG -2.41 0.99 2.66 23.78 -0.52 35.85 #3 TA/TA -2.20 -0.78 3.09 13.17 2.13 31.75 #..... # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # ave. -3.09 -0.34 2.98 14.72 0.94 32.33 # s.d. 1.98 1.76 0.55 9.74 12.96 6.57 for line in splitLines: wordList = line.split() if line.strip() == '': continue # elif '----' in line: # line = line.replace('----', '0.00') #print line elif wordList[0] == 'step': parseLine = True continue elif wordList[0] == '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~': break elif parseLine == True: step = int(wordList[0]) bp = step step_str = wordList[1] bp_str = step2bp(step_str) if wordList[2]=='----': XDISP_A.append(" ") elif wordList[3]=='----': YDISP_A.append(" ") elif wordList[4]=='----': HRISE_A.append(" ") elif wordList[5]=='----': INCL_A.append(" ") elif wordList[6]=='----': TIP_A.append(" ") elif wordList[7]=='----': HTWIST_A.append(" ") else: XDISP_A.append(float(wordList[2])) YDISP_A.append(float(wordList[3])) HRISE_A.append(float(wordList[4])) INCL_A.append(float(wordList[5])) TIP_A.append(float(wordList[6])) HTWIST_A.append(float(wordList[7])) # print i # print "DONE. NEXT",len(TIP_A),len(SHEAR_A) #print "DONE. NEXT",len(TIP_A),len(SHEAR_A) ####OUTPUT############### #initialize output files dotcsv='.csv' shearf=open(file+".shear.csv",'wb') bucklef=open(file+".buckle.csv",'wb') stretchf=open(file+".stretch.csv",'wb') propellerf=open(file+".propeller.csv",'wb') staggerf=open(file+".stagger.csv",'wb') openingf=open(file+".opening.csv",'wb') shiftf=open(file+".shift.csv",'wb') tiltf=open(file+".tilt.csv",'wb') slidef=open(file+".slide.csv",'wb') rollf=open(file+".roll.csv",'wb') risef=open(file+".rise.csv",'wb') twistf=open(file+".twist.csv",'wb') xdispf=open(file+".xdisp.csv",'wb') ydispf=open(file+".ydisp.csv",'wb') inclf=open(file+".incl.csv",'wb') tipf=open(file+".tip.csv",'wb') hrisef=open(file+".hrise.csv",'wb') htwistf=open(file+".htwist.csv",'wb') minppf=open(file+".minpp.csv",'wb') majppf=open(file+".majpp.csv",'wb') def printtofile(fname,outarray,numbp,TotModelNum): table=[] for i in range(TotModelNum): j=numbp tmp2=outarray tmp3=tmp2[0:j] table.append(tmp3) tmp2f=tmp2[0:j] del tmp2[0:j] writer=csv.writer(fname,delimiter=' ') for row in table: writer.writerow(row) def printtofilebpstep(fname,outarray,numbp,TotModelNum): table=[] for i in range(TotModelNum): j=numbp tmp2=outarray tmp3=tmp2[0:j] table.append(tmp3) tmp2f=tmp2[0:j] del tmp2[0:j] writer=csv.writer(fname,delimiter=' ') for row in table: writer.writerow(row) Totbp=BPNUM-1 printtofilebpstep(shiftf,SHIFT_A,Totbp,TotModelNum) printtofilebpstep(tiltf,TILT_A,Totbp,TotModelNum) printtofilebpstep(slidef,SLIDE_A,Totbp,TotModelNum) printtofilebpstep(rollf,ROLL_A,Totbp,TotModelNum) printtofilebpstep(risef,RISE_A,Totbp,TotModelNum) printtofilebpstep(twistf,TWIST_A,Totbp,TotModelNum) printtofilebpstep(xdispf,XDISP_A,Totbp,TotModelNum) printtofilebpstep(inclf,INCL_A,Totbp,TotModelNum) printtofilebpstep(ydispf,YDISP_A,Totbp,TotModelNum) printtofilebpstep(tipf,TIP_A,Totbp,TotModelNum) printtofilebpstep(minppf,MINPP_A,Totbp,TotModelNum) printtofilebpstep(majppf,MAJPP_A,Totbp,TotModelNum) printtofilebpstep(htwistf,HTWIST_A,Totbp,TotModelNum) printtofilebpstep(hrisef,HRISE_A,Totbp,TotModelNum) printtofile(shearf,SHEAR_A,BPNUM,TotModelNum) printtofile(bucklef,BUCKLE_A,BPNUM,TotModelNum) printtofile(stretchf,STRETCH_A,BPNUM,TotModelNum) printtofile(propellerf,PROPELLER_A,BPNUM,TotModelNum) printtofile(staggerf,STAGGER_A,BPNUM,TotModelNum) printtofile(openingf,OPENING_A,BPNUM,TotModelNum) shiftf.close() tiltf.close() slidef.close() rollf.close() risef.close() twistf.close() shearf.close() bucklef.close() stretchf.close() propellerf.close() staggerf.close() openingf.close() xdispf.close() inclf.close() ydispf.close() tipf.close() hrisef.close() htwistf.close() minppf.close() majppf.close()