[adegenet-forum] pairwise.fst

Jombart, Thibaut t.jombart at imperial.ac.uk
Fri Oct 5 13:05:52 CEST 2012


Hi there, 

yes, a month strikes my as a bit long indeed. I get bored after waiting 5 minutes. What is the dimension of you genind object? Is the line:
mat.obs <- pairwise.fst(x, res.type="matrix")
working? 

If so, how long does it take?
To check it out, simply do:

system.time(mat.obs <- pairwise.fst(x, res.type="matrix"))

Feel free to send me a sample of the data off the ML, I can give it a try. Valeria's got a point, hierfstat is back and is probably more computer efficient.

Best

Thibaut

--
######################################
Dr Thibaut JOMBART
MRC Centre for Outbreak Analysis and Modelling
Department of Infectious Disease Epidemiology
Imperial College - School of Public Health
St Mary’s Campus
Norfolk Place
London W2 1PG
United Kingdom
Tel. : 0044 (0)20 7594 3658
t.jombart at imperial.ac.uk
http://sites.google.com/site/thibautjombart/
http://adegenet.r-forge.r-project.org/
________________________________________
From: adegenet-forum-bounces at lists.r-forge.r-project.org [adegenet-forum-bounces at lists.r-forge.r-project.org] on behalf of Valeria Montano [mirainoshojo at gmail.com]
Sent: 04 October 2012 22:50
To: takele taye
Cc: adegenet-forum at lists.r-forge.r-project.org
Subject: Re: [adegenet-forum] pairwise.fst

Dear Takele,

I don't know if this can be the issue in your case, but I experimented problems with the pairwise.fst function as well (and I had forgotten about that). My desktop tried hard though, I remember the CPU running at its maximum for quite a while staring at the script and trying to decide what to do with it... in the end I switched to the dist.genpop function, although in my case a pairwise fst would have been much more appropriate than a molecular distance. After shuffling the inds you will get new pops with different allele frequencies but I am not sure you are interested in doing it... I am not sure either whether there is an issue with the pairwise.fst function, maybe it was just me and my desktop (and I) and it has nothing to do with your case. You can try to run just the pairwise.fst function on a small dataset and see what happens? Uhm, in the meanwhile the package hierfstat was back on CRAN so you could do pairwise fst and bootstrap with it, if I am not wrong.

Are you really waiting since one month? You are patient man.
Ok..nothing else to add

Cheers

Valeria

On 4 October 2012 22:02, takele taye <takele_taye at yahoo.com<mailto:takele_taye at yahoo.com>> wrote:
Dear

I get difficulties to estimate the significance level for pairwise.fst values. It is running on a PC for a month no output yet. Moreover, I  couldn't able to run a parallel computation as this does not allow to perform this. The input data is the one I used for STRUCTURE analysis. My script is

library (adegenet)

library(ade4)

ttt<- read.structure(file="oh.str", n.ind=94, n.loc=47486,  onerowperind=TRUE, col.lab=NULL, col.pop=1, ask=FALSE)
x <- ttt[sample(1:nrow(ttt at tab), )]
mat.obs <- pairwise.fst(x, res.type="matrix")
NBPERM <- 1000
mat.perm <- lapply(1:NBPERM, function(i) pairwise.fst(x, pop=sample(pop(x)), res.type="matrix"))

meanmatobs= mean(c(mat.obs[1,2] < na.omit(sapply(1:NBPERM, function(i) mat.perm[[i]][1,2])), TRUE))


test12 <- as.randtest(na.omit(sapply(1:NBPERM, function(i) mat.perm[[i]][1,2])), mat.obs[1,2], alter="greater")
test12
Monte-Carlo test
Call: as.randtest(sim = na.omit(sapply(1:NBPERM, function(i) mat.perm[[i]][1,
    2])), obs = mat.obs[1, 2], alter = "greater")


allTests <- list()
 for(i in 1:(nrow(mat.obs)-1)){
   for(j in 2:nrow(mat.obs)){
   allTests[[paste(rownames(mat.obs)[i],rownames(mat.obs)[j],sep="-")]] <- as.randtest(na.omit(sapply(1:NBPERM, function(k) mat.perm[[k]][i,j])), mat.obs[i,j], alter="greater")
   }
}

Any help is appreciated


_______________________________________________
adegenet-forum mailing list
adegenet-forum at lists.r-forge.r-project.org<mailto:adegenet-forum at lists.r-forge.r-project.org>
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/adegenet-forum



More information about the adegenet-forum mailing list