!###############################################################################
!                                                                              #
!         CCCCCC  SSSSSS  IIIII  RRRRRR   OOOO          L      IIIII           #
!         C       S         I    R    R  O    O         L        I             #
!         C       SSSSSS    I    RRRRRR  O    O  =====  L        I             #
!         C            S    I    R RR    O    O         L        I             #
!         CCCCCC  SSSSSS  IIIII  R   RR   OOOO          LLLLL  IIIII           #
!                                                                              #
!         Queensland Bioscience Precinct                                       #
!         306 Carmody Road                                                     #
!         St Lucia, QLD 4067                                                   #
!         AUSTRALIA                                                            #
!         Ph. (07) 3214 2392                                                   #
!         Fx. (07) 3214 2900                                                   #
!                                                                              #
!         Objective:      Computes Sensitivity for Microarray data             #
!         STD Input:      A file with 3 Fields: 1. Clone_ID                    #
!                                               2. Avg_Intensity               #
!                                               3. DE_Flag (1/0)               #
!                           ...and sorted in ascending order by Avg_Intensity  #
!                                                                              #
!            Generate it like: join -a1 de.all de.mn | \                       #
!                                 awk 'NF==3 {$3="1"}; \                       #
!                                      NF==2 {$3="0"}; \                       #
!                                      {print $0}' | sort -n +1 > input        #
!                                                                              #
!         STD Output:     A file with 4 fields: 1. Category number             #
!                                               2. Distrib of Total (=MPSS)    #
!                                               3. Distrib of DE               #
!                                               4. % DE                        #
!                                                                              #
!         Author:         Toni Reverter (Tony.Reverter-Gomez@csiro.au)         #
!         Last Revision:  July 2003                                            #
!                                                                              #
!###############################################################################

PROGRAM sensitivity

IMPLICIT NONE

INTEGER :: i,j,k,m
INTEGER :: ntot,nde,deflag

CHARACTER(LEN=3) :: cloneid ! NB. We'll read only 3 char as the name of the
                            ! clone in itself is of no interest

REAL :: cat(9),     & ! % records according to MPSS (and NOISE PNAS) paper
        cat_pde(9), & ! % DE genes in each category (as a % of total genes)
        log_tpm(9), & ! tpm in log10 scale for each category
        val           ! Average intensity real temporarily at std input

INTEGER :: cat_nde(9)    ! N DE genes in each category

TYPE rec_element
     REAL :: intens
     INTEGER :: deflag
END TYPE rec_element
type(rec_element), ALLOCATABLE, DIMENSION(:) ::  gene ! Data records


OPEN(10,file='sensi.tmp',form='formatted',status='unknown')

log_tpm(1)=0.0
log_tpm(2)=0.7
log_tpm(3)=1.0
log_tpm(4)=1.7
log_tpm(5)=2.0
log_tpm(6)=2.7
log_tpm(7)=3.0
log_tpm(8)=3.7
log_tpm(9)=4.0

DO i = 1, 9
   cat(i) = 100.0* exp((-2.0*log_tpm(i)**2)/(log_tpm(i)+1.0))
ENDDO

ntot = 0
nde = 0
DO
   READ(*,*,IOSTAT=k)cloneid,val,deflag
   IF( k .NE. 0 ) EXIT
   WRITE(10,*)val,deflag
   ntot = ntot + 1
   IF( deflag > 0 )nde = nde + 1
ENDDO
WRITE(*,*)'# Total Genes = ',ntot,' Total DE Genes = ',nde

REWIND(10)
ALLOCATE(gene(ntot))

DO i = 1, ntot
   READ(10,*)gene(i)%intens,gene(i)%deflag
ENDDO
!CLOSE(10,status='delete')

cat_nde = 0
cat_pde = 0.0

i=1
cat_nde(i) = nde
cat_pde(i) = 100.0 * nde/ntot
WRITE(*,'(a4,i1,f9.2,f9.3,f9.3,f9.3,f9.1)')'Cat_',i,log_tpm(i),cat(i), &
     100.0*cat_nde(i)/nde,cat_pde(i),gene(i)%intens

DO i = 2, 9
   j = ntot - int(ntot*cat(i)/100.00)
   m = 0 ! Counter for DE genes found so far
   DO k = 1, ntot
      IF( gene(k)%deflag > 0 )THEN
         m = m + 1
         IF( gene(k)%intens > int(gene(j)%intens) )THEN
            cat_nde(i) = nde-m+1
            cat_pde(i) = 100.0*(cat_nde(i)/(ntot*(cat(i)/100.0)))
            EXIT
         ENDIF
      ENDIF
   ENDDO
   WRITE(*,'(a4,i1,f9.2,f9.3,f9.3,f9.3,f9.1)') &
        'Cat_',i,log_tpm(i),cat(i),100.0*cat_nde(i)/nde,cat_pde(i),gene(j)%intens
ENDDO



END PROGRAM sensitivity