yeast quantitative genetics cross study

 

---

title: "yeast power study"

author: "H Qin"

date: "8/31/2021"

output:

  pdf_document: default

  html_document: default

---

```{r simulate genotypes}

rm(list=ls())

N = 150

nuc_means = rpois(10, lambda=10)

mit_means = rpois(15, lambda=10)

summary(nuc_means)

summary(mit_means)

b0= 0

b1= 1 # mito influence on phenotype

b2= 1 # nuclear influence on phenotype

b3 =0.2  # mit X nuc interaction influence on phenotype, p << 0.001

b3 =0.1  # mit X nuc interaction influence on phenotype, p=0.049

#b3 = 0.05 # p = 0.3

```

```{r simulate-phenotype}

debug = 0

phenotype_mit_nuc = function(b0, b1, b2, b3, mit_single_mean, nuc_single_mean, debug){

  y = b0 + b1*mit_single_mean + b2*nuc_single_mean + b3*mit_single_mean * nuc_single_mean

  if (debug > 0) {

    print( paste("pmn:: mit_single_mean =", mit_single_mean, "nuc_single_mean", nuc_single_mean) )

  }

  return (y)

}

nuc_genotypes = sample(1:10, N, replace=TRUE)

mit_genotypes = sample(1:15, N, replace=TRUE)

y = 1:N

for ( i in 1:N ){

  #print(paste("i:", i, "mit_genotypes[i]",mit_genotypes[i] ))

  y[i] = phenotype_mit_nuc(b0, b1, b2, b3, mit_means[mit_genotypes[i]], nuc_means[nuc_genotypes[i]], debug=0) + rnorm(1)

}  

tb = data.frame( cbind( y, mit_genotypes, nuc_genotypes)) 

tb$mit_genotypes = factor( tb$mit_genotypes)

tb$nuc_genotypes = factor( tb$nuc_genotypes)

summary(tb)

```

```{r}

library(nlme);

m1a = glm(y ~ mit_genotypes  , tb, family='gaussian');

m2 = glm(y ~ mit_genotypes + nuc_genotypes , tb, family='gaussian');

m3 =  glm( y ~ mit_genotypes + nuc_genotypes + mit_genotypes:nuc_genotypes, data=tb)

```

```{r}

#summary(m1a)

```

```{r}

anova( m1a, m2, test='F')

```

```{r}

summary(m2)

summary(m3)

anova(m2, m3, test='F')

```