#############
# Read data #
#############

utah<-read.table("utah.txt",sep=" ")
nigeria<-read.table("nigeria.txt",sep=" ")

###############
# Common SNPs #
###############
snps<-intersect(utah[,1],nigeria[,1])

utah2<-utah[is.element(utah[,1],snps),]
utah2<-utah2[order(utah2[,1]),]

nigeria2<-nigeria[is.element(nigeria[,1],snps),]
nigeria2<-nigeria2[order(nigeria2[,1]),]

#####################
# Define data frame #
#####################

my.snps<-rbind(t(nigeria2[,-c(1,2)]),t(utah2[,-c(1,2)]))
names(my.snps)<-nigeria2[,1]


n<-nrow(my.data)
p<-ncol(my.data[,-1])

####################
# Replace NN by NA #
####################

replaceNN<-function(x)
{
NN<-which(x=="NN")
if (length(NN)>0) 
 {
 x[x=="NN"]<-NA
 }
return(x)
}

my.snps<-apply(my.snps,FUN=replaceNN,MARGIN=2)

my.data<-data.frame(y=factor(c(rep(0,(ncol(nigeria2)-2)),rep(1,(ncol(utah2)-2)))),my.snps)

rm(nigeria,nigeria2,utah,utah2)

#########################################################
# Select SNPs with >1 levels and less than 30% missings #
#########################################################

nlevelsSNPs<-sapply(my.data[,-1],FUN=function(x) nlevels(factor(x)))
propNA0<-sapply(my.data[my.data[,1]==0,-1],FUN=function(x) sum(is.na(x)))/sum(my.data[,1]==0)
propNA1<-sapply(my.data[my.data[,1]==1,-1],FUN=function(x) sum(is.na(x)))/sum(my.data[,1]==1)

selected<-which((nlevelsSNPs>1) & (propNA0<0.3) & (propNA1<0.3))
my.data.selected<-my.data[,c(1,selected+1)]

###########################
# Imputing missing values #
###########################

library(randomForest)
my.data.selected.imputed<-rfImpute(y=my.data.selected[,1],x=my.data.selected[,-1])
save(my.data.selected.imputed,file="my.data.selected.imputed.RData")

