!Simulation of the effect of sampling full sibs on tests of HWE Module MyDataType IMPLICIT NONE INTEGER, PARAMETER :: I1B = SELECTED_INT_KIND(2) ! -128 ~ 127 INTEGER, PARAMETER :: SP = KIND(1.0) ! Single Precision INTEGER, PARAMETER :: DP = KIND(1.0D0) ! Double Precision TYPE Element INTEGER(I1B), Pointer :: I END TYPE Element END Module MyDataType ! Program HWEtest USE MyDataType IMPLICIT NONE INTEGER I, J, K, L, M, KK, LL, MM, NN, NPerm, NRepeat INTEGER :: N(0 : 5), Idx(4, 4), GObs(10), SigPm(2, 5) INTEGER(I1B) :: GP(2, 2) INTEGER(I1B), TARGET, Allocatable :: G(:, :) INTEGER, Allocatable :: GT(:) TYPE(Element), Allocatable :: GG(:) REAL(DP), Allocatable :: Factorial(:) REAL(DP) :: CP(4), Q(4), CH, Fac, X, Y REAL :: R EXTERNAL :: DrawAllele !-parameters to be set WRITE(*, '(A)') 'Input number of permutations (Default=10000)?' READ(*, *) NPerm IF(NPerm < 1) THEN WRITE(*, '(A)') 'Number of permutations must be a positive integer!' STOP END IF WRITE(*, '(A)') 'Input number of individuals in the unrelated sample (Default=80)?' READ(*, *) N(0) WRITE(*, '(A)') 'Input 5 alternative fullsib sample sizes (Default=5, 10, 20, 40, 80)?' READ(*, *) N(1 : 5) IF(ANY(N(0:) < 1)) THEN WRITE(*, '(A)') 'The sample sizes must be > 0!' STOP END IF WRITE(*, '(A)') 'Input the frequencies of 4 alleles at a locus (Default=0.1, 0.2, 0.3, 0.4)?' READ(*, *) CP(1 : 4) IF(ANY(CP(:) <= 0.0_DP)) THEN WRITE(*, '(A)') 'Allele frequencies must be > 0!' STOP END IF X = Sum(CP(:)) IF(ABS(X - 1._DP) > 1.E-2_DP) THEN WRITE(*, '(A)') 'Allele frequencies do not sum to 1!' STOP END IF WRITE(*, '(A)') 'Input number of replicates (Default=1000)?' READ(*, *) NRepeat IF(NRepeat < 1) THEN WRITE(*, '(A)') 'Number of replicates must be > 0!' STOP END IF !- K = Sum(N(:)) Allocate(G(2, K), GT(2), GG(1)) ALLOCATE(Factorial(0 : 2 * K)) Factorial(0) = 0.D0 DO L = 1, 2 * K Factorial(L) = Log(DBLE(L)) + Factorial(L - 1) END DO K = 0 DO I = 1, 4 DO J = I, 4 K = K + 1 Idx(I, J) = K Idx(J, I) = K END DO END DO DO I = 2, 3 CP(I) = CP(I - 1) + CP(I) END DO CP(4) = 1.1 DO M = 1, 5 SigPm(:, M) = 0 DO I = 1, NRepeat IF(Mod(I, 10) == 0) WRITE(*, *)'M,I=', M, I DO J = 1, N(0) CALL DrawAllele(CP, G(:, J)) END DO DO J = 1, 2 CALL DrawAllele(CP, GP(:, J)) END DO DO J = 1 + N(0), N(M) + N(0) DO K = 1, 2 CALL RANDOM_NUMBER(R) G(K, J) = GP(Floor(R * 2 + 1.0), K) END DO END DO DO K = 1, 2 !1/2 for unrelated and combined samples GObs(:) = 0 NN = N(M) * (K - 1) + N(0) Q(:) = 0.0 DO J = 1, NN L = Idx(G(1, J), G(2, J)) GObs(L) = GObs(L) + 1 DO L = 1, 2 Q(G(L, J)) = Q(G(L, J)) + 1 END DO END DO Q(:) = Q(:) / 2.D0 / NN !permutation test Fac = Factorial(NN) - Factorial(2 * NN) DO J = 1, 4 Fac = Fac + Factorial(NINT(Q(J) * 2.D0 * NN)) END DO CH = Fac L = 0 MM = 0 DO J = 1, 4 DO LL = J, 4 L = L + 1 IF(GObs(L) < 1) CYCLE IF(J /= LL) MM = MM + GObs(L) CH = CH - Factorial(GObs(L)) END DO END DO CH = CH + MM * Log(2.D0) Y = 0.D0 IF(Ubound(GG, 1) < 2 * NN) THEN Deallocate(GG) Allocate(GG(2 * NN)) END IF DO KK = 1, NPerm FORALL(J = 1 : NN) GG(J)%I => G(1, J) FORALL(J = NN + 1 : NN + NN) GG(J)%I => G(2, J - NN) GObs(:) = 0 MM = 2 * NN DO J = 1, NN DO L = 1, 2 CALL RANDOM_NUMBER(R) LL = Floor(R * MM + 1.0) GT(L) = GG(LL)%I GG(LL)%I => GG(MM)%I MM = MM - 1 END DO L = Idx(GT(1), GT(2)) GObs(L) = GObs(L) + 1 END DO X = Fac L = 0 MM = 0 DO J = 1, 4 DO LL = J, 4 L = L + 1 IF(GObs(L) < 1) CYCLE IF(J /= LL) MM = MM + GObs(L) X = X - Factorial(GObs(L)) END DO END DO X = X + MM * Log(2.D0) IF(X <= CH) Y = Y + 1.D0 END DO IF(Y / NPerm < 0.05D0) SigPm(K, M) = SigPm(K, M) + 1 END DO END DO END DO WRITE(*, '(//3A8)')'#FS', 'Bold', 'Naive' DO I = 1, 5 WRITE(*, '(I8, 2 F8.4)') N(I), SigPm(:, I) / Dble(NRepeat) END DO END Program HWEtest ! Subroutine DrawAllele(CP, G) USE MyDataType IMPLICIT NONE REAL(DP), INTENT(IN) :: CP(4) INTEGER(I1B) :: G(2) INTEGER(I1B) :: I, L REAL :: R DO I = 1_I1B, 2_I1B CALL RANDOM_NUMBER(R) DO L = 1_I1B, 4_I1B IF(R >= CP(L)) CYCLE G(I) = L EXIT END DO END DO END Subroutine DrawAllele