RowNARemove<-function(data,missratio=0.1){
  # remove the rows whose missing ratio is higher than 10%
  threshold<-(missratio)*dim(data)[2]
  NaRaw<-which(apply(data,1,function(x) sum(is.na(x))>threshold))
  zero<-which(apply(data,1,function(x) all(x==0))==T)
  NaRAW<-c(NaRaw,zero)
  if(length(NaRAW)>0){
    data1<-data[-NaRAW,]
  }else{
    data1<-data;
  }
  data1
}   

require(beeswarm)
setwd("/home/sguo/livercancer/methylation")
load("data.RData")
data<-RowNARemove(data)    # time-consuming 
data[,4:dim(data)[2]]<-matrix(rnorm(length(data[,4:dim(data)[2]]),0.001,0.001),dim(data)[1],(dim(data)[2]-3))
# obtain x,y index
output<-matrix(NA,dim(data)[1],3)
for(i in 1:dim(data)[1]){
  a<-shapiro.test(as.numeric(data[i,seq(4,100,by=2)]))$p.value
  b<-shapiro.test(as.numeric(data[i,seq(5,101,by=2)]))$p.value
  if(any(c(a<0.05,b<0.05)==T)){
    tmp1<-try(wilcox.test(as.numeric(data[i,seq(4,100,by=2)]),as.numeric(data[i,seq(5,101,by=2)],paired=T)))
    output[i,1]<-tmp1$p.value
    output[i,2]<-tmp1$statistic
    output[i,3]<-"wilcox"
    if(tmp1$p.value<0.05/dim(data)[1]){
      filename=paste(data[i,3],data[i,2],"pdf",sep=".")
      pdf(filename)
      type=rep(c("cancer","normal"),49)
      beeswarm(as.numeric(data[i,4:101])~ type, data = data, method = 'swarm',pch = 16,xlab = '', ylab = 'Methylation Level (Beta)',labels = c('Cancer', 'Control'),main=paste(data[i,3]," (",data[i,2],")",sep=" "))
      boxplot( as.numeric(data[i,4:101])~ type, data = data, add = T,names = c("",""), col="#0000ff22") 
      dev.off() 
    }
  }else{
    tmp1<-try(t.test(as.numeric(data[i,seq(4,100,by=2)]),as.numeric(data[i,seq(5,101,by=2)],paired=T)))
    output[i,1]<-tmp1$p.value
    output[i,2]<-tmp1$statistic
    output[i,3]<-"ttest"
    if(tmp1$p.value<0.05/dim(data)[1]){
      filename=paste(data[i,3],data[i,2],"pdf",sep=".")
      pdf(filename)
      type=rep(c("cancer","normal"),49)
      beeswarm(as.numeric(data[i,4:101])~ type, data = data, method = 'swarm',pch = 16,xlab = '', ylab = 'Methylation Level',labels = c('Cancer', 'Control'),main=paste(data[i,3]," (",data[i,2],")",sep=" "))
      boxplot( as.numeric(data[i,4:101])~ type, data = data, add = T,names = c("",""), col="#0000ff22") 
      dev.off() 
    }
  }
  print(i)
}

out<-cbind(data[,1:3],output)
write.table(out,file="MethylationTest.txt",sep="\t",quote=F,col.names=NA,row.names=T)