!***************************************************************************

!     Jun 2008

!     Computational identification of Long Homogeneous Genome Regions (LHGRs)

!     IsoFinder

!     Copyright 2003-2008 Pedro Bernaola-Galvan, Pedro Carpena & Jose L. Oliver
!     This program is distributed under the terms of the GNU General Public License (GPL)

!     This file is part of IsoFinder

!     IsoFinder is free software: you can redistribute it and/or modify
!     it under the terms of the GNU General Public License as published by
!     the Free Software Foundation, either version 3 of the License, or
!     (at your option) any later version.

!     IsoFinder is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!     GNU General Public License for more details.

!     You should have received a copy of the GNU General Public License
!     along with IsoFinder.  If not, see <http://www.gnu.org/licenses/>.

!***************************************************************************

!     References:

!     J L Oliver, P Carpena, M Hackenberg, P Bernaola-Galván  (2004)  
!     IsoFinder : computational prediction of isochores in genome sequences
!     Nucleic Acids Research  32:  W287-W292   

!     IsoFinder project: http://isofinder.sourceforge.net/
!     IsoFinder web server: http://bioinfo2.ugr.es/IsoF/isofinder.html

!     Contacts:

!     Pedro Bernaola(1), Pedro Carpena(1) & Jose L. Oliver(2)
!     (1) Dpto. de Física Aplicada II, Universidad de Málaga, Spain
!         rick@ctima.uma.es, pcarpena@ctima.uma.es
!     (2) Dpto. de Genética, Universidad de Granada, Spain
!         oliver@ugr.es

!***************************************************************************

program isofinder
use leeadn
use findenes
use random
use cortes
use special
use diver

implicit none

!**** Secuencia de ADN y lista de cortes *******
character*90  dna_file,fout
integer*1     adn (350000000)  !Dna sequence
integer*4     nwords,enes,haygaps
integer*4     f_words(0:1)
type(corte)  ,pointer::first,last,p_max,p
type(corte)  ,target::pp
character*80  ss
character*10  method,method_des

INTEGER*4     MINSIZE  !minimo tamaño de los segmentos
INTEGER*4     WSIZE    !tamanio de ventana
REAL*8        SIGLEVEL !nivel de signifcación (tanto por 1)
INTEGER*4     NCUTS,NSEGS    !numero de cortes, numero de segmentos
INTEGER*4     IDUM

!******* PROVISIONALES ********
real*8        tmax
integer*4     p_corte,f_corte(0:1)
real*8        mu1,mu2,nv1,nv2,t_filt,sd,n1,n2,p0
integer*4     newcuts,ini,fin,size
real*8        p01,p02,p03,p04,p00,completo


!******* LISTA DE DATOS FINALES ********
integer*4   ,dimension(50000)::l_ini,l_fin,l_size
real*8      ,dimension(50000)::l_gc
integer*4    cont


!**********************
integer       lwin
integer       i,j,l,m,jj,iii,jjj,seq_limpia,N_iso

!********** PARAMETROS *********
call getarg(1,dna_file)

if (trim(dna_file)=='') then
  write(*,*)
  write(*,*) 'Syntax:'
  write(*,*) '  isofinder <DNA file> <Sig.Level> <Sig.Method> <Win.Size> <Output file>'
  write(*,*)
  stop
endif

call getarg(2,ss)
read(ss,*) SIGLEVEL
call getarg(3,method)
  if (method=='r') then
     method_des='   Random.'
  else if (method=='p1') then
     method_des=' t-student'
  else if (method=='p2') then
     method_des=' t dif.var'
  else if (method=='p3') then
     method_des='   Maximum'
  endif

call getarg(4,ss)
read(ss,*) WSIZE
call getarg(5,fout)             
open(2,file=fout)


IDUM=-34241117
MINSIZE=2*WSIZE+2

fin=0
call lee_adn(adn,nwords,enes,dna_file,1,fin,'SW',f_words)

write(2,*)                'DNA file:           ',trim(dna_file)
write(2,'(a17,i10)')      ' Length:          ',nwords
write(2,'(a17,f10.3)')    ' G+C (%):         ',100.d0*f_words(1)/(sum(f_words))
write(2,'(a17,i10)')      ' Undefined        ',enes
write(2,'(a17,f10.5)')    ' Sig.Level:       ',SIGLEVEL
write(2,'(a17,a18,a10)')  ' Method:          ',method,method_des
write(2,'(a17,i10)')      ' Window size:     ',wsize

call init_list(first,last,nwords,f_words)
call insert_enes(adn,first,nwords,haygaps,IDUM)

if (haygaps==1) then
pp=first
do while (associated(pp%next))
  if (pp%intact==1) then
  endif
  pp=pp%next
enddo
endif

!**********************************************
!*************** SEGMENTACION *****************
!**********************************************

completo=0.d0
NCUTS=0
print*
do
  p=>first
  newcuts=0
  do while (associated(p%next))
    if (p%intact==0.and.(p%next%posic-p%posic)>2*MINSIZE) then
      call tstudentmax(adn,p,tmax,p_corte,f_corte)

      call filterout(adn,p%posic+1,p_corte,WSIZE,mu1,nv1,n1)
      call filterout(adn,p_corte+1,p%next%posic,WSIZE,mu2,nv2,n2)
      if ( (n1 <= 0) .or. (n2 <= 0) ) goto 500
	  sd=dsqrt((nv1+nv2)*(1/n1+1/n2)/(n1+n2-2))	 
      t_filt=abs(mu1-mu2)/sd

      if (method=='r') then
         p0=random_ttest(adn,p,p_corte,WSIZE,t_filt,IDUM)          !Randomization
      else if (method=='p1') then
         p0=1.d0-ttest0(t_filt,1.d0*(n1+n2-2))                     !standard t-test
      else if (method=='p2') then
         p0=1.d0-ttest_difvar(mu1,mu2,nv1/(n1-1),nv2/(n2-1),n1,n2) !standard t-test diff. variances
      else if (method=='p3') then
         p0=prob(t_filt,n1+n2)                                     !maximum test
      endif
      if (p0>=SIGLEVEL) then

        call insert(p,p_corte,f_corte)
        p=>p%next
        newcuts=newcuts+1
        NCUTS=NCUTS+1
      else
        p%intact=1
        completo=completo+100.d0*(p%next%posic-p%posic)/(1.d0*nwords)
      endif
    endif
    500 p=>p%next
  enddo
  if (newcuts==0) exit
enddo
NSEGS=NCUTS+1
!**********************************************

!*********************************
!****** Crea  la lista   *********
!*********************************

p=>first
cont=0
do while (associated(p%next))
   cont=cont+1
   l_ini(cont)=p%posic+1
   l_fin(cont)=p%next%posic
   l_size(cont)=l_fin(cont)-l_ini(cont)+1
   l_gc(cont)=100.d0*(p%next%frec(1)-p%frec(1))/l_size(cont)
   p=>p%next
enddo

nsegs=cont   !hay que rectificar el número de segmentos
             !para tener en cuenta las enes.



!***************** INICIOS, FINALES, TAMANIOS                    ******
!***************** Y G+C DE CADA SEGMENTO EN LOS ARRAYS          ******
!***************** l_ini,l_fin,l_size,l_gc                       ******
!***************** que van desde 1 a NSEGS                       ******
!***************** Por ejemplo el listado:

!Imprime isocoras
write(2,*) '___________________________________________'
write(2,*) '       From          To        Size     G+C'
write(2,*) '___________________________________________'
lwin=WSIZE
do cont=1,NSEGS
   iii = l_ini(cont)
   jjj = l_fin(cont)
   if (l_gc(cont) > 0) then
       write(2,'(3i12,f8.2)') l_ini(cont),l_fin(cont),l_size(cont),l_gc(cont)
   endif
enddo

!*******************************************
!*******************************************
!*******************************************

CONTAINS

!*********** PROBABILITY STUFF *********************
!***************************************************

function neff(n)
  real*8  n,neff
  neff=-11.54+4.189*dlog(n)
end function neff


!***************************************************

FUNCTION PROB(x,n)

  real*8   prob,x,n
  real*8   f,expon,nu

  expon=neff(n)

  if (abs(x)<1.0e-5) then
    prob=0
  else
    expon=neff(n)
    f=0.8
    nu=n-2
    prob=(1.0-betai(0.5*f*nu,0.5d0*f,nu/(nu+x**2)))**expon
  endif

end FUNCTION prob


!***************************************************


!*************************************************

SUBROUTINE filterout(adn,ini,fin,wsize,mu,nsigma2,rnwin)
   integer*1   adn(:)
   integer*4   ini,fin,wsize
   real*8      mu,nsigma2,rnwin

   integer*4   i,j,nwin,size
   real*8      sx,sx2,x

   size=fin-ini+1
   nwin=size/wsize
   sx=0.d0
   sx2=0.d0
   do i=1,nwin
     x=0.d0
     do j=(i-1)*wsize,i*wsize
       if (adn(ini+j)==1) x=x+1
     enddo
     sx=sx+x
     sx2=sx2+x**2
   enddo
   sx=sx/wsize
   sx2=sx2/(wsize**2)

   mu=sx/(1.d0*nwin)
   nsigma2=(sx2-nwin*(mu**2))
   rnwin=nwin
end subroutine filterout

!*******************************************************


!*************************************************

function random_ttest(adn,p,p_corte,wsize,t0,idum)
   integer*1   adn(:)
   integer*4   ini,fin,wsize,nexper,p_corte
   real*8      mu1,mu2,nv1,nv2,gc(1000000)
   integer*4   n1,n2,n,idum
   type(corte) ,target::p

   integer*4   i,j,size,jj
   real*8      sd,t,t0,random_ttest,prob,x,df


   nexper=50000

   ini=p%posic+1
   fin=p_corte
   size=fin-ini+1
   n1=size/wsize
   do i=1,n1
     x=0.d0
     do j=(i-1)*wsize,i*wsize
       if (adn(ini+j)==1) x=x+1
     enddo
     gc(i)=x/(wsize)
   enddo

   ini=p_corte+1
   fin=p%next%posic
   size=fin-ini+1
   n2=size/wsize
   n=n1+n2
   do i=1,n2
     x=0.d0
     do j=(i-1)*wsize,i*wsize
       if (adn(ini+j)==1) x=x+1
     enddo
     gc(i+n1)=x/(wsize)
   enddo

   df=n-2.d0
   prob=0.d0
   do i=1,nexper
      do j=n,1,-1
        jj=(ran3(idum)*j)+1
        x=gc(jj)
        gc(jj)=gc(j)
        gc(j)=x
      enddo
      mu1=0.d0
      mu2=0.d0
      nv1=0.d0
      nv2=0.d0
      do j=1,n1
        mu1=mu1+gc(j)
        nv1=nv1+gc(j)**2
      enddo
      do j=n1+1,n
        mu2=mu2+gc(j)
        nv2=nv2+gc(j)**2
      enddo
      mu1=mu1/n1
      mu2=mu2/n2
      nv1=nv1-n1*mu1**2
      nv2=nv2-n2*mu2**2
      sd=dsqrt((nv1+nv2)*(1.d0/n1+1.d0/n2)/df)
      t=abs(mu1-mu2)/sd
      !print*, t,t0
      if (t<=t0) then
         prob=prob+1.d0
      endif
   enddo
   random_ttest=prob/(1.d0*nexper)

end function random_ttest

!******************************************************************

SUBROUTINE tstudentmax(adn,p,tmax,p_corte,f_corte)
  real*8       tmax
  integer*1    adn(:)
  integer*4    p_corte,f_corte(0:1)
  type(corte) ,target::p

  integer*1    b
  integer*4    ini,fin,i
  integer*4   ,dimension(0:1)::f1,f2
  real*8       n1,n2,n,sd,t,mu1,mu2


  ini=p%posic+1
  fin=p%next%posic
  f2=p%next%frec-p%frec
  f1=0
  n=fin-ini-1.d0            !*** resto 2 para la t
  tmax=0.d0
  do i=ini,fin-1
    b=adn(i)
    f1(b)=f1(b)+1
    f2(b)=f2(b)-1
    n1=i-ini+1
    n2=fin-i
    mu1=f1(0)/n1
    mu2=f2(0)/n2
	if ( (mu1*(1-mu1)+mu2*(1-mu2))*(1/n1+1/n2)/(n1+n2-2) <=0 ) goto 600
        sd=dsqrt((mu1*(1-mu1)+mu2*(1-mu2))*(1/n1+1/n2)/(n1+n2-2))
    t=abs(mu1-mu2)/sd
    if (t>tmax.and.(i-ini)>MINSIZE.and.(fin-i-1)>MINSIZE) then
       tmax=t
       p_corte=i
       f_corte=p%frec+f1
    endif
	600 continue
    enddo
END SUBROUTINE tstudentmax

FUNCTION total_divergence(first,f_words)
  !******** PARAMETERS **********
  type(corte),pointer:: first
  integer*4             f_words(:)
  real*8                total_divergence
  !****** LOCAL VARIABLES *******
  type(corte), pointer :: p
  integer*4               size,totsize
  real*8                  totentro,inter

  totentro=entro(f_words)
  totsize=sum(f_words)


  p=>first
  inter=0
  do while (associated(p%next))
    cont=cont+1
    size=sum(p%next%frec-p%frec)

    if (size/=0) then
      inter=inter-size*entro(p%next%frec-p%frec)
    endif
    p=>p%next
  end do


  total_divergence=inter/(1.0*totsize)+totentro
END FUNCTION total_divergence
!***************************************************

end
