[adegenet-commits] r758 - pkg/R

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

Author: jombart
Date: 2011-01-07 17:30:42 +0100 (Fri, 07 Jan 2011)
New Revision: 758

Modified:
   pkg/R/glFunctions.R
Log:
in the middle of glPca


Modified: pkg/R/glFunctions.R
===================================================================
--- pkg/R/glFunctions.R	2011-01-07 14:12:57 UTC (rev 757)
+++ pkg/R/glFunctions.R	2011-01-07 16:30:42 UTC (rev 758)
@@ -72,7 +72,6 @@
 
 
 
-
 ########
 ## glVar
 ########
@@ -101,5 +100,87 @@
 
 
 
+
+
+
+
+
+########
+## glPca
+########
+##
+## PCA for genlight objects
+##
+glPca <- function(x, center=TRUE, scale=FALSE, scannf=TRUE, nf=2, loadings=TRUE){
+    if(!inherits(x, "genlight")) stop("x is not a genlight object")
+
+    ## COMPUTE MEANS AND VARIANCES ##
+    if(center) {
+        vecMeans <- glMean(x)
+        if(any(is.na(vecMeans))) stop("NAs detected in the vector of means")
+    } else {
+        vecMeans <- 0
+    }
+    if(scale){
+        if(any(is.na(vecVar))) stop("NAs detected in the vector of variances")
+        vecSd <- sqrt(glVar(x))
+    } else {
+        vecSd <- 1
+    }
+
+
+    ## COMPUTE DOT PRODUCTS BETWEEN GENOTYPES ##
+    dotProd <- function(a,b){ # a and b are two SNPbin objects
+        res <- sum( ((as.integer(a)-vecMeans)/vecSd) * ((as.integer(b)-vecMeans)/vecSd), na.rm=TRUE )
+        return(res)
+    }
+
+    allComb <- combn(1:nInd(x), 2)
+    allProd <- unlist(lapply(1:ncol(allComb), function(i) dotProd(x at gen[[allComb[1,i]]]), x at gen[[allComb[2,i]]]))
+    allProd <- allProd / nInd(x) # assume uniform weights
+
+    ## shape result as a matrix
+    attr(allProd,"Size") <- nInd(x)
+    attr(allProd,"Diag") <- FALSE
+    attr(allProd,"Upper") <- FALSE
+    class(allProd) <- "dist"
+    allProd <- as.matrix(allProd)
+
+
+    ## PERFORM EIGENANALYSIS ##
+    ## eigenanalysis
+    eigRes <- eigen(allProd, symmetric=TRUE, only.values=FALSE)
+    rank <- sum(eigRes > 1e-10)
+    eigRes$values <- eigRes$value[1:rank]
+    eigRes$vectors <- eigRes$vectors[1:rank]
+
+    ## scan  nb of axes retained
+    if(scannf){
+        barplot(eigRes$values, main="Eigenvalues", col=heat.colors(rank))
+        cat("Select the number of axes: ")
+        nf <- as.integer(readLines(n = 1))
+    }
+
+    ## rescale PCs
+    res <- list()
+    ## use: li = XQU = V\Lambda^(1/2)
+    res$scores <- sweep(eigRes$vectors[, 1:nf],2, sqrt(eigRes$values[1:nf]), FUN="*")
+
+    ## get loadings
+    if(loadings){
+        ## use: c1 = X^TDV
+        ##### TO FINISH !!!
+        res$loadings <- matLoad
+    }
+    ## compute loadings
+
+    names(res) <- locNames(x)
+    return(res)
+
+} # glPca
+
+
+
+
 ## TESTING ##
 ## all.equal(glVar(x), apply(as.matrix(x), 2, function(e) mean((e-mean(e, na.rm=TRUE))^2, na.rm=TRUE)))