###
#'@export
################################################################################
pp.sigma.loc<-function(x,y,dat=dat,diploid=TRUE,...){
if (is.genind(dat)) dat<-genind2hierfstat(dat)
dum<-names(dat)[1]
wpop<-dat[,1]==x | dat[,1]==y
ndat<-dat[wpop,]
ndat[ndat[,1]==x,1]<-1
ndat[ndat[,1]==y,1]<-2
return(wc(ndat,diploid,...)$sigma.loc)
}
###
#'@export
################################################################################
pp.fst<-function(dat=dat,diploid=TRUE,...){
cl<-match.call()
if (is.genind(dat)) dat<-genind2hierfstat(dat)
dat<-dat[order(dat[,1]),]
Pop<-dat[,1]
npop<-length(unique(Pop))
nloc<-dim(dat)[2]-1
if (diploid) nvc<-3 else nvc<-2
ppsl<-array(numeric(npop*npop*nloc*nvc),dim=c(npop,npop,nloc,nvc))
for (i in 1:(npop-1))
for (j in (i+1):npop){
#cat(i," ",j,"\n") #for debugging
# TODO: optimize by using function for 2 pops only. Needs a new function
ppsl[i,j,,]<-as.matrix(pp.sigma.loc(unique(Pop)[i],unique(Pop)[j],dat,diploid,...))
}
ovl<-apply(ppsl,c(1,2,4),sum)
ppfst<-ovl[,,1]/apply(ovl,c(1,2),sum)
res<-list(call=cl,fst.pp=ppfst,vc.per.loc=ppsl)
class(res)<-"pp.fst"
res
}
#' @method print pp.fst
#' @export
print.pp.fst<-function(x,...){
print(x$fst.pp)
invisible(x)
}
###
#'@export
################################################################################
boot.ppfst<-function(dat=dat,nboot=100,quant=c(0.025,0.975),diploid=TRUE,...){
cl<-match.call()
if (is.genind(dat)) dat<-genind2hierfstat(dat)
typ<-dim(dat)
if(length(dim(dat))==2){
#Pop<-dat[,1]
npop<-length(table(dat[,1]))
nloc<-dim(dat)[2]-1
if (diploid) nvc<-3 else nvc<-2
ppsl<-array(numeric(npop*npop*nloc*nvc),dim=c(npop,npop,nloc,nvc))
x<-unique(dat[,1])
if(is.factor(dat[,1])) dat[,1]<-as.integer(dat[,1])
for (i in 1:(npop-1)) {
for (j in (i+1):npop) {
#cat(i," ",j,"\n") #for debugging
ppsl[i,j,,]<-as.matrix(pp.sigma.loc(i,j,dat,diploid,...))
}
}
}
else
{
npop<-typ[1]
nloc<-typ[3]
ppsl<-dat
}
bppfst<-array(numeric(npop*npop*nboot),dim=c(npop,npop,nboot))
for (i in 1:nboot){
dum<-sample(nloc,replace=TRUE)
ovl<-apply(ppsl[,,dum,],c(1,2,4),sum)
bppfst[,,i]<-ovl[,,1]/apply(ovl,c(1,2),sum)
}
#browser()
ll<-apply(bppfst,c(1,2),stats::quantile,quant[1],na.rm=TRUE)
hl<-apply(bppfst,c(1,2),stats::quantile,quant[2],na.rm=TRUE)
dimnames(ll)[[1]]<-dimnames(ll)[[2]]<-sort(x)
dimnames(hl)<-dimnames(ll)
res<-(list(call=cl,ll=ll,ul=hl,vc.per.loc=ppsl))
class(res)<-"boot.ppfst"
res
}
#' @method print boot.ppfst
#' @export
print.boot.ppfst<-function(x,...){
a<-dim(x$ll)[1]
ci<-matrix(nrow=a,ncol=a)
dimnames(ci)<-dimnames(x$ll)
for (i in 2:a)
for (j in 1:(a-1)){
ci[j,i]<-x$ul[j,i]
ci[i,j]<-x$ll[j,i]
}
cat("\n Upper limit above diagonal \n")
cat(" Lower limit below diagonal \n \n")
print(ci)
invisible(x)
}
###
#'@export
################################################################################
boot.ppfis<-function(dat=dat,nboot=100,quant=c(0.025,0.975),diploid=TRUE,dig=4,...){
cl<-match.call()
if (is.genind(dat)) dat<-genind2hierfstat(dat)
bs<-basic.stats(dat)
Ho<-bs$Ho
Hs<-bs$Hs
nloc<-dim(Hs)[1]
npop<-dim(Hs)[2]
my.boot<-matrix(numeric(nboot*npop),ncol=npop)
for (i in 1:nboot){
x<-sample(nloc,replace=TRUE)
my.boot[i,]<-1-colSums(Ho[x,],na.rm=TRUE)/colSums(Hs[x,],na.rm=TRUE)
}
ll<-apply(my.boot,2,stats::quantile,quant[1],na.rm=TRUE)
hl<-apply(my.boot,2,stats::quantile,quant[2],na.rm=TRUE)
res<-data.frame(ll=ll,hl=hl)
return(list(call=cl,fis.ci=round(res,digits=dig)))
}
###########################################################################
#' Estimates bootstrap confidence intervals for pairwise betas FST estimates
#'
#' Estimates bootstrap confidence intervals for pairwise betas FST estimates.
#'
#' @usage boot.ppbetas(dat=dat,nboot=100,quant=c(0.025,0.975),diploid=TRUE,digits=4)
#'
#' @param dat A data frame containing population of origin as the first column
#' and multi-locus genotypes in following columns
#' @param nboot the number of bootstrap samples to draw
#' @param quant the limit of the confidence intervals
#' @param diploid whether the data is from a diploid (default) or haploid organism
#' @param digits how many digits to print out
#'
#' @return a matrix with upper limit of the bootstrap CI
#' above the diagonal and lower limit below the diagonal
#'
#' @seealso \link{betas} \link{pairwise.betas}
#'
#' @examples
#' \dontrun{
#' data(gtrunchier)
#' boot.ppbetas(gtrunchier[,-2])
#' }
#' @export
#####################################################################
boot.ppbetas<-function(dat=dat,nboot=100,quant=c(0.025,0.975),diploid=TRUE,digits=4){
if (is.genind(dat)) dat<-genind2hierfstat(dat)
dat<-dat[order(dat[,1]),]
pops<-unique(dat[,1])
npop<-length(pops)
fstmat <- matrix(nrow=npop,ncol=npop,dimnames=list(pops,pops))
if (is.factor(dat[,1])) {
dat[,1]<-as.numeric(dat[,1])
pops<-as.numeric(pops)
}
for(a in 2:npop){
for(b in 1:(a-1)){
subdat <- dat[dat[,1] == pops[a] | dat[,1]==pops[b],]
tmp<-betas(subdat,nboot=nboot,lim=quant,diploid=diploid)$ci[,3]
fstmat[a,b]<-tmp[1]
fstmat[b,a]<-tmp[2]
}
}
round(fstmat,digits = digits)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.