C
C FUNCTION TO CALCULATE THE EXACT TEST OF HARDY-WEINBERG EQUILIBRIUM AS 
C DESCRIBED IN WIGGINTON JE, CUTLER DJ, AND GR ABECASIS (2005) A NOTE 
C ON EXACT TESTS OF HARDY-WEINBERG EQUILIBRIUM. AJHG: 76:000-000 
C 
C DATA DICTIONARY 
C 
C MAXHET   SETS THE MAXIMUM POSSIBLE NUMBER OF HETEROZYGOTES PER SAMPLE. 
C          IF YOUR DATA SET WILL INCLUDE MORE THAN THE DEFAULT NUMBER OF
C          MINOR ALLELE COPIES AT ANY SINGLE MARKER LOCUS, MAXHET SHOULD BE 
C          INCREMENTED ACCORDINGLY 
C N_ALLR   NUMBER OF MINOR ALLELE COPIES 
C N_ALL    TOTAL NUMBER OF ALLELES 
C N_GENO   TOTAL NUMBER OF GENOTYPES 
C SNPHWE   RETURN VALUE FOR SNPHWE FUNCTION 
C PVALUE   P_VALUE FOR STATISTIC 
C SUM      USED TO ACCUMULATE SUM OF HET PROBABILITIES FOR SCALING IN FINAL STEP 
C PROBS    VECTOR OF HETREOZYGOTE PROBABILITIES 
C MIDPT    MIDPOINT OF THE HETEROZYGOTE PROBABILITY DISTRIBUTION 
C C_HET    NUMBER OF HETEROZYGOTES IN CURRENT STEP 
C C_HOMR   NUMBER OF RARE HOMOZYGOTES IN CURRENT STEP 
C C_HOMC   NUMBER OF COMMON HOMOZYGOTES IN CURRENT STEP 
C 
C WRITTEN BY: JAN WIGGINTON
C
      FUNCTION SNPHWE(N_HET, N_HOM1, N_HOM2, MAXHET) 
C 
C 
      DATA       KIN, KOUT /5,6/ 
      INTEGER    N_ALLR, N_ALL
      INTEGER    N_GENO 
      INTEGER    MIDPT
      INTEGER    C_HET, C_HOMR, C_HOMC 
      REAL*8     SNPHWE, PVALUE, SUM
      REAL*8     PROBS(MAXHET)
      REAL*8     PREV
C
C CALCULATE TOTAL NUMBER OF GENOTYPES TOTAL NUMBER OF MINOR ALLLELES
C
      IF (N_HET .LT. 0 .OR. N_HOM1 .LT. 0 .OR. N_HOM2 .LT. 0) THEN
         WRITE(KOUT, 101)'SNP-HWE: ERROR - GENOTYPE 
     1   CONFIGURATION (HETS = ', N_HETS, ' HOM1 = ', 
     2   N_HOM1, ' HOM2 = ', N_HOM2, ' ) IS 
     3   NOT POSSIBLE'
 101     FORMAT(A48, I5, A8, I5, A8, I5, A18)
      ENDIF

      N_GENO = N_HET + N_HOM1 + N_HOM2
      N_ALL  = N_GENO * 2
C
      N_ALLR = N_HET + (N_HOM1 * 2) 
      IF (N_HOM2 .LT. N_HOM1) THEN
         N_ALLR = N_HET + N_HOM2 * 2
      ENDIF
C
C REINITIALIZE DISTRIUBUTION
C
      DO 10 I=1, MAXHET
         PROBS(I) = 0.0
 10   CONTINUE 
C 
C FIND THE MIDPOINT OF THE DISTRIBUTION
C
      MIDPT = IDINT(N_ALLR * (N_ALL - N_ALLR)/ DBLE(N_ALL)) 
      IF (MOD(N_ALLR, 2) .NE. MOD(MIDPT, 2))  THEN
         MIDPT = MIDPT + 1
      ENDIF 
C
C SET PROBABILITY AT MIDPOINT TO 1.0 FOR SIMPLICITY -> SCALING DONE LATER 
C MAKES SPECIFIC VALUE IRRELEVANT
C
      PROBS(MIDPT + 1) = 1.0
      SUM = 1.0
C 
C CALCULATE PROBABILITIES FOR HETS > # HETS AT MIDPOINT 
C
      C_HET  = MIDPT
      C_HOMR = (N_ALLR - MIDPT) / 2 
      C_HOMC = N_GENO - C_HET - C_HOMR

 20   IF (C_HET .GE. 2) THEN
         PROBS(C_HET - 1) = PROBS(C_HET + 1) * DBLE(C_HET) 
     1   * (DBLE(C_HET) - 1.0) / (4.0 * (DBLE(C_HOMR) + 1.0) 
     2   * (DBLE(C_HOMC) + 1.0 ))
         SUM    = SUM + PROBS(C_HET - 1)
C
         C_HET  = C_HET - 2
         C_HOMR = C_HOMR + 1
         C_HOMC = C_HOMC + 1
         GOTO 20
      ENDIF 
C
C CALCULATE PROBABILITIES FOR HETS < # HETS AT MIDPOINT
C
      C_HET  = MIDPT
      C_HOMR = (N_ALLR - MIDPT) / 2
      C_HOMC = N_GENO - C_HET - C_HOMR 
C
 30   IF (C_HET .LE. N_ALLR - 2) THEN
         PROBS(C_HET + 3) = PROBS(C_HET + 1) * 4.0 
     1    * DBLE(C_HOMR) * DBLE(C_HOMC) / ( (DBLE(C_HET) + 2.0) 
     2    * (DBLE(C_HET) + 1.0) )
         SUM = SUM + PROBS(C_HET + 3)
C
         C_HET  = C_HET + 2
         C_HOMR = C_HOMR - 1
         C_HOMC = C_HOMC - 1
         GOTO 30
      ENDIF 
C 
C SCALING TO ENSURE PROBABILITIES SUM TO 1.0
C
      DO 40 I = 1, N_ALLR + 1
         PREV = PROBS(I)
         PROBS(I) = PROBS(I) / SUM 
 40   CONTINUE 
C
C CALCULATE P_HWE AS SUM OF HETEROZYGOTE PROBABILITIES <= P(N_HET)
C
      PVALUE = 0.0
      C_HET = 1
C 
 50   IF (C_HET .LE. N_ALLR + 1) THEN
         IF (PROBS(C_HET) .LE. PROBS(N_HET + 1)) THEN
            PVALUE = PVALUE + PROBS(C_HET)
         ENDIF
         C_HET = C_HET + 1
         GOTO 50
      ENDIF 
C
      SNPHWE = PVALUE
      IF (SNPHWE .GT. 1.0) THEN
         SNPHWE = 1.0
      ENDIF
C
      RETURN
      END


