Burguete horses and Pirenaica cattle

In Navarre: Burguete mare and foal, and Pirenaica cows.

Chisqueta sheep in Sobrarbe

Still a sheep breed that I heard about but never saw, the Chisqueta or Xisqueta sheep. Here are some references in aragonese and catalan, and here is one in Spanish.

Variance components with high numbers

I have this large data set in which I want to estimate heritability. These are the features for the trait:

n  nmiss   avg         min       max         var      

------------------------------------------------------------------------------

1703515 0 1772.570604    100.470000 7783.789000 524507.442442

So, the variance for the trait is very high. Blupf90 runs nicely. However, both airemlf90 and gibbs2f90 variance component programs in the suite gave me some trouble if I used. The reason is that the programs contain some "tolerance" parameter in the inverses of the updated variance components, and this seems to set some of the numbers to 0.

For instance, if divide the value of the performance by 10, I get this:

 * Start Gibbs iteration04-29-2016  16h 26m 59s 482

 round            1

 G

   768.75    

 G

   157.18    

 R

   0.0000    

 round            2

The reason for this is the scale and the large number of data. The sum of squares of the residual is

 ee=  0.16456E+10

when this is sent to the subroutine, its inverse is lower than 10^-9 and set to 0.

I finally standardized MY to a variance of 1. I put initial variances to $latex \sigma^2_g=\sigma^2_p=\sigma^2_e=1 $latex. After 1000 iterations of Gibbs, I get:

* End of iteration04-29-2016  12h 09m 52s 514

 ave G

  0.11236    

 SD G

  0.31784E-01

 ave G

  0.14050    

 SD G

  0.13210E-01

 ave R

  0.25338    

 SD R

  0.58904E-02

In this case I obtain correct estimates ($latex h^2 =0.22 $latex, repeatability of 0.5).

Standard errors of complex functions of variance components from REML Information matrix

In airemlf90 and other programs, a byproduct of the AIREML is the inverse of the Information Matrix that, asymptotically, describes the sampling covariance of the variance components. It can also be considered as an approximation to the posterior distribution of the parameters.

Getting standard errors of functions of the variance components is rather painful algebraically as one requires the delta method and Taylor series expansions. A much simpler method , implemented in as OPTION se_covar_function in airemlf90, comes through MonteCarlo sampling and was described by Meyer and Houle. Here I will show how this work.
This a model with three variance components, genetic, permanent and residual variance. Estimates from airemlf90 are in vector ahat:

#genetic, permanent, residual

ahat=c(

  0.11478,    

  0.13552,    

  0.25290,

  )

with AI matrix:

# inverse of AI matrix (Sampling Variance)

AI=matrix(c(

  0.16799E-05, -0.96486E-06, -0.82566E-08,

 -0.96486E-06,  0.96167E-06, -0.37113E-07,

 -0.82566E-08, -0.37113E-07,  0.10864E-06)

 ,ncol=3)

asymptotically, the "true" value has a sampling (or posterior) distribution a~N(ahat,AI), and getting samples from this distribution is easy:

require(MASS)

b=mvrnorm(10000,ahat,AI)

> head(b)

          [,1]      [,2]      [,3]

[1,] 0.1146738 0.1357640 0.2529399

[2,] 0.1163889 0.1342926 0.2528479

[3,] 0.1166155 0.1344342 0.2525161

[4,] 0.1142085 0.1358928 0.2534974

[5,] 0.1136835 0.1361108 0.2530133

[6,] 0.1140485 0.1365707 0.2530573

and from here, we can obtain samples of the heritability and its standard deviation:

h2=b[,1]/(b[,1]+b[,2]+b[,3])

sd(h2)

> 0.002318198

This is very useful, i.e. for complex model such as consideration of genomic and dominance variation.

Equations in blogger

http://irrep.blogspot.fr/2011/07/mathjax-in-blogger-ii.html 

Seems to work: $h^2=1$

UPDATE: it does work but then it messes up the $'s corresponding to awk expressions so I desactivate it.

Estimation of allelic base frequencies

Tao Xiang is estimating base frequencies in a pig population. A simple but efficient method is Gengler et al. 2007 . This is actually very similar to GCMTBLUP, nut the latter considers one locus and an associated trait. Thus, Tao implement the model g=1m+Za+e  using blupf90 software, where m=2p. The data file, with genotype coded as 0/1/2 looks like

6051 1 2

6930 1 2

5205 1 2

7478 1 0

3208 1 1

Using direct inversion by fspak, the solution is p=0.88 . However, using the default iterative solver of blupf90 we get 0.44 !! Actually, using a more strict convergence criteria we get the exact solution:

  convcrit iter      20*m         p

1       12    5  9.016301 0.4508151

2       13    5  9.016301 0.4508151

3       14    9  9.016345 0.4508172

4       15   11  9.016358 0.4508179

5       16  103 16.779241 0.8389621

6       17  110 17.503890 0.8751945

7       18  117 17.570641 0.8785320

8       19  124 17.574420 0.8787210

9       20  126 17.574496 0.8787248

where the convergence criteria is actually 10^(-convcrit). So, the iterator is stuck at 0.45. The reason for this is that heritability is 1 so that for a given mu, EBVs for genotyped animals almost do not change and viceversa.

This is how the correct estimate of EBV looks like:

There is no pathology in this graph. The first 7000 individuals have genotype, so the breeding value is simply the deviation of the genotype from the populational mean. The three bars are individuals in genotypes AA, Aa and aa. The rest of individuals are ancestors, vaguely recoded in ascending order, so that their genotype is not neatly estimated.

This is how the wrong estimate of EBV looks like:

With biased EBVs for genotyped animals and, for non genotyped animals, EBVs quite different from the previous graph.

Take home message: if you estimate base allele frequencies using Gengler method by BLUP, use either solutions by inversion or set the convergence criteria to strict values.

Installing Julia and Jupyter in Windows 7

Rohan Fernando came to INRA Toulose last two weeks and he showed us Julia and Jupyter. They look very nice. The instructions to install both are here and here. I installed them in my Mac without much trouble but in the Windows 7 machine of Jean Michel Elsen was more difficult. There were issues related to WinRPM that I did not understand.


To install in JM computer I did

1. Clean up everything, desinstall Julia, Anaconda, etc if already installed; delete folders .Julia, .Jupyter, etc.
2. Install Julia
3. Install Anaconda (Python >3)
4. Start Julia. In its command line, Pkg.add("IJUlia") [ this is not clearly documented]. This seems to install Jupyter automatically. Let Julia open.
5. Open a command line in Windows (in "Accesories"), then write jupyter notebook. The default navigator (Firefox) should open.

It works correctly in the Windows 7 Virtual Machine in my Mac