library(corpcor)
library(nlme)
library(mgcv)
library(sva)
bdata<-read.table("G:/科研/methy-evo/bdata_new.txt",row.names=1)
colnames(bdata)=c("sample","gse","class")
con<-file("G:/gene_avebetas_final.txt","r")
title=readLines(con,1)  #read the title line
f1=readLines(con,1)#read the first data line
k<-0
row_num=500 #assign the row number of the data matrix
number=floor(21231/row_num)#the number of matrix with row_num x 8675 elements
while(length(f1)!=0){
	line=strsplit(f1,"\t");
	if(k==0){                  #when k=0 the first data line 
		mdata<-matrix(nrow=row_num,ncol=8675)
		for(j in (2:length(line[[1]]))){
       mdata[1,j-1]<-as.numeric(line[[1]][j])
    }	#the first line of the matrix	
	}	
	else{
		if(k<number*row_num){
			if(k %% row_num != 0){
				i<-(k %% row_num)+1
			  for(j in (2:length(line[[1]]))){			  	
           mdata[i,j-1]<-as.numeric(line[[1]][j])
        }			
			}
			else{
				mod = model.matrix(~as.factor(class), data=bdata)#create a model according to phenotype
				combat_mdata = ComBat(dat=mdata, batch=bdata$gse, mod=mod, numCovs=NULL, par.prior=TRUE,prior.plots=FALSE)#adjust batch effect 
				write.table(combat_mdata,"G:/科研/methy-evo/combat_mdata.txt",quote=F,sep="\t",row.names = F,col.names = F,append=T)#write the adjusted data to a output file
			  mdata<-matrix(nrow=row_num,ncol=8675)#initial the data matrix
	      for(j in (2:length(line[[1]]))){
           mdata[1,j-1]<-as.numeric(line[[1]][j])
        } 
			}			
		}
		if(k>=number*row_num){
			 if(k==number*row_num){
			 	  mod = model.matrix(~as.factor(class), data=bdata)
			 	  combat_mdata = ComBat(dat=mdata, batch=bdata$gse, mod=mod, numCovs=NULL, par.prior=TRUE,prior.plots=FALSE)
			 	  write.table(combat_mdata,"G:/科研/methy-evo/combat_mdata.txt",quote=F,sep="\t",row.names = F,col.names = F,append=T)
			 	  mdata<-matrix(nrow=21231-number*row_num,ncol=8675)#initial the data matrix
			 	  for(j in (2:length(line[[1]]))){
             mdata[1,j-1]<-as.numeric(line[[1]][j])
          } 
       }
       if(k>number*row_num){
       	  i<-(k %% row_num)+1
			    for(j in (2:length(line[[1]]))){			    	
             mdata[i,j-1]<-as.numeric(line[[1]][j])
          }		
          if(k==21230){
       	     mod = model.matrix(~as.factor(class), data=bdata)#create a model according to phenotype
       	     combat_mdata = ComBat(dat=mdata, batch=bdata$gse, mod=mod, numCovs=NULL, par.prior=TRUE,prior.plots=FALSE)#adjust batch effect 
       	     write.table(combat_mdata,"G:/科研/methy-evo/combat_mdata.txt",quote=F,sep="\t",row.names = F,col.names = F,append=T)#write the adjusted data to a output file
          }   
       }    
		}
	}
  f1=readLines(con,n=1) 
  k=k+1 
}
close(con)