[adegenet-commits] r661 - pkg/R

[noreply at r-forge.r-project.org]

Author: jombart
Date: 2010-09-08 16:47:56 +0200 (Wed, 08 Sep 2010)
New Revision: 661

Modified:
   pkg/R/Hs.R
   pkg/R/fstat.R
   pkg/R/genind2genpop.R
Log:
A few changes. Hs tolerates genind objects. genind2genpop now works without pop factor provided; developping pairwise Fst



Modified: pkg/R/Hs.R
===================================================================
--- pkg/R/Hs.R	2010-09-08 14:12:52 UTC (rev 660)
+++ pkg/R/Hs.R	2010-09-08 14:47:56 UTC (rev 661)
@@ -4,6 +4,9 @@
 Hs <- function(x, truenames=TRUE) {
 
     ## CHECKS
+    if(is.genind(x)){
+        x <- genind2genpop(x, quiet=TRUE)
+    }
     if(!is.genpop(x)) stop("x is not a valid genpop object")
     if(x at type=="PA") stop("not implemented for presence/absence markers")
 
@@ -14,6 +17,6 @@
     lX <- lapply(x.byloc, function(e) makefreq(e, quiet=TRUE, truenames=truenames)$tab)
     lres <- lapply(lX, function(X) 1- apply(X^2,1,sum))
     res <- apply(as.matrix(data.frame(lres)),1,mean)
- 
+
     return(res)
 } # end Hs

Modified: pkg/R/fstat.R
===================================================================
--- pkg/R/fstat.R	2010-09-08 14:12:52 UTC (rev 660)
+++ pkg/R/fstat.R	2010-09-08 14:47:56 UTC (rev 661)
@@ -34,18 +34,41 @@
 #
 # classical fst sensu Nei
 #
-pairwise.fst <- function(x, pop=NULL){
+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)){
+        pop(x) <- pop
+    }
+    temp <- pop(x)
+    if(is.null(temp)) stop("no grouping factor (pop) provided")
+    if(length(levels(temp)) < 2){
+        warning("There is only one pop - returning NULL")
+        return(NULL)
+    }
 
-    if(is.null(pop)) pop <- pop(x)
-    if(is.null(pop)) stop("no pop factor provided")
-    if(length(pop)!=nrow(x at tab)) stop("pop has a wrong length.")
 
-
     ## COMPUTATIONS ##
 
+    ## function to compute one Fst ##
+    f1 <- function(pop1, pop2){ # pop1 and pop2 are genind obj. with a single pop each
+        temp <- repool(pop1,pop2)
+        b <- mean(Hs(temp))
+        pop(temp) <- NULL
+        a <- Hs(temp)
+        return((a-b)/a)
+    }
 
 
+    ## compute pairwise Fst for all pairs
+    lx <- seppop(x)
+    temp <- pop(x)
+    levPop <- levels(temp)
+    allPairs <- combn(1:length(levPop), 2)
+
+    vecRes <- numeric()
+    for(i in 1:ncol(allPairs)){
+        vecRes[i] <- f1(lx[[allPairs[1,i]]], lx[[allPairs[2,i]]])
+    }
     return(res)
 } # end of pairwise.fst

Modified: pkg/R/genind2genpop.R
===================================================================
--- pkg/R/genind2genpop.R	2010-09-08 14:12:52 UTC (rev 660)
+++ pkg/R/genind2genpop.R	2010-09-08 14:47:56 UTC (rev 661)
@@ -7,7 +7,10 @@
   if(!is.genind(x)) stop("x is not a valid genind object")
   checkType(x)
 
-  if(is.null(x at pop) && is.null(pop)) stop("pop is not provided either in x or in pop")
+  if(is.null(x at pop) && is.null(pop)) {
+      if(!quiet) warning("\npop is not provided either in x or in pop - assuming one single group")
+      pop <- factor(rep(1, nrow(x at tab)))
+  }
 
   missing <- match.arg(missing)