[adegenet-forum] Inbreeding calculation
Jombart, Thibaut
t.jombart at imperial.ac.uk
Wed May 14 08:06:10 CEST 2014
Hello,
This is a scaling problem (I think I mentioned it before, now is probably a good time to clarify it). Densities are very small and when using them to get probabilities for 'sample' all probabilities are zero. You need to use functions directly.
Note that I don't get your error, but informative warnings. I suspect your version of adegenet is outdated.
## what you've done ##
U1<- inbreeding(U, N=100)
There were 50 or more warnings (use warnings() to see the first 50)
> warnings()
Warning messages:
1: In FUN(X[[81L]], ...) :
Likelihood uniformly zero likely reflecting precision issue
returning uniformly distributed sample
2: In FUN(X[[81L]], ...) :
Likelihood uniformly zero likely reflecting precision issue
returning uniformly distributed sample
....
## the alternative ##
res <- inbreeding(U, res.type="function")
## look at the first function ##
plot(res[[1]]) # see, the y scale is 1e-245 (!!); still the graph tells you there's no inbreeding here
## look at the first 80 functions ##
par(mfrow=c(8,10), mar=rep(.1,4))
for(i in 1:80) plot(res[[i]], yaxt="n", col="red")
You'll see that a density can be estimated for most individuals. Flat bars or missing curves indicate NA problems.
If you really, really need to derive a sample for one given distribution, you have to scale the y axis:
####
x <- seq(0,1,le=1000)
mySamp <- sample(x, 100, prob=res[[1]](x)*1e250, replace=TRUE)
hist(mySamp, xlim=c(0,1))
####
But still, note that:
##
hist(mySamp, xlim=c(0,1))
##
is just an approximation of:
##
plot(res[[1]])
##
Cheers
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 Tohamy Yousef [tohamyy at yahoo.com]
Sent: 13 May 2014 17:18
To: adegenet-forum at lists.r-forge.r-project.org
Subject: [adegenet-forum] Inbreeding calculation
Dear Dr. Jombart
Firsst of all, I am sorry for disrupting you again. But, really I need your help. I wrote you before and asked about inbreeding calculation. You suggested that I use less number of markers and other advices. SO, I tried to follow all your comments. For example I decided to use only 2247 SNP markers without missing data (before I had 43231 with missing data). But I am trying a lot and still have the same problem.
I am attaching to you the new data (new_inbreeding.stru).
Thanks in advance.
######### my script
> library(adegenet)
> toto<-read.structure(file="new_inbreeding.stru",n.ind=154,n.loc=2247,onerowperind=FALSE,col.lab=1,row.marknames=1,NA.char="-9",missing="mean")
Which column contains the population factor ('0' if absent)? 2
Which other optional columns should be read (press 'return' when done)? 1:
> toto
#####################
### Genind object ###
#####################
- genotypes of individuals -
S4 class: genind
@call: read.structure(file = "inbreeding.stru",
n.ind = 154, n.loc = 2247, onerowperind = FALSE, col.lab = 1,
row.marknames = 1, NA.char = "-9", missing = "mean")
@tab: 154 x 3610 matrix of genotypes
@ind.names: vector of 154 individual names
@loc.names: vector of 2247 locus names
@loc.nall: number of alleles per locus
@loc.fac: locus factor for the 3610 columns of @tab
@all.names: list of 2247 components yielding allele names for each locus
@ploidy: 2
@type: codom
Optionnal contents:
@pop: factor giving the population of each individual
@pop.names: factor giving the population of each individual
@other: - empty -
> toto$pop.names
P1 P2
"US" "IP"
>U<- seppop(toto)$US
> U
#####################
### Genind object ###
#####################
- genotypes of individuals -
S4 class: genind
@call: .local(x = x, i = i, j = j, treatOther = ..1, quiet = ..2, drop = drop)
@tab: 81 x 3610 matrix of genotypes
@ind.names: vector of 81 individual names
@loc.names: vector of 2247 locus names
@loc.nall: number of alleles per locus
@loc.fac: locus factor for the 3610 columns of @tab
@all.names: list of 2247 components yielding allele names for each locus
@ploidy: 2
@type: codom
Optionnal contents:
@pop: factor giving the population of each individual
@pop.names: factor giving the population of each individual
@other: a list containing: elements without names
> U1<- inbreeding(U, N=100)
Error in sample.int(length(x), size, replace, prob) :
NA in probability vector
More information about the adegenet-forum
mailing list