[adegenet-commits] r720 - in pkg: . R

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

Author: jombart
Date: 2010-11-23 21:45:05 +0100 (Tue, 23 Nov 2010)
New Revision: 720

Added:
   pkg/R/inbreeding.R
Modified:
   pkg/DESCRIPTION
Log:
inbreeding stuff


Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION	2010-11-23 15:53:03 UTC (rev 719)
+++ pkg/DESCRIPTION	2010-11-23 20:45:05 UTC (rev 720)
@@ -10,6 +10,6 @@
 Suggests: ade4, genetics, hierfstat, spdep, tripack, ape, pegas, graph, RBGL, seqinr
 Depends: methods, MASS
 Description: Classes and functions for genetic data analysis within the multivariate framework.
-Collate: classes.R basicMethods.R handling.R auxil.R setAs.R find.clust.R hybridize.R scale.R fstat.R import.R seqTrack.R chooseCN.R genind2genpop.R loadingplot.R sequences.R gstat.randtest.R makefreq.R  colorplot.R  monmonier.R spca.R coords.monmonier.R haploGen.R old2new.R spca.rtests.R dapc.R haploPop.R PCtest.R dist.genpop.R Hs.R propShared.R export.R HWE.R propTyped.R  zzz.R
+Collate: classes.R basicMethods.R handling.R auxil.R setAs.R find.clust.R hybridize.R scale.R fstat.R import.R seqTrack.R chooseCN.R genind2genpop.R loadingplot.R sequences.R gstat.randtest.R makefreq.R  colorplot.R  monmonier.R spca.R coords.monmonier.R haploGen.R old2new.R spca.rtests.R dapc.R haploPop.R PCtest.R dist.genpop.R Hs.R propShared.R export.R HWE.R propTyped.R  inbreeding.R zzz.R
 License: GPL (>=2)
 LazyLoad: yes

Added: pkg/R/inbreeding.R
===================================================================
--- pkg/R/inbreeding.R	                        (rev 0)
+++ pkg/R/inbreeding.R	2010-11-23 20:45:05 UTC (rev 720)
@@ -0,0 +1,69 @@
+#############
+## inbreeding
+#############
+inbreeding <- function(x, pop=NULL, truenames=TRUE, res.type=c("mean","byloc"), plot=TRUE, ...){
+    ## CHECKS ##
+    if(!is.genind(x)) stop("x is not a valid genind object")
+    checkType(x)
+    res.type <- match.arg(res.type)
+
+    if(x$ploidy != 2) stop("this inbreeding coefficient is designed for diploid genotypes only")
+
+    if(!is.null(pop)) pop(x) <- 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)))
+    }
+
+
+    ## COMPUTATIONS ##
+
+    ## get allele frequencies and \sum p_i^2 by pop and loc ##
+    tabfreq2 <- (makefreq(x = genind2genpop(x, quiet = TRUE), quiet=TRUE, truenames=truenames)$tab) ^2
+    sumpi2 <- t(apply(tabfreq2, 1, tapply, x$loc.fac, sum))
+
+    ## function to check a 1-locus genotype for heterozigocity
+    ## returns 1 if heteroz, 0 otherwise
+    f1 <- function(gen){
+        if(any(is.na(gen))) return(NA)
+        if(sum(abs(gen-0.5) < 1e-10)==2) return(1)
+        return(0)
+    }
+
+    ## get the table of binary hetero/homo data
+    if(truenames) {
+        X <- truenames(x)$tab
+    } else
+    X <- x$tab
+
+    heterotab <- t(apply(X, 1, tapply, x at loc.fac, f1))
+
+
+    ## get pi2 for the appropriate pop
+    if(truenames){
+    popx <- pop(x)
+    } else {
+        popx <- x$pop
+    }
+
+    popx <- as.character(popx)
+    tabpi2 <- sumpi2[popx, , drop=FALSE]
+
+
+    ## COMPUTE FINAL RESULT ##
+    num <- heterotab - tabpi2
+    denom <- tabpi2 * (1 - tabpi2)
+    res <- num / denom
+    if(res.type=="byloc") return(res)
+
+    res <- apply(res, 1, mean, na.rm=TRUE)
+    if(plot){
+        par(bg="grey")
+        nPop <- length(unique(popx))
+        myCol <- rainbow(nPop)[as.integer(pop(x))]
+        plot(res, col=myCol, type="h", ylab="Inbreeding", xlab="Individuals", ...)
+    }
+
+    return(res)
+
+} # end inbreeding