GROUP3

## Haplotype counts
AB <- 40
Ab <- 10
aB <- 20
ab <- 30
# Total sample size
n <- AB + Ab + aB + ab
n

## [1] 100

# Allele frequencies
pA <- pAB + pAb
pa <- paB + pab
pB <- pAB + paB
pb <- pAb + pab
c(pA = pA, pa = pa, pB = pB, pb = pb)

##  pA  pa  pB  pb 
## 0.5 0.5 0.6 0.4

# LD coefficient D
D <- pAB - pA * pB
D

## [1] 0.1

# Maximum possible D for positive D
Dmax <- min(pA * pb, pa * pB)
# Standardized D and r-squared
Dprime <- D / Dmax
r2 <- D^2 / (pA * pa * pB * pb)
round(c(D = D, Dprime = Dprime, r2 = r2), 4)

##      D Dprime     r2 
## 0.1000 0.5000 0.1667

## A reusable function

ld_measures <- function(AB, Ab, aB, ab) { 
  n <- AB + Ab + aB + ab
pAB <- AB/n
pAb <- Ab/n
paB <- aB/n
pab <- ab/n
pA <- pAB + pAb
pa <- paB + pab
pB <- pAB + paB
pb <- pAb + pab 
D <- pAB - pA * pB
if(
D >= 0) {
Dmax <- min(pA* pb, pa* pB)
} else{
Dmax <- min(pA* pB, pa* pb)
}
Dprime <-D/ Dmax
r2 <-D^2/ (pA* pa* pB*pb)

return(c(D=D, Dprime = Dprime, r2 = r2))
}

### ASSIGNMENT
##TASK 1
##Haplotype counts
AB <- 25
Ab <- 15
aB <- 35
ab <- 25
#TOTAL
n <- AB + Ab + aB + ab
n

## [1] 100

##
ld_measures <- function(AB, Ab, aB, ab) {
  N <- AB + Ab + aB + ab
  
  pA <- (AB + Ab) / N
  pa <- (aB + ab) / N
  pB <- (AB + aB) / N
  pb <- (Ab + ab) / N
  
  P_AB <- AB / N
  D <- P_AB - (pA * pB)
  
  return(list(
    pA = pA,
    pa = pa,
    pB = pB,
    pb = pb,
    D = D
  ))
}

ld_measures(AB = 25, Ab = 15, aB = 35, ab = 25)

## $pA
## [1] 0.4
## 
## $pa
## [1] 0.6
## 
## $pB
## [1] 0.6
## 
## $pb
## [1] 0.4
## 
## $D
## [1] 0.01

##TASK 2
#Haplotype frequencies
pAB<- 0.30
pAb<- 0.20
paB<- 0.10
pab<- 0.40

#Allele frequencies
pA <- pAB + pAb
pa <- paB + pab
pB <- pAB + paB
pb <- pAb + pab
c(pA = pA, pa = pa, pB = pB, pb = pb)

##  pA  pa  pB  pb 
## 0.5 0.5 0.4 0.6

#Expected pAB under independence
pAB_exp <- pA * pB
pAB_exp

## [1] 0.2

#D
D <- pAB - pAB_exp
D

## [1] 0.1

#Interpretation
if (D > 0) {
  "LD is positive"
} else if (D < 0) {
  "LD is negative"
} else {
  "LD is zero"
}

## [1] "LD is positive"

##TASK 3
#Values
pA <- 0.7
pa <- 0.3
pB <- 0.6
pb <- 0.4
pAB <- 0.50

# Expected frequency
pAB_exp <- pA * pB
pAB_exp

## [1] 0.42

# D
D <- pAB - pAB_exp
D

## [1] 0.08

#D'
Dmax <- min(pA * pb, pa * pB)
Dprime <- D / Dmax
Dprime

## [1] 0.4444444

# r^2
r2 <- D^2 / (pA * pa * pB * pb)
r2 

## [1] 0.1269841