R/thetaWC.pair.R

In FinePop: Fine-Scale Population Analysis

thetaWC.pair <-
function(popdata){
numpop <- popdata$npops
numlocus <- popdata$nloci
rpop <- 2

cfstmat <- matrix(0, nrow=numpop, ncol=numpop)
dimnames(cfstmat) <- list(popdata$pop_names,popdata$pop_names)

LFx <- intToUtf8(0x0A)
BSx <- intToUtf8(0x08)

message("Calculating population ", appendLF=F)
cstep.pop <- ""
for(cpop1 in 1:(numpop-1)){
for(cpop2 in (cpop1+1):numpop){
  message(paste0(rep(BSx, nchar(cstep.pop)), collapse=""), appendLF=F)
  cstep.pop <- paste0(cpop1, ":", cpop2, " ")
  message(cstep.pop, appendLF=F); flush.console()

  sum_a <- 0
  sum_abc <- 0
  for(cloc in 1:numlocus){
    nsamples12 <- popdata$indtyp[[cloc]][c(cpop1, cpop2)]

    genotype1A <- popdata$pop_allele[[cpop1]][[1]][,cloc]
    genotype1a <- popdata$pop_allele[[cpop1]][[2]][,cloc]
    genotype2A <- popdata$pop_allele[[cpop2]][[1]][,cloc]
    genotype2a <- popdata$pop_allele[[cpop2]][[2]][,cloc]
    allele_names <- rownames(popdata$allele_freq[[cloc]])

    for(callele in 1:((popdata$nalleles)[cloc])){
      n_bar <- mean(nsamples12)
      n_c <- (n_bar * rpop - sum(nsamples12^2)/(n_bar*rpop)) / (rpop-1)

      freqA <- popdata$allele_freq[[cloc]][callele,c(cpop1,cpop2)]
      p_bar <- sum(freqA * nsamples12) / (rpop*n_bar)

      s2 <- sum(nsamples12 * (freqA - p_bar)^2) / ((rpop-1)*n_bar)
      #s2 <- sum(nsamples12 * (freqA - p_bar)^2) / ((rpop-1)*n_bar) * (rpop-1)/rpop

      callelename <- allele_names[callele]
      hetero1 <- ((genotype1A==callelename)&(genotype1a!=callelename)) |
               ((genotype1A!=callelename)&(genotype1a==callelename))
      hetero1 <- hetero1[!is.na(hetero1)]
      hetero2 <- ((genotype2A==callelename)&(genotype2a!=callelename)) |
               ((genotype2A!=callelename)&(genotype2a==callelename))
      hetero2 <- hetero2[!is.na(hetero2)]
      freqAH <- c(sum(hetero1),sum(hetero2))/nsamples12
      #freqAH <- 2 * freqA * (1-freqA)
      #freqAH <- 1 - sum(freqA^2)
      h_bar <- sum(freqAH * nsamples12) / (rpop*n_bar)
             #0

      #WCa <- n_bar/n_c * (s2 - 1/(n_bar-1) * (p_bar*(1-p_bar) - (rpop-1)/rpop * s2 - h_bar/4 ))
      WCa <- (rpop-1)/rpop * n_bar/n_c * (s2 - 1/(n_bar-1) * (p_bar*(1-p_bar) - (rpop-1)/rpop * s2 - h_bar/4 ))
      #WCa <- n_bar/n_c * (s2 - 1/(n_bar-1) * (p_bar*(1-p_bar) - s2 - h_bar/4 ))
      #WCa <- n_bar/n_c * (s2*(rpop-1)/rpop - 1/(n_bar-1) * (p_bar*(1-p_bar) - (rpop-1)/rpop * s2 - h_bar/4 ))

      WCb <- n_bar/(n_bar-1) * (p_bar*(1-p_bar) - (rpop-1)/rpop * s2 - (2*n_bar-1)/(4*n_bar) * h_bar)
      #WCb <- n_bar/(n_bar-1) * (p_bar*(1-p_bar) - s2 - (2*n_bar-1)/(4*n_bar) * h_bar)
      #WCb <- p_bar*(1-p_bar)# - (rpop-1)/rpop * s2 - (2*n_bar-1)/(4*n_bar) * h_bar)

      WCc <- h_bar/2
      #WCc <- 0

      if(is.finite(WCa)){
         sum_a <- sum_a + WCa;
         sum_abc <- sum_abc + WCa + WCb + WCc
      }
    }#allele
  }#loc

  cfstmat[cpop1, cpop2] <- cfstmat[cpop2, cpop1] <-(sum_a / sum_abc)
}}#pop1 pop2
message(" done.")

return(cfstmat)
}