#args<-commandArgs(TRUE)
#chrid<-toString(args[1]) 

library(nlme)
dir <- '/home/hum/LCLHiCNov18/110918_FIREQTL_search/ann/' # FIRE score

for(chrid in 1:22)
{
F    <- read.table( paste(dir, 'testable_YRI11com_FIREscore_chr',chrid,'_110918.txt',sep=''), head=T)
Frep <- read.table( paste(dir, 'testable_YRI11rep_FIREscore_chr',chrid,'_110918.txt',sep=''), head=T)
G    <- read.table( paste(dir, 'testable_FIRE_SNP_chr',chrid,'_110918.txt',sep=''), head=T)

final<-matrix(-1, nrow=nrow(G), ncol=5)
for(i in 1:nrow(G))
{
# combined data
v0<-as.numeric( G[i,-c(1:2)] )
u0<-as.numeric( F[ F[,2]+1<=G[i,2] & G[i,2]<=F[,3] , -c(1:3) ] )
fitlm <- lm( u0 ~ v0 ) #simple linear regression
final[i,1] <- summary(fitlm)$coefficients[2,1]    
final[i,2] <- summary(fitlm)$coefficients[2,4] 

# two biological replicates
v0rep<-rep(v0, each=2)
u0rep<-as.numeric( Frep[ Frep[,2]+1<=G[i,2] & G[i,2]<=Frep[,3] , -c(1:3) ] )
subject <- rep( seq(1,11,1), each=2 )

fitlmm <- gls( u0rep ~ v0rep, correlation=corAR1(form=~1|subject), method="REML") #linear mixed effect model
cslmm <- as.data.frame(summary(fitlmm)$tTable)  
 
final[i,3] <- cslmm[2,1]    
final[i,4] <- cslmm[2,4]   
final[i,5] <- coef(fitlmm$modelStruct$corStruct, uncons = FALSE, allCoef = TRUE)

 if(i%%1000==0)
    {
     print(c(chrid, i, round(i/nrow(G),2)))
    }
}

final<-as.data.frame(final)
out<-cbind(G[,1:2], final)
colnames(out)<-c('chr', 'pos', 'lm.est', 'lm.p', 'lmm.est', 'lmm.p', 'lmm.phi')

write.table(out, file=paste('FIREQTL_SLR_LMM_pvalue_chr', chrid, '_110918.txt', sep=''), 
row.names=FALSE, col.names=TRUE, sep='\t', quote=FALSE)
}