#Let's load library SNPassoc
library(SNPassoc)
#get the data example:
#both data.frames SNPs and SNPs.info.pos are loaded typing data(SNPs)
data(SNPs)
#look at the data (only first four SNPs)
SNPs[1:10,1:9]
table(SNPs[,2])
mySNP<-snp(SNPs$snp10001,sep="")
mySNP
summary(mySNP)
# left figure
plot(mySNP,label="snp10001",col="darkgreen")
# right figure
plot(mySNP,type=pie,label="snp10001",col=c("darkgreen","yellow","red"))
reorder(mySNP,ref="minor")
gg<-c("het","hom1","hom1","hom1","hom1","hom1","het","het","het","hom1","hom2","hom1","hom2")
snp(gg,name.genotypes=c("hom1","het","hom2"))
myData<-setupSNP(data=SNPs,colSNPs=6:40,sep="")
myData.o<-setupSNP(SNPs, colSNPs=6:40, sort=TRUE,info=SNPs.info.pos, sep="")
labels(myData)
summary(myData)
plot(myData,which=20)
plotMissing(myData)
res<-tableHWE(myData)
res
res<- tableHWE(myData,strata=myData$sex)
res

data(HapMap)
myDat.HapMap<-setupSNP(HapMap, colSNPs=3:9307, sort = TRUE,info=HapMap.SNPs.pos, sep="")
resHapMap<-WGassociation(group, data=myDat.HapMap, model="log-add")
plot(resHapMap, whole=FALSE, print.label.SNPs = FALSE)
summary(resHapMap)
plot(resHapMap, whole=TRUE, print.label.SNPs = FALSE)
resHapMap.scan<-scanWGassociation(group, data=myDat.HapMap, model="log-add")
resHapMap.perm<-scanWGassociation(group, data=myDat.HapMap,model="log-add", nperm=1000)
res.perm<- permTest(resHapMap.perm)
print(resHapMap.scan[1:5,])
print(resHapMap.perm[1:5,])
perms <- attr(resHapMap.perm, "pvalPerm")
print(res.perm)
plot(res.perm)
res.perm.rtp<- permTest(resHapMap.perm,method="rtp",K=20)
print(res.perm.rtp)

getSignificantSNPs(resHapMap,chromosome=5)
association(casco~snp(snp10001,sep=""), data=SNPs)
myData<-setupSNP(data=SNPs,colSNPs=6:40,sep="")
association(casco~snp10001, data=myData)
association(casco~snp10001, data=myData, model=c("cod","log"))
association(casco~sex+snp10001+blood.pre, data=myData)
association(casco~snp10001+blood.pre+strata(sex), data=myData)
association(casco~snp10001+blood.pre, data=myData,subset=sex=="Male")
association(log(protein)~snp100029+blood.pre+strata(sex), data=myData)
ans<-association(log(protein)~snp10001*sex+blood.pre, data=myData,model="codominant")
print(ans,dig=2)
ans<-association(log(protein)~snp10001*factor(recessive(snp100019))+blood.pre, data=myData, model="codominant")
print(ans,dig=2)
sigSNPs<-getSignificantSNPs(resHapMap,chromosome=5,sig=5e-8)$column
myDat2<-setupSNP(HapMap, colSNPs=sigSNPs, sep="")
resHapMap2<-WGassociation(group~1, data=myDat2)
plot(resHapMap2,cex=0.8)

myData<-setupSNP(SNPs, colSNPs=6:40, sep="")
myData.o<-setupSNP(SNPs, colSNPs=6:40, sort=TRUE,info=SNPs.info.pos, sep="")
ans<-WGassociation(protein~1,data=myData.o)
library(Hmisc)
SNP<-pvalues(ans)
out<-latex(SNP,file="c:/temp/ans1.tex", where="'h",caption="Summary of case-control study for SNPs data set.",center="centering", longtable=TRUE, na.blank=TRUE,
size="scriptsize", collabel.just=c("c"), lines.page=50,rownamesTexCmd="bfseries")
WGstats(ans,dig=5)
plot(ans)
Bonferroni.sig(ans, model="log-add", alpha=0.05,include.all.SNPs=FALSE)
pvalAdd<-additive(resHapMap)
pval<-pval[!is.na(pval)]
library(qvalue)
qobj<-qvalue(pval)
max(qobj$qvalues[qobj$pvalues <= 0.001])
procs<-c("Bonferroni","Holm","Hochberg","SidakSS","SidakSD","BH","BY")
res2<-mt.rawp2adjp(rawp,procs)
mt.reject(cbind(res$rawp,res$adjp),seq(0,0.1,0.001))$r

datSNP<-setupSNP(SNPs,6:40,sep="")
tag.SNPs<-c("snp100019", "snp10001", "snp100029")
geno<-make.geno(datSNP,tag.SNPs)
mod<-haplo.glm(log(protein)~geno,data=SNPs,family=gaussian,locus.label=tag.SNPs,allele.lev=attributes(geno)$unique.alleles,
control = haplo.glm.control(haplo.freq.min=0.05))
mod
intervals(mod)
ansCod<-interactionPval(log(protein)~sex, data=myData.o,model="codominant")
plot(ansCod)