Gregor Gorjanc (gg)

Compress PDF

noreply@blogger.com (Gregor Gorjanc) — Tue, 03 Jan 2012 09:58:00 +0000

See nice description here. In short using the command bellow on my Linux box (Windows users see PDFCreator) allowed me to reduce scanned color PDF for a factor of 4 to 5, e.g., from ~14 MB to ~3 MB. Supposedly PDFCreator also

gs -sDEVICE=pdfwrite -dMaxSubsetPct=100 -dPDFSETTINGS=/ebook -sOutputFile="out.pdf" -dNOPAUSE -dBATCH "in.pdf"

Update (2012-01-16): you can get a nice shell script called compressPDF.sh that saves some typing and can process multiple files, e.g.,

compressPDF.sh file.pdf

VIPER - cool pedigree&genotype data view

noreply@blogger.com (Gregor Gorjanc) — Fri, 16 Dec 2011 19:47:00 +0000

Interested in genetics - in looking genotype data on pedigrees? Check out VIPER (Visual Interactive Pedigree ExploreR). Looks really cool!!!

Ocenjevanje parametrov v Bayesovski statistiki (IBMI)

noreply@blogger.com (Gregor Gorjanc) — Mon, 24 Oct 2011 22:17:00 +0000

Prosojnice z današnje predstavitve o ocenjevanju parametrov v Bayesovski statistiki na IBMI so na voljo tukaj.

Ocenjevanje parametrov v Bayesovski statistiki

Recently discovered study resources

noreply@blogger.com (Gregor Gorjanc) — Thu, 13 Oct 2011 20:18:00 +0000

Khan Academy - a set of 2600 free videos from basic to more advances quantitative skills
Introduction to Artifical Inteligence - Stanford on-line course for masses
Machine Learning - Stanford on-line course for masses

Ovčereja na Norveškem

noreply@blogger.com (Gregor Gorjanc) — Wed, 05 Oct 2011 11:43:00 +0000

Ovčereja na Norveškem

ASD 2011 - Primošten

noreply@blogger.com (Gregor Gorjanc) — Tue, 04 Oct 2011 14:22:00 +0000

Animal Science Days (ASD) conference was held in Primošten this year. I participated with the following contributions as a co-author:

Lambing Interval in Jezersko-Solčava and Improved Jezersko-Solčava Breed. PDF
Partitioning of Genetic Trends by Originin Croatian Simmental Cattle. PDF
Estimation of Variance Components for Litter Size in the First and Later Parities in Improved Jezersko-Solcava Sheep. PDF
The Effect of Coenzyme Q10 and Lipoic Acid Added to the Feed of Hens on Physical Characteristics of Eggs. PDF

InterBull: Partitioning of international genetic trends by origin in Brown Swiss bulls

noreply@blogger.com (Gregor Gorjanc) — Thu, 29 Sep 2011 22:59:00 +0000

I attended the InterBull meeting this year in Stavanger (Norway), which was jointly organised with the EAAP conference in the same place. I participated with contribution titled "Partitioning of international genetic trends by origin in Brown Swiss bulls" co-authored with colleagues. In essence we partitioned breeding values of Brown Swiss animals by origin of selection and summarized those partitions as origin specific genetic trends. This gives us an opportunity to evaluate the effect of selection performed in different countries and how this affects the global genetic trends. Results for this breed are quite shocking! Look into the paper and talk (see bellow) for more ;)

We were looking forward for comments and/or critiques about the applied method from the audience, but did not get any direct questions. Several people approached me after the talk and one of the comments was that our method is nothing more than multiplying breeding values by the share of genes coming from different origin. This is not true and I will demonstrate this in with a simple example.

Let us assume that we have a simple small pedigree as shown bellow with R code, where column names are: id =individual code, fid = father code, mid = mother code, ori = origin, and bv = breeding value. This pedigree could represent situation where we constantly use foreign sires - here only two generations are being shown. Origin represents the country of registering/selecting aninmal.

## Simple example
example <- data.frame( id=c(  1,   2,   3,   4,   5),
                      fid=c( NA,  NA,   1,  NA,   3),
                      mid=c( NA,  NA,   2,  NA,   4),
                      ori=c("A", "B", "A", "B", "A"),
                       bv=c(100, 106, 104, 106, 103))

Now we would like to perform gene proportion analysis and partitioning of supplied breeding values, both according to origin. This is easy to achieve "by hand" for this simple example (see paper bellow for maths), but tedious for bigger examples. I have wrote an R package partAGV (not yet publicly available, but you can contact me) that can be used for such analyses.

## For gene proportion analysis
example$gp <- 1
 
library(partAGV)
partAGV(example, colAGV=6:5)
 
## Gene proportions
##    id  fid  mid ori gp gp_pa gp_w gp_A gp_B
## 1   1 <NA> <NA>   A  1     0    1 1.00 0.00
## 2   2 <NA> <NA>   B  1     0    1 0.00 1.00
## 3   3    1    2   A  1     1    0 0.50 0.50
## 4   4 <NA> <NA>   B  1     0    1 0.00 1.00
## 5   5    3    4   A  1     1    0 0.25 0.75
 
## Partitions of breeding values 
##    id  fid  mid ori  bv bv_pa bv_w  bv_A  bv_B
## 1   1 <NA> <NA>   A 100     0  100 100.0   0.0
## 2   2 <NA> <NA>   B 106     0  106   0.0 106.0
## 3   3    1    2   A 104   103    1  51.0  53.0
## 4   4 <NA> <NA>   B 106     0  106   0.0 106.0
## 5   5    3    4   A 103   105   -2  23.5  79.5

As we can see gene proportions are as expected: 1/2 for each origin in individual 3 and 1/4 vs. 3/4 for individual 5. Partitions of breeding values show that in animal 5 we 23.5 (out of 105) is attributed to selection work done in country A and 79.5 is attributed to selection work done in country B. Now, if we multiply breeding value of individual 5 (105) by origin specific gene proportions (1/4 and 3/4) we get 25.75 and 77.5, which is similar to partitions, but not the same. This shows that our method enables separation of gene flow and selection work preformed by particular country. In order to see this algebraically we can write breeding values of this individual as:

a_5 = 1/2 a_3 + 1/2 a_4 + w_5

= 1/2 (1/2 a_1 + 1/2 a_2 + w_3) + 1/2 w_4 + w_5

= 1/4 a_1 + 1/4 a_2 + 1/2 w_3 + 1/2 w_4 + w_5

= 1/4 w_1 + 1/4 w_2 + 1/2 w_3 + 1/2 w_4 + w_5

= 1/4 100 + 1/4 106 + 1/2 1 + 1/2 106 + -2

= 25 + 26.5 + 0.5 + 53 + -2

If we now collect terms specific to each origin we get:

a_5_A = 25 + 0.5 + + -2 = 23.5

a_5_B = 26.5 + 53 = 79.5

Partitioning of international genetic trends by origin in Brown Swiss bulls

Polyploidy in sugarcane

noreply@blogger.com (Gregor Gorjanc) — Wed, 28 Sep 2011 12:50:00 +0000

While reading UseR conference abstracts I came across this sentence: "Sugarcane is polypoid, i.e., has 8 to 14 copies of every chromosome, with individual alleles in varying numbers." Vau! This generates really complex genotype system. Say we have biallelic gene with alleles being A and B. In diploids the possible genotypes are AA, AB, and BB. Given the above sentence in sugarcane possible genotypes are any permutation of A's and B's in a series of 8 to 14 alleles. I am not sure if 9, 11, and 13 are also allowed, that is having even number of chromosomes. In any case such permutations result in really large numbers!

Thinking about this a bit further it appears that the whole system is not that complex once we realize that genotyping does not tell as about the order of alleles (we can not distinguish between AB and BA), which simplifies from all possible permutations to all possible combinations, e.g., for biallelic gene in tetraploids this would correspond to 5 combinations and 16 permutations.

Bellow is an R snippet that shows how to enumerate all possible combinations or permutations

## Load package having nice combinatorial functions
library(package="gtools")
 
## Specify alleles - just two for simplicity
alleles <- c("A", "B")
 
## Possible genotypes for diploids
combinations(n=length(alleles), r=2, v=alleles, repeats.allowed=TRUE)
##      [,1] [,2]
## [1,] "A"  "A" 
## [2,] "A"  "B" 
## [3,] "B"  "B" 
 
## Possible genotypes for tetraploids
combinations(n=length(alleles), r=4, v=alleles, repeats.allowed=TRUE)ž
##      [,1] [,2] [,3] [,4]
## [1,] "A"  "A"  "A"  "A" 
## [2,] "A"  "A"  "A"  "B" 
## [3,] "A"  "A"  "B"  "B" 
## [4,] "A"  "B"  "B"  "B" 
## [5,] "B"  "B"  "B"  "B" 
 
permutations(n=length(alleles), r=4, v=alleles, repeats.allowed=TRUE)
##       [,1] [,2] [,3] [,4]
##  [1,] "A"  "A"  "A"  "A" 
##  [2,] "A"  "A"  "A"  "B" 
##  [3,] "A"  "A"  "B"  "A" 
##  [4,] "A"  "A"  "B"  "B" 
##  [5,] "A"  "B"  "A"  "A" 
##  [6,] "A"  "B"  "A"  "B" 
##  [7,] "A"  "B"  "B"  "A" 
##  [8,] "A"  "B"  "B"  "B" 
##  [9,] "B"  "A"  "A"  "A" 
## [10,] "B"  "A"  "A"  "B" 
## [11,] "B"  "A"  "B"  "A" 
## [12,] "B"  "A"  "B"  "B" 
## [13,] "B"  "B"  "A"  "A" 
## [14,] "B"  "B"  "A"  "B" 
## [15,] "B"  "B"  "B"  "A" 
## [16,] "B"  "B"  "B"  "B" 
 
## Possible genotypes for 8-14 ploids
spectrum <- seq(from=8, to=14, by=2)
nS <- length(spectrum)
retC <- vector(mode="list", length=nS)
retP <- vector(mode="list", length=nS)
for(i in 1:nS) {
  retC[[i]] <- combinations(n=length(alleles), r=spectrum[i], v=alleles, repeats.allowed=TRUE)
  retP[[i]] <- permutations(n=length(alleles), r=spectrum[i], v=alleles, repeats.allowed=TRUE)
}
combC <- sapply(retC, nrow)
combP <- sapply(retP, nrow)
cbind(spectrum, combC, combP)
##      spectrum combC combP
## [1,]        8     9   256
## [2,]       10    11  1024
## [3,]       12    13  4096
## [4,]       14    15 16384

Created by Pretty R at inside-R.org

Setup up the inverse of additive relationship matrix in R

noreply@blogger.com (Gregor Gorjanc) — Thu, 11 Aug 2011 22:58:00 +0000

Additive genetic covariance between individuals is one of the key concepts in (quantitative) genetics. When doing the prediction of additive genetic values for pedigree members, we need the inverse of the so called numerator relationship matrix (NRM) or simply A. Matrix A has off-diagonal entries equal to numerator of Wright's relationship coefficient and diagonal elements equal to 1 + inbreeding coefficient. I have blogged before about setting up such inverse in R using routine from the ASReml-R program or importing the inverse from the CFC program. However, this is not the only way to "skin this cat" in R. I am aware of the following attempts to provide this feature in R for various things (the list is probably incomplete and I would grateful if you point me to other implementations):

pedigree R package has function makeA() and makeAinv() with obvious meanings; there is also calcG() if you have a lot of marker data instead of pedigree information; there are also some other very handy functions calcInbreeding(), orderPed(), trimPed(), etc.
pedigreemm R package does not have direct implementation to get A inverse, but has all the needed ingredients, which makes the package even more interesting
MCMCglmm R package has function inverseA() which works with pedigree or phlyo objects; there are also handy functions such as prunePed(), rbv(), sm2asreml(), etc.
kinship and kinship2 R packages have function kinship() to setup kinship matrix, which is equal to the half of A; there are also nice functions for plotting pedigrees etc. (see also here)
see also a series of R scripts for relationship matrices

As I described before, the interesting thing is that setting up inverse of A is easier and cheaper than setting up A and inverting it. This is very important for large applications. This is an old result using the following matrix theory. We can decompose symmetric positive definite matrix as A = LU = LL' (Cholesky decomposition) or as A = LDU = LDL' (Generalized Cholesky decomposition), where L (U) is lower (upper) triangular, and D is diagonal matrix. Note that L and U in previous two equations are not the same thing (L from Cholesky is not equal to L from Generalized Cholesky decomposition)! Sorry for sloppy notation. In order to confuse you even more note that Henderson usually wrote A = TDT'. We can even do A = LSSU, where S diagonal is equal to the square root of D diagonal. This can get us back to A = LU = LL' as LSSU = LSSL' = LSS'L' = LS(LS)' = L'L (again be ware of sloppy notation)! The inverse rule says that inv(A) = inv(LDU) = inv(U) inv(D) inv(L) = inv(L)' inv(D) inv(L) = inv(L)' inv(S) inv(S) inv(L). I thank to Martin Maechler for pointing out to the last (obviously) bit to me. In Henderson's notation this would be inv(A) = inv(T)' inv(D) inv(T) = inv(T)' inv(S) inv(S) inv(T) Uf ... The important bit is that with NRM (aka A) inv(L) has nice simple structure - it shows the directed graph of additive genetic values in pedigree, while inv(D) tells us about the precision (inverse variance) of additive genetic values given the additive genetic values of parents and therefore depends on knowledge of parents and their inbreeding (the more they are inbred less variation can we expect in their progeny). Both inv(L) and inv(D) are easier to setup.

Packages MCMCglmm and pedigree give us inv(A) directly (we can also get inv(D) in MCMCglmm), but pedigreemm enables us to play around with the above matrix algebra and graph theory. First we need a small example pedigree. Bellow is an example with 10 members and there is also some inbreeding and some individuals have both, one, or no parents known. It is hard to see inbreeding directly from the table, but we will improve that later (see also here).

ped <- data.frame( id=c(  1,   2,   3,   4,   5,   6,   7,   8,   9,  10),
                  fid=c( NA,  NA,   2,   2,   4,   2,   5,   5,  NA,   8),
                  mid=c( NA,  NA,   1,  NA,   3,   3,   6,   6,  NA,   9))

Now we will create an object of a pedigree class and show the A = U'U stuff:

## install.packages(pkgs="pedigreemm")
library(package="pedigreemm")
 
ped2 <- with(ped, pedigree(sire=fid, dam=mid, label=id))
 
U <- relfactor(ped2)
A <- crossprod(U)
 
round(U, digits=2)
## 10 x 10 sparse Matrix of class "dtCMatrix"                                     
##  [1,] 1 . 0.50 .    0.25 0.25 0.25 0.25 . 0.12
##  [2,] . 1 0.50 0.50 0.50 0.75 0.62 0.62 . 0.31
##  [3,] . . 0.71 .    0.35 0.35 0.35 0.35 . 0.18
##  [4,] . . .    0.87 0.43 .    0.22 0.22 . 0.11
##  [5,] . . .    .    0.71 .    0.35 0.35 . 0.18
##  [6,] . . .    .    .    0.71 0.35 0.35 . 0.18
##  [7,] . . .    .    .    .    0.64 .    . .   
##  [8,] . . .    .    .    .    .    0.64 . 0.32
##  [9,] . . .    .    .    .    .    .    1 0.50
## [10,] . . .    .    .    .    .    .    . 0.66

## To check

U - chol(A)

round(A, digits=2)
## 10 x 10 sparse Matrix of class "dsCMatrix"                                     
##  [1,] 1.00 .    0.50 .    0.25 0.25 0.25 0.25 .   0.12
##  [2,] .    1.00 0.50 0.50 0.50 0.75 0.62 0.62 .   0.31
##  [3,] 0.50 0.50 1.00 0.25 0.62 0.75 0.69 0.69 .   0.34
##  [4,] .    0.50 0.25 1.00 0.62 0.38 0.50 0.50 .   0.25
##  [5,] 0.25 0.50 0.62 0.62 1.12 0.56 0.84 0.84 .   0.42
##  [6,] 0.25 0.75 0.75 0.38 0.56 1.25 0.91 0.91 .   0.45
##  [7,] 0.25 0.62 0.69 0.50 0.84 0.91 1.28 0.88 .   0.44
##  [8,] 0.25 0.62 0.69 0.50 0.84 0.91 0.88 1.28 .   0.64
##  [9,] .    .    .    .    .    .    .    .    1.0 0.50
## [10,] 0.12 0.31 0.34 0.25 0.42 0.45 0.44 0.64 0.5 1.00

Note that pedigreemm package uses Matrix classes in order to store only what we need to store, e.g., matrix U is triangular (t in "dtCMatrix") and matrix A is symmetric (s in "dsCMatrix"). To show the generalized Cholesky A = LDU (or using Henderson notation A = TDT') we use gchol() from the bdsmatrix R package. Matrix T shows the "flow" of genes in pedigree.

## install.packages(pkgs="bdsmatrix")
library(package="bdsmatrix")
tmp <- gchol(as.matrix(A))
D <- diag(tmp)
(T <- as(as.matrix(tmp), "dtCMatrix"))
## 10 x 10 sparse Matrix of class "dtCMatrix"                                     
##  [1,] 1.000 .      .    .     .    .    . .   .   .
##  [2,] .     1.0000 .    .     .    .    . .   .   .
##  [3,] 0.500 0.5000 1.00 .     .    .    . .   .   .
##  [4,] .     0.5000 .    1.000 .    .    . .   .   .
##  [5,] 0.250 0.5000 0.50 0.500 1.00 .    . .   .   .
##  [6,] 0.250 0.7500 0.50 .     .    1.00 . .   .   .
##  [7,] 0.250 0.6250 0.50 0.250 0.50 0.50 1 .   .   .
##  [8,] 0.250 0.6250 0.50 0.250 0.50 0.50 . 1.0 .   .
##  [9,] .     .      .    .     .    .    . .   1.0 .
## [10,] 0.125 0.3125 0.25 0.125 0.25 0.25 . 0.5 0.5 1

## To check

L <- T %*% diag(sqrt(DInv))

L - t(U)

Now the A inverse part (inv(A) = inv(T)' inv(D) inv(T) = inv(T)' inv(S) inv(S) inv(T) using Henderson's notation, note that ). The nice thing is that pedigreemm authors provided functions to get inv(T) and D.

(TInv <- as(ped2, "sparseMatrix"))
## 10 x 10 sparse Matrix of class "dtCMatrix" (unitriangular)            
## 1   1.0  .    .    .    .    .   .  .    .   .
## 2   .    1.0  .    .    .    .   .  .    .   .
## 3  -0.5 -0.5  1.0  .    .    .   .  .    .   .
## 4   .   -0.5  .    1.0  .    .   .  .    .   .
## 5   .    .   -0.5 -0.5  1.0  .   .  .    .   .
## 6   .   -0.5 -0.5  .    .    1.0 .  .    .   .
## 7   .    .    .    .   -0.5 -0.5 1  .    .   .
## 8   .    .    .    .   -0.5 -0.5 .  1.0  .   .
## 9   .    .    .    .    .    .   .  .    1.0 .
## 10  .    .    .    .    .    .   . -0.5 -0.5 1
round(DInv <- Diagonal(x=1/Dmat(ped2)), digits=2)
## 10 x 10 diagonal matrix of class "ddiMatrix"
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
##  [1,]    1    .    .    .    .    .    .    .    .     .
##  [2,]    .    1    .    .    .    .    .    .    .     .
##  [3,]    .    .    2    .    .    .    .    .    .     .
##  [4,]    .    .    . 1.33    .    .    .    .    .     .
##  [5,]    .    .    .    .    2    .    .    .    .     .
##  [6,]    .    .    .    .    .    2    .    .    .     .
##  [7,]    .    .    .    .    .    . 2.46    .    .     .
##  [8,]    .    .    .    .    .    .    . 2.46    .     .
##  [9,]    .    .    .    .    .    .    .    .    1     .
## [10,]    .    .    .    .    .    .    .    .    .  2.33
 
round(t(TInv) %*% DInv %*% TInv, digits=2)
## 10 x 10 sparse Matrix of class "dgCMatrix"
## ...
round(crossprod(sqrt(DInv) %*% TInv), digits=2)

## 10 x 10 sparse Matrix of class "dsCMatrix"

##  [1,]  1.5  0.50 -1.0  .     .     .     .     .     .     .   
##  [2,]  0.5  2.33 -0.5 -0.67  .    -1.00  .     .     .     .   
##  [3,] -1.0 -0.50  3.0  0.50 -1.00 -1.00  .     .     .     .   
##  [4,]  .   -0.67  0.5  1.83 -1.00  .     .     .     .     .   
##  [5,]  .    .    -1.0 -1.00  3.23  1.23 -1.23 -1.23  .     .   
##  [6,]  .   -1.00 -1.0  .     1.23  3.23 -1.23 -1.23  .     .   
##  [7,]  .    .     .    .    -1.23 -1.23  2.46  .     .     .   
##  [8,]  .    .     .    .    -1.23 -1.23  .     3.04  0.58 -1.16
##  [9,]  .    .     .    .     .     .     .     0.58  1.58 -1.16
## [10,]  .    .     .    .     .     .     .    -1.16 -1.16  2.33

## To check

solve(A) - crossprod(sqrt(DInv) %*% TInv)

The second method (using crossprod) is preferred as it leads directly to symmetric matrix (dsCMatrix), which stores only upper or lower triangle. And make sure you do not do crossprod(TInv %*% sqrt(DInv)) as it is the wrong order of matrices.

As promised we will display (plot) pedigree by use of conversion functions of matrix objects to graph objects using the following code. Two examples are provided using the graph and igraph packages. The former does a very good job on this example, but otherwise igraph seems to have much nicer support for editing etc.

## source("http://www.bioconductor.org/biocLite.R")
## biocLite(pkgs=c("graph", "Rgraphviz"))
library(package="graph")
library(package="Rgraphviz")
g <- as(t(TInv), "graph")
plot(g)

## install.packages(pkgs="igraph")
library(package="igraph")
i <- igraph.from.graphNEL(graphNEL=g)
V(i)$label <- 1:10
plot(i, layout=layout.kamada.kawai)
## tkplot(i)

Kermode Bear Genetics Graphic From National Geographic Magazine

noreply@blogger.com (Gregor Gorjanc) — Mon, 01 Aug 2011 19:27:00 +0000

Kermode bear is a black bear from British Columbia, which is not black, but white as demonstrated in this photo (linked from National Geographic (NG) story about these bears).

The story also explains the genetic mechanism behind this coloring. It is a type of autosomal recessive gene, where white color is seen only in animals carrying two copies of the white allele (gene type). NG story tried to demonstrate this with the following graphic (again linked from the National Geographic (NG) story about these bears). However, they did one pretty annoying flaw. Why did not they draw two circles in each bear representing two gene copies each animal carries? This would be the most natural way to boost this graphics! They should have used two circles with black and white color. Bears would be white only if they would carry two white circles and black otherwise. In addition it would be nice to add also a mating between two animals carrying each two black genes.

Genomic selection FAQ

noreply@blogger.com (Gregor Gorjanc) — Mon, 18 Jul 2011 19:37:00 +0000

Was compiled by Ignacy Misztal and is available here.

Agilent's Integrated Biology Grant Program

noreply@blogger.com (Gregor Gorjanc) — Fri, 03 Jun 2011 12:29:00 +0000

Fostering integrated, whole-systems approaches to biological research with two $75,000US grants for open source data-integration tool development. See more here.

LaTeXsearch.com

noreply@blogger.com (Gregor Gorjanc) — Mon, 18 Apr 2011 09:54:00 +0000

Want to search over published equations? Try out LaTeXsearch.com. I wonder how useful this might be in practice?

gRain: genetics example

noreply@blogger.com (Gregor Gorjanc) — Thu, 14 Apr 2011 22:33:00 +0000

gRain is an R package for "probability propagation in graphical independence networks, also known as probabilistic expert systems (which includes Bayesian networks as a special case)". This package caught my attention because some genetic models can be seen as graphical models, which is a broader class of models encompassing Bayesian networks or sometimes also called Directed Acyclic graphs (DAGs).

A very simple example of DAG in genetics is to draw a pedigree of a small family (say father, mother, and child) with three nodes representing their genotype (g01, g02, and g03). Nodes are logically connected, father to child (g01 --> g03) and mother to child (g02 --> g03). Assume we observe some phenotype condition in a child and if we postulate that this phenotype is affected by genotype, we draw another connection between genotype and phenotype of a child (g03 --> p03).

This kind of problems/models are often postulated and one might want to do some queries about the genotype of individuals given some data (observed phenotype and/or genotype). Can we do this in R? Yes, this can be neatly done with the gRain package.

Before digging into the code we need to pause to see how can we describe the above model? If you have read at least some parts of Wikipedia links provided above than it should be straightforward to understand the following factorization of this simple model: p(g01, g02, g03, p03) = p(g01) p(g02) p(g03|g01, g02) p(p03|g03). The terms p(g01) and p(g02) are often called prior distributions (probabilities) of genotypes in some population. Term p(g03|g01, g02) describes probability of genotype in child (g03) given that we know the genotype of parents (g01 and g02). This is actually very simple part and often referred as transmission probability If there are only two alleles (gene variants) in population (A and a) and father is AA and mother is AA then child must also be AA or p(g03|g01=AA, g02=AA) = (1, 0, 0), where number show probability for genotype AA, Aa, and aa. The last term p(p03|g03) is often called penetrance function and tells us the conditional probability of some phenotype given the genotype data. Let us assume that we are interested in some binary condition, which is under very very large influence of genotype, but expressed primarily in individuals with genotype aa. Bellow is R code for specifying such a model with some additional details skipped in the text.

## Set of genotype values
geno <- c("AA", "Aa", "aa")

## Prior allele frequencies in (founder) population
PrA <- 2/3
(Pra <- 1 - PrA)

## Prior genotype frequencies in (founder) population - according to Hardy-Weinberg's law
(gP <- c(PrA*PrA, 2*PrA*Pra, Pra*Pra))

## Genotype frequencies in children given the genotype of parents - according to Mendel's law 
gM_AA_AA <- c(  1,   0,   0)
gM_AA_Aa <- c(1/2, 1/2,   0)
gM_AA_aa <- c(  0,   1,   0)

gM_Aa_AA <- c(1/2, 1/2,   0)
gM_Aa_Aa <- c(1/4, 1/2, 1/4)
gM_Aa_aa <- c(  0, 1/2, 1/2)

gM_aa_AA <- c(  0,   1,   0)
gM_aa_Aa <- c(  0, 1/2, 1/2)
gM_aa_aa <- c(  0,   0,   1)

gM <- c(gM_AA_AA, gM_AA_Aa, gM_AA_aa,
gM_Aa_AA, gM_Aa_Aa, gM_Aa_aa,
gM_aa_AA, gM_aa_Aa, gM_aa_aa)

## Set of phenotype values
phen <- c("OK", "NOK")

## Error rate (= 1 - penetrance)
e <- 0.01

## Phenotype frequencies in individual given the genotype of individual - penetrance
##      OK  NOK
p <- c(1-e, 0+e, ## AA
1-e, 0+e, ## Aa
0+e, 1-e) ## aa

## Graph - joint distribution
g01 <- cptable(~g01,             values=gP, levels=geno)
g02 <- cptable(~g02,             values=gP, levels=geno)
g03 <- cptable(~g03 + g01 + g02, values=gM, levels=geno)
p03 <- cptable(~p03 + g03,       values=p,  levels=phen)

## Compile model
plist <- compileCPT(list(g01, g02, g03, p03))
g <- grain(plist)
summary(g)

## Plot the model
iplot(g)

Now if we want to query our model without any observed data, but just using parameter values for this population we use the following code:

## Query model given parameter values (=prior)
querygrain(g, nodes=c("g01", "g02", "g03"))

which returns probabilities shown bellow. Not really interesting, but we should not expect much as we did not observe any data. Genotype probabilities are equal for all individuals and they are completely defined by population parameters.

$g01
g01
   AA        Aa        aa
0.4444444 0.4444444 0.1111111

$g02
g02
   AA        Aa        aa
0.4444444 0.4444444 0.1111111

$g03
g03
   AA        Aa        aa
0.4444444 0.4444444 0.1111111

Now we can add some data and repeat the query again.

g <- setFinding(g, nodes="p03", states="NOK")
summary(g)
## Query model given parameter values (=prior) + some data
querygrain(g, nodes=c("g01", "g02", "g03"))

Results show probability of such finding (that phenotype condition was observed) given the defined model and parameters. Probability is not really high. Further on we get updated genotype probabilities and we can see the change due to introduction of phenotype information. Now we can claim with high confidence that child has genotype aa, while there is still quite some ambiguity about genotypes of parents as there are several possible combinations of parent genotypes that would give aa genotype in child (Aa x Aa, Aa x aa, aa x Aa, and aa x aa).

Independence network: Compiled: TRUE Propagated: TRUE
Nodes : chr [1:4] "g01" "g02" "g03" "p03"
Number of cliques:                 2
Maximal clique size:               3
Maximal state space in cliques:   27
Finding:
   variable state
[1,] p03      NOK
Pr(Finding)= 0.1188889

$g01
g01
      AA         Aa         aa
0.03738318 0.64797508 0.31464174

$g02
g02
      AA         Aa         aa
0.03738318 0.64797508 0.31464174

$g03
g03
      AA         Aa         aa
0.03738318 0.03738318 0.92523364

This is all for the purpose of this blog post, but it should be clear that it is very easy to extend our model either with model individuals (just defined them with cptable and compile in joint model) or with other stuff (more alleles leads to more genotypes, other penetrance function, sex-linked genes, ...).

Redish color of screen for evening or night work

noreply@blogger.com (Gregor Gorjanc) — Fri, 01 Apr 2011 07:23:00 +0000

While reading local computer magazine (Monitor) I find out about f.lux, which is program that changes the color of screen to more redish for evening or night work. They claim that this protects our eyes and sleep pattern. I like the idea. Linux support has some problems, but their is alternative that works with Linux - redshift (see here and here for installation).

Different view options for my blog

noreply@blogger.com (Gregor Gorjanc) — Fri, 01 Apr 2011 07:11:00 +0000

Blogger now offers several dynamic views for public blogs. Click on the following links to see them in action for my blog:

Flipcard (quite good overview of posts - reminds me of the memory game where cards are flipped to find a pair)
Mosaic (similar as Flipcard, but cards have different sizes)
Sidebar (more or less like standard view, but with bar on the left side)
Snapshot (shows only preview-snapshots of posts with images stored in Picasa and linked with Blogger)
Timeslide (a bit confused layout by time)

In my view these are cool views, but the help page is much nicer due to use of graphics. If you do not have much graphics in your posts you will not be that much impressed.

Plot large networks

noreply@blogger.com (Gregor Gorjanc) — Wed, 30 Mar 2011 21:45:00 +0000

After fiding out Gephi I came across hive plots, which seem to be genuine way to plot large networks with information gain in mind. The tool is called linnet, which probably relates to "linearizing the network" and was developed in the group of Martin Krzywinski. Bellow is a plot from their webpage, while you can find more detailes in talks and posters that are linked on their website. Cool stuff indeed!!! They also have a tool called circo for circular display of relations among many nodes, something like correlations.

Exploring graphs with Gephi

noreply@blogger.com (Gregor Gorjanc) — Thu, 17 Mar 2011 10:01:00 +0000

This tool (Gephi) looks awesome. See it in action at Diseasome!

Nice simple notes on running R in parallel by Geyer

noreply@blogger.com (Gregor Gorjanc) — Tue, 15 Mar 2011 22:32:00 +0000

Here, by Charles J. Geyer.

AB Chapman Lectures

noreply@blogger.com (Gregor Gorjanc) — Tue, 15 Mar 2011 09:02:00 +0000

There is a nice set of videos on animal breeding and genetics topics from the AB Chapman Lectures series available at the University of Wisconsin site. This is now also included in the ABG-Hub site.

Prevod: Bayesian statistics in Bayesian updating

noreply@blogger.com (Gregor Gorjanc) — Fri, 04 Mar 2011 00:53:00 +0000

Prevajanje strokovnih izrazov iz tujega jezika (primarno angleščine) v slovenski jezik je zadeva, ki je ne maram. Najbrž iz razloga, da nisem ravno ne vem kako "dober v slovenščini". Poleg tega je prevajanje zamudna reč in vsak ima svojo varianto prevoda. Pred nekaj leti sem ob druženju s študenti statistike predstavil "Bayesian statistics". Takrat smo se ubadali glede prevoda in prišli do sklepa, da je "Bayesovska statistika" primeren prevod in ne "Bayesova statistika". Danes sem bil ponovno vključen v podobno "debato", kjer smo prevedli "Bayesian updating" v "Bayesovsko posodabljanje". Pomen tega izraza tiči v tem, da lahko pri Bayesovski statistiki rezultat (posteriorna porazdelitev) iz ene analize uporabimo kot predhodno znanje/vedenje (apriorno porazdelitev) pri drugi/novi analizi. Ima mogoče kdo kakšen boljši predlog prevoda?

ABG-HUB - Animal Breeding and Genetics Hub

noreply@blogger.com (Gregor Gorjanc) — Thu, 03 Mar 2011 10:56:00 +0000

Gabor Meszaros sent out today this message

Dear colleagues,

I would like to officially announce the website “Animal Breeding and Genetics Hub” we created with Gregor Gorjanc (Slovenia). The site supposed to serve as a place holding links to free animal breeding relatedsites throughout the Web.

At this point we have links to some online course materials,conference proceedings, links to ABG software sites, mailing lists, biodiversity and agricultural development resources. We also make use of the Web 2.0features such as RSS feeds, which show you the links to papers from various scientific journals in real time on one page meaning that if anew paper comesout it shows up in the RSS feed). Also there are some pictures and videos fromsome organizations, and more…

The link to the website is: http://www.netvibes.com/ABG-Hub#General

We tried to include most of our links to have some contentright at the start. But it is likely that we missed some good resources/journals out there. We would like to ask you that if you know about something that is not yet on the site, please let us know via the email: ABGHub@gmail.com

Also we would be happy to get any kind of feedback on ourwork.

onlinecourses.net

noreply@blogger.com (Gregor Gorjanc) — Thu, 17 Feb 2011 19:12:00 +0000

The website http://www.onlinecourses.net/ offers a nice index to online courses. The section with free courses is particularly rich.

Claim: I was asked to spread the word about this website and I believe it is worth doing it;)

Use of INLA for animal breeding applications

noreply@blogger.com (Gregor Gorjanc) — Fri, 04 Feb 2011 20:35:00 +0000

I am working with Ingelin Steinsland and Sara Martino on the use of approximate Bayesian method called Integrated Nested Laplace Approximation (INLA). We presented some of our preliminary results at Workshop in Bayesian Inference for Latent Gaussian Models with Applications in Zurich (Switzerland). We had two contributions:

Converting strsplit() output to a data.frame

noreply@blogger.com (Gregor Gorjanc) — Fri, 28 Jan 2011 17:44:00 +0000

R has a nice set of utilities to work with strings. Function paste is surely one among these. It can be used to "glue" several strings with optional separator. The following example shows how paste can be used to create a new variable in a dataset:

dat <- data.frame(x=1:5, y=letters[1:5])

(dat$z <- with(dat, paste(x, y, sep="-")))

Today I was in a situation where I only had column z and wanted to reverse the action of paste. Is there a way to do it? Not directly (AFAIK), but strsplit seems to be quite useful for this:

(tmp <- strsplit(x=dat$z, split="-"))

However, the output of strsplit is a list object with elements (vectors) by the elements of my column z and not by split components. Consequently one can not convert strsplit output easily back to a data.frame as you can test yourself with:

as.data.frame(tmp)

Argh. I understand that strsplit is meant to be very general (say we could have unequal number of components in one element, e.g., c("1-a-0", "1-a")), but its output is really inconvenient for transformation to a data.frame. I came up with the following solution, which seems to work nicely and is quite fast.

tmp <- unlist(strsplit(dat$z, split="-"))
cols <- c("x2", "y2")
nC <- length(cols)
ind <- seq(from=1, by=nC, length=nrow(dat))
for(i in 1:nC) {
  dat[, cols[i]] <- tmp[ind + i - 1]
}

Does anyone have a better (more obvious) solution?