gde_linux/ZUKER/rna.f

354 lines
12 KiB
FortranFixed
Raw Permalink Normal View History

2022-03-08 04:43:05 +08:00
c MFOLD - Prediction of RNA secondary structure by free energy
c minimization.
c - Version 2.0
c - Michael Zuker and John Jaeger
c - LRNA : folds linear RNA sequences
c - CRNA : folds circular RNA sequences
c
c The original version (1.0) was designed by Michael Zuker and
c programmed by Eric Nelson in the summer of 1987 in the Division of
c Biological Sciences at the National Research Council of Canada. John
c Jaeger added the tetraloop bonus energy feature and created the
c BATGEN program for batch file generation.
c
c Version 2.0 corrects a number of small bugs from the original
c program. These were added to version 1 and itemized in the
c ERRATA.LIST file that was distributed along with version 1. The major
c improvements of version 2 are : 1. During the generation of
c suboptimal foldings, the number of new base pairs that are
c sufficiently different from base pairs that have already been found
c must be greater than the WINDOW parameter. This feature was added
c during the summer of 1989, and was made part of version 1 (item 11 in
c the ERRATA.LIST file distributed with this version). The effect is to
c eliminate structures that contain just a few new base pairs. 2.
c Temperature dependent folding. This was added in the fall of 1989 and
c was never a feature of version 1.
c
c METHOD : A dynamic programming algorithm is used to find optimal and
c suboptimal foldings of an RNA molecule starting from linear sequence
c data. Auxiliary information can be used to constrain the folding.
c
c NB : Base pairs are forced by giving them a bonus energy (EPARAM(9)
c in the program code). These energies are subtracted during the
c traceback algorithm so that the computed sturctures have the correct
c energies. Unfortunately, there is no way to subtract the bonus
c energies from the energy dot plots. Moreover, each forced base pair
c contain two bonus energies because of the nature of the algorithm.
c For example, suppose that an optimal folding of an RNA contains 3
c forced base pairs ( default bonus energy is 50.0 kcal per forced base
c pair ) and that the correct folding energy is -180.0 kcal/mole.
c Internally, the energy will be -180.0 - (3+1) x 50.0 = -380.0
c kcal/mole. To find foldings within 10% of the correct energy, one
c needs to compute foldings to within 18.0 kcal of -180.0 - 3 x 50.0 =
c -330.0 kcal/mole. This comes out to -312.0 kcal/mole. The ratio of
c -312.0 to -380.0 is 82%, so that one would request the 18% level of
c suboptimality! This confustion only exists when base pairs are
c forced. Each closed excision counts as one forced base pair.
c
c Energy data from :
c S.M. Freier et al., Proc. Natl. Acad. Sci. USA, 83, 9373-9377, 1986.
c D.H. Turner et al., Cold Spring Harbor Symposia on Quantitative Biology,
c 52, 123-133, 1987.
c D.H. Turner et al., Annu. Rev. Biophys. Biophys. Chem 17, 167-192 (1988).
c This last reference has all the dangling end and terminal mismatch data.
c
c References :
c M. Zuker
c On Finding All Suboptimal Foldings of an RNA Molecule.
c Science, 244, 48-52, (1989)
c
c J. A. Jaeger, D. H. Turner and M. Zuker
c Improved Predictions of Secondary Structures for RNA.
c Proc. Natl. Acad. Sci. USA, BIOCHEMISTRY, 86, 7706-7710, (1989)
c
c J. A. Jaeger, D. H. Turner and M. Zuker
c Predicting Optimal and Suboptimal Secondary Structure for RNA.
c in "Molecular Evolution: Computer Analysis of Protein and
c Nucleic Acid Sequences", R. F. Doolittle ed.
c Methods in Enzymology, 183, 281-306 (1989)
c
include 'rfd.inc'
real energy
logical flag,mark
c Fill screen with author and reference data.
call begin
c Initial setup for run.
5 call device
if (cntrl(1).ne.2) then
c Read energy information if this is not a continuation run.
call enefiles
call erg(1,0,0,0,0)
else
c dump out information read in from continuation file
call cdump
endif
c Determine output specifications if this is not a save run.
if (cntrl(1).ne.1) call outputs
c Call the menu if this is not a continuation run.
if (cntrl(1).ne.2) call menu
mrep = 1
c CNTRL(7) = 0 - suboptimal dot plot
c 1 - N best sorted by energy
c 2 - best folding for all sequences in a file
10 if (cntrl(7).eq.2) call mseq(mrep)
c Process sequence before folding.
do i = 1,mxbits
marks(i) = 0
force2(i) = 0
enddo
call process
if (n*2.gt.fldmax) then
tt = fldmax/2
c Fragment is too long. Try again.
write (6,*) 'Segment larger than ',tt
stop
endif
if (cntrl(1).ne.2) then
c Fill the optimal energy arrays except in a continuation run.
call fill
endif
if (cntrl(1).eq.1) then
c Save the results from FILL in a SAVE run and then stop.
call putcont
stop
endif
rep = 1
jump = 1
flag = .true.
err = 0
do while (flag)
if (cntrl(7).eq.0) then
c Interactive dot plot returns IRET, JRET (new numbering).
c Zuker comments out call to dotplt : do not choose this option
C If you do, the program stops dead here.
c call dotplt(iret,jret,jump)
if(1.eq.1) stop ' Energy dot plot disabled.'
jump = 2
else
c Automatic sort returns IRET,JRET (new numbering).
print *,'traceback'
call sortout(iret,jret,rep,err)
if (err.eq.30) then
flag = .false.
call errmsg(err,rep-1,0)
err = 0
endif
rep = rep + 1
endif
c First traceback yields the best structure on the included fragment
c from IRET to JRET.
if (flag) call trace(iret,jret,nforc1,err)
if (err.ne.0) call errmsg(err,iret,jret)
if (flag) then
it = iret+n
c Second traceback yields the best structure on the excluded fragment
c from IRET to JRET.
call trace(jret,it,nforc2,err)
if (err.ne.0) then
call errmsg(err,jret,it)
else
c The energy of the best structure containing the base-pair IRET,
c JRET is the sum of the energies of the optimal foldings on
c the included and excluded fragments. A correction is made for
c forced base-pairs.
c* CALL EFN(ENE,1,N)
c* WRITE (6,*) 'NEW ENERGY ',ENE
ene = v(iret,jret) + v(jret,iret+n) - eparam(9) * (nforc1+nforc2)
energy = float(ene) / 10.0
c Count the number of new base pairs not within WINDOW
c of existing base pairs.
numbp = 0
do k = 1,n
if(k.lt.basepr(k)) then
if(.not.mark(k,basepr(k))) numbp = numbp + 1
endif
enddo
do k = 1,n
if(k.lt.basepr(k)) then
c Mark "traced-back" base pairs and also base-pairs
c which are close (within WINDOW = CNTRL(9) ).
call smark(k,basepr(k))
if(cntrl(9).gt.0) then
do k1 = -cntrl(9),cntrl(9)
do k2 = -cntrl(9),cntrl(9)
if(k+k1.gt.0.and.k+k1.lt.basepr(k)+k2.and.
1 basepr(k)+k2.le.n) call smark(k+k1,basepr(k)+k2)
enddo
enddo
endif
endif
enddo
if(numbp.le.cntrl(9)) then
rep = rep - 1
go to 900
endif
write (6,1010) rep - 1
c 1010 format('+',i5$)
1010 format('+',i5)
if (cntrl(2).ne.2.and.cntrl(2).ne.3.and.cntrl(2).ne.6) then
c Line printer output.
call linout(1,n,energy,iret,jret,err)
endif
if (err.ne.0) then
call errmsg(err,iret,jret)
else
if (cntrl(2).ge.3.and.cntrl(2).ne.4) then
c Region table output.
call regtab
endif
if (mod(cntrl(2),2).eq.0.or.cntrl(2).eq.7) then
c CT file output.
call ct(energy)
endif
endif
endif
endif
if (cntrl(7).eq.1.and.rep.gt.cntrl(6)) flag = .false.
900 continue
enddo
c
c Multiple sequence option (CNTRL(7) = 2)
c If sequence number (MREP) is < total number of sequences
c (CNTRL(5)), go get another sequence.
c
if (cntrl(7).eq.2.and.mrep.lt.cntrl(5)) then
mrep = mrep + 1
goto 10
endif
stop
end
c Marks a base-pair I,J.
c Assumes that 1 <= I <= J <= N.
c The information is stored in a single bit in the MARKS
c array.
c The conversion from double dimension to single is through the
c transformation I,J ==> (J-1)*J/2 + I .
subroutine smark(i,j)
include 'rfd.inc'
integer*2 bit
posn = (((j-1)*j)/2) + i
word = (posn+15) / 16
bit = mod(posn,16)
c marks(word) = iibset(marks(word),bit)
marks(word) = ibset(marks(word),bit)
return
end
c Marks a forced base-pair I,J.
c The incoming base-pair II,JI is processed to an I,J
c base-pair satisfying 1 <= I <= J <= N.
c The information is stored in a single bit in the FORCE2
c array.
c The conversion from double dimension to single is through the
c transformation I,J ==> (J-1)*J/2 + I .
subroutine sfce(ii,ji)
include 'rfd.inc'
integer*2 bit
if (ii.gt.n) then
i = ii - n
j = ji - n
elseif (ji.gt.n) then
i = ji-n
j = ii
else
i = ii
j = ji
endif
posn = (((j-1)*j)/2) + i
word = (posn+15) / 16
bit = mod(posn,16)
c force2(word) = iibset(force2(word),bit)
force2(word) = ibset(force2(word),bit)
return
end
c Retrieves information on whether or not the base-pair I,J
c has been marked by a traceback passing through or close to
c this pair.
logical function mark(i,j)
include 'rfd.inc'
integer*2 bit,one
one = 1
posn = (((j-1)*j)/2) + i
word = (posn+15) / 16
bit = mod(posn,16)
c set = iibits(marks(word),bit,one)
set = ibits(marks(word),bit,one)
mark = .false.
if (set.ne.0) mark = .true.
return
end
c Retrieves information on whether or not the base-pair I,J
c has been forced.
logical function fce(ii,ji)
include 'rfd.inc'
integer*2 bit,one
if (ii.gt.n) then
i = ii - n
j = ji - n
elseif (ji.gt.n) then
i = ji-n
j = ii
else
i = ii
j = ji
endif
one = 1
posn = (((j-1)*j)/2) + i
word = (posn+15) / 16
bit = mod(posn,16)
c set = iibits(force2(word),bit,one)
set = ibits(force2(word),bit,one)
fce = .false.
if (set.ne.0) fce = .true.
return
end
c fills screen with author and reference information
subroutine begin
character*1 ans
character*80 record
open(3,file='begin.dat',status='old',err=5)
write(6,1010)
1010 format(' ')
1 read(3,1020,end=2) record
1020 format(a80)
write(6,1030) record
c1030 format(' ',a80)
1030 format(a80)
go to 1
2 write(6,1040)
c 1040 format(' Press <return> to continue ...'$)
c 1040 format(' Press <return> to continue ...'$)
1040 format(' Press <return> to continue ...',$)
read(5,1050) ans
1050 format(a1)
return
5 write(6,1060)
1060 format(//' Author and reference file not available.'//)
c
c C.Wang copied the next 2 lines here. 6/12/91
c to keep consistency.
3 write(6,1040)
read(5,1050) ans
return
end