Skip to content

Commit

Permalink
fixed fstat things coming from adegenet to work with current version …
Browse files Browse the repository at this point in the history
…2.0.0 of adegenet; package installs; check fails but unrelated errors
  • Loading branch information
thibautjombart committed Jun 22, 2015
1 parent fa89487 commit b63d8dc
Show file tree
Hide file tree
Showing 8 changed files with 82 additions and 22 deletions.
19 changes: 6 additions & 13 deletions R/fstat.from.adegenet.R
Original file line number Diff line number Diff line change
Expand Up @@ -28,8 +28,6 @@
#' to any grouping of individuals'.
#' @param res.type the type of result to be returned: a \code{dist} object, or
#' a symmetric matrix
#' @param truenames a logical indicating whether true labels (as opposed to
#' generic labels) should be used to name the output.
#' @param fstonly a logical stating whether only the Fst should be returned.
#' @return A vector, a matrix, or a dist object containing F statistics.
#' @author Thibaut Jombart \email{t.jombart@@imperial.ac.uk}
Expand Down Expand Up @@ -84,15 +82,14 @@ fstat <- function(x, pop=NULL, fstonly=FALSE){
## misc checks
if(!is.genind(x)) stop("x is not a valid genind object")
## if(!require(hierfstat)) stop("hierfstat package is required. Please install it.")
if(x@ploidy != as.integer(2)) stop("not implemented for non-diploid genotypes")
checkType(x)
if(any(x@ploidy != 2L)) stop("not implemented for non-diploid genotypes")

if(is.null(pop)) pop <- x@pop
if(is.null(pop)) pop <- pop(x)
if(is.null(pop)) stop("no pop factor provided")
if(length(pop)!=nrow(x@tab)) stop("pop has a wrong length.")
if(length(pop)!=nInd(x)) stop("pop has a wrong length.")

## computations
dat <- genind2hierfstat(x)[,-1]
dat <- .genind2hierfstat(x)[,-1]
res <- varcomp.glob(levels=data.frame(pop), loci=dat)$F

if(fstonly) {res <- res[1,1]}
Expand All @@ -108,7 +105,7 @@ fstat <- function(x, pop=NULL, fstonly=FALSE){
#' @export
#' @aliases pairwise.fst
#'
pairwise.fst <- function(x, pop=NULL, res.type=c("dist","matrix"), truenames=TRUE){
pairwise.fst <- function(x, pop=NULL, res.type=c("dist","matrix")){
## MISC CHECKS ##
if(!is.genind(x)) stop("x is not a valid genind object")
if(!is.null(pop)){
Expand Down Expand Up @@ -159,11 +156,7 @@ pairwise.fst <- function(x, pop=NULL, res.type=c("dist","matrix"), truenames=TRU

if(res.type=="matrix"){
res <- as.matrix(res)
if(truenames){
lab <- x@pop.names
} else {
lab <- names(x@pop.names)
}
lab <- popNames(x)

colnames(res) <- rownames(res) <- lab
}
Expand Down
73 changes: 70 additions & 3 deletions R/gstat.randtest.R
Original file line number Diff line number Diff line change
@@ -1,3 +1,70 @@


###################
# Function .genlab
###################
# recursive function to have labels of constant length
# base = a character string
# n = number of labels
.genlab <- function(base, n) {
f1 <- function(cha,n){
if(nchar(cha)<n){
cha <- paste("0",cha,sep="")
return(f1(cha,n))
} else {return(cha)}
}
w <- as.character(1:n)
max0 <- max(nchar(w))
w <- sapply(w, function(cha) f1(cha,max0))
return(paste(base,w,sep=""))
}



############################
# Function genind2hierfstat - not exported
############################
.genind2hierfstat <- function(x, pop=NULL){
## if(!inherits(x,"genind")) stop("x must be a genind object (see ?genind)")
## invisible(validObject(x))
if(!is.genind(x)) stop("x is not a valid genind object")
if(any(ploidy(x) != 2L)) stop("not implemented for non-diploid genotypes")
## checkType(x)

if(is.null(pop)) pop <- pop(x)
if(is.null(pop)) pop <- as.factor(rep("P1",nrow(x@tab)))

## ## NOTES ON THE CODING IN HIERFSTAT ##
## - interpreting function is genot2al
## - same coding has to be used for all loci
## (i.e., all based on the maximum number of digits to be used)
## - alleles have to be coded as integers
## - alleles have to be sorted by increasing order when coding a genotype
## - for instance, 121 is 1/21, 101 is 1/1, 11 is 1/1

## find max number of alleles ##
max.nall <- max(nAll(x))
x@all.names <- lapply(alleles(x), function(e) .genlab("",max.nall)[1:length(e)])


## VERSION USING GENIND2DF ##
gen <- genind2df(x, sep="", usepop=FALSE)
gen <- as.matrix(data.frame(lapply(gen, as.numeric)))
res <- cbind(as.numeric(pop),as.data.frame(gen))
colnames(res) <- c("pop", locNames(x))
if(!any(table(indNames(x))>1)){
rownames(res) <- indNames(x)
} else {
warning("non-unique labels for individuals; using generic labels")
rownames(res) <- 1:nrow(res)
}

return(res)
} # end .genind2hierfstat




##########################
## Function gstat.randtest
##########################
Expand Down Expand Up @@ -78,8 +145,8 @@ gstat.randtest <- function(x,pop=NULL, method=c("global","within","between"),
## return(invisible())

if(!is.genind(x)) stop("x is not a valid genind object")
if(x@ploidy != as.integer(2)) stop("not implemented for non-diploid genotypes")
checkType(x)
if(any(ploidy(x) != 2L)) stop("not implemented for non-diploid genotypes")
## checkType(x)
## if(!require(hierfstat)) stop("hierfstat package is required. Please install it.")

if(is.null(pop)) pop <- x@pop
Expand All @@ -91,7 +158,7 @@ gstat.randtest <- function(x,pop=NULL, method=c("global","within","between"),
if(met=="between" && is.null(sub.pop)) stop("Method 'between' chosen but 'sub.pop' is not provided.")

## make data for hierfstat
X <- genind2hierfstat(x=x,pop=pop)
X <- .genind2hierfstat(x=x,pop=pop)

## compute obs gstat
obs <- g.stats.glob(X)$g.stats
Expand Down
2 changes: 1 addition & 1 deletion man/betai.unw.Rd
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
% Generated by roxygen2 (4.1.0): do not edit by hand
% Generated by roxygen2 (4.1.1): do not edit by hand
% Please edit documentation in R/my.betai.R
\name{betai.unw}
\alias{betai.unw}
Expand Down
2 changes: 1 addition & 1 deletion man/betas.Rd
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
% Generated by roxygen2 (4.1.0): do not edit by hand
% Generated by roxygen2 (4.1.1): do not edit by hand
% Please edit documentation in R/my.betai.R
\name{betas}
\alias{betas}
Expand Down
2 changes: 1 addition & 1 deletion man/fstat.from.adegenet.Rd
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
% Generated by roxygen2 (4.1.0): do not edit by hand
% Generated by roxygen2 (4.1.1): do not edit by hand
% Please edit documentation in R/fstat.from.adegenet.R
\name{fstat}
\alias{FST}
Expand Down
2 changes: 1 addition & 1 deletion man/gstat.randtest.Rd
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
% Generated by roxygen2 (4.1.0): do not edit by hand
% Generated by roxygen2 (4.1.1): do not edit by hand
% Please edit documentation in R/gstat.randtest.R
\name{gstat.randtest}
\alias{gstat.randtest}
Expand Down
2 changes: 1 addition & 1 deletion man/sim.genot.t.Rd
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
% Generated by roxygen2 (4.1.0): do not edit by hand
% Generated by roxygen2 (4.1.1): do not edit by hand
% Please edit documentation in R/simgenot_new.R
\name{sim.genot.t}
\alias{sim.genot.t}
Expand Down
2 changes: 1 addition & 1 deletion man/write.bayescan.Rd
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
% Generated by roxygen2 (4.1.0): do not edit by hand
% Generated by roxygen2 (4.1.1): do not edit by hand
% Please edit documentation in R/writebayescan.R
\name{write.bayescan}
\alias{write.bayescan}
Expand Down

0 comments on commit b63d8dc

Please sign in to comment.