[adegenet-commits] r633 - in pkg: R man

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

Author: jombart
Date: 2010-05-18 13:18:18 +0200 (Tue, 18 May 2010)
New Revision: 633

Modified:
   pkg/R/seqTrack.R
   pkg/man/seqTrack.Rd
Log:
H1N1 example


Modified: pkg/R/seqTrack.R
===================================================================
--- pkg/R/seqTrack.R	2010-05-18 10:05:19 UTC (rev 632)
+++ pkg/R/seqTrack.R	2010-05-18 11:18:18 UTC (rev 633)
@@ -307,10 +307,9 @@
     anc.dates <- as.POSIXct(x$ances.date)
     nb.days <- abs(as.integer(anc.dates-dates))
     nb.mut <- x$weight
-    ##mu <- mu/365
-    ##mu <- mu*nb.days
 
-    res <- dbinom(nb.mut, size=seq.length*nb.days, prob=mu)
+    ##    res <- dbinom(nb.mut, size=haplo.length*nb.days, prob=mu)
+    res <- dpois(x=nb.mut, lambda=mu*haplo.length*nb.days)
 
     return(res)
 } # end get.likelihood.seqTrack

Modified: pkg/man/seqTrack.Rd
===================================================================
--- pkg/man/seqTrack.Rd	2010-05-18 10:05:19 UTC (rev 632)
+++ pkg/man/seqTrack.Rd	2010-05-18 11:18:18 UTC (rev 633)
@@ -145,8 +145,119 @@
 }
 \examples{
 if(require(ape)){
+## ANALYSIS OF SIMULATED DATA ##
+## SIMULATE A GENEALOGY
+dat <- haploGen(seq.l=1e4, repro=function(){sample(1:4,1)}, gen.time=1, t.max=3)
+
+
+## SEQTRACK ANALYSIS
+res <- seqTrack(dat, mu=0.0001, haplo.length=1e4) 
+
+
+## PROPORTION OF CORRECT RECONSTRUCTION
+mean(dat$ances==res$ances,na.rm=TRUE)
+
+
+## PLOT RESULTS
+if(require(graph) && require(Rgraphviz)){
+dat.g <- as(dat, "graphNEL")
+res.g <- as(res, "graphNEL")
+
+
+## ORIGINAL DATA
+dat.annot <- as.character(unlist(edgeWeights(dat.g)))
+names(dat.annot) <- edgeNames(dat.g)
+renderGraph(layoutGraph(dat.g, edgeAttrs = list(label = dat.annot)))
+
+
+## SEQTRACK RESULTS
+res.annot <- as.character(unlist(edgeWeights(res.g)))
+names(res.annot) <- edgeNames(res.g)
+renderGraph(layoutGraph(res.g, edgeAttrs = list(label = res.annot)))
+}
+
+
+
+
+## ANALYSIS OF PANDEMIC A/H1N1 INFLUENZA DATA ##
 dat <- read.csv(system.file("files/pdH1N1-data.csv",package="adegenet"))
 ha <-  read.dna(system.file("files/pdH1N1-HA.fasta",package="adegenet"), format="fa")
 na <- read.dna(system.file("files/pdH1N1-NA.fasta",package="adegenet"), format="fa")
+
+
+## COMPUTE NUCLEOTIDIC DISTANCES
+nbNucl <- ncol(as.matrix(ha)) + ncol(as.matrix(na))
+D <- dist.dna(ha,model="raw")*ncol(as.matrix(ha)) + dist.dna(na,model="raw")*ncol(as.matrix(na))
+D <- round(as.matrix(D))
+
+
+## MATRIX OF SPATIAL CONNECTIVITY
+## (to promote local transmissions)
+xy <- cbind(dat$lon, dat$lat)
+temp <- as.matrix(dist(xy))
+M <- 1* (temp < 1e-10)
+
+
+## SEQTRACK ANALYSIS
+dat$date <- as.POSIXct(dat$date)
+res <- seqTrack(D, rownames(dat), dat$date, prox.mat=M, mu=.00502/365, haplo.le=nbNucl)
+
+
+## COMPUTE GENETIC LIKELIHOOD
+p <- get.likelihood(res, mu=.00502/365, haplo.length=nbNucl)
+# (these could be shown as colors when plotting results)
+# (but mutations will be used instead)
+
+
+## EXAMINE RESULTS
+head(res)
+tail(res)
+range(res$weight, na.rm=TRUE)
+barplot(table(res$weight)/sum(!is.na(res$weight)), ylab="Frequency", xlab="Mutations between ancestor and descendent")
+
+
+## DISPLAY SPATIO-TEMPORAL DYNAMICS 
+if(require(maps)){
+myDates <- as.integer(difftime(dat$date, as.POSIXct("2009-01-21"), unit="day"))
+myMonth <- as.POSIXct(c("2009-02-01", "2009-03-01","2009-04-01","2009-05-01","2009-06-01","2009-07-01"))
+x.month <-  as.integer(difftime(myMonth, as.POSIXct("2009-01-21"), unit="day"))
+
+
+## FIRST STAGE:
+## SPREAD TO THE USA AND CANADA
+curRange <- as.POSIXct(c("2009-03-29","2009-04-25"))
+par(bg="deepskyblue")
+map("world", fill=TRUE, col="grey")
+opal <- palette()
+palette(rev(heat.colors(10)))
+plotSeqTrack(res, round(xy), add=TRUE,annot=FALSE,lwd=2, date.range=curRange, col=res$weight+1)
+title(paste(curRange, collapse=" to "))
+legend("bottomright", lty=1, leg=0:8, title="number of mutations", col=1:9, lwd=2)
+
+
+## SECOND STAGE
+## SPREAD WITHIN AMERICA, FIRST SEEDING OUTSIDE AMERICA
+curRange <- as.POSIXct(c("2009-04-30","2009-05-07"))
+par(bg="deepskyblue")
+map("world", fill=TRUE, col="grey")
+opal <- palette()
+palette(rev(heat.colors(10)))
+plotSeqTrack(res, round(xy), add=TRUE,annot=FALSE,lwd=2, date.range=curRange, col=res$weight+1)
+title(paste(curRange, collapse=" to "))
+legend("bottom", lty=1, leg=0:8, title="number of mutations", col=1:9,lwd=2, horiz=TRUE)
+
+
+## THIRD STAGE
+## PANDEMIC
+curRange <- as.POSIXct(c("2009-05-15","2009-05-25"))
+par(bg="deepskyblue")
+map("world", fill=TRUE, col="grey")
+opal <- palette()
+palette(rev(heat.colors(10)))
+plotSeqTrack(res, round(xy), add=TRUE,annot=FALSE,lwd=2, date.range=curRange, col=res$weight+1)
+title(paste(curRange, collapse=" to "))
+legend("bottom", lty=1, leg=0:8, title="number of mutations", col=1:9,lwd=2, horiz=TRUE)
+
 }
+}
 }
\ No newline at end of file