# Exercise 1 N <- 10 sim.genotypes <- matrix(NA, nrow = N, ncol = 2) for( i in 1:N){ sim.genotypes[i,] <- sample(nucleotides, 2, replace = TRUE, prob = c(0.3,0.7)) } min(as.matrix(table(sim.genotypes))) / sum(as.matrix(table(sim.genotypes))) # Exercise 2 Xd <- 1 - 2*abs(Xa) # Exercise 3 X.mx <- cbind(1, Xa,Xd) beta.vec <- c(0.3,1.5,0.05) phenotypes <- X.mx %*% beta.vec + rnorm(10,0,1) # Exercise 4 MLE.beta <- solve(t(X.mx) %*% X.mx) %*% t(X.mx) %*% phenotypes y.hat <- X.mx %*% MLE.beta SSM <- sum((y.hat - mean(phenotypes))^2) SSE <- sum((phenotypes - y.hat)^2) df.M <- 2 df.E <- 10 - 3 MSM <- SSM / df.M MSE <- SSE / df.E Fstatistic <- MSM / MSE # to check if it is correct pf(Fstatistic, 2, 7,lower.tail = FALSE) summary(lm(phenotypes ~ X.mx)) # Exercise 5 N <- 200 n.geno <- 50 sim.genotypes <- matrix(NA, nrow = N, ncol = 2 * n.geno) for( i in 1:n.geno){ nucleotides <- sample(c("a","t","g","c"), 2,replace = FALSE) maf <- runif(1,min = 0.1,max =0.5) major <- 1 - maf for (j in 1:N){ sim.genotypes[j,c(2*i-1,2*i)] <- sample(nucleotides, 2, replace = TRUE, prob = c(maf,major)) } } Xa <- matrix(NA, nrow = N, ncol = n.geno) for ( i in 1: n.geno){ minor.allele <- names(which.min(table(sim.genotypes[,c(2*i -1, 2*i)]))) temp <- 1 * (sim.genotypes[,c(2*i -1, 2*i)] == minor.allele) Xa[,i] <- temp[,1] + temp[,2] - 1 } Xd <- 1 - 2*abs(Xa) # Lets make genotype # 20 the significant one X.mx <- cbind(1, Xa[,20],Xd[,20]) beta.vec <- c(0.3,0.8,0.05) phenotypes <- X.mx %*% beta.vec + rnorm(10,0,1) pval <- vector(length = n.geno) for (i in 1 : n.geno){ N.samples <- dim(Xa)[1] X.mx <- cbind(1,Xa[,i],Xd[,i]) MLE.beta <- solve(t(X.mx) %*% X.mx) %*% t(X.mx) %*% phenotypes y.hat <- X.mx %*% MLE.beta SSM <- sum((y.hat - mean(phenotypes))^2) SSE <- sum((phenotypes - y.hat)^2) df.M <- 2 df.E <- N.samples - 3 MSM <- SSM / df.M MSE <- SSE / df.E Fstatistic <- MSM / MSE # to check if it is correct pval[i] <- pf(Fstatistic, df.M, df.E,lower.tail = FALSE) summary(lm(phenotypes ~ X.mx)) } plot(-log10(pval), col = "blue", main = "Example manhattan plot")