install.packages("gdata")
library("gdata")
packageDescription("gdata")
source("http://bioconductor.org/biocLite.R")
biocLite("ShortRead")
int=20
b <- rep(0, 2001)
b[1] = 3.4
i = 1
while (i <= int) {
bv = rnorm(1, 0, 30)
index = sample(seq(2:2001), 1)
b[index] = bv
i = i + 1
}
ppa2<-runif(2000)
nsb=200
re_b=b
ppa_start=ppa2
testPi=1
mcmc=2000
BayesInREx<-function(nsb,re_b,ppa_start,testPi, mcmc)
{
nsubjects=nsb
nmarkers = 2000
numiter = mcmc
pi = testPi
vara = 1/20
logPi = log(pi)
logPiComp = log(1 - pi)
mean2pq = 0.5
nua = 4
cat("Simulation of Bayes\n")
data<-matrix(sample(c(0,1,2),nsubjects*nmarkers,replace=T),nrow=nsubjects,ncol = nmarkers , byrow = TRUE)
nrecords = dim(data)[1]
startMarker = 1
x = cbind(1, data[, startMarker:nmarkers])
b=re_b
y = x %*% b + rnorm(nsubjects)
oldy = y
oldb = b
storePi = array(0, numiter)
nmarkers = nmarkers - startMarker + 1
varEffects = vara/(nmarkers * (1 - pi) * mean2pq)
scalec = varEffects * (nua - 2)/nua
meanb = b
b[1] = mean(y)
var = array(0, nmarkers)
ppa = ppa_start
piMean = 0
meanVar = 0
ycorr = y - x[, 1] * b[1]
acf_mu <- array()
bf <- array(0, dim = c(mcmc, 2000))
b_effect <- array(0, dim = c(mcmc, 2001))
for (iter in 1:numiter) {
vare = (t(ycorr) %*% ycorr)/rchisq(1, nrecords + 3)
ycorr = ycorr + x[, 1] * b[1]
rhs = sum(ycorr)/vare
invLhs = 1/(nrecords/vare)
mean = rhs * invLhs
b[1] = rnorm(1, mean, sqrt(invLhs))
ycorr = ycorr - x[, 1] * b[1]
meanb[1] = meanb[1] + b[1]
acf_mu[iter] = b[1]
nLoci = 0
b_effect[iter, 1] = b[1]
for (locus in 1:nmarkers) {
ycorr = ycorr + x[, 1 + locus] * b[1 + locus]
rhs = t(x[, 1 + locus]) %*% ycorr
xpx = t(x[, 1 + locus]) %*% x[, 1 + locus]
v0 = xpx * vare
v1 = (xpx^2 * varEffects + xpx * vare)
logDelta0 = -0.5 * (log(v0) + rhs^2/v0) + logPi
logDelta1 = -0.5 * (log(v1) + rhs^2/v1) + logPiComp
probDelta1 = 1/(1 + exp(logDelta0 - logDelta1))
bf[iter, locus] = probDelta1
u = runif(1)
if (u < probDelta1) {
nLoci = nLoci + 1
lhs = xpx/vare + 1/varEffects
invLhs = 1/lhs
mean = invLhs * rhs/vare
b[1 + locus] = rnorm(1, mean, sqrt(invLhs))
ycorr = ycorr - x[, 1 + locus] * b[1 + locus]
meanb[1 + locus] = meanb[1 + locus] + b[1 + locus]
ppa[locus] = ppa[locus] + 1
var[locus] = varEffects
}
else {
b[1 + locus] = 0
var[locus] = 0
}
b_effect[iter, locus] = b[1 + locus]
}
if (iter%%100 == 0)
cat("iteration ", iter, " number of loci in model = ",
nLoci, "\n")
countLoci = 0
sum = 0
for (locus in 1:nmarkers) {
if (var[locus] > 0) {
countLoci = countLoci + 1
sum = sum + b[1 + locus]^2
}
}
cat(countLoci, nLoci, "\n")
varEffects = (scalec * nua + countLoci * sum)/rchisq(1,
nua + countLoci)
meanVar = meanVar + varEffects
aa = nmarkers - countLoci + 1
bb = countLoci + 1
pi = rbeta(1, aa, bb)
storePi[iter] = pi
piMean = piMean + pi
scalec = (nua - 2)/nua * vara/((1 - pi) * nmarkers *
mean2pq)
logPi = log(pi)
logPiComp = log(1 - pi)
}
meanb = 0.5*meanb/numiter
ppa = ppa/numiter
piMean = piMean/numiter
meanVar = meanVar/numiter
aHat_tr = x %*% meanb
results <- list()
results[[1]] <- meanb
results[[2]] <- ppa
results[[3]] <- storePi
results[[4]] <- oldb
results[[5]] <- bf
results[[6]] <- dim(x)
results[[7]] <- y
results[[8]] <- b_effect
results[[9]] <- x
corr_realb_meanb = cor(oldb, meanb)
mse_realb_meanb = sum((oldb - meanb)^2)/(length(meanb))
cat("corr_realb_meanb=", corr_realb_meanb, "\n")
cat("mse_realb_meanb=", int, mse_realb_meanb, "\n")
corrY = cor(y, aHat_tr)
cat("corrY = ", corrY, "\n")
return(results)
}
result.ex<-BayesInREx(nsb,re_b,ppa_start,testPi,mcmc)
plot(result.ex[[4]],col="blue",type="b",xlab="Marker index",ylab="Effect")
lines(result.ex[[1]],col="red",type="b")
source("http://bioconductor.org/biocLite.R")
biocLite("gdsfmt")
biocLite("SNPRelate")
library(gdsfmt)
library(SNPRelate)
snpgdsSummary(snpgdsExampleFileName())
vcf.fn <- system.file("extdata", "sequence.vcf", package="SNPRelate")
snpgdsVCF2GDS(vcf.fn, "test.gds", method="biallelic.only")
quit("yes)
quit("yes")
first = list(a = 1, b = 2, c = 3)
second = list(a = 2, b = 3, c = 4)
first
second
mapply(c, first, second, SIMPLIFY=FALSE)
load(system.file("extdata/addTracks.RData",package="LDheatmap"))
ll<-LDheatmap(GIMAP5.CEU$snp.data,GIMAP5.CEU$snp.support$Position,flip=TRUE)
library(LDheatmap)
library(ggplot2)
require(chopsticks)
ll<-LDheatmap(GIMAP5.CEU$snp.data,GIMAP5.CEU$snp.support$Position,flip=TRUE)
ll<-LDheatmap(GIMAP5.CEU$snp.data,GIMAP5.CEU$snp.support$Position,flip=TRUE)
data(GIMAP5.CEU)
ll<-LDheatmap(GIMAP5.CEU$snp.data,GIMAP5.CEU$snp.support$Position,flip=TRUE)
grid.newpage()
grid.draw(llGenes$LDheatmapGrob)
llGenesRecomb <- LDheatmap.addRecombRate(llGenes, chr="chr7", genome="hg18")
grid.newpage()
grid.draw(llGenes$LDheatmapGrob)
atests<-runif(nrow(GIMAP5.CEU$snp.support))
names(atests)<-rownames(GIMAP5.CEU$snp.support)
atests["rs6598"]<-1e-5
llGenesRecombScatter<-LDheatmap.addScatterplot(llGenesRecomb,-log10(atests),ylab="-log10(p-values)")
posn<-GIMAP5.CEU$snp.support$Position
manhattan2<-ggplotGrob(
{
qplot(posn,-log10(atests),xlab="", xlim=range(posn),asp=1/10)
last_plot() + theme(axis.text.x=element_blank(),
axis.title.y = element_text(size = rel(0.75)))
})
llQplot<-LDheatmap.addGrob(ll,manhattan2,height=.7)
llQplot2<-LDheatmap.addGrob(ll,rectGrob(gp=gpar(col="white")),height=.2)
pushViewport(viewport(x=.48,y=.76,width=.99,height=.2))
grid.draw(manhattan2)
popViewport(1)
dev.off()
llImage<-LDheatmap.addGrob(ll,rasterGrob(GIMAP5ideo))
names(ll$LDheatmapGrob$children)
names(llGenesRecombScatter$LDheatmapGrob$children)
names(llQplot$LDheatmapGrob$children)
names(llImage$LDheatmapGrob$children)
setwd("/Users/guillaume/Desktop/Rqtl_DataSets/rqtl_Carotenoids")
install.packages("qtl")
install.packages("devtools")
library(qtl)
library(qtlbim)
pop11 = read.cross(format='csv', file='ASCarotFamily091014.csv')
pop11 <- jittermap(pop11)
??sim.geno
pop11 <- sim.geno(pop11)
pop11
summary(pop11)
nind(pop11)
nphe(pop11)
nchr(pop11)
nmar(pop11)
totmar(pop11)
plot(pop11)
plot(pop11)
??calc.genoprob
smout.hk1 <- scanone(pop11, pheno.col=1, method="hk")
smout.hk2 <- scanone(pop11, pheno.col=2, method="hk")
smout.hk6 <- scanone(pop11, pheno.col=6, method="hk")
pop11 <- calc.genoprob(pop11, step=2.5, error.prob=0.01)
smout.hk1 <- scanone(pop11, pheno.col=1, method="hk")
smout.hk2 <- scanone(pop11, pheno.col=2, method="hk")
smout.hk6 <- scanone(pop11, pheno.col=6, method="hk")
plot(smout.em, smout.hk2, lty=1, col=c("blue","red"))
smout.em <- scanone(pop11)
plot(smout.em, smout.hk2, lty=1, col=c("blue","red"))
