#Let's load library SNPassoc
library(SNPassoc)

#get the data example
data(SNPs)

#data description
?SNPs


#look at the data (only first five SNPs)
SNPs[1:10,1:9]

#the snp function
mySNP<-snp(SNPs$snp10001,sep="")

# print, summary
print(mySNP)
summary(mySNP)

# plot
plot(mySNP,label="snp10001",col="darkgreen")
plot(mySNP,type=pie,label="snp10001",col=c("darkgreen","yellow","red"))

# reorder 
reorder(mySNP,"minor")

# other codes
gg<-c("het","hom1","hom1","hom1","hom1","hom1","het","het","het",
  "hom1","hom2","hom1","hom2")
snp(gg,name.genotypes=c("hom1","het","hom2"))

# setup SNPs
myData<-setupSNP(data=SNPs,colSNPs=6:40,sep="")

# setup SNPs with genomic order
myData.o<-setupSNP(SNPs,6:40,sort=TRUE,info=SNPs.info.pos,sep="")


# summary information for an object of class setupSNP
summary(myData)


#plot a given SNP
plot(myData,which=20)

# Missing information in SNPs
plotMissing(myData)

# the same with SNPs ordered
plotMissing(myData.o)
     


#H-W equilibrium
?tableHWE
res<-tableHWE(myData)
res
print(res,sig=0.0001)

res<-tableHWE(myData,strata=myData$sex)
res



#Whole genome
data(HapMap) 
ptm <- proc.time()
myDat<-setupSNP(HapMap, colSNPs=3:9307, sort = TRUE,
   info=HapMap.SNPs.pos, sep="")
proc.time() - ptm


args(WGassociation)
?WGassociation
resHapMap<-WGassociation(group~1, data=myDat, model="log")

print(resHapMap)
WGstats(resHapMap)
plot(resHapMap)

# By default if there are more than 10 chrs this kind of plot.
# To change the plot, just whole=FALSE

plot(resHapMap, whole=FALSE, print.label.SNPs = FALSE)


# Short summary of this plot
summary(resHapMap)


# Faster analysis, only p values
resHapMap<-scanWGassociation(group~1, data=myDat, model="log")



# Identify SNPs statistically significant
getSignificantSNPs(resHapMap,chromosome=5)

sigSNPs<-getSignificantSNPs(resHapMap,chrom=5,sig=5e-8)$column


# Permutation test
resHapMap.perm<-scanWGassociation(group, data=myDat, nperm=1000, model="log-add")
res.perm<- permTest(resHapMap.perm)
plot(res.perm)
print(res.perm)

# rank truncated product
res.perm.rtp<- permTest(resHapMap.perm,method="rtp",K=20)
print(res.perm.rtp)


# Analysis with those SNPs that are statistically significant
myDat2<-setupSNP(HapMap, colSNPs=sigSNPs, sep="")
resHapMap2<-scanWGassociation(group~1, data=myDat2)
plot(resHapMap2,cex=0.8)


# association study. case-control
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 study. case-control adjusted by sex and blood.pre
association(casco~sex+snp10001+blood.pre, data=myData)

#stratified analysis
association(casco~snp10001+blood.pre+strata(sex), data=myData)

#subset analysis
association(casco~snp10001+blood.pre, data=myData, 
      subset=sex=="Male")


#quantitative trait. Adjusted by blood.pre and stratified by sex
association(log(protein)~snp100029+blood.pre+strata(sex), data=myData)

#interaction
ans<-association(log(protein)~snp10001*sex+blood.pre, data=myData, model="codominant")
print(ans,dig=2)

ans<-association(log(protein)~snp10001*recessive(snp100019)+blood.pre, data=myData, model="codominant")


# more than one SNP
myData.o<-setupSNP(SNPs,6:40,sort=TRUE,info=SNPs.info.pos,sep="")
ans<-WGassociation(protein~1,data=myData.o)
ans
print(ans)
WGstats(ans,dig=5)
plot(ans)

ans<-WGassociation(protein~1,data=myData.o,model=c("log","co"))
plot(ans)

# export to LaTeX (needs Hmisc)
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")



#significant SNPs after Bonferroni correction
Bonferroni.sig(ans, model="log-add", alpha=0.05, include.all.SNPs=FALSE)


# Get p values

pvalAdd<-additive(resHapMap)
pval<-pvalAdd[!is.na(pvalAdd)]


# FDR
library(qvalue)
qobj<-qvalue(pval)
max(qobj$qvalues[qobj$pvalues <= 0.001])*100


# Multiple test
library(multtest)
procs<-c("Bonferroni","Holm","Hochberg","SidakSS","SidakSD","BH","BY")
res<-mt.rawp2adjp(pval,procs)

qobj<-qvalue(pval)
max(qobj$qvalues[qobj$pvalues <= 0.001])*100




# Some plots


# Haplotype analysis

# an object of class 'setupSNP' is required
datSNP<-setupSNP(SNPs,6:40,sep="")

# sigle analysis
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)


# adjusted analysis
mod<-haplo.glm(log(protein)~geno+sex,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)




# interaction
ans<-haplo.interaction(log(protein)~int(sex),data=datSNP,
    SNPs.sel=c("snp100019","snp10001","snp100029"))
ans


ans<-haplo.interaction(casco~int(sex),data=datSNP,
    SNPs.sel=c("snp100019","snp10001","snp100029"))
ans


# adjusted by blood.pre
ans<-haplo.interaction(log(protein)~blood.pre+int(sex),data=datSNP,
    SNPs.sel=c("snp100019","snp10001","snp100029"))
ans


# SNPs interactions
ansCod<-interactionPval(log(protein)~sex, data=myData.o, 
          model="codominant")
print(ansCod)
plot(ansCod)