##Data structure
##RILineIndicator marker1-markerM trait1 ....traitK 

genopheno.dat<- read.table("sampleRI.dat",header=T)

##Marker Regression analysis
scanoneMR<- function(dat=genopheno.dat, n.marker = 1, 
		   pheno.col=1, n.perm=-1) {

       ##sort the data first to make the later permutation easier 
       genopheno.dat<- genopheno.dat[order(genopheno.dat[,1]),]
       
       ##check number of replicates within each line
       xx<- table(genopheno.dat[,1])
       n.line<- length(xx)
       nobs.line<- as.vector(xx)
       
       pheno.true.col<- n.marker+ pheno.col +1

      ##error message of the phenotype column is out the data range
      if(pheno.true.col> length(genopheno.dat[1,]))
         warning("Phenotype is not available!")

      else{
        ##analyze raw data
        geno.dat<- genopheno.dat[,(2:(n.marker+1))]
        trait.dat<- genopheno.dat[,pheno.true.col]
        lrt<- getFstat(geno.dat, trait.dat)

        ##get the reduced line genotype
        geno.line<- geno.dat[c(1, cumsum(nobs.line[-n.line])+1),]

    
       ##do permutation
       if(n.perm>0){
         k<- 0
         permuted.geno<- geno.dat
 
        
        while(k<n.perm){
           permuteLine<- sample(1:n.line)

           permuted.geno[1:nobs.line[1],]<- geno.line[permuteLine[1],]

           for(j in 2:n.line){
             permuted.geno[(sum(nobs.line[1:(j-1)])+1):sum(nobs.line[1:j]),]<- 
               geno.line[permuteLine[j],]  
            }

           lrt<- data.frame(lrt,  getFstat(permuted.geno, trait.dat))

           k<- k+1     
         }      
       }

      return(lrt)
     }
  }

## Function to call regression analysis and return Fstatistics
getFstat<- function(geno=geno.dat, ytrait = trait.dat){
    nMarker<- length(geno[1,])
    
    Fstat<- rep(-1, nMarker)
    
    for(i in 1:n.marker){
	   flm<- lm(ytrait~as.factor(geno[,i])) 
           Fstat[i]<- anova(flm)[1,4]
    }

 return(Fstat)
}


lrt.all<- scanoneMR(genopheno.dat,10,1,10)

##permutation based P-value
permutedP<- sum(apply(lrt.all[,-1],1, max)>= max(lrt.all[,1]))/length(lrt.all[1,])