subroutine sbspot(m,ndim,nc,nr,mm,mp,nmed,nopen,nseg,n5,n6,nest, * ma,mma,mw,z) c c Subroutine to implement methods in: "Combinatorial image analysis of c DNA microarray features" by C.A. Glasbey and P. Ghazal, for data used c in this paper. c c input: c m integer array, ndim by nr, first nc by nr elements are pixel c values in channel 1 of microarray c ndim integer first dimension of m and mm c nc integer number of columns in m and mm containing microarray values c nr integer number of rows in m and mm containing microarray values c mm as for m, but for channel 2 of microarray c mp integer array, 2 by 1080, of coordinates of approximate centres c of 1080 microarray spots c nmed integer size of median filter (0 for no filter) c nopen integer size of top-hat filter (>=15 for no filter) c nseg integer choice of segmentation method c 1 fixed-radius c 2 fixed-proportion c 3 k-means clustering c n5 integer, specifying RS if nseg=1, pS if nseg=2 c n6 integer, specifying RB if nseg=1, pB if nseg=2 c nest integer choice of estimation method c 2 square-root transformation, no background subtraction c 4 original scale, with background subtraction c 6 square-root transformation, no background subtraction c 8 original scale, with background subtraction c (options 1,3,5,7 use sum of pixels instead of mean of pixels, c and are not reported in the paper) c c workspace: c ma integer array, ndim by nr c mma integer array, ndim by nr c mw integer array, ndim by nr c c output: c z real array, 3 by 1080, of estimated log-transformed spot c intensities in channel 1, channel 2, and ratio of channels c dimension m(ndim,nr),ma(ndim,nr), * mm(ndim,nr),mma(ndim,nr),mw(ndim,nr),mp(2,1080),z(3,1080) c c median filter c call sbmedi(m,ndim,nc,nr,nmed,nmed,ma,ndim,mw,65535) call sbmedi(mm,ndim,nc,nr,nmed,nmed,mma,ndim,mw,65535) c c top-hat filter c if (nopen.lt.15) then call sbmin(ma,ndim,nc,nr,nopen,nopen,m,ndim,mw) call sbmax(m,ndim,nc,nr,nopen,nopen,mm,ndim,mw) do 11 i=1,nr do 11 j=1,nc 11 ma(j,i)=ma(j,i)-mm(j,i) call sbmin(mma,ndim,nc,nr,nopen,nopen,m,ndim,mw) call sbmax(m,ndim,nc,nr,nopen,nopen,mm,ndim,mw) do 14 i=1,nr do 14 j=1,nc 14 mma(j,i)=mma(j,i)-mm(j,i) endif c c segmentation and estimation c if (nseg.eq.1) then call sbrad(ma,ndim,nc,nr,mma,m,mm,mw,mp,n5,n6,nest,z) elseif (nseg.eq.2) then call sbprop(ma,ndim,nc,nr,mma,mp,n5,n6,nest,z) elseif (nseg.eq.3) then call sbclus(ma,ndim,nc,nr,mma,mp,nest,z) endif c return end ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*cc subroutine sbrad(ma,ndim,nc,nr,mma,m,mm,mw,mp,nradf,nradb,nest,z) dimension m(ndim,nr),ma(ndim,nr), * mm(ndim,nr),mma(ndim,nr),mw(ndim,nr),mp(2,1080),zz(2,5,1080), * mc(2,1080),z(3,1080) c c fixed-radius segmentation c do 13 i=1,nr do 13 j=1,nc 13 m(j,i)=ma(j,i)+mma(j,i) nmav=10 call sbmav(m,ndim,nc,nr,nmav,nmav,mm,ndim,mw) c do 6 kk=1,1080 nmax=0 nbox=5 do 6 i=mp(1,kk)-nbox,mp(1,kk)+nbox do 6 j=mp(2,kk)-nbox,mp(2,kk)+nbox if (mm(j,i).gt.nmax) then nmax=mm(j,i) mc(1,kk)=i mc(2,kk)=j endif 6 continue c nradf2=nradf*(nradf+1) nradb2=nradb*(nradb+1) c do 2 kk=1,1080 bb1=0 bb2=0 bbb1=0 bbb2=0 nbb=0 ff1=0 ff2=0 fff1=0 fff2=0 nff=0 do 3 i=mp(1,kk)-20,mp(1,kk)+20 do 3 j=mp(2,kk)-20,mp(2,kk)+20 if ((i-mc(1,kk))**2+(j-mc(2,kk))**2.gt.nradb2) then bb1=bb1+sqrt(float(max(ma(j,i),0))) bb2=bb2+sqrt(float(max(mma(j,i),0))) bbb1=bbb1+ma(j,i) bbb2=bbb2+mma(j,i) nbb=nbb+1 elseif ((i-mc(1,kk))**2+(j-mc(2,kk))**2.le.nradf2) then ff1=ff1+sqrt(float(max(ma(j,i),0))) ff2=ff2+sqrt(float(max(mma(j,i),0))) fff1=fff1+ma(j,i) fff2=fff2+mma(j,i) nff=nff+1 endif 3 continue c zz(1,1,kk)=nff zz(1,2,kk)=ff1 zz(1,3,kk)=ff2 zz(1,4,kk)=fff1 zz(1,5,kk)=fff2 zz(2,1,kk)=nbb zz(2,2,kk)=bb1 zz(2,3,kk)=bb2 zz(2,4,kk)=bbb1 zz(2,5,kk)=bbb2 2 continue c 12 call sbestf(zz,nest,z) return end ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*cc subroutine sbprop(ma,ndim,nc,nr,mma,mp,npropf,npropb,nest,z) dimension ma(ndim,nr),mma(ndim,nr),mp(2,1080),zz(2,5,1080), * m1(1681),m2(1681),zzz(35,2,5,1080),z(3,1080) c c fixed-proportion segmentation c n=1681 kkk=0 ntf=int(1.5+n*float(npropf)/100) ntb=int(1.5+n*float(npropb)/100) kkk=kkk+1 c do 2 kk=1,1080 c k=0 do 6 i=mp(1,kk)-20,mp(1,kk)+20 do 6 j=mp(2,kk)-20,mp(2,kk)+20 k=k+1 m1(k)=ma(j,i) m2(k)=mma(j,i) 6 continue c call sbsrti(m1,n,1,n) call sbsrti(m2,n,1,n) c bb1=0 bb2=0 bbb1=0 bbb2=0 nbb=0 ff1=0 ff2=0 fff1=0 fff2=0 nff=0 do 3 k=1,n if (k.lt.ntb) then bb1=bb1+sqrt(float(max(m1(k),0))) bb2=bb2+sqrt(float(max(m2(k),0))) bbb1=bbb1+m1(k) bbb2=bbb2+m2(k) nbb=nbb+1 elseif (k.ge.ntf) then ff1=ff1+sqrt(float(max(m1(k),0))) ff2=ff2+sqrt(float(max(m2(k),0))) fff1=fff1+m1(k) fff2=fff2+m2(k) nff=nff+1 endif 3 continue c zzz(kkk,1,1,kk)=nff zzz(kkk,1,2,kk)=ff1 zzz(kkk,1,3,kk)=ff2 zzz(kkk,1,4,kk)=fff1 zzz(kkk,1,5,kk)=fff2 zzz(kkk,2,1,kk)=nbb zzz(kkk,2,2,kk)=bb1 zzz(kkk,2,3,kk)=bb2 zzz(kkk,2,4,kk)=bbb1 zzz(kkk,2,5,kk)=bbb2 2 continue c kkk=0 kkk=kkk+1 do 14 i=1,2 do 14 j=1,5 do 14 k=1,1080 14 zz(i,j,k)=zzz(kkk,i,j,k) 13 call sbestf(zz,nest,z) return end ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*cc subroutine sbclus(ma,ndim,nc,nr,mma,mp,nest,z) dimension ma(ndim,nr),mma(ndim,nr),mp(2,1080),zz(2,5,1080), * y1(1681),y2(1681),z(3,1080) c c k-means clustering c n=1681 nb=nint(n*0.4) nf=nint(n*0.9) do 2 kk=1,1080 c k=0 do 6 i=mp(1,kk)-20,mp(1,kk)+20 do 6 j=mp(2,kk)-20,mp(2,kk)+20 k=k+1 y1(k)=sqrt(float(max(ma(j,i),0))) y2(k)=sqrt(float(max(mma(j,i),0))) 6 continue c call sbsrtr(y1,n,1,n) call sbsrtr(y2,n,1,n) c b1=y1(nb) b2=y2(nb) f1=y1(nf) f2=y2(nf) do 7 iter=1,100 c bb1=0 bb2=0 nbb=0 ff1=0 ff2=0 nff=0 do 4 k=1,n if ((y1(k)-b1)**2+(y2(k)-b2)**2.lt. * (y1(k)-f1)**2+(y2(k)-f2)**2) then bb1=bb1+y1(k) bb2=bb2+y2(k) nbb=nbb+1 else ff1=ff1+y1(k) ff2=ff2+y2(k) nff=nff+1 endif 4 continue c if (abs(b1-bb1/nbb)+abs(b2-bb2/nbb)+ * abs(f1-ff1/nff)+abs(f2-ff2/nff).lt.1.0d-5) go to 5 b1=bb1/nbb b2=bb2/nbb f1=ff1/nff f2=ff2/nff 7 continue write (6,*) ' *** k-means failed to converge, for spot ',kk 5 continue c if (b1*b2.gt.f1*f2) then x=f1 f1=b1 b1=x x=f2 f2=b2 b2=x endif c bb1=0 bb2=0 bbb1=0 bbb2=0 nbb=0 ff1=0 ff2=0 fff1=0 fff2=0 nff=0 do 8 k=1,n if ((y1(k)-b1)**2+(y2(k)-b2)**2.lt. * (y1(k)-f1)**2+(y2(k)-f2)**2) then bb1=bb1+y1(k) bb2=bb2+y2(k) bbb1=bbb1+y1(k)**2 bbb2=bbb2+y2(k)**2 nbb=nbb+1 else ff1=ff1+y1(k) ff2=ff2+y2(k) fff1=fff1+y1(k)**2 fff2=fff2+y2(k)**2 nff=nff+1 endif 8 continue c zz(1,1,kk)=nff zz(1,2,kk)=ff1 zz(1,3,kk)=ff2 zz(1,4,kk)=fff1 zz(1,5,kk)=fff2 zz(2,1,kk)=nbb zz(2,2,kk)=bb1 zz(2,3,kk)=bb2 zz(2,4,kk)=bbb1 zz(2,5,kk)=bbb2 2 continue c call sbestf(zz,nest,z) return end ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*cc subroutine sbestf(zz,nest,z) dimension zz(2,5,1080),z(3,1080) c c estimation c do 13 k=1,1080 if (nest.eq.1) then z(1,k)=funlog(zz(1,2,k)) z(2,k)=funlog(zz(1,3,k)) z(3,k)=funlog(zz(1,2,k)/zz(1,3,k)) elseif (nest.eq.2) then z(1,k)=funlog(zz(1,2,k)/zz(1,1,k)) z(2,k)=funlog(zz(1,3,k)/zz(1,1,k)) z(3,k)=funlog(zz(1,2,k)/zz(1,3,k)) elseif (nest.eq.3) then z(1,k)=funlog(zz(1,4,k)) z(2,k)=funlog(zz(1,5,k)) z(3,k)=funlog(zz(1,4,k)/zz(1,5,k)) elseif (nest.eq.4) then z(1,k)=funlog(zz(1,4,k)/zz(1,1,k)) z(2,k)=funlog(zz(1,5,k)/zz(1,1,k)) z(3,k)=funlog(zz(1,4,k)/zz(1,5,k)) elseif (nest.eq.5) then z(1,k)=funlog(zz(1,2,k)-zz(1,1,k)*zz(2,2,k)/zz(2,1,k)) z(2,k)=funlog(zz(1,3,k)-zz(1,1,k)*zz(2,3,k)/zz(2,1,k)) z(3,k)=funlog((zz(1,2,k)-zz(1,1,k)*zz(2,2,k)/zz(2,1,k))/ * (zz(1,3,k)-zz(1,1,k)*zz(2,3,k)/zz(2,1,k))) elseif (nest.eq.6) then z(1,k)=funlog(zz(1,2,k)/zz(1,1,k)-zz(2,2,k)/zz(2,1,k)) z(2,k)=funlog(zz(1,3,k)/zz(1,1,k)-zz(2,3,k)/zz(2,1,k)) z(3,k)=funlog((zz(1,2,k)/zz(1,1,k)-zz(2,2,k)/zz(2,1,k))/ * (zz(1,3,k)/zz(1,1,k)-zz(2,3,k)/zz(2,1,k))) elseif (nest.eq.7) then z(1,k)=funlog(zz(1,4,k)-zz(1,1,k)*zz(2,4,k)/zz(2,1,k)) z(2,k)=funlog(zz(1,5,k)-zz(1,1,k)*zz(2,5,k)/zz(2,1,k)) z(3,k)=funlog((zz(1,4,k)-zz(1,1,k)*zz(2,4,k)/zz(2,1,k))/ * (zz(1,5,k)-zz(1,1,k)*zz(2,5,k)/zz(2,1,k))) elseif (nest.eq.8) then z(1,k)=funlog(zz(1,4,k)/zz(1,1,k)-zz(2,4,k)/zz(2,1,k)) z(2,k)=funlog(zz(1,5,k)/zz(1,1,k)-zz(2,5,k)/zz(2,1,k)) z(3,k)=funlog((zz(1,4,k)/zz(1,1,k)-zz(2,4,k)/zz(2,1,k))/ * (zz(1,5,k)/zz(1,1,k)-zz(2,5,k)/zz(2,1,k))) endif 13 continue c return end ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*cc real function funlog(x) if (x.lt.5.0e-5) then funlog=-10.0 else funlog=log(x) endif return end ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*cc subroutine sbmav(m,ndim,nc,nr,ncc,nrr,mm,nmdim,mw) dimension m(ndim,nr),mm(nmdim,nr),mw(nc,nr) c c moving average filter of size (2 ncc + 1) x (2 nrr + 1) for integer imwge c nk=(2*ncc+1)*(2*nrr+1) do 1 i=1,nr do 1 j=1,nc 1 mm(j,i)=m(j,i) c do 31 j=1,nc nav=0.0 do 32 i=1,2*nrr+1 32 nav=nav+m(j,i) iz=0 mw(j,2*nrr+1)=nav do 31 i=2*nrr+2,nr iz=iz+1 nav=nav-m(j,iz)+m(j,i) 31 mw(j,i)=nav c do 33 i=2*nrr+1,nr nav=0.0 do 34 j=1,2*ncc+1 34 nav=nav+mw(j,i) jz=0 mm(ncc+1,i-nrr)=(nav+nk/2)/nk do 33 j=2*ncc+2,nc jz=jz+1 nav=nav-mw(jz,i)+mw(j,i) 33 mm(j-ncc,i-nrr)=(nav+nk/2)/nk c return end ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*cc subroutine sbmin(m,ndim,nc,nr,ncc,nrr,mm,nmdim,mw) dimension m(ndim,nr),mm(nmdim,nr),mw(nc,nr) c c min filter of size (2 ncc + 1) x (2 nrr + 1) for real image c c along rows c do 2 i=1,nr c nmin=1000000 do 3 j=1,2*ncc+1 nmin=min(nmin,m(j,i)) if (j.gt.ncc) mw(j-ncc,i)=nmin 3 continue c do 2 jj=2,nc-ncc c if (jj+2*ncc.le.nc) nmin=min(nmin,m(jj+2*ncc,i)) if (nmin.eq.m(jj-1,i)) then nmin=1000000 do 4 j=jj,min(jj+2*ncc,nc) 4 nmin=min(nmin,m(j,i)) endif c 2 mw(jj+ncc,i)=nmin c c down columns c do 5 j=1,nc c nmin=1000000 do 6 i=1,2*nrr+1 nmin=min(nmin,mw(j,i)) if (i.gt.nrr) mm(j,i-nrr)=nmin 6 continue c do 5 ii=2,nr-nrr c if (ii+2*nrr.le.nr) nmin=min(nmin,mw(j,ii+2*nrr)) if (nmin.eq.mw(j,ii-1)) then nmin=1000000 do 7 i=ii,min(ii+2*nrr,nr) 7 nmin=min(nmin,mw(j,i)) endif c 5 mm(j,ii+nrr)=nmin c return end ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*cc subroutine sbmax(m,ndim,nc,nr,ncc,nrr,mm,nmdim,mw) dimension m(ndim,nr),mm(nmdim,nr),mw(nc,nr) c c max filter of size (2 ncc + 1) x (2 nrr + 1) for real image c c along rows c do 2 i=1,nr c nmax=-1000000 do 3 j=1,2*ncc+1 nmax=max(nmax,m(j,i)) if (j.gt.ncc) mw(j-ncc,i)=nmax 3 continue c do 2 jj=2,nc-ncc c if (jj+2*ncc.le.nc) nmax=max(nmax,m(jj+2*ncc,i)) if (nmax.eq.m(jj-1,i)) then nmax=-1000000 do 4 j=jj,min(jj+2*ncc,nc) 4 nmax=max(nmax,m(j,i)) endif c 2 mw(jj+ncc,i)=nmax c c down columns c do 5 j=1,nc c nmax=-1000000 do 6 i=1,2*nrr+1 nmax=max(nmax,mw(j,i)) if (i.gt.nrr) mm(j,i-nrr)=nmax 6 continue c do 5 ii=2,nr-nrr c if (ii+2*nrr.le.nr) nmax=max(nmax,mw(j,ii+2*nrr)) if (nmax.eq.mw(j,ii-1)) then nmax=-1000000 do 7 i=ii,min(ii+2*nrr,nr) 7 nmax=max(nmax,mw(j,i)) endif c 5 mm(j,ii+nrr)=nmax c return end ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*cc subroutine sbmedi(m,ndim,nc,nr,ncc,nrr,mm,nmdim,mh,nhdim) dimension m(ndim,nr),mm(nmdim,nr),mh(0:nhdim) c c median filter of size (2 ncc + 1) x (2 nrr + 1) for integer image c do 1 i=1,nr do 1 j=1,nc 1 mm(j,i)=m(j,i) if ((ncc.eq.0).and.(nrr.eq.0)) return c nt=(2*ncc+1)*(2*nrr+1)/2 do 2 ii=1,nr-2*nrr c do 3 k=0,nhdim 3 mh(k)=0 c do 4 i=ii,ii+2*nrr do 4 j=1,2*ncc+1 k=m(j,i) 4 mh(k)=mh(k)+1 c nsum=0 do 5 k=0,nhdim nsum=nsum+mh(k) if (nsum.gt.nt) go to 6 5 continue k=nhdim 6 nmed=k mm(ncc+1,ii+nrr)=nmed nsum=nsum-mh(k) c do 2 jj=2,nc-2*ncc c j=jj-1 ja=jj+2*ncc do 8 i=ii,ii+2*nrr k=m(j,i) mh(k)=mh(k)-1 if (k.lt.nmed) nsum=nsum-1 k=m(ja,i) mh(k)=mh(k)+1 if (k.lt.nmed) nsum=nsum+1 8 continue c if (nsum.le.nt) go to 9 c 10 nmed=nmed-1 nsum=nsum-mh(nmed) if (nsum.le.nt) go to 2 go to 10 c 12 nsum=nsum+mh(nmed) nmed=nmed+1 9 if (nsum+mh(nmed).le.nt) go to 12 c 2 mm(jj+ncc,ii+nrr)=nmed c return end ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*cc RECURSIVE SUBROUTINE sbsrti(my,n,ns,nf) ! Quick sort routine from: ! Brainerd, W.S., Goldberg, C.H. & Adams, J.C. (1990) "Programmer's Guide to ! Fortran 90", McGraw-Hill ISBN 0-07-000248-7, pages 149-150. ! Modified. dimension my(n) IF (nf < ns + 6) THEN DO i = ns, nf - 1 DO j = i+1, nf IF (my(i) > my(j)) THEN mxx = my(i); my(i) = my(j); my(j) = mxx END IF END DO END DO ELSE myy = my((ns + nf)/2) i = ns - 1; j = nf + 1 DO DO i = i + 1 IF (my(i) >= myy) EXIT END DO DO j = j - 1 IF (my(j) <= myy) EXIT END DO IF (i < j) THEN mxx = my(i); my(i) = my(j); my(j) = mxx ELSE IF (i == j) THEN i = i + 1 EXIT ELSE EXIT END IF END DO IF (ns < j) CALL sbsrti(my,n,ns,j) IF (i < nf) CALL sbsrti(my,n,i,nf) END IF END SUBROUTINE sbsrti ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*ccccccccc*cc RECURSIVE SUBROUTINE sbsrtr(y,n,ns,nf) ! Quick sort routine from: ! Brainerd, W.S., Goldberg, C.H. & Adams, J.C. (1990) "Programmer's Guide to ! Fortran 90", McGraw-Hill ISBN 0-07-000248-7, pages 149-150. ! Modified. dimension y(n) IF (nf < ns + 6) THEN DO i = ns, nf - 1 DO j = i+1, nf IF (y(i) > y(j)) THEN xx = y(i); y(i) = y(j); y(j) = xx END IF END DO END DO ELSE yy = y((ns + nf)/2) i = ns - 1; j = nf + 1 DO DO i = i + 1 IF (y(i) >= yy) EXIT END DO DO j = j - 1 IF (y(j) <= yy) EXIT END DO IF (i < j) THEN xx = y(i); y(i) = y(j); y(j) = xx ELSE IF (i == j) THEN i = i + 1 EXIT ELSE EXIT END IF END DO IF (ns < j) CALL sbsrtr(y,n,ns,j) IF (i < nf) CALL sbsrtr(y,n,i,nf) END IF END SUBROUTINE sbsrtr