Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- fastAssoc <- function(y, x)
- {
- index <- is.finite(y) & is.finite(x)
- n <- sum(index)
- y <- y[index]
- x <- x[index]
- vx <- var(x)
- bhat <- cov(y, x) / vx
- ahat <- mean(y) - bhat * mean(x)
- # fitted <- ahat + x * bhat
- # residuals <- y - fitted
- # SSR <- sum((residuals - mean(residuals))^2)
- # SSF <- sum((fitted - mean(fitted))^2)
- rsq <- (bhat * vx)^2 / (vx * var(y))
- fval <- rsq * (n-2) / (1-rsq)
- tval <- sqrt(fval)
- se <- abs(bhat / tval)
- # Fval <- (SSF) / (SSR/(n-2))
- # pval <- pf(Fval, 1, n-2, lowe=F)
- p <- pf(fval, 1, n-2, lowe=F)
- return(list(
- ahat=ahat, bhat=bhat, se=se, fval=fval, pval=p, n=n
- ))
- }
- getFittedVals <- function(y, x)
- {
- n <- length(x)
- bhat <- cov(y, x) / var(x)
- ahat <- mean(y) - bhat * mean(x)
- fitted <- ahat + x * bhat
- return(fitted)
- }
- runGwas <- function(x, y)
- {
- require(dplyr)
- y <- as.matrix(y)
- x <- as.matrix(x)
- ntrait <- ncol(y)
- nsnp <- ncol(x)
- l <- list()
- for(i in 1:ntrait)
- {
- message(i)
- res <- array(0, c(nsnp, 7))
- res[,1] <- 1:nsnp
- for(j in 1:nsnp)
- {
- res[j,2:7] <- unlist(fastAssoc(y[,i], x[,j]))
- }
- l[[i]] <- data.frame(res)
- names(l[[i]]) <- c("snp", "a", "b", "se", "fval", "pval", "n")
- l[[i]]$trait <- paste("trait", i)
- }
- return(bind_rows(l))
- }
- makeGeno <- function(nid, nsnp, maf=0.5)
- {
- matrix(rbinom(nid * nsnp, 2, maf), nid, nsnp)
- }
- makePhen <- function(effs, indep, vy=1, vx=rep(1, length(effs)))
- {
- if(is.null(dim(indep))) indep <- cbind(indep)
- stopifnot(ncol(indep) == length(effs))
- stopifnot(length(vx) == length(effs))
- cors <- effs * vx / sqrt(vx) / sqrt(vy)
- stopifnot(sum(cors^2) <= 1)
- cors <- c(cors, sqrt(1-sum(cors^2)))
- indep <- t(t(scale(cbind(indep, rnorm(nrow(indep))))) * cors * c(vx, 1))
- y <- drop(scale(rowSums(indep)) * sqrt(vy))
- return(y)
- }
- spreadGwas <- function(gwas, col)
- {
- require(tidyr)
- l <- list()
- stopifnot(all(col %in% names(gwas)))
- for(i in 1:length(col))
- {
- x <- select(gwas, snp, trait, get(col[i]))
- l[[i]] <- spread_(x, "trait", col[i])
- }
- if(length(col)==1)
- {
- return(l[[i]])
- } else {
- names(l) <- col
- return(l)
- }
- }
- chooseEffects <- function(nsnp, totvar, sqrt=TRUE)
- {
- eff <- rnorm(nsnp)
- aeff <- abs(eff)
- sc <- sum(aeff) / totvar
- out <- eff / sc
- if(sqrt)
- {
- out <- sqrt(abs(out)) * sign(out)
- }
- return(out)
- }
- # g <- makeGeno(1000, 20)
- # u <- makePhen(chooseEffects(20, 0.3), g)
- # x <- makePhen(sqrt(0.6), u)
- # y <- makePhen(sqrt(0.2), u)
- # gwas <- runGwas(g, cbind(x,y))
- # out <- spreadGwas(gwas, c("b"))
- # plot(out[,3] ~ out[,2])
- # g <- makeGeno(10000, 20)
- # x <- makePhen(chooseEffects(20, 0.3), g)
- # u <- rnorm(10000)
- # y <- makePhen(sqrt(c(0.4, 0.02)), cbind(x,u))
- # gwas <- runGwas(g, cbind(x,y))
- # out <- spreadGwas(gwas, c("b"))
- # plot(out[,3] ~ out[,2])
- # gwas
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement