PROGRAM DIRECTED ! READS THE DIRECTED NETWORK FILE !! DOES ANNEALING TO PUT INTO nlev LAYERS MINIMIZING THE NUMBER OF !! COUNTER-HIERARCHICAL ('BAD') LINKS !! needs net.dat file which has the number of interacting pairs on the first !! and then each line consiting of an interecting pair A->B of characters !! distinuishable by the first 7 elements (this could be adjusted by changing !! all character lenght announcement. Examlpe of net.dat file !! 2 !! trg 444 !! trg edf !! The output consists of the file order.txt with a list of nodes and their !! level assignments and a list of 'bad' links on the screen. (could be !! modified). !! For beginning, only number of layers (or levels) NLEV needs to set to the !! desired one. AC could be set larger than 1 if several distinct annealing !! runs are necessary for statistics (for example, to determine the minimum !! number of bad links among several realizations, or how often a particular !! lins is 'bad'). Nrec could be set larger than 0 to randomize the network, !! keeping the degree of each vertex unchanged. !! MAximum number of iterations itmax and number of iterations between !! annealing or measurement steps is set conservatively and may be reduced if !! faster performance is desired. Initial annealing temperature TE is also !! selected impirically after many tests on different network and hardly !! needs to be adjusted IMPLICIT NONE INTEGER :: IWARM,NB,KB,IP,JB,NREC,key,NP,J1,J2,NLEV,ITMAX,IMEAS,IT INTEGER :: nbadmin,NBAD,NBADN,p1,p2,vp,TL,IC,AC,TTL,TTLN,SL INTEGER, ALLOCATABLE :: I1(:),I2(:),UP(:),NNU(:),NND(:),NEIGHU(:,:) INTEGER, ALLOCATABLE :: NEIGHD(:,:),BADL(:) CHARACTER Q1*7,Q2*7 CHARACTER, ALLOCATABLE :: M1*7(:),M2*7(:),PROT*7(:) double precision ran3,te,prob,ETL logical yes ! setting the parameters NLEV=10 ! NUMBER OF HIERARCHICAL LEVELS AC=1 ! NUMBER OF repeated ANNEALING CYCLES FOR STATISTICS TE=1.5d0 ! Initial temperature for annealing (in the number of bad links) ETL=2.D-2 ! ENERGETIC COST OF LONG good LINKS (IN THE NUMBER OF BAD LINKS) ITMAX=40000 ! MAXIMUM NUMBER OF ITERATIONS PER PROTEIN IMEAS=ITMAX/200 ! ANNEALING AND OUTPUT ITMEAS nrec=0 ! number of reconnecions of each bond in trunk reconnection INQUIRE(FILE='rand.dat', EXIST=yes) ! seed for Random number generator if(yes) then open(40,file='rand.dat') read(40,*)iwarm close (40) else iwarm=-9846 end if Open (10, file='net.dat',status='old') ! READING PAIRS OF INTERACTIONS READ(10,*) nb !first line in net.dat file gives the max. number of ALLOCATE(i1(nb),i2(nb)) ! pairs for memory allocation ALLOCATE(UP(nb),M1(NB),M2(NB),prot(2*NB)) jb=0 do kb=1,nb read (10,*,ERR=44,END=55) Q1, Q2 ! reads characters of the length 7 if(Q1.eq.Q2) GO TO 44 key=1 do ip=1,jb ! Purifying: throwing out identical pairs if(Q1.eq.M1(ip).and.Q2.eq.M2(ip)) key=0 end do if(key.eq.1) then ! adding new bond jb=jb+1 M1(jb)=Q1 M2(jb)=Q2 end if 44 end do 55 NB=JB ! TOTAL NUMBER OF DISTINCT EDGES WITHOUT DIMERS close(10) NP=0 do kb=1,nb ! RECOGNIZING LINKED PROTEINS AND LABELING THEM BY NUMBERS Q1=M1(kb) Q2=M2(kb) key=0 ! first vertex of the edge do ip=1,np if (Q1.eq.prot(ip)) then ! the protein on the list key=1 j1=ip end if end do if(key.eq.0) then ! no protein on the list np=np+1 PROT(NP)=Q1 j1=np end if key=0 ! second vertex of the edge do ip=1,np if (Q2.eq.prot(ip)) then ! the protein on the list key=1 j2=ip end if end do if(key.eq.0) then ! no protein on the list np=np+1 PROT(NP)=Q2 j2=np end if i1(kb)=j1 i2(kb)=j2 end do IF(nrec.gt.0) then ! RANDOMIZING NETWORK PRESERVING DEGREES CALL trunk(nb,i1,i2,nrec,iwarm) ! FOR TRUNK RECONNECTION open(40,file='reconnected.dat') do kb=1,nb write(40,*) i1(kb),i2(kb) end do close (40) END IF open(40,file='rand.dat') write(40,*) -nint(1.e6*ran3(iwarm)) close (40) ! finding the nearest neighbours of each protein ALLOCATE(neighd(np,np),nnd(np)) ALLOCATE(neighu(np,np),nnu(np)) NNU=0 NND=0 do kb=1,nb p1=i1(kb) p2=i2(kb) nnd(p1)=nnd(p1)+1 nnu(p2)=nnu(p2)+1 neighd(p1,nnd(p1))=p2 neighu(p2,nnu(p2))=p1 end do write(*,*) 'end of sorting', ' nlinks',nb,'nprot',np open(15, file='netn.dat') ! writing out the list of links in write(15,*)nb,'links' ! protein label notation do kb=1,nb write(15,*)i1(kb),i2(kb) end do close(15) ALLOCATE(BADL(NB)) ! BAD (COUNTER-HIERARCHICAL) LINKS BADL=0 DO IC=1,AC ! CYCLE ON ANNEALINGS (IF WE NEED STATISTICS FOR, SAY, NUMBER OF ! BAD LINKS TO CHOOSE THE MINIMUN DO IP=1,NP ! RANDOM SEEDING ON LEVELS UP(IP)=RAN3(IWARM)*NLEV+1 END DO IT=0 !!! ANNEALING !!!!! nbadmin=nb 10 IP=int(ran3(iwarm)*np)+1 VP=UP(IP) IF (VP.NE.0.AND.VP.NE.NLEV+1) THEN ! ATTEMPTING MOVE NBAD=0 NBADN=0 TTL=0 TTLN=0 TL=int(ran3(iwarm)*NLEV)+1 ! NEW LEVEL DO J1=1,NNU(IP) SL=UP(NEIGHU(IP,J1)) IF(SL.GE.VP) then NBAD=NBAD+1 ELSE IF(SL.LT.VP) then TTL=TTL+VP-SL END IF IF(SL.GE.TL) then NBADN=NBADN+1 ELSE IF(SL.LT.TL) then TTLN=TTLN+TL-SL END IF END DO DO J1=1,NND(IP) SL=UP(NEIGHD(IP,J1)) IF(SL.LE.VP) THEN NBAD=NBAD+1 ELSE IF(SL.GT.VP) THEN TTL=TTL-VP+SL END IF IF(SL.LE.TL) THEN NBADN=NBADN+1 ELSE IF(SL.GT.TL) THEN TTLN=TTLN-TL+SL END IF END DO if(dexp(((nbad-nbadn) + (TTL-TTLN)*ETL)/te).GE.ran3(iwarm)) UP(IP)=TL IT=IT+1 ! ITERATION NUMBER INCREMENTED IF(IT/(IMEAS*NP)*(IMEAS*NP).EQ.IT) THEN ! MEASUREMENT TE=TE*0.95 ! ANNEALING NBAD=0 TTL=0 DO KB=1,NB IF(UP(I1(KB)).GE.UP(I2(KB))) THEN NBAD=NBAD+1 ELSE IF(UP(I1(KB)).LT.UP(I2(KB))) THEN TTL=TTL+UP(I2(KB))-UP(I1(KB)) END IF END DO nbadmin=min(nbadmin,nbad) WRITE(*,*)'ITERATION',IT,' BAD LINKS',NBAD,'LENGTH',TTL END IF END IF IF(IT.LT.ITMAX*NP) GO TO 10 !!!!!!!!!! OUTPUT FOR EACH RUN !!!!!!!!!!!!!!!!!!!!!!!!!!!! open(unit=10,file='order.txt') do ip=1,np write(10,*)ip,prot(ip),up(ip) !Prot. index,name,level(between 1 and NLEV) end do write(*,*)'N_badlinks', nbadmin, 'LENGTH',TTL! THE NUMBER OF BAD LINKS close(10) DO KB=1,NB IF(UP(I1(KB)).GE.UP(I2(KB))) THEN write(*,*)kb,prot(i1(kb)),prot(i2(kb)) ! WRITING EACH BAD LINK BADL(KB)=BADL(KB)+1 ! HOW MANY TIMES A LINK IS BAD END IF END DO END DO DO KB=1,NB ! OUTPUTTING THE STATISTICS ON HOW OFTEN A LINK IS BAD. !IF(BADL(KB).GE.AC/4) write(*,*)FLOAT(BADL(KB))/AC,kb,prot(i1(kb)),prot(i2(kb)) END DO END PROGRAM DIRECTED SUBROUTINE trunk(nb,i1,i2,nrec,iwarm) ! NETWORK RANDOMIZATION IMPLICIT NONE INTEGER NB,I1(NB),I2(NB),nrec,k1,k2,i,iwarm,temp,j,key double precision ran3 do i=1,nrec*nb k1=int(ran3(iwarm)*(nb+1)) k2=int(ran3(iwarm)*(nb+1)) if(k1*k2.eq.0.or.k1.eq.k2) go to 10 if(i1(k1).eq.i2(k2).or.i2(k1).eq.i1(k2)) go to 10 ! checking if we create a doublelink key=0 do j=1,nb if(i1(k1).eq.i1(j).and.i2(k2).eq.i2(j))key=1 if(i2(k1).eq.i2(j).and.i1(k2).eq.i1(j))key=1 end do if(key.eq.1) go to 10 temp=i1(k1) i1(k1)=i1(k2) i1(k2)=temp 10 continue end do ! check for the number of doublelinks END SUBROUTINE trunk ! ran3 from Numerical Recipes, 2nd edition !-------------------------------------------------------------------------- function ran3(idum) !========================================================================== ! implicit real*4(m) integer idum integer mbig,seed,mz double precision ran3,fac parameter (mbig=1000000000,mseed=161803398,mz=0,fac=1./mbig) integer i,iff,ii,inext,inextp,k integer mj,mk,ma(55) save iff,inext,inextp,ma data iff /0/ 1 if(idum.lt.0.or.iff.eq.0)then iff=1 mj=mseed-iabs(idum) mj=mod(mj,mbig) ma(55)=mj mk=1 do 11 i=1,54 ii=mod(21*i,55) ma(ii)=mk mk=mj-mk if(mk.lt.mz)mk=mk+mbig mj=ma(ii) 11 continue do 13 k=1,4 do 12 i=1,55 ma(i)=ma(i)-ma(1+mod(i+30,55)) if(ma(i).lt.mz)ma(i)=ma(i)+mbig 12 continue 13 continue inext=0 inextp=31 idum=1 endif inext=inext+1 if(inext.eq.56)inext=1 inextp=inextp+1 if(inextp.eq.56)inextp=1 mj=ma(inext)-ma(inextp) if(mj.lt.mz)mj=mj+mbig ma(inext)=mj ran3=mj*fac if(ran3.le.0.or.ran3.ge.1) goto 1 return end