Plotting histograms and line plots

[nicolas]

Hi,

I managed to run simulations and produce rejctions files for my
dataset (3 populations). Designed a prior file accordingly (3 pops,
3Ne, 3 AncNe, 3 mig rates,etc) but when I try and plot histograms or
plot lines in R, it only returns the following plots:

- mut rate
- recomb rate
- tev
- Ne1
- Ne2
- NeA
- mig1
- mig2

I guess I have to manually change the R script and add in it commands
for missing plots such as:

# Plot histograms for effective size of population 3 and save it in
a .eps file
hist(abc.rej[,12],freq=F,main="prior (black) and posterior (blue)
distributions",xlab="Ne3",col="blue",density=9)
hist(priors[,12],freq=F,col="white",add=T)
dev.copy2eps(file="Ne3.eps", horizontal=F)
print("Ne3 done.")

Right? I guess I also have to change that bit "abc.rej[,12]" ? Does
that 12 stands for parameter number 12? if yes, how do I know what is
the number of all my missing parameters?

Cheers

Nic

[Joao Sollari Lopes]

Hi Nic,

Good to see that you solved your previous problems.
You are absolutely right regarding the usage of the R scripts.

The number of the parameters correspond to the columns in the .dat
files. These can be obtained from the the .txt files produced by
running the simulations.

Note that you should also change the columns for the summary
statistics when performing the regression step.

Hope that helps,
Joao

[nicolas]

Hi Joao,

Thanks heaps for that!

Nic

[nicolas]

Hi Joao,

I tried to modify the Plot_line file provided with popABC but I get
the following message error I get after typing the command line:

> plot_line(rej="NZFS5rej.rej",pri="NZFS5rej.pri")


Erreur dans plot.window(...) : valeurs finies requises pour 'ylim'
De plus : Messages d'avis :
1: In lfproc(x, y, weights = weights, cens = cens, base = base, geth =
geth, :
procv: density estimate, empty integration region
2: In lfproc(x, y, weights = weights, cens = cens, base = base, geth =
geth, :
procv: density estimate, empty integration region
3: In min(x) : aucun argument trouvé pour min ; Inf est renvoyé
4: In max(x) : aucun argument pour max ; -Inf est renvoyé


Here's the first bit of the modified script:

"plot_line <- function(rej_file , pri_file){

#import locfit library
library(locfit)

#demographic parameters' priors (minimum and maximum values)
mintev1 <- 0
maxtev1 <- 10000
mintev2 <- 0
maxtev2 <- 10000
minNe1 <- 0
maxNe1 <- 100000
minNe2 <- 0
maxNe2 <- 100000
minNe3 <- 0
maxNe3 <- 100000
minNeAnc1 <- 0
maxNeAnc1 <- 100000
minNeAnc2 <- 0
maxNeAnc2 <- 100000
minmig1 <- 0
maxmig1 <- 0.0001
minmig2 <- 0
maxmig2 <- 0.0001
minmig3 <- 0
maxmig3 <- 0.0001
minmigAnc1 <- 0
maxmigAnc1 <- 0.0001

#import the .rej file
abc.rej <- data.matrix(read.table(rej_file))

#import the .pri files
priors <- data.matrix(read.table(pri_file))

# Plot line for sequence DNA mutation rate and save it in a .eps file
plot(locfit(~abc.rej[,2]),main="prior (black) and posterior (blue)
distributions",xlab="mut rate",col="blue")
plot(locfit(~priors[,2]),col="black",add=T)
dev.copy2eps(file="mutrate_line.eps", horizontal=F)
print("mut rate done.")

# Plot line for sequence DNA recombination rate and save it in a .eps
file
plot(locfit(~abc.rej[,4]),main="prior (black) and posterior (blue)
distributions",xlab="rec rate",col="blue")
plot(locfit(~priors[,4]),col="black",add=T)
dev.copy2eps(file="recrate_line.eps", horizontal=F)
print("rec rate done.")

# Plot line for splitting time 1 and save it in a .eps file
plot(locfit(~abc.rej[,6],xlim=c(mintev1,maxtev1)),main="prior (black)
and posterior (blue) distributions",xlab="tev1",col="blue")
plot(locfit(~priors[,6],xlim=c(mintev1,maxtev1)),col="black",add=T)
dev.copy2eps(file="tev1_line.eps", horizontal=F)
print("tev1 done.")

# Plot line for splitting time 2 and save it in a .eps file
plot(locfit(~abc.rej[,7],xlim=c(mintev2,maxtev2)),main="prior (black)
and posterior (blue) distributions",xlab="tev2",col="blue")
plot(locfit(~priors[,7],xlim=c(mintev2,maxtev2)),col="black",add=T)
dev.copy2eps(file="tev2_line.eps", horizontal=F)
print("tev2 done.")
......"

Is the problem due to the priors range in the plot_line script?

Thanks!

[Joao Sollari Lopes]

Hi Nic,

The error message comes from lfproc( ), a function from the locfit
library.
My guess is that the prior distribution for the mutation rate was set
to be constant, so column 2 in your .dat file is just a constant
number. That presents a problem when using locfit.
Anyway, if that is the case there is no reason to plot the posterior
distribution for mutation. Just delete or comment those lines in the
script.

Let me know if that helped,
Joao

[nicolas]

Hi Joao,

Yes it worked fine! :-)

Cheers

Nic

On 14 avr, 04:31, Joao Sollari Lopes <j.sollari.lo...@gmail.com>

> > > > > # Plothistogramsfor effective size of population 3 and save it in