load("Hapmap/my.data.selected.imputed.RData")
load("Hapmap/my.gini.perm.RData")
load("Hapmap/my.perm.perm.RData")

nlevelsSNPs<-sapply(my.data.selected.imputed[,-1],FUN=function(x) nlevels(factor(x)))


##############################
# select  3-categorical SNPs #
##############################

my.snps3<-my.data.selected.imputed[,1+which(nlevelsSNPs==3)]
my.gini.perm3<-my.gini.perm[which(nlevelsSNPs==3),]
my.perm.perm3<-my.perm.perm[which(nlevelsSNPs==3),]

######################
# estimate their MAF #
######################

maf3<-function(x)
{
tab<-table(x)
return((tab[2]+2*min(tab[1],tab[3]))/(2*sum(tab)))
}

my.maf3<-sapply(my.snps3,FUN=maf3)

#################
# draw boxplots #
#################

pdf("Figure3.pdf",width=12,height=6)
par(mfrow=c(1,2),cex.main=1.5,cex.axis=1.2,cex.lab=1.2,las=1)

# Boxplot Gini
boxplot(as.numeric(my.gini.perm3[my.maf3<0.1,]),as.numeric(my.gini.perm3[my.maf3>=0.1&my.maf3<0.2,]),as.numeric(my.gini.perm3[my.maf3>=0.2&my.maf3<0.3,]),as.numeric(my.gini.perm3[my.maf3>=0.3& my.maf3<0.4,]),as.numeric(my.gini.perm3[my.maf3>=0.4,]),ylim=c(0,0.1),main="randomForest Gini VIM\n HapMap data",ylab="VIM",xlab="MAF",names=c("0-0.1","0.1-0.2","0.2-0.3","0.3-0.4","0.4-0.5"))

# Boxplot perm
boxplot(as.numeric(my.perm.perm3[my.maf3<0.1,]),as.numeric(my.perm.perm3[my.maf3>=0.1&my.maf3<0.2,]),as.numeric(my.perm.perm3[my.maf3>=0.2&my.maf3<0.3,]),as.numeric(my.perm.perm3[my.maf3>=0.3& my.maf3<0.4,]),as.numeric(my.perm.perm3[my.maf3>=0.4,]),ylim=c(-5e-5,1.5e-4),main="randomForest permutation VIM\n HapMap data",ylab="",xlab="MAF",names=c("0-0.1","0.1-0.2","0.2-0.3","0.3-0.4","0.4-0.5"))

dev.off()

