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

substitute number with sed

I want to automatically automatically substitute the "number of levels" of the pedigree in a blupf90 parameter file. For that, I keep a basic.par:

DATAFILE

data

NUMBER_OF_TRAITS

           1

NUMBER_OF_EFFECTS

           2

OBSERVATION(S)

3

WEIGHT(S)

EFFECTS: POSITIONS_IN_DATAFILE NUMBER_OF_LEVELS TYPE_OF_EFFECT[EFFECT NESTED]

1 1 cross

2 30000 cross

RANDOM_RESIDUAL VALUES

0.7

 RANDOM_GROUP

           2

 RANDOM_TYPE

add_animal

 FILE

ped                                                

(CO)VARIANCES

0.3

then I read the actual number of levels in ped:

# get number in pedigree

n=$(wc -l ped | awk '{print $1}')

in my example, 2800

and I use sed to modify the 30000 to 'n':

# change 2nd number in line 13 to n 

sed  '13 s/[0-9]* /'$n' /2' < basic.par > new.par

I am not good with regular expressions, but I believe this means "in line 13, any number ([0-9] repeated from 0 to infinite times -that's what * means' and only substitute the 2nd occurrence".  The spaces in the regular expression are by trial and error :-(  
The little that I know about sed comes from the Grymoire .

split and substr in awk to handle genotypes

Genotype handling often involves loops across genotypes. In awk, I typically use loops. To pick one locus from the long string described in the previous post, we have two possibilities, substr() or split()

QmSim2uga.awk changes formats from output of QMSim to input of blupf90. Here is one such loop using substr() , which selects a substring from a string

for (i=1; i<=nsnp; i++){

                        out=substr($2,i,1)

                        # QmSim -> UGA

                        # 0 -> AA -> 0

                        # 2 -> aa -> 2

                        # 3 -> Aa -> 1

                        # 4 -> aA -> 1

                        # (missing does not exist) -> 5

                        if(out>2) {out=1}

                        genotype=genotype out

                }

For 9 individuals, in my mac

time ./QmSim2uga.awk p1_mrk_001.txt >/dev/null

nsnp, nanim     78519         9

user 1m47.467s

this takes 2 min .

Alternatively, we split() the string into an array, then we loop through the array:

split($2,geno,"")

for (i=1; i<=nsnp; i++){

                        out=geno[i]

                        #substr($2,i,1)

                        # QmSim -> UGA

                        # 0 -> AA -> 0

                        # 2 -> aa -> 2

                        # 3 -> Aa -> 1

                        # 4 -> aA -> 1

                        # (missing does not exist) -> 5

                        if(out>2) {out=1}

                        genotype=genotype out

                }

This form is waaaay faster:

user 0m0.285s

And with all individuals:

nsnp, nanim     78519      6037

real 3m16.279s

user 3m13.201s

format genotypes for blupf90, GS3

Blupf90 and GS3 require genotypes to be in this form:

       345 1111212111212112

       346 1121111211211021

       347 2022222220202022

       348 1111111211211021

      1349 2022222220202022

     12350 1111212111212112

       351 1121111211211021

       352 1121111211211021

       353 2022222220202022

this also works

     345   1111212111212112

       346 1121111211211021

      347  2022222220202022

       348 1111111211211021

      1349 2022222220202022

     12350 1111212111212112

       349 2022222220202022

       350 1111212111212112

       351 1121111211211021

       352 1121111211211021

       353 2022222220202022

or this

345   1111212111212112

346   1121111211211021

347   2022222220202022

348   1111111211211021

1349  2022222220202022

12350 1111212111212112

 this will be read erroneously and it will give wrong results:

345 1111212111212112

346 1121111211211021

347 2022222220202022

348 1111111211211021

1349 2022222220202022

12350 1111212111212112

Id's and genotypes (coded as 0/1/2) need to be separated by 1 or several spaces (not tabs) and genotypes need to start at exactly the same column.

A simple fix is to use awk.
Imagine that your genotype file is gene.txt :

$ cat gene.txt

345 1111212111212112

346 1121111211211021

347 2022222220202022

348 1111111211211021

1349 2022222220202022

12350 1111212111212112

then you can do

awk 'printf("%10s%1s%" length($2) "s\n",$1," ",$2) gene.txt >gene2.txt

On gene2.txt, things are formatted:

awk '{printf("%10s%1s%" length($2) "s\n",$1," ",$2)}' gene.txt >gene2.txt

cat gene2.txt 

       345 1111212111212112

       346 1121111211211021

       347 2022222220202022

       348 1111111211211021

      1349 2022222220202022

     12350 1111212111212112

Bivariate plots in R with many many pairs

Bivariate plots with many many pairs

Think that we need to plot many pairs, e.g. comparing A vs. G. What we can do in R:

#generate two correlated variables
a=rnorm(1000000)
b=a+rnorm(a)
plot is rather horrible

plot(a,b)

We can draw a map with level curves
 require(MASS)
 contour(kde2d(a,b))


But I find that hexagonal binning is even nicer (and very fast !!)
require(hexbin)
tt=hexbin(a,b)
plot(tt)

Finally, we have persp:
persp(kde2d(a,b))