R/update_slim.R In BlockFeST: Bayesian Calculation of Region-Specific Fixation Index to Detect Local Adaptation

```#source("math.R")

#constants
#e_f 		       <-  0.05
#sd_prior_beta          <-  1 #3
#mean_prior_beta        <- -1  #-5.15
#lambda 	       <-  1

X_update_d_alphaco_i <- function(id){ # id

if(id[1]==FALSE){ # Pilot runs
id      <- 1:BBB2\$locnum
locnum2 <- BBB2\$locnum
}else{                # Main Loop
locnum2 <- length(id)
}
##

sub_freq_locus    <- BBB2\$freq_locus[id,]
sub_freq_pop      <- lapply(BBB2\$freq_pop,function(yy){return(yy[id,])})
old_alpha         <- BBB2\$d_alpha[id]
sub_sample_size   <- BBB2\$sample_size[id,]
sub_GROUP         <- BBB2\$GROUP[id]

##

old_theta     <- outer(old_alpha,BBB2\$d_beta,"+")
old_theta     <- exp(-(old_theta))
old1          <- lgamma(old_theta)-lgamma(sub_sample_size + old_theta)

val <- matrix(0,locnum2,dim(BBB2\$freq_locus)[2])
for(xx in 1:BBB2\$popnum){
xyz <- old_theta[,xx]*sub_freq_locus
val <- val + lgamma(sub_freq_pop[[xx]] + xyz) - lgamma(xyz)
}

old1 <- rowSums(old1)
oldL <- old1 + rowSums(val,na.rm=TRUE)

## modification ------------------------------------
oldL <- tapply(oldL,sub_GROUP,sum)
##--------------------------------------------------

# e             <- rnorm(locnum2,0,sqrt(BBB2\$var_prop_alpha))
# new_alpha     <- old_alpha + e

### modification -----------------------------------------------

if(!BBB2\$PILOT){

whatgroup       <- unique(sub_GROUP)
xxx             <- BASIX.match(whatgroup,BBB2\$GROUP)
e_group         <- rnorm(length(whatgroup),0,sqrt(BBB2\$var_prop_alpha[xxx]))

#
BBB2\$new_alpha   <- old_alpha
BBB2\$jjj <- 1
sapply(whatgroup,function(x){
what            <- x==sub_GROUP
BBB2\$new_alpha[what] <- old_alpha[what] + e_group[BBB2\$jjj]
BBB2\$jjj <- BBB2\$jjj + 1
})

new_alpha <- BBB2\$new_alpha
#

}else{

e_group    <- rnorm(BBB2\$N.REGIONS,0,sqrt(BBB2\$var_prop_alpha))
new_alpha  <- old_alpha + e_group[BBB2\$GROUP]
whatgroup  <- 1:BBB2\$N.REGIONS

}

#
##--------------------------------------------------------------

new_theta <- outer(new_alpha,BBB2\$d_beta,"+")
new_theta <- exp(-(new_theta))
new1      <- lgamma(new_theta)-lgamma(sub_sample_size+new_theta)

val <- matrix(0,locnum2,dim(BBB2\$freq_locus)[2])
for(xx in 1:BBB2\$popnum){
xyz <- new_theta[,xx]*sub_freq_locus
val <- val + lgamma(sub_freq_pop[[xx]] + xyz) - lgamma(xyz)
}

new1 <- rowSums(new1)
newL <- new1 + rowSums(val,na.rm=TRUE)

## modification ----------------------------
newL            <- tapply(newL,sub_GROUP,sum)

old_alpha_group <- tapply(old_alpha,sub_GROUP,function(x){return(x[1])})
new_alpha_group <- tapply(new_alpha,sub_GROUP,function(x){return(x[1])})

# print(old_alpha)
#  print(new_alpha)

##-----------------------------------------

#A        	       <-  -oldL + newL + (old_alpha^2-new_alpha^2)/(2*BBB2\$sd_prior_alpha^2)

## mod ------------------------
A         	       <-  -oldL + newL + (old_alpha_group^2-new_alpha_group^2)/(2*BBB2\$sd_prior_alpha^2)
##-----------------------------

# r        	      <-  runif(locnum2)

## mod ---------------------------
r                    <-  runif(length(whatgroup))
## mod ---------------------------

check    	       <-  log(r) < A

# BBB2\$d_alpha[id[check]]      <- new_alpha[check]
# BBB2\$acc_alpha[id[check]]    <- BBB2\$acc_alpha[id[check]] + 1

# modification -------------------------------------
# update groups

group_idX   <- which(check)
group_id    <- whatgroup[group_idX]

#  make this also in C

if(length(group_idX)>0){

#res <- .Call("update_group2",group_id,BBB2\$GROUP,BBB2\$d_alpha,new_alpha,BBB2\$acc_alpha)
#BBB2\$d_alpha   <- res[[1]]
#BBB2\$acc_alpha <- res[[2]]

mm         <- sapply(group_id,function(xx){xx==BBB2\$GROUP})
BBB2\$iii   <- 1
apply(mm,2,function(x){
BBB2\$d_alpha[x]   <- new_alpha_group[group_idX[BBB2\$iii]] # new_alpha[x]
BBB2\$acc_alpha[x] <- BBB2\$acc_alpha[x] + 1
BBB2\$iii <- BBB2\$iii + 1

})
}

# ---------------------------------------------------

}# End of Function

update_d_betaco <- function(){

# Init
#val <- rep(NA,2)
val <- BBB2\$GLOBAL_INIT2
# ---------------

old_beta    <- BBB2\$d_beta
old_theta   <- outer(BBB2\$d_alpha,old_beta,"+")
old_theta   <- exp(-(old_theta))

old1 <- lgamma(old_theta)-lgamma(BBB2\$sample_size+old_theta)

for(xx in 1:BBB2\$popnum){
xyz     <- old_theta[,xx]*BBB2\$freq_locus
val[xx] <- sum(lgamma(BBB2\$freq_pop[[xx]]+xyz)-lgamma(xyz),na.rm=TRUE)
}

old1 <- colSums(old1)
oldL <- old1 + val
#--------- calculate new likelihood
e          <- rnorm(BBB2\$popnum,0,BBB2\$var_prop_beta)
BBB2\$d_beta <- BBB2\$d_beta + e

new_theta <- outer(BBB2\$d_alpha,BBB2\$d_beta,"+")
new_theta <- exp(-(new_theta))
new1      <- lgamma(new_theta)-lgamma(BBB2\$sample_size+new_theta)
for(xx in 1:BBB2\$popnum){
xyz     <- new_theta[,xx]*BBB2\$freq_locus
val[xx] <- sum(lgamma(BBB2\$freq_pop[[xx]]+xyz)-lgamma(xyz),na.rm=TRUE)
}
new1 <- colSums(new1)
newL <- new1 + val

A                <-  -oldL + newL #+ (old_beta^2-BBB2\$d_beta^2)/(2*BBB2\$sd_prior_beta^2)
#+ (old_alpha_group^2-new_alpha_group^2)/(2*BBB2\$sd_prior_alpha^2)
#+ (old_beta-BBB2\$d_beta)*(old_beta+BBB2\$d_beta-2*BBB2\$mean_prior_beta)/(2*BBB2\$sd_prior_beta^2)
r                <-  runif(BBB2\$popnum)
check            <-  log(r)>A
BBB2\$d_beta[check]        <- old_beta[check]
BBB2\$acc_beta[!check]     <- BBB2\$acc_beta[!check] + 1

}# End of Function

update_freq_codominant <- function(){

lambda    <- 1
# choose randomly two haplotype of each locus
# samp      <-  sapply(BBB2\$hapcount,function(x){return(sample(1:x,size=2)) })

theta     <-  outer(BBB2\$d_alpha,BBB2\$d_beta,"+")
theta     <-  exp(-theta)
# ok

#BBB2\$iii       <- 1
# frequencies of haplotypes of the two random haps
#samp_freq <- apply(samp,2,function(x){
#	      r   <- BBB2\$freq_locus[BBB2\$iii,x]
#	      BBB2\$iii <- BBB2\$iii+1
#	      return(r)
#          })
#ok

## mod
samp_freq <- t(BBB2\$freq_locus)
## end mod

#print(dim(samp_freq))
#print(dim(BBB2\$freq_locus))

# haplotype distribution for each pop of the two random haps
# BBB2\$iii           <- 1
samp_pop_freq <- vector("list",BBB2\$popnum)
for(yy in 1:BBB2\$popnum){

samp_pop_freq[[yy]] <- t(BBB2\$freq_pop[[yy]])

#samp_pop  <- apply(samp,2,function(x){
#	      r       <- BBB2\$freq_pop[[yy]][BBB2\$iii,x]
#	      BBB2\$iii <- BBB2\$iii+1
#	      return(r)
#	   })
#samp_pop_freq[[yy]] <- samp_pop
#BBB2\$iii <- 1
}
#ok
#print(dim(samp_pop_freq[[1]]))

old_L <- BBB2\$GLOBAL_INIT3
# 1 <- m
# 2 <- n
for(xx in 1:BBB2\$popnum){

old_L <- old_L + lgamma(samp_pop_freq[[xx]][1,] + theta[,xx] * samp_freq[1,]) -
lfactorial (samp_pop_freq[[xx]][1,]) - lgamma(theta[,xx]*samp_freq[1,])+
lgamma(samp_pop_freq[[xx]][2,]+theta[,xx]*samp_freq[2,]) -
lfactorial(samp_pop_freq[[xx]][2,]) - lgamma(theta[,xx]*samp_freq[2,])

}

old_m <- samp_freq[1,]
old_n <- samp_freq[2,]

#### BwPrp
old_m_n 	 <- old_m + old_n
old_m_ancestral  <- old_m + BBB2\$e_ancestral
old_m_ancestral2 <- old_m - BBB2\$e_ancestral

checkmin <- cbind(old_m_n,old_m_ancestral)
checkmax <- cbind(BBB2\$GLOBAL_INIT3,old_m_ancestral2)
minn     <- apply(checkmin,1,min)
maxx     <- apply(checkmax,1,max)
BwPrp    <- minn - maxx
###################

# change frequencies
u     <- runif(BBB2\$locnum,0,BwPrp) + maxx
# samp_freq_new
samp_freq_new     <- samp_freq
samp_freq_new[1,] <- u
samp_freq_new[2,] <- old_n + old_m - u

### FwPrp
ss1   <- colSums(samp_freq_new)
ss2   <- samp_freq_new[1,] + BBB2\$e_ancestral
ss3   <- samp_freq_new[1,] - BBB2\$e_ancestral
minn  <- apply(cbind(ss1,ss2),1,min)
maxx  <- apply(cbind(BBB2\$GLOBAL_INIT3,ss3),1,max)
FwPrp <- minn - maxx
###################

new_L <- BBB2\$GLOBAL_INIT3
for(xx in 1:BBB2\$popnum){

new_L <- new_L + lgamma(samp_pop_freq[[xx]][1,] + theta[,xx] * samp_freq_new[1,]) -
lfactorial (samp_pop_freq[[xx]][1,]) - lgamma(theta[,xx]*samp_freq_new[1,])+
lgamma(samp_pop_freq[[xx]][2,]+theta[,xx]*samp_freq_new[2,]) -
lfactorial(samp_pop_freq[[xx]][2,]) - lgamma(theta[,xx]*samp_freq_new[2,])

}

A     <- -old_L + new_L + log(BwPrp) - log(FwPrp) +
(lambda-1)*(log(samp_freq_new[1,]) + log(samp_freq_new[2,]) - log(old_m) - log(old_n) )
r     <- runif(BBB2\$locnum)
check <- log(r) < A

#for(xx in which(check)){
# BBB2\$freq_locus[xx,samp[1,xx]]    <- samp_freq_new[1,xx]
# BBB2\$freq_locus[xx,samp[2,xx]]    <- samp_freq_new[2,xx]
# BBB2\$acc_freq_ancestral[xx]       <- BBB2\$acc_freq_ancestral[xx] + 1
#}

for(xx in which(check)){
BBB2\$freq_locus[xx,1]             <- samp_freq_new[1,xx]
BBB2\$freq_locus[xx,2]             <- samp_freq_new[2,xx]
BBB2\$acc_freq_ancestral[xx]       <- BBB2\$acc_freq_ancestral[xx] + 1
}

}# End of Function

X_jump_model_codominant <- function(){

old_alpha 	   <- BBB2\$d_alpha
excluded   	   <- which(!BBB2\$alpha_included)

BBB2\$d_alpha[BBB2\$alpha_included]  <-  0 #rnorm(1,0,1)

#cat(1)

if(length(excluded)!=0){

### mod -------------------------------------------
whatgroups    <- unique(BBB2\$GROUP[excluded])
group_id      <- match(whatgroups,BBB2\$GROUP)
e_group       <- rnorm(length(whatgroups),mean=BBB2\$mean_alpha[group_id],sd=sqrt(BBB2\$var_alpha[group_id]))
#e_group      <- rnorm(length(whatgroups),mean=BBB2\$d_alpha[group_id],1)
#e_group      <- rlaplace(length(whatgroups),BBB2\$mean_alpha[group_id],BBB2\$var_alpha[group_id])

.Call("update_group_Jump",
BBB2\$d_alpha,
as.integer(BBB2\$GROUP),
as.integer(whatgroups),
e_group
)

}

BBB2\$alpha_included     <- !BBB2\$alpha_included

## old_theta
old_theta          <- outer(old_alpha,BBB2\$d_beta,"+")
old_theta          <- exp(-(old_theta))
old1               <- lgamma(old_theta)-lgamma(BBB2\$sample_size+old_theta)

val  <- BBB2\$GLOBAL_INIT1

for(xx in 1:BBB2\$popnum){
xyz <- old_theta[,xx]*BBB2\$freq_locus
val <- val + lgamma(BBB2\$freq_pop[[xx]]+xyz) - lgamma(xyz)
}

old1 <- rowSums(old1)
oldL <- old1 + rowSums(val,na.rm=TRUE)
oldL <- tapply(oldL,BBB2\$GROUP,sum)

#+ tapply(oldL,BBB2\$GROUP,min)
#- log(1-exp(tapply(oldL,BBB2\$GROUP,max)))

## new_theta
new_theta          <- outer(BBB2\$d_alpha,BBB2\$d_beta,"+")
new_theta          <- exp(-(new_theta))
new1               <- lgamma(new_theta)-lgamma(BBB2\$sample_size+new_theta)

val  <- BBB2\$GLOBAL_INIT1
for(xx in 1:BBB2\$popnum){
xyz <- new_theta[,xx]*BBB2\$freq_locus
val <- val + lgamma(BBB2\$freq_pop[[xx]]+xyz) - lgamma(xyz)
}

new1 <- rowSums(new1)
newL <- new1 + rowSums(val,na.rm=TRUE)
newL <- tapply(newL,BBB2\$GROUP,sum)

#+ tapply(newL,BBB2\$GROUP,min)
#- log(1-exp(tapply(newL,BBB2\$GROUP,max)))

A_group <- -oldL + newL
A 	<- A_group[BBB2\$GROUP]

included       		<-  BBB2\$alpha_included #==TRUE
excluded       		<-  !included

A_included <- 0
A_excluded <- 0

# include in the selection model
if(any(included)){
A_included     		<- - log(BBB2\$tempting)
- log(dnorm(BBB2\$d_alpha[included],BBB2\$mean_alpha[included],sqrt(BBB2\$var_alpha[included])))
+ log(dnorm(BBB2\$d_alpha[included],0,1))
}

# exclude from the selection model
if(any(excluded)){
A_excluded     		<-   log(BBB2\$tempting)
+ log(dnorm(BBB2\$d_alpha[excluded],BBB2\$mean_alpha[excluded],sqrt(BBB2\$var_alpha[excluded])))
- log(dnorm(BBB2\$d_alpha[excluded],0,1))
}

A[included] 	           <-  A[included] + A_included
A[excluded]                <-  A[excluded] + A_excluded

A                       <-  tapply(A,BBB2\$GROUP,unique)

write.table(A,"Jump_A")

r                       <-  runif(BBB2\$N.REGIONS)

check            	   <-  log(r) > A

group_id       <- which(check)

# tempting
#  BBB2\$acc.ratio <- BBB2\$acc.ratio + length(which(!check))/length(A)
#

if(length(group_id)>0){

.Call("update_group_Jump2",
BBB2\$d_alpha,
old_alpha,
as.integer(BBB2\$GROUP),
as.integer(group_id),
BBB2\$alpha_included
)
}

}# End of Function

X_update_d_alphaco_i_new <- function(){

##

sub_freq_locus    <- BBB2\$freq_locus
sub_freq_pop      <- BBB2\$freq_pop
old_alpha         <- BBB2\$d_alpha
sub_sample_size   <- BBB2\$sample_size
sub_GROUP         <- BBB2\$GROUP
locnum2           <- length(BBB2\$d_alpha)

##

old_theta     <- outer(old_alpha,BBB2\$d_beta,"+")
old_theta     <- exp(-(old_theta))
old1          <- lgamma(old_theta)-lgamma(sub_sample_size + old_theta)

val <- matrix(0,locnum2,dim(BBB2\$freq_locus)[2])
for(xx in 1:BBB2\$popnum){
xyz <- old_theta[,xx]*sub_freq_locus
val <- val + lgamma(sub_freq_pop[[xx]] + xyz) - lgamma(xyz)
}

old1 <- rowSums(old1)
oldL <- old1 + rowSums(val,na.rm=TRUE)
oldL <- tapply(oldL,sub_GROUP,sum)

#+ tapply(oldL,sub_GROUP,min)
#- log(1-exp(tapply(oldL,sub_GROUP,max)))

e_group   <- rnorm(BBB2\$N.REGIONS,0,tapply(BBB2\$var_prop_alpha,BBB2\$GROUP,unique)) # BUG FIXED
#e_group   <- rnorm(BBB2\$N.REGIONS,0,1)

new_alpha  <- old_alpha + e_group[BBB2\$GROUP]
whatgroup  <- 1:BBB2\$N.REGIONS

new_theta <- outer(new_alpha,BBB2\$d_beta,"+")
new_theta <- exp(-(new_theta))
new1      <- lgamma(new_theta)-lgamma(sub_sample_size+new_theta)

val <- matrix(0,locnum2,dim(BBB2\$freq_locus)[2])
for(xx in 1:BBB2\$popnum){
xyz <- new_theta[,xx]*sub_freq_locus
val <- val + lgamma(sub_freq_pop[[xx]] + xyz) - lgamma(xyz)
}

new1 <- rowSums(new1)
newL <- new1 + rowSums(val,na.rm=TRUE)
newL <- tapply(newL,sub_GROUP,sum)

#+ tapply(newL,sub_GROUP,min)
#- log(1-exp(tapply(newL,sub_GROUP,max)))

old_alpha_group <- tapply(old_alpha,sub_GROUP,function(x){return(x[1])})
new_alpha_group <- tapply(new_alpha,sub_GROUP,function(x){return(x[1])})

A         	     <-  -oldL + newL + (old_alpha_group^2-new_alpha_group^2)/(2*BBB2\$sd_prior_alpha^2)
r                    <-  runif(length(whatgroup))

check    	     <-  log(r) < A

group_idX   <- which(check)
group_id    <- whatgroup[group_idX]

if(length(group_idX)>0){

#CCode
.Call("update_group2",
as.integer(group_id),
as.integer(BBB2\$GROUP),
BBB2\$d_alpha,
new_alpha,
BBB2\$acc_alpha
)

}

}# End of function
X_update_d_alphaco_i_new_id <- function(id){

sub_freq_locus    <- BBB2\$freq_locus[id,]
sub_freq_pop      <- lapply(BBB2\$freq_pop,function(yy){return(yy[id,])})
old_alpha         <- BBB2\$d_alpha[id]
new_alpha         <- old_alpha
sub_sample_size   <- BBB2\$sample_size[id,]
sub_GROUP         <- BBB2\$GROUP[id]
locnum2           <- length(id)

old_theta     <- outer(old_alpha,BBB2\$d_beta,"+")
old_theta     <- exp(-(old_theta))
old1          <- lgamma(old_theta)-lgamma(sub_sample_size + old_theta)

val <- matrix(0,locnum2,dim(BBB2\$freq_locus)[2])
for(xx in 1:BBB2\$popnum){
xyz <- old_theta[,xx]*sub_freq_locus
val <- val + lgamma(sub_freq_pop[[xx]] + xyz) - lgamma(xyz)
}

old1 <- rowSums(old1)
oldL <- old1 + rowSums(val,na.rm=TRUE)
oldL <- tapply(oldL,sub_GROUP,sum)

# create new alpha
group_id   <- unique(sub_GROUP) # what groups have to be updated
xxx        <- BASIX.match(group_id,BBB2\$GROUP)
e_group    <- rnorm(length(group_id),0,sqrt(BBB2\$var_prop_alpha[xxx]))

# update the whole groups
.Call("update_groupJumpX",
new_alpha,
as.integer(sub_GROUP),
as.integer(group_id),
e_group
)

#whatgroup  <- 1:BBB2\$N.REGIONS

new_theta <- outer(new_alpha,BBB2\$d_beta,"+")
new_theta <- exp(-(new_theta))
new1      <- lgamma(new_theta)-lgamma(sub_sample_size+new_theta)

val <- matrix(0,locnum2,dim(BBB2\$freq_locus)[2])
for(xx in 1:BBB2\$popnum){
xyz <- new_theta[,xx]*sub_freq_locus
val <- val + lgamma(sub_freq_pop[[xx]] + xyz) - lgamma(xyz)
}

new1 <- rowSums(new1)
newL <- new1 + rowSums(val,na.rm=TRUE)
newL <- tapply(newL,sub_GROUP,sum)

old_alpha_group <- tapply(old_alpha,sub_GROUP,function(x){return(x[1])})
new_alpha_group <- tapply(new_alpha,sub_GROUP,function(x){return(x[1])})

A         	     <-  -oldL + newL + (old_alpha_group^2-new_alpha_group^2)/(2*BBB2\$sd_prior_alpha^2)
r                    <-  runif(length(group_id))

check    	     <-  log(r) < A

group_idX   <- which(check)
group_id    <- group_id[group_idX]

if(length(group_idX)>0){

#CCode
.Call("update_group2",
as.integer(group_id),
as.integer(BBB2\$GROUP),
BBB2\$d_alpha,
new_alpha,
BBB2\$acc_alpha
)

}

}# End of function update alpha new id
```

Try the BlockFeST package in your browser

Any scripts or data that you put into this service are public.

BlockFeST documentation built on May 2, 2019, 3:25 a.m.