R/COR_Fstd.R

In xinghuq/HierD: R package for calculating hierarchical diversity across metrixs and scales, from genes to ecosystem

###### updated on 09-08-2018, the function for correlation of pirwise Fst and distance 

COR_Fstd= function(x,d,ncode) { 
require(adegenet)
require(hierfstat)
g=read.genepop(x,ncode, quiet = FALSE)
PFst=pairwise.fst(g) # genind file
npops=length(levels(g$pop))
## if M is matrix, then
if (d==TRUE) {
if (is.matrix(d)==TRUE) { 
Dgeo=as.dist(d, diag = FALSE, upper = FALSE)
COR_Fstd=cor(PFst,Dgeo,method = "pearson")
return(COR_Fstd)
}
  else {
    print("d must be a matrix")
  }}
else { d==FALSE
M=matrix(data=0,nrow =npops ,ncol =npops) ## nrow=npops, ncol=npops, which is 16 here
colnames(M)=levels(g$pop)
rownames(M)=levels(g$pop)
for (i in 1:npops){ 
  for (j in 1:npops){
    M[i,j]=abs(i-j)
    }
}

Dgeo=as.dist(M, diag = FALSE, upper = FALSE)
# if we want test IBD by mantel test
# ibd <- mantel.randtest(Dgen2,Dgeo1)
# plot(ibd)
# plot(Dgeo1, Dgen2)
# abline(lm(Dgen2~Dgeo1), col="red",lty=2)
# if we want use the euclidean distance rather the real matrix distance
# Dgeo1=dist(M, method = "euclidean", diag = FALSE, upper = FALSE, p = 2)### this is euclidean distance
# COR_Fsteud= cor(Dgen2,Dgeo1,method = "pearson") #cor(c(matrix1), c(matrix2)).

if(class(PFst)!="matrix"& class(PFst)!="dist") stop("PFst has to be a matrix")
if(class(Dgeo)!="matrix"& class(Dgeo)!="dist") stop("Dgeo has to be a matrix")

if(sum(is.na(PFst))!=0|sum(is.na(Dgeo))!=0) stop("Missing data in the dataset")
if(length(PFst)!=length(Dgeo)) stop("Numbers of rows in PFst and Dgeo are not equal")

COR_Fstd=cor(PFst,Dgeo,method = "pearson")
return(list(pwFst=PFst,Dgeo=Dgeo,COR_Fstd=COR_Fstd))
}
}