[adegenet-forum] Inbreeding coefficient calculation by adegenet

[Jombart, Thibaut]

Dear Tohamy

sorry, I am currently teaching at a workshop for the week and won't have time to look into this before April. Meanwhile I would suggest trying on a sample of SNPs with growing size (e.g. 10, 50, 100, 500, etc) to ensure that nothing fishy is happening here.

Cheers
Thibaut

________________________________________
From: Tohamy Yousef [tohamyy at yahoo.com]
Sent: 21 March 2014 13:50
To: Jombart, Thibaut; adegenet-forum at lists.r-forge.r-project.org
Subject: Re: [adegenet-forum] Inbreeding coefficient calculation by adegenet

Dear Dr. Jombart
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 still have the same problem. I attached 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



On Thursday, December 19, 2013 6:32 AM, "Jombart, Thibaut" <t.jombart at imperial.ac.uk> wrote:

Hello,

this was a tricky one. There is a problem in your command line, and then there is a problem inherent to likelihood and large datasets.

First, the problem is you replace missing data by "zero"s  allele frequencies (missing=0), which means that for some loci and some individuals, you have no allele, but still not treated as missing data. This screws the definition of homozygotes and thus the computations of inbreeding.

Second, for large datasets, the sum of log-likelihoods is so low that reverting them back to likelihoods gives 0. You end up with a density distribution that is approximately zero everywhere. You may still be able to visualize it as a density, but when deriving samples, this results in probabilities of zero and 'sample' complains about it. Note that this is not a theoretical issue, only a numerical precision problem.

I have just committed a patch so that now, a meaningful warning will be issued.

In principle, one could just add a constant to the sum of log-likelihood values as a workaround. But which value to add is not a trivial choice, and may vary from one individual to another. I'll pass on that for now.

As a workaround for your problem:
- don't use "missing=0"
- ask for the distributions, not the samples: "res.type="function""
- alternatively, use a smaller subset of loci
- use the new patch (attached) to replace the error with a warning - you'll get flat (uniform) distributions for the problematic individuals.

Cheers
Thibaut
________________________________________

From: adegenet-forum-bounces at lists.r-forge.r-project.org<mailto:adegenet-forum-bounces at lists.r-forge.r-project.org> [adegenet-forum-bounces at lists.r-forge.r-project.org<mailto:adegenet-forum-bounces at lists.r-forge.r-project.org>] on behalf of Tohamy Yousef [tohamyy at yahoo.com<mailto:tohamyy at yahoo.com>]
Sent: 18 December 2013 08:44
To: adegenet-forum at lists.r-forge.r-project.org<mailto:adegenet-forum at lists.r-forge.r-project.org>
Subject: [adegenet-forum] Inbreeding coefficient calculation by adegenet

Dear Dr. Jombart,
I am trying to calculate the inbreeding coefficient for two populations (US=81 and IP=73)  with 43231 SNPs. I am using your your manual, an introduction to adgenet 1.4-0, but I have a problem in the calculation. It gives me an error message:
Error in sample.int(length(x), size, replace, prob) :
  NA in probability vector
I do not know why? I did it as following:

> toto<-read.structure(file="finalFiltered_noLowCov_e0_LOCUS_POP.stru",n.ind=154,n.loc=43231,onerowperind=FALSE,col.lab=1,row.marknames=1,NA.char="-9",missing=0)
Which column contains the population factor ('0' if absent)? 2

Which other optional columns should be read (press 'return' when done)? 1:
Converting data from a STRUCTURE .stru file to a genind object...
> is.genind(toto)
[1] TRUE
> toto$pop.names
  P1  P2
"US" "IP"
> sa1<- seppop(toto)$US
> sa1
  #####################
  ### 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 61169 matrix of genotypes
@ind.names: vector of  81 individual names
@loc.names: vector of  43231 locus names
@loc.nall: number of alleles per locus
@loc.fac: locus factor for the  61169 columns of @tab
@all.names: list of  43231 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

> temp1<- inbreeding(sa1, N=81)
Error in sample.int(length(x), size, replace, prob) :
  NA in probability vector

> sa2<- seppop(toto)$IP
> sa2
  #####################
  ### Genind object ###
  #####################
- genotypes of individuals -
S4 class:  genind
@call: .local(x = x, i = i, j = j, treatOther = ..1, quiet = ..2, drop = drop)
@tab:  73 x 61169 matrix of genotypes
@ind.names: vector of  73 individual names
@loc.names: vector of  43231 locus names
@loc.nall: number of alleles per locus
@loc.fac: locus factor for the  61169 columns of @tab
@all.names: list of  43231 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

> temp2<- inbreeding(sa2, N=73)
Error in sample.int(length(x), size, replace, prob) :
  NA in probability vector


Could you please tell me how can I fix this problem?
I attached you a part of my data. and I hope that  help me.

Thank you advance
Best regards,
Tohamy
-------------------------------------------
Eltohamy Yousef  M.Sc.
Crop Biodiversity and Breeding Informatics
Institute of Plant Breeding, Seed Science and Population Genetics (350)
University of Hohenheim
Fruwirtstrasse 21, 70599 Stuttgart
Office phone:0049711 459-24437
Email: tohamyy at yahoo.com<mailto:tohamyy at yahoo.com>
            E.yousef at uni-hohenheim.de<mailto:E.yousef at uni-hohenheim.de>
            tohamy_yousef at agr.suez.edu.eg<mailto:tohamy_yousef at agr.suez.edu.eg>
Web.www.evoplant.uni-hohenheim.de