[adegenet-commits] r612 - in www: . files/patches
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue May 4 14:40:24 CEST 2010
Author: jombart
Date: 2010-05-04 14:40:24 +0200 (Tue, 04 May 2010)
New Revision: 612
Added:
www/files/patches/loadingplot.R
Removed:
www/files/patches/classes.R
www/files/patches/genind2df.R
www/files/patches/makefreq.R
www/files/patches/propShared.R
www/files/patches/spca.R
Modified:
www/download.html
Log:
Added a patch for loadingplot; removed old patches.
Modified: www/download.html
===================================================================
--- www/download.html 2010-05-04 12:37:28 UTC (rev 611)
+++ www/download.html 2010-05-04 12:40:24 UTC (rev 612)
@@ -44,8 +44,11 @@
as') in your working directory and
type <span style="font-family: monospace;">source("[your-patch-file.R]")</span>
to
-use a patch.<br>
-- no patch for now<br>
+use
+a patch.<br>
+- <a href="files/patches/loadingplot.R"><span
+ style="font-family: monospace;">loadingplot.R</span></a>: a few checks
+added when the 'at' argument the specified<br>
<br>
<img alt="" src="images/bullet.png" style="width: 10px; height: 10px;">
<span style="font-weight: bold;">Older
Deleted: www/files/patches/classes.R
===================================================================
--- www/files/patches/classes.R 2010-05-04 12:37:28 UTC (rev 611)
+++ www/files/patches/classes.R 2010-05-04 12:40:24 UTC (rev 612)
@@ -1,457 +0,0 @@
-########################################################################
-# adegenet classes definitions. All classes are S4.
-#
-# Thibaut Jombart, November 2007
-# t.jombart at imperial.ac.uk
-########################################################################
-
-###############################
-# Two classes of R object are
-# defined :
-# gen - common part to genind and genpop
-# genind - genotypes of individuals
-# genpop - allelic frequencies of populations
-###############################
-
-
-###############################################################
-###############################################################
-# AUXILIARY FUNCTIONS
-###############################################################
-###############################################################
-
-
-
-
-###############################################################
-###############################################################
-# CLASSES DEFINITION
-###############################################################
-###############################################################
-
-##.initAdegenetClasses <- function(){
-
-
-####################
-# Unions of classes
-####################
-setClassUnion("listOrNULL", c("list","NULL"))
-setClassUnion("factorOrNULL", c("factor","NULL"))
-setClassUnion("charOrNULL", c("character","NULL"))
-setClassUnion("callOrNULL", c("call","NULL"))
-setClassUnion("intOrNum", c("integer","numeric","NULL"))
-
-
-
-####################
-# virtual class gen
-####################
-.gen.valid <- function(object){
- # this function tests only the consistency
- # of the length of each component
- p <- ncol(object at tab)
- k <- length(unique(object at loc.names))
-
-
- if(!is.null(object at loc.fac)){
- if(length(object at loc.fac) != p) {
- cat("\ninvalid length for loc.fac\n")
- return(FALSE)
- }
-
- if(length(levels(object at loc.fac)) != k) {
- cat("\ninvalid number of levels in loc.fac\n")
- return(FALSE)
- }
- }
-
- if(!is.null(object at loc.nall)){
- if(length(object at loc.nall) != k) {
- cat("\ninvalid length in loc.nall\n")
- return(FALSE)
- }
- }
-
- temp <- table(object at loc.names[object at loc.names!=""])
- if(any(temp>1)) {
- warning("\nduplicate names in loc.names:\n")
- print(temp[temp>1])
- }
-
- if(!is.null(object at all.names)){
- if(length(unlist(object at all.names)) != p) {
- cat("\ninvalid length in all.names\n")
- return(FALSE)
- }
- }
-
- return(TRUE)
-
-}# end .gen.valid
-
-
-setClass("gen", representation(tab = "matrix",
- loc.names = "character",
- loc.fac = "factorOrNULL",
- loc.nall = "intOrNum",
- all.names = "listOrNULL",
- call = "callOrNULL",
- "VIRTUAL"),
- prototype(tab=matrix(ncol=0,nrow=0), loc.nall=integer(0), call=NULL))
-
-setValidity("gen", .gen.valid)
-
-
-
-
-
-########################
-# virtual class indInfo
-########################
-setClass("indInfo", representation(ind.names = "character",
- pop = "factorOrNULL",
- pop.names = "charOrNULL",
- ploidy = "integer",
- type = "character",
- other = "listOrNULL", "VIRTUAL"),
- prototype(pop=NULL, pop.names = NULL, type = "codom", ploidy = as.integer(2), other = NULL))
-
-
-
-
-
-###############
-# Class genind
-###############
-.genind.valid <- function(object){
- if(!.gen.valid(object)) return(FALSE)
-
- if(length(object at ind.names) != nrow(object at tab)) {
- cat("\ninvalid length in ind.names\n")
- return(FALSE)
- }
-
- temp <- table(object at ind.names[object at ind.names!=""])
- if(any(temp>1)) {
- warning("\nduplicate names in ind.names:\n")
- print(temp[temp>1])
- }
-
- if(!is.null(object at pop)){ # check pop
-
- if(length(object at pop) != nrow(object at tab)) {
- cat("\npop is given but has invalid length\n")
- return(FALSE)
- }
-
- if(is.null(object at pop.names)) {
- cat("\npop is provided without pop.names")
- }
-
-
- if(length(object at pop.names) != length(levels(object at pop))) {
- cat("\npop.names has invalid length\n")
- return(FALSE)
- }
-
- temp <- table(object at pop.names[object at pop.names!=""])
- if(any(temp>1)) {
- warning("\nduplicate names in pop.names:\n")
- print(temp[temp>1])
- }
-
- } # end check pop
-
- ## check ploidy
- if(object at ploidy < as.integer(1)){
- cat("\nploidy inferior to 1\n")
- return(FALSE)
- }
-
- ## check type of marker
- if(!object at type %in% c("codom","PA") ){
- cat("\nunknown type of marker\n")
- return(FALSE)
- }
-
-
- return(TRUE)
-} #end .genind.valid
-
-setClass("genind", contains=c("gen", "indInfo"))
-setValidity("genind", .genind.valid)
-
-
-
-########################
-# virtual class popInfo
-########################
-setClass("popInfo", representation(pop.names = "character", ploidy = "integer",
- type = "character", other = "listOrNULL", "VIRTUAL"),
- prototype(type = "codom", ploidy = as.integer(2), other = NULL))
-
-
-
-###############
-# Class genpop
-###############
-.genpop.valid <- function(object){
- if(!.gen.valid(object)) return(FALSE)
- if(length(object at pop.names) != nrow(object at tab)) {
- cat("\ninvalid length in pop.names\n")
- return(FALSE)
- }
-
- temp <- table(object at pop.names[object at pop.names!=""])
- if(any(temp>1)) {
- warning("\nduplicate names in pop.names:\n")
- print(temp[temp>1])
- }
-
- ## check ploidy
- if(object at ploidy < as.integer(1)){
- cat("\nploidy inferior to 1\n")
- return(FALSE)
- }
-
- ## check type of marker
- if(!object at type %in% c("codom","PA") ){
- cat("\nunknown type of marker\n")
- return(FALSE)
- }
-
- return(TRUE)
-} #end .genpop.valid
-
-setClass("genpop", contains=c("gen", "popInfo"))
-setValidity("genpop", .genpop.valid)
-
-
-
-
-
-
-
-###############################################################
-###############################################################
-# MAIN CLASS METHODS
-###############################################################
-###############################################################
-
-
-
-#################
-# Function names
-#################
-setMethod("names", signature(x = "genind"), function(x){
- return(slotNames(x))
-})
-
-setMethod("names", signature(x = "genpop"), function(x){
- return(slotNames(x))
-})
-
-
-
-
-
-##################
-# Function genind
-##################
-## constructor of a genind object
-genind <- function(tab,pop=NULL,prevcall=NULL,ploidy=2,type=c("codom","PA")){
- ## handle arguments
- X <- as.matrix(tab)
- if(is.null(colnames(X))) stop("tab columns have no name.")
- if(is.null(rownames(X))) {rownames(X) <- 1:nrow(X)}
-
- type <- match.arg(type)
- ploidy <- as.integer(ploidy)
- nind <- nrow(X)
-
-
- ## HANDLE LABELS ##
-
- ## loc names is not type-dependent
- temp <- colnames(X)
- ## temp <- gsub("[.].*$","",temp)
- temp <- gsub("[.][^.]*$", "", temp)
- temp <- .rmspaces(temp)
- loc.names <- unique(temp)
- nloc <- length(loc.names)
- loc.codes <- .genlab("L",nloc)
- names(loc.names) <- loc.codes
-
- ## ind names is not type-dependent either
- ind.codes <- .genlab("", nind)
- ind.names <- .rmspaces(rownames(X))
- names(ind.names) <- ind.codes
- rownames(X) <- ind.codes
-
-
- if(type=="codom"){
- ## loc.nall
- loc.nall <- table(temp)[match(loc.names,names(table(temp)))]
- loc.nall <- as.integer(loc.nall)
- names(loc.nall) <- loc.codes
-
- ## loc.fac
- loc.fac <- rep(loc.codes,loc.nall)
-
- ## alleles name
- temp <- colnames(X)
- temp <- gsub("^.*[.]","",temp)
- temp <- .rmspaces(temp)
- all.names <- split(temp,loc.fac)
- all.codes <- lapply(all.names,function(e) .genlab("",length(e)))
- for(i in 1:length(all.names)){
- names(all.names[[i]]) <- all.codes[[i]]
- }
-
- colnames(X) <- paste(loc.fac,unlist(all.codes),sep=".")
- loc.fac <- as.factor(loc.fac)
- } else { # end if type=="codom" <=> if type=="PA"
- colnames(X) <- loc.codes
- loc.fac <- NULL
- all.names <- NULL
- loc.nall <- NULL
- }
-
- ## Ideally I should use an 'initialize' method here
- res <- new("genind")
- res at tab <- X
- res at ind.names <- ind.names
- res at loc.names <- loc.names
- res at loc.nall <- loc.nall
- res at loc.fac <- loc.fac
- res at all.names <- all.names
-
- ## populations name (optional)
- ## beware, keep levels of pop sorted in
- ## there order of appearance
- if(!is.null(pop)) {
- # convert pop to a factor if it is not
- if(!is.factor(pop)) {pop <- factor(pop)}
- pop.lab <- .genlab("P",length(levels(pop)) )
- # put pop levels in appearance order
- pop <- as.character(pop)
- pop <- factor(pop, levels=unique(pop))
- temp <- pop
- # now levels are correctly ordered
- levels(pop) <- pop.lab
- res at pop <- pop
- pop.names <- as.character(levels(temp))
- names(pop.names) <- as.character(levels(res at pop))
- res at pop.names <- pop.names
- }
-
- ## ploidy
- plo <- as.integer(ploidy)
- if(plo < as.integer(1)) stop("ploidy inferior to 1")
- res at ploidy <- plo
-
- ## type of marker
- res at type <- as.character(type)
-
- if(is.null(prevcall)) {prevcall <- match.call()}
- res at call <- prevcall
-
- return(res)
-
-} # end genind
-
-######################
-# alias for as.genind
-######################
-as.genind <- genind
-
-
-
-##################
-# Function genpop
-##################
-genpop <- function(tab,prevcall=NULL,ploidy=as.integer(2),type=c("codom","PA")){
-
- ## handle args
- X <- as.matrix(tab)
- if(is.null(colnames(X))) stop("tab columns have no name.")
- if(is.null(rownames(X))) {rownames(X) <- 1:nrow(X)}
-
- type <- match.arg(type)
- ploidy <- as.integer(ploidy)
- npop <- nrow(X)
-
-
- ## HANDLE LABELS ##
-
- ## loc names is not type-dependent
- temp <- colnames(X)
- ## temp <- gsub("[.].*$","",temp)
- temp <- gsub("[.][^.]*$", "", temp)
- temp <- .rmspaces(temp)
- loc.names <- unique(temp)
- nloc <- length(loc.names)
- loc.codes <- .genlab("L",nloc)
- names(loc.names) <- loc.codes
-
- ## pop names is not type-dependent either
- pop.codes <- .genlab("", npop)
- pop.names <- .rmspaces(rownames(X))
- names(pop.names) <- pop.codes
- rownames(X) <- pop.codes
-
- ## type-dependent stuff
- if(type=="codom"){
- ## loc.nall
- loc.nall <- table(temp)[match(loc.names,names(table(temp)))]
- loc.nall <- as.integer(loc.nall)
- names(loc.nall) <- loc.codes
-
- ## loc.fac
- loc.fac <- rep(loc.codes,loc.nall)
-
- ## alleles name
- temp <- colnames(X)
- temp <- gsub("^.*[.]","",temp)
- temp <- .rmspaces(temp)
- all.names <- split(temp,loc.fac)
- all.codes <- lapply(all.names,function(e) .genlab("",length(e)))
- for(i in 1:length(all.names)){
- names(all.names[[i]]) <- all.codes[[i]]
- }
-
- rownames(X) <- pop.codes
- colnames(X) <- paste(loc.fac,unlist(all.codes),sep=".")
- loc.fac <- as.factor(loc.fac)
- } else { # end if type=="codom" <=> if type=="PA"
- colnames(X) <- loc.codes
- loc.fac <- NULL
- all.names <- NULL
- loc.nall <- NULL
- }
-
- res <- new("genpop")
-
- res at tab <- X
- res at pop.names <- pop.names
- res at loc.names <- loc.names
- res at loc.nall <- loc.nall
- res at loc.fac <- loc.fac
- res at all.names <- all.names
- res at ploidy <- ploidy
- res at type <- as.character(type)
-
- if(is.null(prevcall)) {prevcall <- match.call()}
- res at call <- prevcall
-
- return(res)
-
-} # end genpop
-
-
-
-######################
-# alias for as.genpop
-######################
-as.genpop <- genpop
-
Deleted: www/files/patches/genind2df.R
===================================================================
--- www/files/patches/genind2df.R 2010-05-04 12:37:28 UTC (rev 611)
+++ www/files/patches/genind2df.R 2010-05-04 12:40:24 UTC (rev 612)
@@ -1,78 +0,0 @@
-#####################
-# Function genind2df
-#####################
-genind2df <- function(x, pop=NULL, sep="", usepop=TRUE, oneColPerAll=FALSE){
-
- if(!is.genind(x)) stop("x is not a valid genind object")
- ## checkType(x)
-
- if(is.null(pop)) {
- pop <- x at pop
- levels(pop) <- x at pop.names
- }
-
- if(oneColPerAll){
- sep <- "/"
- }
-
- ## PA case ##
- if(x at type=="PA"){
- temp <- truenames(x)
- if(is.list(temp) & usepop){
- res <- cbind.data.frame(pop=temp[[2]],temp[[1]])
- } else{
- if(is.list(temp)) {
- res <- temp[[1]]
- } else{
- res <- temp
- }
- }
-
- return(res) # exit here
- }
-
- ## codom case ##
- # make one table by locus from x at tab
- kX <- seploc(x,res.type="matrix")
- kX <- lapply(kX, function(X) round(X*x at ploidy)) # take data as numbers of alleles
- ## (kX is a list of nloc tables)
-
- ## function to recode a genotype in form "A1[sep]...[sep]Ak" from frequencies
- recod <- function(vec,lab){
- if(any(is.na(vec))) return(NA)
- res <- paste( rep(lab,vec), collapse=sep)
- return(res)
- }
-
-
- # kGen is a list of nloc vectors of genotypes
- kGen <- lapply(1:length(kX), function(i) apply(kX[[i]],1,recod,x at all.names[[i]]))
- names(kGen) <- x at loc.names
-
- ## if use one column per allele
- if(oneColPerAll){
- f1 <- function(vec){ # to repeat NA with seperators
- vec[is.na(vec)] <- paste(rep("NA",x at ploidy), collapse=sep)
- return(vec)
- }
- temp <- lapply(kGen, f1)
- temp <- lapply(temp, strsplit,sep)
-
- res <- lapply(temp, function(e) matrix(unlist(e), ncol=x at ploidy, byrow=TRUE))
- res <- data.frame(res,stringsAsFactors=FALSE)
- names(res) <- paste(rep(locNames(x),each=x at ploidy), 1:x at ploidy, sep=".")
-
- ## handle pop here
- if(!is.null(pop) & usepop) res <- cbind.data.frame(pop,res,stringsAsFactors=FALSE)
-
- return(res) # exit here
- } # end if oneColPerAll
-
- ## build the final data.frame
- res <- cbind.data.frame(kGen,stringsAsFactors=FALSE)
-
- ## handle pop here
- if(!is.null(pop) & usepop) res <- cbind.data.frame(pop,res,stringsAsFactors=FALSE)
-
- return(res)
-}
Added: www/files/patches/loadingplot.R
===================================================================
--- www/files/patches/loadingplot.R (rev 0)
+++ www/files/patches/loadingplot.R 2010-05-04 12:40:24 UTC (rev 612)
@@ -0,0 +1,76 @@
+##############
+# loadingplot
+##############
+loadingplot <- function(x, at=NULL, threshold=quantile(x,0.75), axis=1, fac=NULL, byfac=FALSE,
+ lab=names(x), cex.lab=0.7, cex.fac=1, lab.jitter=0,
+ main="Loading plot", xlab="Variables", ylab="Loadings",...){
+ ## some checks
+ if(is.data.frame(x) || is.matrix(x)){
+ temp <- rownames(x)
+ x <- x[,axis]
+ names(x) <- temp
+ }
+ if(!is.numeric(x)) stop("x is not numeric")
+ if(any(is.na(x))) stop("NA entries in x")
+ if(any(x<0)) {
+ warning("Some values in x are less than 0\n Using abs(x) instead, but this might not be optimal.")
+ x <- abs(x)
+ }
+ if(is.null(at)){
+ at <- 1:length(x)
+ } else {
+ if(length(at) != length(x)) stop("x and at do not have the same length.")
+ }
+
+ ## preliminary computations
+ y.min <- min(min(x),0)
+ y.max <- max(max(x),0)
+ y.offset <- (y.max-y.min)*0.02
+ if(is.null(lab)) {lab <- 1:length(x)}
+
+ if(!is.null(fac)){
+ if(byfac){
+ x <- tapply(x, fac, mean)
+ if(length(lab) != length(x)) lab <- names(x)
+ } else {
+ fac <- factor(fac, levels=unique(fac))
+ grp.idx <- cumsum(table(fac)) + 0.5
+ grp.lab.idx <- tapply(1:length(x), fac, mean)
+ grp.lab <- names(grp.idx)
+ grp.idx <- grp.idx[-length(grp.idx)]
+ }
+ } # end fac handling
+
+
+ ## start the plot
+ dat <- cbind(at, x)
+ plot(dat, type="h", xlab=xlab, ylab=ylab,
+ main=main, xaxt="n", ylim=c(y.min,y.max*1.2), ...)
+
+ ## add groups of variables (optional)
+ if(!is.null(fac) & !byfac) {
+ abline(v=grp.idx,lty=2) # split groups of variables
+ text(x=grp.lab.idx,y=y.max*1.15, labels=grp.lab, cex=cex.fac) # annotate groups
+ }
+
+ ## annotate variables that are above the threshold
+ x.ann <- at[x > threshold]
+ x.ann <- jitter(x.ann,fac=lab.jitter)
+
+ y.ann <- x[x > threshold] + y.offset
+ y.ann <- jitter(y.ann,fac=lab.jitter)
+
+ txt.ann <- lab[x > threshold]
+ text(x=x.ann, y=y.ann, label=txt.ann, cex=cex.lab)
+
+ ## indicate the threshold
+ abline(h=threshold, col="grey")
+
+ ## build the result
+ res <- list(threshold=threshold,
+ var.names=txt.ann,
+ var.idx=which(x > threshold),
+ var.values=x[x > threshold])
+ return(invisible(res))
+
+} # end loadingplot
Deleted: www/files/patches/makefreq.R
===================================================================
--- www/files/patches/makefreq.R 2010-05-04 12:37:28 UTC (rev 611)
+++ www/files/patches/makefreq.R 2010-05-04 12:40:24 UTC (rev 612)
@@ -1,58 +0,0 @@
-####################
-# Function makefreq
-####################
-makefreq <- function(x,quiet=FALSE,missing=NA,truenames=TRUE){
-
- if(!is.genpop(x)) stop("x is not a valid genpop object")
- ##if(x at type=="PA") stop("frequencies not computable for presence/asbence data")
- checkType(x)
-
- if(!quiet) cat("\n Finding allelic frequencies from a genpop object... \n")
-
- f1 <- function(v){
- if(all(is.na(v)) || sum(v,na.rm=TRUE)==0) return(rep(NA,length(v)))
- return(v/(sum(v,na.rm=TRUE)))
- }
-
- res <- list()
-
- tabcount <- x at tab
-
- eff.pop <- t(apply(tabcount,1,function(r) tapply(r,x at loc.fac,sum,na.rm=TRUE)))
- if(nLoc(x)==1){ # fix for nloc==1
- eff.pop <- t(eff.pop)
- }
-
- # tabfreq is a pop x loci table of allelic frequencies
- tabfreq <- t(apply(tabcount,1,function(r) unlist(tapply(r,x at loc.fac,f1))))
- colnames(tabfreq) <- colnames(x at tab)
-
- # NA treatment
- # NA can be kept as is, or replaced 0 or by the mean frequency of the allele.
- if(!is.na(missing)){
- if(missing==0) tabfreq[is.na(tabfreq)] <- 0
- if(toupper(missing)=="MEAN") {
- moy <- apply(tabfreq,2,function(c) mean(c,na.rm=TRUE))
- for(j in 1:ncol(tabfreq)) {tabfreq[,j][is.na(tabfreq[,j])] <- moy[j]}
- }
- }
-
- if(!quiet) cat("\n...done.\n\n")
-
- res$tab <- tabfreq
- res$nobs <- eff.pop
- res$call <- match.call()
-
- ## handle truenames
- if(truenames){
- temp <- rep(x at loc.names,x at loc.nall)
- colnames(res$tab) <- paste(temp,unlist(x at all.names),sep=".")
- rownames(res$tab) <- x at pop.names
-
- colnames(res$nobs) <- x at loc.names
- rownames(res$nobs) <- x at pop.names
- }
-
- return(res)
-} #end makefreq
-
Deleted: www/files/patches/propShared.R
===================================================================
--- www/files/patches/propShared.R 2010-05-04 12:37:28 UTC (rev 611)
+++ www/files/patches/propShared.R 2010-05-04 12:40:24 UTC (rev 612)
@@ -1,127 +0,0 @@
-## propShared computes the proportion of shared alleles
-## in a genind object
-
-
-######################
-# Function propShared
-######################
-propShared <- function(obj){
- x <- obj
-
- ## convert alleles to integers (alleles may be characters)
- x at all.names <- lapply(x at all.names, function(v) 1:length(v))
-
- ## check that this is a valid genind
- if(!inherits(x,"genind")) stop("obj must be a genind object.")
- invisible(validObject(x))
-
- ## check ploidy level
- if(x$ploidy > 2) stop("not implemented for ploidy > 2")
- checkType(x)
-
-
- ## if ploidy = 1
- if(x$ploidy == as.integer(1)){
- ## stop("not implemented for ploidy = 1")
- ## compute numbers of alleles used in each comparison
- nAllByInd <- propTyped(x,"both")
- nAll <- nAllByInd %*% t(nAllByInd)
-
- ## compute numbers of common alleles
- X <- x at tab
- X[is.na(X)] <- 0
- M <- X %*% t(X)
-
- ## result
- res <- M / nAll
- res[is.nan(res)] <- NA # as 0/0 is NaN (when no common locus typed)
- colnames(res) <- rownames(res) <- x$ind.names
- return(res)
- }
-
- ## if ploidy = 2
- if(x$ploidy == as.integer(2)){
- ## build a matrix of genotypes (in rows) coded by integers
- ## NAs are coded by 0
- ## The matrix is a cbind of two matrices, storing respectively the
- ## first and the second allele.
- temp <- genind2df(x,usepop=FALSE)
- alleleSize <- max(apply(temp,1:2,nchar))/2
- mat1 <- apply(temp, 1:2, substr, 1, alleleSize)
- mat2 <- apply(temp, 1:2, substr, alleleSize+1, alleleSize*2)
-
- matAll <- cbind(mat1,mat2)
- matAll <- apply(matAll,1:2,as.integer)
- matAll[is.na(matAll)] <- 0
-
- n <- nrow(matAll)
- resVec <- double(n*(n-1)/2)
- res <- .C("sharedAll", as.integer(as.matrix(matAll)),
- n, ncol(matAll), resVec, PACKAGE="adegenet")[[4]]
-
- attr(res,"Size") <- n
- attr(res,"Diag") <- FALSE
- attr(res,"Upper") <- FALSE
- class(res) <- "dist"
- res <- as.matrix(res)
- } # end if ploidy = 2
-
- diag(res) <- 1
- rownames(res) <- x at ind.names
- colnames(res) <- x at ind.names
- return(res)
-}
-
-
-
-
-## ######################
-## # Function propShared (old, pure-R version)
-## ######################
-## propShared <- function(obj){
-
-## x <- obj
-## ## check that this is a valid genind
-## if(!inherits(x,"genind")) stop("obj must be a genind object.")
-## invisible(validObject(x))
-
-## ## replace NAs
-## x <- na.replace(x, method="0")
-
-## ## some useful variables
-## nloc <- length(x at loc.names)
-
-## ## fnorm: auxiliary function for scaling
-## fnorm <- function(vec){
-## norm <- sqrt(sum(vec*vec))
-## if(length(norm) > 0 && norm > 0) {return(vec/norm)}
-## return(vec)
-## }
-
-## ## auxiliary function f1
-## ## computes the proportion of shared alleles in one locus
-## f1 <- function(X){
-## X <- t(apply(X, 1, fnorm))
-## res <- X %*% t(X)
-## res[res>0.51 & res<0.9] <- 0.5 # remap case one heteroZ shares the allele of on homoZ
-## return(res)
-## }
-
-## ## separate data per locus
-## temp <- seploc(x)
-## listProp <- lapply(temp, function(e) f1(e at tab))
-
-## ## produce the final result
-## res <- listProp[[1]]
-## if(nloc>2){
-## for(i in 2:nloc){
-## res <- res + listProp[[i]]
-## }
-## }
-
-## res <- res/nloc
-## rownames(res) <- x at ind.names
-## colnames(res) <- x at ind.names
-
-## return(res)
-## }
Deleted: www/files/patches/spca.R
===================================================================
--- www/files/patches/spca.R 2010-05-04 12:37:28 UTC (rev 611)
+++ www/files/patches/spca.R 2010-05-04 12:40:24 UTC (rev 612)
@@ -1,451 +0,0 @@
-##############################################
-#
-# spatial Principal Components Analysis
-#
-# require ade4, spdep and eventually tripack
-#
-# generic functions were derived from
-# those of multispati class (ade4)
-#
-# T. Jombart (t.jombart at imperial.ac.uk)
-# 31 may 2007
-##############################################
-
-
-################
-# spca genind
-################
-spca <- function(obj, xy=NULL, cn=NULL, matWeight=NULL,
- scale=FALSE, scale.method=c("sigma","binom"),
- scannf=TRUE, nfposi=1, nfnega=1,
- type=NULL, ask=TRUE, plot.nb=TRUE, edit.nb=FALSE,
- truenames=TRUE, d1=NULL, d2=NULL, k=NULL, a=NULL, dmin=NULL){
-
- ## first checks
- if(!any(inherits(obj,c("genind","genpop")))) stop("obj must be a genind or genpop object.")
- invisible(validObject(obj))
- ## checkType(obj)
- if(!require(ade4, quiet=TRUE)) stop("ade4 library is required.")
- if(!require(spdep, quiet=TRUE)) stop("spdep library is required.")
-
-
- ## handle xy coordinates
- if(is.null(xy) & (inherits(cn,"nb") & !inherits(cn,"listw")) ){
- xy <- attr(cn,"xy") # xy can be retrieved from a nb object (not from listw)
- }
- if(is.null(xy) & !is.null(obj$other$xy)) {
- xy <- obj$other$xy # xy from @other$xy if it exists
- }
- if(is.null(xy)) stop("xy coordinates are not provided")
- if(is.data.frame(xy)) {
- xy <- as.matrix(xy)
- }
- if(!is.matrix(xy)) stop("wrong 'xy' provided")
- if(ncol(xy) != 2) stop("xy does not have two columns.")
- if(nrow(xy) != nrow(obj at tab)) stop("obj at tab and xy must have the same row numbers.")
-
-
- appel <- match.call()
-
- ## == To find the spatial weights ==
- resCN <- NULL
-
- ## connection network from matWeight
- if(!is.null(matWeight)){
- if(!is.matrix(matWeight)) stop("matWeight is not a matrix")
- if(!is.numeric(matWeight)) stop("matWeight is not numeric")
- if(nrow(matWeight) != ncol(matWeight)) stop("matWeight is not square")
- if(nrow(matWeight) != nrow(obj at tab)) stop("dimension of datWeight does not match genetic data")
- diag(matWeight) <- 0
- matWeight <- prop.table(matWeight, 1)
- resCN <- mat2listw(matWeight)
- resCN$style <- "W"
-
- }
-
- ## connection network from cn argument
- if(is.null(resCN) & !is.null(cn)) {
- if(inherits(cn,"nb")) {
- if(!inherits(cn,"listw")){ # cn is a 'pure' nb object (i.e., nb but not listw)
- cn <- nb2listw(cn, style="W", zero.policy=TRUE)
- }
- resCN <- cn
- } else stop("cn does not have a recognized class")
- }
-
-
- ## connection network from xy coordinates
- if(is.null(resCN)) {
- resCN <- chooseCN(xy=xy, ask=ask, type=type, plot.nb=plot.nb, edit.nb=edit.nb,
- result.type="listw", d1=d1, d2=d2, k=k, a=a, dmin=dmin)
- }
-
- ## == spatial weights are done ==
-
-
- ## handle NAs warning
- if(any(is.na(obj at tab))){
- warning("NAs in data are automatically replaced (to mean allele frequency)")
- }
-
- ## handle NAs, centring and scaling
- X <- scaleGen(obj, center=TRUE, scale=scale, method=scale.method, missing="mean", truenames=truenames)
-
- ## perform analyses
- pcaX <- dudi.pca(X, center=FALSE, scale=FALSE, scannf=FALSE)
-
- spcaX <- multispati(dudi=pcaX, listw=resCN, scannf=scannf, nfposi=nfposi, nfnega=nfnega)
-
- nfposi <- spcaX$nfposi
- nfnega <- spcaX$nfnega
-
- spcaX$xy <- xy
- rownames(spcaX$xy) <- rownames(spcaX$li)
- colnames(spcaX$xy) <- c("x","y")
-
- spcaX$lw <- resCN
-
- spcaX$call <- appel
-
- posaxes <- if(nfposi>0) {1:nfposi} else NULL
- negaxes <- if(nfnega>0) {(length(spcaX$eig)-nfnega+1):length(spcaX$eig)} else NULL
- keptaxes <- c(posaxes,negaxes)
-
- ## set names of different components
- colnames(spcaX$c1) <- paste("Axis",keptaxes)
- colnames(spcaX$li) <- paste("Axis",keptaxes)
- colnames(spcaX$ls) <- paste("Axis",keptaxes)
- row.names(spcaX$c1) <- colnames(X)
- colnames(spcaX$as) <- colnames(spcaX$c1)
- temp <- row.names(spcaX$as)
- row.names(spcaX$as) <- paste("PCA",temp)
-
- class(spcaX) <- "spca"
-
- return(spcaX)
-
-} # end spca
-
-
-
-
-
-######################
-# Function print.spca
-######################
-print.spca <- function(x, ...){
- cat("\t########################################\n")
- cat("\t# spatial Principal Component Analysis #\n")
- cat("\t########################################\n")
- cat("class: ")
- cat(class(x))
- cat("\n$call: ")
- print(x$call)
- cat("\n$nfposi:", x$nfposi, "axis-components saved")
- cat("\n$nfnega:", x$nfnega, "axis-components saved")
-
- cat("\nPositive eigenvalues: ")
- l0 <- sum(x$eig >= 0)
- cat(signif(x$eig, 4)[1:(min(5, l0))])
- if (l0 > 5)
- cat(" ...\n")
- else cat("\n")
- cat("Negative eigenvalues: ")
- l0 <- sum(x$eig <= 0)
- cat(sort(signif(x$eig, 4))[1:(min(5, l0))])
- if (l0 > 5)
- cat(" ...\n")
- else cat("\n")
- cat('\n')
- sumry <- array("", c(1, 4), list(1, c("vector", "length",
- "mode", "content")))
- sumry[1, ] <- c('$eig', length(x$eig), mode(x$eig), 'eigenvalues')
- class(sumry) <- "table"
- print(sumry)
- cat("\n")
- sumry <- array("", c(4, 4), list(1:4, c("data.frame", "nrow", "ncol", "content")))
- sumry[1, ] <- c("$c1", nrow(x$c1), ncol(x$c1), "principal axes: scaled vectors of alleles loadings")
- sumry[2, ] <- c("$li", nrow(x$li), ncol(x$li), "principal components: coordinates of entities ('scores')")
- sumry[3, ] <- c("$ls", nrow(x$ls), ncol(x$ls), 'lag vector of principal components')
- sumry[4, ] <- c("$as", nrow(x$as), ncol(x$as), 'pca axes onto spca axes')
-
- class(sumry) <- "table"
- print(sumry)
-
- cat("\n$xy: matrix of spatial coordinates")
- cat("\n$lw: a list of spatial weights (class 'listw')")
-
- cat("\n\nother elements: ")
- if (length(names(x)) > 10)
- cat(names(x)[11:(length(names(x)))], "\n")
- else cat("NULL\n")
-}
-
-
-
-
-
-########################
-# Function summary.spca
-########################
-summary.spca <- function (object, ..., printres=TRUE) {
- if (!inherits(object, "spca"))stop("to be used with 'spca' object")
- if(!require(ade4,quietly=TRUE)) stop("the library ade4 is required; please install this package")
- if(!require(spdep,quietly=TRUE)) stop("the library spdep is required; please install this package")
-
- #util <- function(n) { ## no longer used
- # x <- "1"
- # for (i in 2:n) x[i] <- paste(x[i - 1], i, sep = "+")
- # return(x)
- #}
- norm.w <- function(X, w) {
- f2 <- function(v) sum(v * v * w)/sum(w)
- norm <- apply(X, 2, f2)
- return(norm)
- }
-
- resfin <- list()
-
- if(printres) {
- cat("\nSpatial principal component analysis\n")
- cat("\nCall: ")
- print(object$call)
- }
-
- appel <- as.list(object$call)
- ## compute original pca
- # prepare data
- obj <- eval(appel$obj)
- if(is.null(appel$truenames)) appel$truenames <- FALSE
-
- f1 <- function(vec){
- m <- mean(vec,na.rm=TRUE)
- vec[is.na(vec)] <- m
- return(vec)
- }
-
- if(is.genind(obj)) { X <- obj at tab }
- if(is.genpop(obj)) { X <- makefreq(obj, quiet=TRUE)$tab }
-
- X <- apply(X,2,f1)
-
- if(appel$truenames){
- rownames(X) <- rownames(truenames(obj))
- colnames(X) <- colnames(truenames(obj))
- }
-
- nfposi <- object$nfposi
- nfnega <- object$nfnega
-
- dudi <- dudi.pca(X, center=TRUE, scale=FALSE, scannf=FALSE, nf=nfposi+nfnega)
- ## end of pca
-
- lw <- object$lw
-
- # I0, Imin, Imax
- n <- nrow(X)
- I0 <- -1/(n-1)
- L <- listw2mat(lw)
- eigL <- eigen(0.5*(L+t(L)))$values
- Imin <- min(eigL)
- Imax <- max(eigL)
- Ival <- data.frame(I0=I0,Imin=Imin,Imax=Imax)
- row.names(Ival) <- ""
- if(printres) {
- cat("\nConnection network statistics:\n")
- print(Ival)
- }
-
- Istat <- c(I0,Imin,Imax)
- names(Istat) <- c("I0","Imin","Imax")
- resfin$Istat <- Istat
-
-
- # les scores de l'analyse de base
- nf <- dudi$nf
- eig <- dudi$eig[1:nf]
- cum <- cumsum(dudi$eig)[1:nf]
- ratio <- cum/sum(dudi$eig)
- w <- apply(dudi$l1,2,lag.listw,x=lw)
- moran <- apply(w*as.matrix(dudi$l1)*dudi$lw,2,sum)
- res <- data.frame(var=eig,cum=cum,ratio=ratio, moran=moran)
- row.names(res) <- paste("Axis",1:nf)
- if(printres) {
- cat("\nScores from the centred PCA\n")
- print(res)
- }
-
- resfin$pca <- res
-
-
- # les scores de l'analyse spatiale
- # on recalcule l'objet en gardant tous les axes
- eig <- object$eig
- nfposimax <- sum(eig > 0)
- nfnegamax <- sum(eig < 0)
-
- ms <- multispati(dudi=dudi, listw=lw, scannf=FALSE,
- nfposi=nfposimax, nfnega=nfnegamax)
-
- ndim <- dudi$rank
- nf <- nfposi + nfnega
- agarder <- c(1:nfposi,if (nfnega>0) (ndim-nfnega+1):ndim)
- varspa <- norm.w(ms$li,dudi$lw)
- moran <- apply(as.matrix(ms$li)*as.matrix(ms$ls)*dudi$lw,2,sum)
- res <- data.frame(eig=eig,var=varspa,moran=moran/varspa)
- row.names(res) <- paste("Axis",1:length(eig))
-
- if(printres) {
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/adegenet -r 612
More information about the adegenet-commits
mailing list