[adegenet-commits] r1044 - in pkg: R man src vignettes vignettes/Robjects
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Nov 2 18:02:00 CET 2012
Author: jombart
Date: 2012-11-02 18:02:00 +0100 (Fri, 02 Nov 2012)
New Revision: 1044
Added:
pkg/vignettes/Robjects/x.dapc1.RData
pkg/vignettes/Robjects/x.genlight.RData
pkg/vignettes/Robjects/x.snpbin.RData
Modified:
pkg/R/glHandle.R
pkg/R/glSim.R
pkg/R/import.R
pkg/man/dapcGraphics.Rd
pkg/man/fasta2genlight.Rd
pkg/man/genlight.Rd
pkg/man/glSim.Rd
pkg/man/seppop.Rd
pkg/src/snpbin.c
pkg/vignettes/adegenet-genomics.Rnw
pkg/vignettes/loadingplotGL.png
Log:
Changes for the fixed release of 1.3-5
Modified: pkg/R/glHandle.R
===================================================================
--- pkg/R/glHandle.R 2012-10-31 15:55:45 UTC (rev 1043)
+++ pkg/R/glHandle.R 2012-11-02 17:02:00 UTC (rev 1044)
@@ -81,7 +81,7 @@
##x <- as.matrix(x)[, j, drop=FALSE] # maybe need to process one row at a time
x <- new("genlight", gen=new.gen, pop=ori.pop, ploidy=ori.ploidy,
ind.names=old.ind.names, loc.names=new.loc.names,
- chromosome=new.chr, position=new.position, alleles=new.alleles, other=old.other)
+ chromosome=new.chr, position=new.position, alleles=new.alleles, other=old.other, multicore=FALSE,...)
}
return(x)
@@ -168,7 +168,7 @@
res[[i]] <- Reduce(function(a,b) {cbind(a,b,checkPloidy=FALSE)}, lapply(myList, function(e) e at gen[[i]]) )
}
- res <- new("genlight",res)
+ res <- new("genlight",res,...)
## handle loc.names, alleles, etc. ##
indNames(res) <- indNames(myList[[1]])
@@ -205,7 +205,7 @@
## build output
- res <- new("genlight", Reduce(c, lapply(myList, function(e) e at gen)))
+ res <- new("genlight", Reduce(c, lapply(myList, function(e) e at gen)), ...)
locNames(res) <- locNames(myList[[1]])
alleles(res) <- alleles(myList[[1]])
indNames(res) <- unlist(lapply(myList, indNames))
@@ -224,7 +224,7 @@
##########
## seppop
##########
-setMethod("seppop", signature(x="genlight"), function(x, pop=NULL, treatOther=TRUE, quiet=TRUE){
+setMethod("seppop", signature(x="genlight"), function(x, pop=NULL, treatOther=TRUE, quiet=TRUE, ...){
## HANDLE POP ARGUMENT ##
if(!is.null(pop)) {
pop(x) <- pop
@@ -233,7 +233,7 @@
if(is.null(pop(x))) stop("pop not provided and pop(x) is NULL")
## PERFORM SUBSETTING ##
- kObj <- lapply(levels(pop(x)), function(lev) x[pop(x)==lev, , treatOther=treatOther, quiet=quiet])
+ kObj <- lapply(levels(pop(x)), function(lev) x[pop(x)==lev, , treatOther=treatOther, quiet=quiet, ...])
names(kObj) <- levels(pop(x))
return(kObj)
Modified: pkg/R/glSim.R
===================================================================
--- pkg/R/glSim.R 2012-10-31 15:55:45 UTC (rev 1043)
+++ pkg/R/glSim.R 2012-11-02 17:02:00 UTC (rev 1044)
@@ -2,8 +2,7 @@
##########
## glSim
##########
-glSim <- function(n.ind, n.snp.nonstruc, n.snp.struc=0, grp.size=round(n.ind/2), ploidy=1, alpha=0,
- block.size=NULL, LD=FALSE){
+glSim <- function(n.ind, n.snp.nonstruc, n.snp.struc=0, grp.size=round(n.ind/2), ploidy=1, alpha=0, block.size=NULL, LD=FALSE, ...){
## BASIC CHECKS ##
if( any(c(n.ind, n.snp.nonstruc+n.snp.struc) <1)) stop("null numbers of individuals and/or SNPs requested")
@@ -38,7 +37,7 @@
Sig <- t(Sig) %*% Sig/2 # get the covariance matrix
temp <- mvrnorm(n, rep(0,p), Sig) # continuous data
temp <- matrix(as.integer(cut(temp, breaks=QUANT))-1, nrow=n, ncol=p)
- return(new("genlight", temp, ploidy=ploidy))
+ return(new("genlight", temp, ploidy=ploidy, ...))
}
if(n.snp.struc > 0){
@@ -54,7 +53,7 @@
f1 <- function(n,p){
temp <- sapply(1:p, function(i) rbinom(n, ploidy, runif(1)))
if(n==1) {temp <- matrix(temp,nrow=1)}
- return(new("genlight", temp, ploidy=ploidy))
+ return(new("genlight", temp, ploidy=ploidy, ...))
}
## draw p snp for i indiv and convert into a genlight - differences between 2 groups
@@ -66,7 +65,7 @@
if(grpA.size==1) {tempA <- matrix(tempA,nrow=1)}
tempB <- sapply(probB, function(i) rbinom(grpB.size, ploidy, i) )
if(grpB.size==1) {tempB <- matrix(tempB,nrow=1)}
- return(new("genlight", rbind(tempA,tempB), ploidy=ploidy))
+ return(new("genlight", rbind(tempA,tempB), ploidy=ploidy, ...))
}
}
}
Modified: pkg/R/import.R
===================================================================
--- pkg/R/import.R 2012-10-31 15:55:45 UTC (rev 1043)
+++ pkg/R/import.R 2012-11-02 17:02:00 UTC (rev 1044)
@@ -861,7 +861,7 @@
other <- list(chromosome = misc.info$chromosome)
}
- res <- new("genlight", gen=res, ind.names=ind.names, position=misc.info$position, loc.all=misc.info$allele, ploidy=misc.info$ploidy, pop=misc.info$population, other=other)
+ res <- new("genlight", gen=res, ind.names=ind.names, position=misc.info$position, loc.all=misc.info$allele, ploidy=misc.info$ploidy, pop=misc.info$population, other=other, multicore=multicore)
if(!quiet) cat("\n...done.\n\n")
@@ -1000,7 +1000,7 @@
## BUILD FINAL OBJECT ##
if(!quiet) cat("\n Building final object... \n")
- res <- new("genlight",res, ploidy=2)
+ res <- new("genlight",res, ploidy=2, multicore=multicore)
indNames(res) <- misc.info$IID
pop(res) <- misc.info$FID
locNames(res) <- loc.names
@@ -1163,7 +1163,7 @@
## BUILD FINAL OBJECT ##
if(!quiet) cat("\n Building final object... \n")
- res <- new("genlight",res, ploidy=1)
+ res <- new("genlight",res, ploidy=1, multicore=multicore)
indNames(res) <- IND.LAB
alleles(res) <- sapply(POOL[snp.posi], paste, collapse="/")
position(res) <- which(snp.posi)
Modified: pkg/man/dapcGraphics.Rd
===================================================================
--- pkg/man/dapcGraphics.Rd 2012-10-31 15:55:45 UTC (rev 1043)
+++ pkg/man/dapcGraphics.Rd 2012-11-02 17:02:00 UTC (rev 1044)
@@ -140,6 +140,7 @@
}
\author{ Thibaut Jombart \email{t.jombart at imperial.ac.uk} }
\examples{
+\dontrun{
data(H3N2)
dapc1 <- dapc(H3N2, pop=H3N2$other$epid, n.pca=30,n.da=6)
@@ -187,6 +188,6 @@
data(microbov)
dapc3 <- dapc(microbov, n.pca=20, n.da=15)
compoplot(dapc3, lab="")
-
}
+}
\keyword{multivariate}
\ No newline at end of file
Modified: pkg/man/fasta2genlight.Rd
===================================================================
--- pkg/man/fasta2genlight.Rd 2012-10-31 15:55:45 UTC (rev 1043)
+++ pkg/man/fasta2genlight.Rd 2012-11-02 17:02:00 UTC (rev 1044)
@@ -75,6 +75,7 @@
}
\author{Thibaut Jombart \email{t.jombart at imperial.ac.uk} }
\examples{
+\dontrun{
## show the example file ##
## this is the path to the file:
myPath <- system.file("files/usflu.fasta",package="adegenet")
@@ -93,6 +94,6 @@
temp <- density(position(obj), bw=10)
plot(temp, xlab="Position in the alignment", lwd=2, main="Location of the SNPs")
points(position(obj), rep(0, nLoc(obj)), pch="|", col="red")
-
}
+}
\keyword{manip}
Modified: pkg/man/genlight.Rd
===================================================================
--- pkg/man/genlight.Rd 2012-10-31 15:55:45 UTC (rev 1043)
+++ pkg/man/genlight.Rd 2012-11-02 17:02:00 UTC (rev 1044)
@@ -243,7 +243,9 @@
warning unless the argument \code{quiet} is set to TRUE.}
\item{quiet}{a logical indicating whether warnings should be issued
when trying to subset components of the \code{@other} slot which
- do not match the number of individuals (TRUE), or not (FALSE, default). }
+ do not match the number of individuals (TRUE), or not (FALSE,
+ default). }
+ \item{\dots}{further arguments passed to the genlight constructor.}
}
}
\author{Thibaut Jombart (\email{t.jombart at imperial.ac.uk})}
Modified: pkg/man/glSim.Rd
===================================================================
--- pkg/man/glSim.Rd 2012-10-31 15:55:45 UTC (rev 1043)
+++ pkg/man/glSim.Rd 2012-11-02 17:02:00 UTC (rev 1044)
@@ -9,7 +9,7 @@
}
\usage{
glSim(n.ind, n.snp.nonstruc, n.snp.struc = 0, grp.size = round(n.ind/2),
- ploidy = 1, alpha = 0, block.size = NULL, LD = FALSE)
+ ploidy = 1, alpha = 0, block.size = NULL, LD = FALSE, \dots)
}
\arguments{
\item{n.ind}{an integer indicating the number of individuals to be simulated.}
@@ -41,6 +41,7 @@
linkage disequilibrium (TRUE) or be generated independently
(FALSE, default). When set to TRUE, data are generated by blocks of
correlated SNPs (see details).}
+ \item{\dots}{arguments to be passed to the genlight constructor.}
}
\details{
=== Allele frequencies in contrasted groups ===
@@ -89,6 +90,7 @@
}
\author{ Thibaut Jombart \email{t.jombart at imperial.ac.uk} }
\examples{
+\dontrun{
## no structure
x <- glSim(100, 1e3, ploid=2)
plot(x)
@@ -97,7 +99,6 @@
x <- glSim(100, 1e3, n.snp.struc=100, ploid=2)
plot(x)
-\dontrun{
## 1,000 non structured SNPs, 100 structured SNPs, ploidy=4
x <- glSim(100, 1e3, n.snp.struc=100, ploid=4)
plot(x)
Modified: pkg/man/seppop.Rd
===================================================================
--- pkg/man/seppop.Rd 2012-10-31 15:55:45 UTC (rev 1043)
+++ pkg/man/seppop.Rd 2012-11-02 17:02:00 UTC (rev 1044)
@@ -20,7 +20,7 @@
\S4method{seppop}{genind}(x,pop=NULL,truenames=TRUE,res.type=c("genind","matrix"),
drop=FALSE, treatOther=TRUE, quiet=TRUE)
-\S4method{seppop}{genlight}(x,pop=NULL, treatOther=TRUE, quiet=TRUE)
+\S4method{seppop}{genlight}(x,pop=NULL, treatOther=TRUE, quiet=TRUE, \dots)
}
\arguments{
\item{x}{a \linkS4class{genind} object}
@@ -42,6 +42,7 @@
\item{quiet}{a logical indicating whether warnings should be issued
when trying to subset components of the \code{@other} slot (TRUE),
or not (FALSE, default). }
+ \item{\dots}{further arguments passed to the genlight constructor.}
}
\value{
According to 'res.type': a list of \linkS4class{genind} object
Modified: pkg/src/snpbin.c
===================================================================
--- pkg/src/snpbin.c 2012-10-31 15:55:45 UTC (rev 1043)
+++ pkg/src/snpbin.c 2012-11-02 17:02:00 UTC (rev 1044)
@@ -331,7 +331,7 @@
*temp=nLoc(x);
bytesToDouble(x->bytevec, x->byteveclength, x->bytevecnb, out, temp);
int i;
-
+
for(i=0; i < nLoc(x); i++){
out[i] = out[i] / ploid;
}
@@ -354,13 +354,13 @@
Rprintf("%i ", (int) (x->bytevec)[i]);
/* printf("%i ", (int) (x->bytevec)[i]); */
}
- printf(" ");
+ Rprintf(" ");
for(i=0;i<nLoc(x);i++){
Rprintf("%i ", temp[i]);
/* printf("%i ", temp[i]); */
}
- printf("NA posi: ");
+ Rprintf("NA posi: ");
for(i=0;i< *(x->nanb);i++){
Rprintf("%i ", (x->naposi)[i]);
/* printf("%i ", (x->naposi)[i]); */
Added: pkg/vignettes/Robjects/x.dapc1.RData
===================================================================
(Binary files differ)
Property changes on: pkg/vignettes/Robjects/x.dapc1.RData
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: pkg/vignettes/Robjects/x.genlight.RData
===================================================================
(Binary files differ)
Property changes on: pkg/vignettes/Robjects/x.genlight.RData
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: pkg/vignettes/Robjects/x.snpbin.RData
===================================================================
(Binary files differ)
Property changes on: pkg/vignettes/Robjects/x.snpbin.RData
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Modified: pkg/vignettes/adegenet-genomics.Rnw
===================================================================
--- pkg/vignettes/adegenet-genomics.Rnw 2012-10-31 15:55:45 UTC (rev 1043)
+++ pkg/vignettes/adegenet-genomics.Rnw 2012-11-02 17:02:00 UTC (rev 1044)
@@ -118,17 +118,24 @@
New objects are created using \texttt{new}, with these slots as arguments.
If no argument is provided, an empty object is created:
-<<>>=
+<<eval=FALSE>>=
new("SNPbin")
@
+<<echo=FALSE>>=
+new("SNPbin", multicore=FALSE)
+@
In practice, only the \texttt{snp} information and possibly the ploidy has to be provided; various
formats are accepted for the \texttt{snp} component, but the simplest is a vector of integers (or
numeric) indicating the number of second allele at each locus.
The argument \texttt{snp}, if provided alone, does not have to be named:
-<<>>=
+<<eval=FALSE>>=
x <- new("SNPbin", c(0,1,1,2,0,0,1))
x
@
+<<echo=FALSE>>=
+x <- new("SNPbin", c(0,1,1,2,0,0,1), multicore=FALSE)
+x
+@
If not provided, the ploidy is detected from the data and determined as the largest number in the
input vector. Obviously, in many cases this will not be adequate, but ploidy can always be rectified
@@ -155,7 +162,14 @@
<<>>=
dat <- sample(0:1, 1e6, replace=TRUE)
print(object.size(dat),unit="auto")
-x <- new("SNPbin", dat)
+@
+<<eval=FALSE>>=
+x <- new("SNPbin", dat, multicore=FALSE)
+@
+<<echo=FALSE>>=
+load("Robjects/x.snpbin.RData")
+@
+<<>>=
print(object.size(x),unit="auto")
@
here, we converted a million SNPs into a \texttt{SNPbin} object, which turns out to be
@@ -210,9 +224,12 @@
\noindent Like \texttt{SNPbin} object, \texttt{genlight} object are created using the constructor \texttt{new},
providing content for the slots above as arguments.
When none is provided, an empty object is created:
-<<>>=
+<<eval=FALSE>>=
new("genlight")
@
+<<echo=FALSE>>=
+new("genlight", multicore=FALSE)
+@
The most important information to provide is obviously the genotypes (argument \texttt{gen}); these
can be provided as:
\begin{itemize}
@@ -224,8 +241,13 @@
Ploidy has to be consistent across loci for a given individual, but individuals do not have to have
the same ploidy, so that it is possible to have hapoid,
diploid, and tetraploid individuals in the same dataset; for instance:
+<<eval=FALSE>>=
+x <- new("genlight", list(indiv1=c(1,1,0,1,1,0), indiv2=c(2,1,1,0,0,0), toto=c(2,2,0,0,4,4)))
+@
+<<echo=FALSE>>=
+x <- new("genlight", list(indiv1=c(1,1,0,1,1,0), indiv2=c(2,1,1,0,0,0), toto=c(2,2,0,0,4,4)), multicore=FALSE)
+@
<<>>=
-x <- new("genlight", list(indiv1=c(1,1,0,1,1,0), indiv2=c(2,1,1,0,0,0), toto=c(2,2,0,0,4,4)))
x
ploidy(x)
@
@@ -240,15 +262,20 @@
\noindent In practice, \texttt{genlight} objects can be handled as if they were matrices of integers
as the one above returned by \texttt{as.matrix}.
However, they offer the advantage of efficient storage of the information; for instance, we can
-simulate 50 individuals typed for 1,00,000 SNPs each (including occasional NAs):
+simulate 50 individuals typed for 100,000 SNPs each (including occasional NAs):
<<>>=
-dat <- lapply(1:50, function(i) sample(c(0,1,NA), 1e6, prob=c(.5, .499, .001), replace=TRUE))
+dat <- lapply(1:50, function(i) sample(c(0,1,NA), 1e5, prob=c(.5, .499, .001), replace=TRUE))
names(dat) <- paste("indiv", 1:length(dat))
print(object.size(dat),unit="auto")
-
+@
+<<eval=FALSE>>=
x <- new("genlight", dat)
+@
+<<echo=FALSE>>=
+load("Robjects/x.genlight.RData")
+@
+<<>>=
print(object.size(x),unit="auto")
-
object.size(dat)/object.size(x)
@
here again, the storage if the data is much more efficient in \texttt{genlight} than using integers: converted data occupy
@@ -316,7 +343,7 @@
<<>>=
dat <- lapply(1:3, function(i) sample(0:2, 10, replace=TRUE))
dat
-x <- new("genlight", dat)
+x <- new("genlight", dat, multicore=FALSE)
x
indNames(x)
indNames(x) <- paste("individual", 1:3)
@@ -364,7 +391,7 @@
Therefore, subsetting can be achieved using $[$ \texttt{idx.row , idx.col} $]$ where \texttt{idx.row}
and \texttt{idx.col} are indices for rows (individuals) and columns (SNPs).
For instance, using the previous toy dataset, we try a few classical subsetting of rows and columns:
-<<>>=
+<<asmatrix>>=
x
as.matrix(x)
as.matrix(x[c(1,3),])
@@ -377,17 +404,18 @@
This can be achieved by specifying either a number of blocks (argument \texttt{n.block}) or the size
of the blocks (argument \texttt{block.size}). The function also allows for randomizing the
distribution of the SNPs in the blocks (argument \texttt{random=TRUE}), which is especially useful
-to replace computations that cannot be achieved on the whole dataset with parallelized computations performed on random blocks.
+to replace computations that cannot be achieved on the whole dataset with parallelized computations
+performed on random blocks (for parallelization, remove the argument \texttt{multicore=FALSE}).
For instance:
-<<>>=
+<<seploc>>=
x
as.matrix(x)
-seploc(x, n.block=2)
-lapply(seploc(x, n.block=2),as.matrix)
+seploc(x, n.block=2, multicore=FALSE)
+lapply(seploc(x, n.block=2, multicore=FALSE),as.matrix)
@
splits the data into two blocks of contiguous SNPs, while:
<<>>=
-lapply(seploc(x, n.block=2, random=TRUE),as.matrix)
+lapply(seploc(x, n.block=2, random=TRUE, multicore=FALSE),as.matrix)
@
generates blocks of randomly selected SNPs.
@@ -490,10 +518,10 @@
1100
\end{verbatim}
We read the file in using:
-<<>>=
-obj <- read.snp(system.file("files/exampleSnpDat.snp",package="adegenet"), chunk=2)
+<<readsnp>>=
+obj <- read.snp(system.file("files/exampleSnpDat.snp",package="adegenet"), chunk=2, multicore=FALSE)
obj
-as.matrix(obj)
+as.matrix(obj, multicore=FALSE)
alleles(obj)
pop(obj)
indNames(obj)
@@ -554,7 +582,7 @@
The file is first located using \texttt{system.file}, and then processed using \texttt{fasta2genlight}:
<<>>=
myPath <- system.file("files/usflu.fasta",package="adegenet")
-flu <- fasta2genlight(myPath, chunk=10)
+flu <- fasta2genlight(myPath, chunk=10, multicore=FALSE)
flu
@
@@ -571,9 +599,10 @@
\noindent It is usually informative to assess the position of the polymorphic sites within the
genome; this is very easily done in R, using \texttt{density} with an appropriate bandewidth:
-<<fig=TRUE>>=
+<<fig=TRUE,keep.source=TRUE>>=
temp <- density(position(flu), bw=10)
-plot(temp, type="n", xlab="Position in the alignment", main="Location of the SNPs", xlim=c(0,1701))
+plot(temp, type="n", xlab="Position in the alignment", main="Location of the SNPs",
+ xlim=c(0,1701))
polygon(c(temp$x,rev(temp$x)), c(temp$y, rep(0,length(temp$x))), col=transp("blue",.3))
points(position(flu), rep(0, nLoc(flu)), pch="|", col="blue")
@
@@ -588,7 +617,7 @@
It is however possible to ask \texttt{fasta2genlight} to keep track of the number of alleles for
each site of the original alignment, by specifying:
<<>>=
-flu <- fasta2genlight(myPath, chunk=10,saveNbAlleles=TRUE, quiet=TRUE)
+flu <- fasta2genlight(myPath, chunk=10,saveNbAlleles=TRUE, quiet=TRUE, multicore=FALSE)
flu
@
@@ -716,9 +745,10 @@
possible centring and scaling.
\end{itemize}
For instance, one can easily derive the distributiong of allele frequencies using:
-<<fig=TRUE>>=
+<<fig=TRUE,keep.source=TRUE>>=
myFreq <- glMean(flu)
-hist(myFreq, proba=TRUE, col="gold", xlab="Allele frequencies", main="Distribution of (second) allele frequencies")
+hist(myFreq, proba=TRUE, col="gold", xlab="Allele frequencies",
+ main="Distribution of (second) allele frequencies")
temp <- density(myFreq)
lines(temp$x, temp$y*1.8,lwd=3)
@
@@ -726,10 +756,11 @@
\noindent In biallelic loci, one allele is always entirely redundant with the other, so it is generally
sufficient to analyse a single allele per loci.
However, the distribution of allele frequencies may be more interpretable by restoring its native symmetry:
-<<fig=TRUE>>=
+<<fig=TRUE,keep.source=TRUE>>=
myFreq <- glMean(flu)
myFreq <- c(myFreq, 1-myFreq)
-hist(myFreq, proba=TRUE, col="darkseagreen3", xlab="Allele frequencies", main="Distribution of allele frequencies", nclass=20)
+hist(myFreq, proba=TRUE, col="darkseagreen3", xlab="Allele frequencies",
+ main="Distribution of allele frequencies", nclass=20)
temp <- density(myFreq, bw=.05)
lines(temp$x, temp$y*2,lwd=3)
@
@@ -751,10 +782,11 @@
Besides these considerations, we can also map missing data across loci as we have done for SNP positions in the
US influenza dataset (see previous section) using \texttt{glNA} and \texttt{density}:
-<<>>=
+<<keep.source=TRUE>>=
head(glNA(flu),20)
temp <- density(glNA(flu), bw=10)
-plot(temp, type="n", xlab="Position in the alignment", main="Location of the missing values (NAs)", xlim=c(0,1701))
+plot(temp, type="n", xlab="Position in the alignment", main="Location of the missing values (NAs)",
+ xlim=c(0,1701))
polygon(c(temp$x,rev(temp$x)), c(temp$y, rep(0,length(temp$x))), col=transp("blue",.3))
points(glNA(flu), rep(0, nLoc(flu)), pch="|", col="blue")
@
@@ -783,14 +815,14 @@
dataset, and results are finally reunited.
\\
-Let us illustrate this procedure using 40 simulated individuals with 100,000 SNPs each:
+Let us illustrate this procedure using 40 simulated individuals with 10,000 SNPs each:
<<>>=
-x <- glSim(40, 1e5, LD=FALSE)
+x <- glSim(40, 1e4, LD=FALSE, multicore=FALSE)
x
@
\texttt{seploc} is used to create a list of smaller objects (here, 10 blocks of 10,000 SNPs):
<<>>=
-x <- seploc(x, n.block=10)
+x <- seploc(x, n.block=10, multicore=FALSE)
class(x)
names(x)
x[1:2]
@@ -835,7 +867,7 @@
As a matter of fact, it does not matter when all individuals have the same ploidy.
For instance, if we take the following data:
<<>>=
-x <- new("genlight", list(a=c(0,0,1,1), b=c(1,1,0,0), c=c(1,1,1,1)))
+x <- new("genlight", list(a=c(0,0,1,1), b=c(1,1,0,0), c=c(1,1,1,1)), multicore=FALSE)
locNames(x) <- 1:4
x
as.matrix(x)
@@ -847,7 +879,7 @@
@
Let us no consider a sightly different case:
<<>>=
-x <- new("genlight", list(a=c(0,0,2,2), b=c(1,1,0,0), c=c(1,1,1,1)))
+x <- new("genlight", list(a=c(0,0,2,2), b=c(1,1,0,0), c=c(1,1,1,1)), multicore=FALSE)
locNames(x) <- 1:4
x
as.matrix(x)
@@ -898,7 +930,6 @@
It also uses compiled C code and possibly multicore ressources if available to speed up computations.
We illustrate the method on the previously introduced influenza dataset (object \texttt{flu}):
<<eval=FALSE>>=
-flu
pca1 <- glPca(flu)
@
<<echo=FALSE>>=
@@ -923,7 +954,8 @@
scatter(pca1, posi="bottomright")
title("PCA of the US influenza data\n axes 1-2")
@
-The first PC suggests the existence of two clades in the data, while the second one shows groups of
+
+\noindent The first PC suggests the existence of two clades in the data, while the second one shows groups of
closely related isolates arranged along a cline of genetic differentiation.
This structure is confirmed by a simple neighbour-joining (NJ) tree:
<<>>=
@@ -966,19 +998,23 @@
allele frequencies at a time, thereby minimizing RAM requirements.
\\
-Below, we illustrate DAPC on a \texttt{genlight} including only 50 structured SNPs out of 50,000
+Below, we illustrate DAPC on a \texttt{genlight} with 100 individuals, including only 50 structured SNPs out of 10,000
non-structured SNPs:
-<<>>=
-x <- glSim(100, 5e4, 50)
-x
-dapc1 <- dapc(x, n.pca=20, n.da=1)
+<<eval=FALSE>>=
+x <- glSim(100, 1e4, 50)
+dapc1 <- dapc(x, n.pca=10, n.da=1)
@
-For the last 50 structured SNPs (located at the end of the alignment), the two groups of individuals have different (random) distribution of
+<<echo=FALSE>>=
+load("Robjects/x.dapc1.RData")
+@
+
+For the last 10 structured SNPs (located at the end of the alignment), the two groups of individuals have different (random) distribution of
allele frequencies, while they share the same distributions in other loci.
DAPC can still make some decent discrimination:
-<<fig=TRUE>>=
-scatter(dapc1,scree.da=FALSE, bg="white", posi.pca="topright", legen=TRUE, txt.leg=paste("group", 1:2), col=c("red","blue"))
+<<fig=TRUE,keep.source=TRUE>>=
+scatter(dapc1,scree.da=FALSE, bg="white", posi.pca="topright", legend=TRUE,
+ txt.leg=paste("group", 1:2), col=c("red","blue"))
@
\noindent While the composition plot confirms that groups are not entirely disentangled...
@@ -988,7 +1024,7 @@
\noindent ... the loading plot identifies pretty well the most discriminating alleles:
<<eval=FALSE>>=
-loadingplot(dapc1$var.contr, thres=4e-4)
+loadingplot(dapc1$var.contr, thres=1e-3)
@
\begin{center}
@@ -997,9 +1033,12 @@
\noindent And we can zoom in to the contributions of the last 100 SNPs to make sure that the tail indeed
corresponds to the 50 last structured loci:
-<<fig=TRUE>>=
-loadingplot(tail(dapc1$var.contr[,1],100))
+<<eval=FALSE>>=
+loadingplot(tail(dapc1$var.contr[,1],100), thres=1e-3)
@
+\begin{center}
+ \includegraphics{loadingplotGLzoom}
+\end{center}
\noindent Here, we indeed identified the structured region of the genome fairly well.
Modified: pkg/vignettes/loadingplotGL.png
===================================================================
(Binary files differ)
More information about the adegenet-commits
mailing list