fastsimcoal2
Christine Dudgeon
Bernd.Gruber
Hi Chrisitine,
I recently played with dartR and fastsimcoal2 and have some snippets lying around. Can you be more specific. I guess you want a function that generates a multidimensional sfs? If so happy to provide a function.
Also out of interest I use strataG to interface between R and fastsimcoal2. What is your approach.
Cheers, Bernd
==============================================================================
Dr Bernd Gruber )/_
_.--..---"-,--c_
Professor Ecological Modelling \|..' ._O__)_
Tel: (02) 6206 3804 ,=. _.+ _ \..--( /
Fax: (02) 6201 2328 \\.-''_.-' \ ( \_
Institute for Applied Ecology `''' `\__ /\
Faculty of Science and Technology ')
University of Canberra ACT 2601 AUSTRALIA
Email: bernd....@canberra.edu.au
WWW: bernd-gruber
Australian Government Higher Education Provider Number CRICOS #00212K
NOTICE & DISCLAIMER: This email and any files transmitted with it may contain
confidential or copyright material and are for the attention of the addressee
only. If you have received this email in error please notify us by email
reply and delete it from your system. The University of Canberra accepts
no liability for any damage caused by any virus transmitted by this email.
==============================================================================
--
You received this message because you are subscribed to the Google Groups "dartR" group.
To unsubscribe from this group and stop receiving emails from it, send an email to
dartr+un...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dartr/ab8b50c2-5e9d-4dd3-aec2-3ef6cb21077cn%40googlegroups.com.
Christine Dudgeon
You received this message because you are subscribed to a topic in the Google Groups "dartR" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dartr/2SzqUIhkwPw/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dartr+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dartr/ME3PR01MB7944B0369781CAE00AD2A216D45D9%40ME3PR01MB7944.ausprd01.prod.outlook.com.
Christine Dudgeon
You received this message because you are subscribed to a topic in the Google Groups "dartR" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dartr/2SzqUIhkwPw/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dartr+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dartr/ME3PR01MB7944B0369781CAE00AD2A216D45D9%40ME3PR01MB7944.ausprd01.prod.outlook.com.
Bernd.Gruber
Hi Chris,
I quickly wrote the msfs function (untested, but it seems to work).
#function multidimensional sfs (untested)
gl.msfs<- function(x) {
pp <- seppop(x)
mi <- max(sapply(pp,nInd))
mm <- sapply(pp, function(x) table((0.5-abs(0.5-gl.alf(x)[,1]))*nInd(x)*2), simplify = FALSE)
# cut off zero if necessary....
mm2 <- sapply(mm, function(x) x[names(x)!="0"], simplify = FALSE)
#max number of individuals (fill with zeros)
mm3 <- t(sapply(mm2, function(x) c(x, rep(0, max(0,mi-length(x))))))
colnames(mm3)<- paste0("d",1:mi)
return(mm3)
}
msfs <- gl.msfs(possums.gl)
If you test if can you let me know if it also works for you, then I can add it to the next dartR version.
Also in case you are interested here is my workflow for a parameter estimation in fastsimcoal2 in R (converting the arp file in a genlight object is fairly simple).
This is an example mainly based on the fscTutorial using simulated SNPs from fastsimcoal and try to estimate parameters from the simulated data set.
In the interim I load the simulated data set, convert to dartR (genlight), calculate the SFS and then run an estimation in fastsimcoal (though the repeats [num.sims] are not high enough here)
Cheers, Bernd
#Code to run fastsimcoal from R (based mainly on fscTutorial)
#Simulate some data with a bottleneck 50 generations ago...
#coalescent -> SFS -> estimate parameters
library(strataG)
library(dartR)
setwd("~/fsc/fsc27_linux64/")
demes <- fscSettingsDemes(fscDeme(deme.size = 50, sample.size = 50))
genetics <- fscSettingsGenetics(fscBlock_snp(sequence.length = 10, mut.rate = 1e-5), num.chrom = 5000)
events <- fscSettingsEvents(
fscEvent(
event.time = 50,
source = 0,
sink = 0,
prop.migrants = 0,
new.size = 0.5,
new.growth = 0,
migr.mat = 0
),
fscEvent(
event.time = 100,
source = 0,
sink = 0,
prop.migrants = 0,
new.size = 10,
new.growth = 0,
migr.mat = 0
))
library(dartR)
library(strataG)
p <- fscWrite(demes = demes, genetics = genetics, events = events, label = "habitat_destroy")
obs.p <- fscRun(p, dna.to.snp = TRUE, no.arl.output = FALSE, num.cores = 3, sfs.type = "maf" , all.sites = F)
#convert to genlight
snp.df <- fscReadArp(obs.p)
gl <- gtypes2genlight(df2gtypes(snp.df, ploidy = 2))
gl.smearplot(gl)
#not filtered for multiple snp on the 10 base sequence!
#folded sfs
sfs <- table((0.5-abs(0.5-gl.alf(gl)[,1]))*nInd(gl)*2)
#we do not need/know/use the first frequency [number of unaltered loci]
sfs_fsc <- as.numeric(sfs[-1])
barplot(sfs_fsc)[2:10]
#for a comparison I look at the fsc SFS
#sfs from fastsimcoal simulation
obs.sfs <- fscReadSFS(obs.p)
barplot((obs.sfs$sfs$marginal[[1]])[2:10])
#test if parameters can be reconstructed....(not run long enought so fairly poor)
demes <- fscSettingsDemes(fscDeme(deme.size = "NCUR", sample.size = 50))
events <- fscSettingsEvents(
fscEvent(event.time = "TBOTSTART", 0, 0, 0, new.size = "NBOT"),
fscEvent(event.time = "TBOTEND", 0, 0, 0, new.size =10)
)
est <- fscSettingsEst(
fscEstParam("NCUR", is.int = TRUE, distr = "unif", min = 10, max = 1000),
# default for is.int = TRUE and distr = "unif"
# fscEstParam("NANC", min = 10, max = 1000),
fscEstParam("NBOT", min = 5, max = 200),
fscEstParam("TBOTSTART", min = 1, max = 200),
# fscEstParam("RESEND", is.int = FALSE, value = "NBOT/NCUR", output = FALSE),
# fscEstParam("RESBEG", is.int = FALSE, value = "NANC/NBOT", output = FALSE),
fscEstParam("TBOTEND", value = "TBOTSTART+50", output = FALSE),
obs.sfs = sfs_fsc
)
#est
est.p <- fscWrite(
demes = demes,
events = events,
genetics = fscSettingsGenetics(fscBlock_freq(1e-5)),
est = est,
label = "est.habdestroy"
)
system.time(est.p <- fscRun(est.p, num.sims = 200000, num.cores = 10))
#not sure what this function does as it does not read all outputs....
param.est <- fscReadParamEst(est.p )
#check bestlhoods
(readLines(file.path( est.p$folder, "est.habdestroy/est.habdestroy.bestlhoods")))
To view this discussion on the web visit https://groups.google.com/d/msgid/dartr/CADZtJzHmaBNCqmhetLU-2njuff3EQF7FkeseCedFgsV%2B-1QhRA%40mail.gmail.com.