ENORA and multi-state structure calculations: Difference between revisions
Line 797: | Line 797: | ||
With a text editor, edit the 'CALC_multistate.cya' macro to activate the inactive commands (by deleting the preceeding hashtag #) necessary to perform the grouping of states and splitting of the conformers. | With a text editor, edit the 'CALC_multistate.cya' macro to activate the inactive commands (by deleting the preceeding hashtag #) necessary to perform the grouping of states and splitting of the conformers. | ||
With a text editor, change the number of states from one to two in the PREP.cya macro. | |||
With a text editor, change the number of states (nbundle) from one to two in the 'PREP.cya' macro: | |||
syntax nbundle=@i=1 togetherweight=@r=0.1 multitensor | |||
When you are done preparing the macros as outlined perform the calculation by running the 'CALC_multistate.cya' macro: | When you are done preparing the macros as outlined perform the calculation by running the 'CALC_multistate.cya' macro: |
Revision as of 16:31, 4 April 2019
In this tutorial we will provide you with a guided example for calculating eNOEs and multi-state structure calculations.
To this end we will first run the modules of eNORA and then use the obtained eNOEs to calculate a single state and a two-state structure model using automated sorting to group the states. Along the way you will learn some additional CYANA skills useful for other purposes as well.
The eNORA module offers in principle two methods to calculate spin diffusion, FRM and TSS. FRM is the generally applicable way to do these calculations and we will set a main focus on that method, will however have one section which explains the principles at play for TSS and how to set it up.
To finalize you will .... And ultimately you can try to improve ....
CYANA setup
Obtaining and installing the CYANA demo version and data
Please follow the following steps carefully (exact Linux commands are given below; you may copy them to a terminal):
- Go to your home directory (or data directory).
- Get the demo data from the server.
- Unpack the input data for the practical.
- Get the demo version of CYANA.
- Unpack CYANA.
- Setup the CYANA environment variables.
- Change into the newly created directory 'eNORA'.
- Copy the demo_data directory to 'enoe'.
- Change into the subdirectory 'enoe'.
- Test whether CYANA can be started by typing its name, 'cyana'.
- Exit from CYANA by typing 'q' or 'quit'.
cd ~ wget 'http://www.cyana.org/wiki/images/6/64/eNORA_multiState.tar.gz'
on mac OS X (use curl instead of wget):
curl http://www.cyana.org/wiki/images/6/64/eNORA_multiState.tar.gz -o eNORA_multiState.tar.gz
then
tar zxf eNORA_multiState.tar.gz wget 'http://www.cyana.org/wiki/images/6/64/Cyana-3.98.9_Demo.tgz'
again, on mac OS X:
curl http://www.cyana.org/wiki/images/6/64/Cyana-3.98.9_Demo.tgz -o Cyana-3.98.9_Demo.tgz
then
tar zxf Cyana-3.98.9_Demo.tgz cd cyana-3.98.9/ ./setup cd ~ cd eNORA cp -r demo_data enoe cd enoe
Try to run CYANA by entering 'cyana' at the command prompt of your terminal (q to quit cyana):
cyana ___________________________________________________________________ CYANA 3.98 (mac-intel) Copyright (c) 2002-17 Peter Guentert. All rights reserved. ___________________________________________________________________ Demo license valid for specific sequences until 2018-12-31 cyana> q
If all worked, you are ready to go in terms of everything related to CYANA!
Execution scripts or "macros" in CYANA
For more complex task within CYANA, rather than to enter the execution commands line by line at the CYANA prompt, the necessary commands are collected in a file named '*.cya'. Collecting the commands in macros has the added advantage, that the macros serve as a record allowing to reconstruct previous calculations.
Hint: For comprehensive information on the CYANA commands etc. is in the CYANA 3.0 Reference Manual.
eNOE calculations
All eNOE related calculations within cyana are carried out using the eNORA modules.
NOESY experiment measured at different mixing times (keeping the mixing times as much as possible within the linear regime of NOE buildup) supply very precise distance restraints used for a structure calculation. In addition other restraints such as backbone angles from chemical shifts and scalar couplings for backbone and aromatic side chains are also used.
Experimental input data
Peak lists in XEASY format are prepared by automatic peak picking with a visualization program such as CcpNmr Analysis, NMRdraw or NMRview and saved as XXX.peaks, where XXX denotes the name of the xeasy peak list file. Since NMRdraw peak lists are of different file type, cyana provides the command read tab to convert the files to XEASY format.
# Number of dimensions 3 #FORMAT xeasy3D #INAME 1 H #INAME 2 HN #INAME 3 N #SPECTRUM N15NOESY H HN N 17086 4.098 4.099 57.441 1 U 6.990943E+08 0.000000E+00 e 0 HA.5 HA.5 CA.5 89532 4.355 1.829 33.507 1 U 1.720779E+06 0.000000E+00 e 0 HA.6 HB2.6 CB.6 89544 4.353 1.757 33.513 1 U 2.939628E+06 0.000000E+00 e 0 HA.6 HB3.6 CB.6
The first line specifies the number of dimensions (3 in this case). The '#SPECTRUM' (no space between characters) lines gives the experiment type (N15NOESY, which refers to the corresponding experiment definition in the CYANA library), followed by an identifier for each dimension of the peak list (H HN N) that specifies which chemical shift is stored in the corresponding dimension of the peak list. The experiment type and identifiers must correspond to an experiment definition in the general CYANA library (see below) in most uses of the definition, here however we cheat slightly and get away with it. We are cheating, because for eNOE calculations we record our NOESY spectra with simultanous evolution of 13C and 15N dimensions, since we require 15N and 13C bound spins within the same spectrum for purposes of normalization (see...).
After the '#SPECTRUM' line follows one line for every peak. For example, the first peak in the 'HNCA.peaks' list has
- Peak number 17086
- H chemical shift 4.098 ppm
- ("HN") chemical shift 4.099 ppm (in this case 13C bound)
- Heavy atom chemical shift 57.441 ppm (in this case 13C labeled)
The other data are relevant entry for the eNOE mudules is the peak volume or intensity (6.990943E+08).
Assignment dimensions have to be arranged with the flow of magnetization, first row spin i, second row spin j
The protein sequence is supplied by three-letter code in the demo.seq file.
SPECTRUM definitions in the CYANA library
When you start CYANA, the program reads the library and displays the full path name of the library file. You can open the standard library file to inspect, for example, the NMR experiment definitions . For instance, the definition for the N15NOESY spectrum (search for 'N15NOESY' in the library file 'cyana.lib') is
SPECTRUM N15NOESY H HN N 0.900 N:N_AM* HN:H_AMI ~4.0 H:H_* 0.800 N:N_AM* HN:H_AMI ~4.5 H:H_* 0.700 N:N_AM* HN:H_AMI ~5.0 H:H_* 0.600 N:N_AM* HN:H_AMI ~5.5 H:H_* 0.500 N:N_AM* HN:H_AMI ~6.0 H:H_*
The first line corresponds to the '#SPECTRUM' line in the peak list. It specifies the experiment name and identifies the atoms that are detected in each dimension of the spectrum. The number of identifiers defines the dimensionality of the experiment (3 in case of N15NOESY).
Each line below defines a (formal) magnetization transfer pathway that gives rise to an expected peak. in the case of N15NOESY there are five lines, corresponding to the through space magnetization transfer by dipol-dipol mechanism. The peak definition starts with the probability to observe the peak (0.900), followed by a series of atom types, e.g. H_AMI for amide proton etc. The atoms whose chemical shifts appear in the spectrum are identified by their labels followed by ':', e.g. for N15NOESY 'H:', 'HN:', and 'N:'. If you were to use the CYANA functions to simulate peaks, expected peaks are generated for each molecular fragment in which these atom types occur.
You may have realized that our peak list contains peaks that are 13C bound, therefore the spectrum definition is wrong, since we are only reading the peak lists and not generating any, this is not a problem.
eNORA
- work in the copy of the data directory ('cd enoe')
Using the text editor of your choice, create your 'init.cya' macro as outlined below (The init macro) and also your 'CALC_enoe.cya' macro. Be extra careful to avoid typos and unwanted spaces in coma lists etc.
The init macro
The initialization macro file has the fixed name 'init.cya' and is executed automatically each time CYANA is started. It can also be called any time one wants to reinitialize the program by typing 'init'. It contains normally at least two commands that read the CYANA library and the protein sequence:
rmsdrange:=8-33 cyanalib read seq demo.seq swap=0 expand=1
The first line sets the appropriate rmsdrange, and the command 'cyanalib' reads the standard CYANA library. The next command reads the protein sequence.
The protein sequence is stored in three-letter code in the file 'demo.seq'.
Stereo-specificity of dia-stereocenters
atoms stereo "HA? 10" atoms stereo "HB? 7 8 11 13 14 21 23 24 25 26 27 33 34 35 37 38" atoms stereo "HG? 8 12 14 17 36 37" atoms stereo "HD2? 18 26 30" atoms stereo "QG? 22" atoms stereo "QD? 7" atoms stereo "HE2? 33" atoms stereo "HD? 8 14 21 37" atoms stereo "HG1? 28"
However, one may do the following to supply all atoms as stereo specific:
atoms select atoms stereo
or to supply all atoms as non stereo specific, use:
atoms select atoms stereo delete
To get a feedback of the supplied stereo specific assignment add to your 'init.cya' the command:
atoms stereo list
D20 exchange
With 3% D2O in the nmr buffer for exchange of backbone amide atoms
atoms set "H" protlev=0.97
The eNORA CALC macro
The 'CALC_enoe.cya' starts with the following:
- Conversion of the nmrDraw peak file to XEASY format with the atom assignments in the file.
- Sorting of the peaks in each peaks list per mixing time.
- Writing out the peaks in XEASY format with atom names contained in the file, one peak list for each mixing time.
- The parameter definitions, their function is explained later with respect to the functions that use them.
- The structure input and which conformers to use for spin-diffusion calculations. If multiple conformers are used to average the spin-diffusion values, the command 'structure select' may be used to select multiple conformers of the pdb.
- Reading the XEASY peak lists in the order of increasing mixing time.
# ------------------------------- nmrDraw peak file conversion ------------------------------- # read nmrDraw tab file read tab demo.tab # sort the peaks peaks sort # write out peak lists do i 1 npkl peaks select "** list=${i}" write peaks $i.peaks names end do # ---------------------------------- eNORA routine ---------------------------------- # ---------------------------------- basic parameter definitions ---------------------------------- echo:= on mixingtimes = '0.02,0.03,0.04,0.05,0.06' protlevel = 0.97 avgrho = 5.3 tauc = 4.25 b0field = 700 maxdistance = 6.5 rhofile = 'rhoInApo.rho' normed = 0 normspin = 2 mode = 1 bname = 'bupplots' dname = 'decplots' # ---------------------------------- structure input ---------------------------------- # specify the conformers for calculations read pdb demo rigid #structure select 1-40 structure select 1 # ---------------------------------- peak file reading ---------------------------------- # read in the peak lists do i 1 npkl read peaks $i.peaks $if(i.eq.1,' ','append') end do # supply averaged rho or izero values if (existfile('$rhofile')) then read rho $rhofile end if # ---------------------------------- initializing the routine ---------------------------------- # initialize the routine, fit experimental decays and buildups enoe init normalize=normspin normed=normed rhoavg=avgrho time=$mixingtimes # print average experimental auto-relaxation and I(0) values to screen enoe avgExpVal # plot the diagonal decay's enoe decay plot=$dname graf $dname.pdf # write the auto-relaxation values to file write rhoOut.rho exit # ---------------------------------- spin-diffusion ---------------------------------- # calculate the spin-diffusion correction enoe spindiff b0field=b0field tauc=tauc maxdist=maxdistance mode=mode # ---------------------------------- apply spin-diffusion ---------------------------------- # apply spin-diffusion correction to experimental buildups (if calculated) and calculate reff enoe reff b0field=b0field tauc=tauc # prepare the cyana restraints (scaling, error margins) enoe restraint errStereoFlag=1 errStereo=-1 chiN=-1 b0field=b0field tauc=tauc # ---------------------------------- write restraints ---------------------------------- # delete fixed distances from upl/lol output distance delete fixed # write the distance restraints to file write upl enoe.upl write lol enoe.lol # --------------------------------- plotting buildups -------------------------------- enoe buildup b0field=b0field tauc=tauc plot=$bname graf $bname.pdf # ---------------------------------- overview ---------------------------------- enoe overview
When you have prepared the 'init.cya' and the 'CALC_enoe.cya' run the eNOE calculation; one could start CYANA and execute the 'CALC_enoe.cya' macro from the CYANA prompt as such:
cyana CALC_enoe.cya
The program if all is setup properly the program will run and display averages per spin type of autorelaxation and I (0) values and then stop as soon as it finished plotting the diagonal decay's and has written the 'rhoOut.rho' file.
The command
exit
following the command 'write rhoOut.rho' stops the routine from running to completion. Before commenting out (#) the exit command and running the routine to its end, do the following exercise.
Exercise x: Compiling the autorelaxation file
Before compiling the rhoIn.rho file, check the diagonal decays for their quality. Edit out any diagonal peaks from the tab file that give bad decays, then run the routine again the same way to get averages of rho and I(0) values. Then compile the rhoIn.rho file by copying the 'rhoOut.rho' file and adding the lines from 'missRhoIzero.out', filling in the rho values as needed.
When compiling the rhoIn.rho file see * auto-relaxation and I(0) values.
eNORA output files
The FLYA algorithm will produce the following output files:
- missRhoIzero: List of spins that lack a diagonal decay, and do not provide a rho or I(0) value (if you set normed=0).
- nonNormalizableNOEs.out: List of NOEs that lack a diagonal decay and are not normalizable (if you set the parameter normed=0 to normed=1, all non normalizable NOES are left out).
- rhoOut.rho: Experimentally fitted rho and I(0) values values form diagonal decays.
- enoe.upl and enoe.lol: Upper limit and lower limit restraint files with tolerances applied.
- enoe.ovw: Collated results file.
- bupplots:
- decplots:
The missRhoIzero file
List of spins that lack a diagonal decay and therefore have no experimentally determined rho and I(0) values. The file is formatted to be used as basis for a rhoIn.rho file derived from average values. The file may also be used to setup generic normalized eNOE calcuations, see .....
The rhoOut.rho file
The enoe.upl/lol distance restrains
Tolerances Comments
The enoe.ovw file
Collated file that may be generated at any time during the routine and will be populated with the values available at the momentary progress of calculations.
- ASSIGNMENT(i->j): Assignment arranged with the flow of magnetization.
- REFFixEXP: The experimental sigma of spin x, where x=i,j
- REFFxCORR: The reff of spin x, where x=i,j
- SIGxEXP: The experimental sigma of spin x, where x=i,j
- SIGxCORR: The spin-diffusion corrected experimental sigma of spin x, where x=i,j
- SDCx: The spin-diffusion correction of spin x, where x=i,j
- SDPx: The number of other partner spins involved in spin-diffusion for spin x, where x=i,j
- IZEROx: The back calculated I(0) value of spin x, where x=i,j
- exp_chiNx: The goodness of fit of the experimental data, where x=i,j
- corr_chiN x: The goodness of fit of the spin-diffusion corrected experimental data, where x=i,j
Considerations regarding the obtained eNOE restraints
Exercise xx: Preparing an xray structure to use within CYANA
Deposited structures often lack specific features. i.e. Xray structures often lack proton coordinates or contain mutations and ligands.
Copy your noecc results to a new directory call regulabb, then delete all the previous, unnecessary output files to reduce clutter and have better oversight.
cp -r noecc regulabb cd regulabb
After reading the sequence file, the pdb file can be read with the option unknown=warn or unknown=skip, this will then skip the parts of the molecule not specified in the sequence file.
read pdb xxxx.pdb unknown=warn
Other options to read pdb's:
read xxxx.pdb unknown=warn hetatm new
where the option 'hetatm' allows for reading of coordinate labeled HETATM, rather than ATOM in the pdb. 'new' will read the sequence from the pdb.
To write back out pdb's and sequences:
write pdb XXX.pdb write seq XXX.seq
Inspect the pdb using chimera: Now, there are several issues besides HETATM, that make the comparison to the calculated NMR structure not possible within CYANA before you fix them.
Using the regularize command one can get a structure calculated within CYANA that has these issues fixed but still is very close to the input structure of your choice. Download the Xray structure used for this excercise:
wget 'https://files.rcsb.org/view/1PIN.pdb'
or for mac os x:
curl https://files.rcsb.org/view/1PIN.pdb -o 1PIN.pdb
Create an 'init.cya' macro with:
cyanalib
Then create a 'CALC_reg.cya' macro with:
read 1PIN.pdb unknown=warn hetatm new write 1PIN.seq write 1PIN_2.pdb
In the provided data, there are two mutations to the sequence (S18N and W34F), furthermore we do not have a ligand and the expressed protein for NMR is truncated. To have cyana truncate the xray structure and build in the mutations:
read demo.seq read 1PIN_2.pdb rigid unknown=warn write 1PIN_ed.pdb ./init read pdb 1PIN_ed.pdb regularize steps=20000 link=LL keep
Execute the 'CALC_reg.cya' macro in the CYANA shell (or use only one processor, do not distribute the job):
cyana CALC_reg.cya
Mapping calculated eNOE restraints onto a known structure
One can map the calculated restraints, such as distance restraints (upl/lol) onto a known structure (in the example here the modified xray structure). This is an approach to analyze restraints and their influence on the results.
Below you find the commands to accomplish this. You see by studying the commands, which files are needed to execute the macro. Therefore, create a new directory ('mkdir') or copy a directory containing the respective files. Delete what you do not need. Use the regularized xray structure from exercise xxxx.
You need an init file:
rmsdrange:=8-33 cyanalib read seq demo.seq
And the main macro (name it 'CALC_xraymap.cya'):
read upl enoe.upl read lol enoe.lol read regula.pdb unknown=warn weight_vdw=0 overview enoe_xray.ovw
- If the restraints do not match with the xray structure, does it mean they are wrong?
- If you tried the two options, what is (are) the difference(s)?
eNOE calculations and the TSS approach to spin-diffusion calculations
Labeling Schemes
Deuterated labeling schemes often involve methyl labeling with 3% D2O in the nmr buffer, i.e.
atoms set "H" protlev=0.97 # labeling scheme: VAL_G1 0% LEU_D1 0% ILE_D1 0% atoms set "QG1 @VAL + QD1 @LEU + QD1 @ILE" protlev=0.0
Using Talos to generate torsion angle restraints
Torsion angle restraints from the backbone chemical shifts help restrict angular conformation space. We wish to use only "strong assignments" to generate these restraints.
If you do not have TALOS installed get it from here. It is part of the nmrpipe software package.
Exercise x: Calculate backbone torsion angle restraints using Talos
cp -r enoe acoPREP cd acoPREP rm *.tab enoe* *.out *.rho
Hint: You only need the 'init.cya' and the '1.peaks' file to calculate the torsion angle restraints.
Use a text editor of your choice to create a 'CALC_talos.cya' file with the commands to calculate the talos angle restraints:
# convert read 1.peaks shifts initialize shifts adapt atom set "* shift=990.0.." shift=none write demo.prot read prot demo.prot unknown=skip talos talos=talos+ talosaco pred.tab write aco talos.aco
This will call the program TALOS+ and store the resulting torsion angle restraints in the file 'talos.aco'.
Since this is not a calculation suited for the MPI scheduler, start CYANA first, then call the 'CALC_talos.cya' macro from the prompt.
Hint: change to a cshell before running cyana (since talos needs a cshell to run):
csh
Multi-state structure calculation
We will perform calculations based on eNOEs by using torsion angle dynamics in order to compute the three-dimensional structure of the protein.
The 'enoe.upl' and 'enoe.lol' files will be used together with the aco based on chemical shifts of the backbone and scalar couplings from backbone, Ha-HB and aromatic residues determined by experiment.
Calculating a single state structure
In principle you could easily do the single-state structure calculation as a regular calculation, using your upl/lol, aco and cco files as input. This would look something like this:
syntax inputseed=@i=3771 # ------ Structure calculation ------ read upl XXX.upl read lol XXX.lol read aco XXX.aco read cco XXX.cco anneal_weight_aco := 1.0, 1.0, 0.0, 0.0 anneal_weight_cco := 0.0, 0.5 seed=inputseed calc_all 100 steps=50000 if (master) then cut_cco=1.0 cut_rdc=3.0 weight_aco = 0.0 rmsdrange:=8-33,108-133 overview sstate structures=20 pdb end if
However, since we end up calculating also multi-state structures later on, it makes sense to setup the single-state calculation exactly the same way as the multi-state calculations, and only edit as few parameters as possible. As soon as you see the PREP script below, you will realize why this makes sense.
Exercise: The single state CALC macro
Copy the 'enoe' directory and give it the name 'sstate', then delete all the files and data we do not need to reduce clutter and have better oversight.
cp -r enoe sstate cd sstate rm *.peaks *.out *.job rm -rf details
From the directory 'acoPREP' copy the calculated talos restraints ('talos.aco').
The init macro
The init macro needs to be updated to read the multi-state sequence in order to prepare the restraints and run the structure calculation:
cyanalib if (existfile('bundle.seq')) then read seq bundle.seq molecules define * atom set * vdwgroup=bundle else read seq demo.seq end if rmsdrange:=8-33 swap=0 expand=1 atoms stereo "HA? 10" atoms stereo "HB? 7 8 11 13 14 21 23 24 25 26 27 33 34 35 37 38" atoms stereo "HG? 8 12 14 17 36 37" atoms stereo "HD2? 18 26 30" atoms stereo "QG? 22" atoms stereo "QD? 7" atoms stereo "HE2? 33" atoms stereo "HD? 8 14 21 37" atoms stereo "HG1? 28" atoms set "H" protlev=0.97
The PREP macro
syntax nbundle=@i=1 togetherweight=@r=0.1 multitensor #multitensor=.true. together=.true. moloffset=100 # ------ Sequence file ------ read seq demo.seq print "\# Bundle of $nbundle conformers" >bundle.seq do j 1 nbundle do i 1 nr if (j.lt.nbundle .and. rnam(i).eq.'PL') break print "$rnam(i) ${$rnum(i)+moloffset*(j-1)}" >> end do if (j.lt.nbundle) print "PL LL2 LL2 LL2 LL2 LL2 LL2 LL2 LL2 LL2 LL2 LL2 LL2 LL2 LL2 LP" >> end do print >>. # ------ Make bundle angle restraints ------ read aco demo.aco write aco bundle.aco do j 2 nbundle atom set * residue=residue+moloffset write aco bundle.aco append end do # ------ Make bundle coupling constant restraints ------ read demo.seq read cco demo_backbone.cco read cco demo_aro.cco append read cco demo_JHaHb.cco append print "\# Coupling constant restraint file" >bundle.cco do i 1 ncco i1=iccoa(1,i); i2=iccoa(2,i) do j 1 nbundle m=moloffset*(j-1) print "${$rnum(iar(i1))+m} $rnam(iar(i1)) $anam(i1) ${$rnum(iar(i2))+m} $rnam(iar(i2)) $anam(i2) $cco(i) $tolcco(i) 1.0 $karplus(1,i) $karplus(2,i) $karplus(3,i) $if(j.ne.nbundle,'&',' ')" >>bundle.cco end do end do print >>. read seq bundle.seq read cco bundle.cco write cco bundle.cco karplus # ------ Make bundle RDC restraints ------ #read bundle.seq #read rdc demo.rdc #print "\# RDC restraint file" >bundle.rdc #do i 1 orientations # if (multitensor) then # do j 1 nbundle # print "${i+orientations*(j-1)} $magnitude(i) $rhombicity(i) ${$rnum(iar(irtena(4,i)))+moloffset*(j-1)}" >> # end do # else # print "$i $magnitude(i) $rhombicity(i) $rnum(iar(irtena(4,i)))" >> # end if #end do #do i 1 nrdc # i1=irdca(1,i); i2=irdca(2,i) # do j 1 nbundle # m=moloffset*(j-1) # iten=irdct(i); if (multitensor) iten=iten+orientations*(j-1) # print "${$rnum(iar(i1))+m} $rnam(iar(i1)) $anam(i1) ${$rnum(iar(i2))+m} $rnam(iar(i2)) $anam(i2) $rdc(i) $tolrdc(i) $weirdc(i) $iten $rdcsca(i) $if(j.lt.nbundle,'&',' ')" >>bundle.rdc # end do #end do #print >>. #read seq bundle.seq #read rdc bundle.rdc #write rdc bundle.rdc # ------ Make ambiguous bundle distance restraints ------ subroutine PURGE #distance delete "HA 9, HB2 9" end init read upl enoe.upl PURGE distance modify info=full molecules symmetrize if (nbundle.gt.1) distances set "$moloffset.., $moloffset.." bound=0.0 distances set "*, *" bound=bound*(1.0*nbundle)**(-1.0/6.0) write upl bundle.upl init read lol enoe.lol PURGE distance modify info=full molecules symmetrize if (nbundle.gt.1) distances set "$moloffset.., $moloffset.." bound=0.0001 distances set "*, *" bound=bound*(1.0*nbundle)**(-1.0/6.0) write lol bundle.lol # ------ Make restraints to keep corresponding atoms together ------ if (together .and. nbundle.gt.1 .and. togetherweight.gt.0.0) then read seq bundle.seq molecules define * atom set * vdwgroup=bundle atom select "N C*" do i 1 na if (iamol(i).ne.1) break if (asel(i)) then distance make "$atom(i)" "$anam(i) ${$rnum(iar(i))+moloffset}" upl=1.2 weight=$togetherweight info=none end if end do distances set "* - N CA C CB, * - N CA C CB" weight=weight*0.1 molecules symmetrize distances unique write upl together.upl end if
The single-state CALC macro
Inside the 'sstate' directory, use a text editor to create the 'CALC_sstate.cya' file for structure calculation as outlined below.
syntax inputseed=@i=3771 if (master) then PREP end if # ------ Structure calculation ------ read upl bundle.upl read lol bundle.lol read aco bundle.aco read cco bundle.cco #read rdc bundle.rdc if (existfile('together.upl')) read upl together.upl append anneal_weight_rdc := 0.0, 0.5 anneal_weight_aco := 1.0, 1.0, 0.0, 0.0 anneal_weight_cco := 0.0, 0.5 seed=inputseed calc_all 100 steps=50000 if (master) then cut_cco=1.0 cut_rdc=3.0 weight_aco = 0.0 rmsdrange:=8-33,108-133 overview bundle structures=20 pdb #read pdb bundle.pdb #rmsdrange:=11-16,21-26,31-34,111-116,121-126,131-134 #overview bundleSec structures=20 pdb # reference=xxx.pdb #molecules sort "BACKBONE 26-30" base=1 #write sortStates.pdb all #SPLIT end if
These commands tell the program to calculate 100 conformers, and to analyze the best 20 of them. 50000 torsion angle dynamics steps will be applied per conformer.
When you are done preparing the macros as outlined run the calculation.
The structure calculation will be performed by running the 'CALC_sState.cya' macro:
cyana -n 33 CALC_sState.cya
Doing this, basically means each processor will calculate 100/33=3 conformers. If you changed the setup to calculate 50 structures, you would start the calculation with 'cyana -n 25 CALC_sState.cya'.
Carefully analyze the WARNING and ERROR messages if any.
Statistics on the the structure calculation will be displayed to screen. The final structure will be 'bundle.pdb'.
Multi-state structure calculation
The SPLIT macro
Each conformer now contains to sets of coordinates corresponding to two states, and the
moloffset=100 read pdb sortStates.pdb n=nstruct write_all split nbundle=$rnum(nr)/moloffset+1 show nbundle do i 1 nstruct read seq bundle.seq read pdb split$i(I3.3).pdb do j 1 nbundle atoms select 1-80 write pdb split$i(I3.3)-$j.pdb selected atoms set * residue=residue-moloffset end do read seq demo.seq read_all split$i(I3.3)-*.pdb unknown=skip write split$i(I3.3).pdb all do j 1 nbundle remove split$i(I3.3)-$j.pdb end do end do read seq demo.seq do i 1 n read pdb split$i(I3.3).pdb append end do write pdb splitall.pdb all rmsd
Grouping the coordinates of multi-state calculations
Exercise x: A two-state structure calculation
Copy the 'sstate' directory and give it the name 'twostate', then delete all the previous, unnecessary output files to reduce clutter and have better oversight. Copy the 'CALC_sstate.cya' and rename it 'CALC_multistate.cya'.
cp -r sstate twostate cd twostate mv CALC_sstate.cya CALC_multistate.cya rm *.out *.job final* rama*
With a text editor, edit the 'CALC_multistate.cya' macro to activate the inactive commands (by deleting the preceeding hashtag #) necessary to perform the grouping of states and splitting of the conformers.
With a text editor, change the number of states (nbundle) from one to two in the 'PREP.cya' macro:
syntax nbundle=@i=1 togetherweight=@r=0.1 multitensor
When you are done preparing the macros as outlined perform the calculation by running the 'CALC_multistate.cya' macro:
cyana -n 33 CALC_multistate.cya.cya
Carefully analyze the WARNING and ERROR messages if any.
Results: analysis
Download and install the molecular viewer Chimera
- Download Chimera (to your personal laptop) from: Chimera
Exercise xx: Single state structure analysis
The final structure will be 'final.pdb'. You can visualize it, for example, with the command
chimera final.pdb
Analyze the result, the bundle seems unnatuarly tight for an NMR structure bundle. Why?
Exercise xx: Calclulate the RMSD of NMR vs. xray structure using a CYANA macro
Using the INCLAN language of CYANA (Writing and using INCLAN macros,Using INCLAN variables,Using INCLAN control statements) it is possible to write complex macros that interact with the FORTRAN code of CYANA. Reading internal variables and manipulating them to achieves custom task.
- save the manually edited xray structure (exercise 11) or the the regularized xray structure (containing the ligand and called 'regula.pdb') as 'reg_xray.pdb' to use the macro below (or change the name in the macro accordingly).
- what do you think about the RMSD, does the value make sense? Does the range make sense?
Below you find the commands for a macro (call it 'CALC_RMSD.cya') that will read the regularized xray structure and the calculated nmr structure, then calculating the rmsd of both the protein and ligand parts of the complex:
read demoLong.seq rmsd range=8-33 structure=final.pdb reference=reg_xray.pdb atom select "BACKBONE 8-33" t=rmsdmean j=rindex('333') n=0 s=0.0 do i ifira(j) ifira(j+1)-1 if (element(i).gt.1) then n=n+1 s=s+displacement(i) end if end do print "RMSD of the LIG: ${s/n} ($n atoms)" read pdb final.pdb structure mean write pdb mean.pdb read pdb mean.pdb read pdb reg_xray.pdb append atom select "BACKBONE 15-111" t=rmsdmean atom select "WITHCOORDALL" j=rindex('333') n=0 s=0.0 do i ifira(j) ifira(j+1)-1 if (element(i).gt.1.and.asel(i)) then n=n+1 s=s+displacement(i)*2 end if end do print "Displacement of the LIG (to ref xray): ${s/n} ($n atoms)"
Exercise xx: Analysis of multi-state calculations
On improving the final structure
Using what you have learned so far, employing some of the options
General questions to answer regarding this task:
- Name additional experimental restraints (or inputs) you could use for structure calculation.
- Name additional NMR experiments you could measure, to acquire experimental data that are not supplied with the demo_data.
eNORA extensions and options
There are a variety of commands to modify eNORA runs to accommodate experimental labeling schemes or etc...
Averaging Of Spin-Diffusion Over Multiple Conformers
Generating XEASY peak list with expected FRM or two spin intensities
# -------------------------- get the shifts from a XEASY peaks list --------------------------- ./init # convert read 1.peaks shifts adapt contribution=0.0 shifts renumber atoms set "* shift=900.0.." shift=none write demo.prot # -------------------------------- basic parameter definitions -------------------------------- ./init echo:= on mixingtimes:= 0.02,0.03,0.04,0.05,0.06 tauc = 4.25 b0field = 700 maxdistance = 6.5 # ---------------------------------- structure input ----------------------------------- # specify the conformers for calculations read pdb demo rigid structure select 1 # ---------------------------------- peak input ---------------------------------- loadspectra structure=demo.pdb peaks=N15NOESY,C13NOESY prot=demo.prot simulate # ---------------------------------- run eNORA elements and write peaks ---------------------------------- do n 1 length('mixingtimes') # FM enoe spindiff b0field=b0field tauc=tauc time=$mixingtimes(n) mode=3 labilatom='NONE' read peaks N15NOESY_exp.peaks enoe intensities write peaks N15NOESY_FM_$n.peaks names read peaks C13NOESY_exp.peaks enoe intensities write peaks C13NOESY_FM_$n.peaks names # 2 spin enoe twospin b0field=b0field tauc=tauc time=$mixingtimes(n) mode=3 labilatom='NONE' read peaks N15NOESY_exp.peaks enoe intensities mode=2 write peaks N15NOESY_2spin_$n.peaks names read peaks C13NOESY_exp.peaks enoe intensities mode=2 write peaks C13NOESY_2spin_$n.peaks names end do read peaks C13NOESY_FM_1.peaks peaks2dplot dimensions=12 read peaks C13NOESY_FM_1.peaks read peaks N15NOESY_FM_1.peaks append shifts initialize shifts adapt atom set "* shift=990.0.." shift=none write prot NOESY_1.prot write peaks NOESY_1.peaks