LOD.QTLdetect=function(LOD,chro,thre=3,gape=20){
  nt=nrow(LOD)
  ns=ncol(LOD)
  nc=nrow(chro)
  
  cr0=c()
  for(i in 1:nc){
    cr0=c(cr0,rep(i,chro[i,2]))
  }
  
  lcr=chro[,2]
  ncr=c()
  for(i in 1:nc){
    ncr[i]=sum(chro[1:i,2])
  }
  colnames(chro)=c("chromosome","# of bins")
  
  LOD[is.na(LOD)]=0
  LOD[LOD<thre]=0
  
  det=matrix(0,nt,ns)
  cat("step","\t", "process", "\n")
  for(j in 1:nt){
    if(j%%50==0 | j==nt){
      cat("detect", paste(j, nt, sep = "/"), "\n", sep = "\t")
    }
    k1=c()
    for(i in 1:nc){
      x=LOD[j,cr0==i]
      l=lcr2[i]
      k0=rep(0,l)
      for(m in 1:l){
        x0=x[m]
        if(x0==0){det[m]=0
        } else { 
          if(m%in%c(1:gape)){
            if(x0<max(x[1:(m+gape)])){
              k0[m]=NA
            } else {
              k0[m]=1
            }
          } else if(m%in%c((l-gape+1):l)){
            if(x0<max(x[(m-gape):l])){
              k0[m]=NA
            } else {
              k0[m]=1
            }
          } else {
            if(x0<max(x[(m-gape):(m+gape)])){
              k0[m]=NA
            } else {
              k0[m]=1
            }
          }
        }
      }
      k1=c(k1,k0)
    }
    det[j,]=k1
  }
  
  summ=function(x){sum(x,na.rm=T)}
  qtldetect=apply(det,1,summ)
  
  link=matrix(0,nt,nc)
  cat("step","\t", "process", "\n")
  for(j in 1:nt){
    if(j%%50==0 | j==nt){
      cat("linkage", paste(j, nt, sep = "/"), "\n", sep = "\t")
    }
    for(i in 1:nc){
      x=det[j,cr0==i]
      a=table(x)
      if(length(a)==2){link[j,i]=a[2]}
    }
  }
  link=table(link)
  
  EQF=matrix(0,nt,ns)
  cat("step","\t", "process", "\n")
  for(i in 1:nt){
    if(i%%50==0 | i==nt){
      cat("linkage", paste(i, nt, sep = "/"), "\n", sep = "\t")
    }
    QTL=det[i,]
    if(!1%in%QTL){next}
    lod=LOD[i,]
    ka=c()
    for(j in 1:nc){
      QTLcr=QTL[cr0==j]
      n=length(QTLcr)
      k0=rep(0,n)
      if(!1%in%QTLcr){
        ka=c(ka,k0)
        next
      }
      QTLlo=which(QTLcr==1 & !is.na(QTLcr))
      QTLci=lod[cr0==j]
      for(k in 1:length(QTLlo)){
        pii=QTLlo[k]
        QTLci1=rep(-1,length(QTLci))
        QTLci1[QTLci>(QTLci[pii]-1) & QTLci>thre]=1
        QTLci2=QTLci1*c(-1,(QTLci1[-n]))
        QTLci3=which(QTLci2==-1)
        QTLci4=sort(c(QTLci3,pii))
        QTLci5=QTLci4-c(QTLci3,(n+1))
        ci0=min(which(QTLci5<0))
        ci1=QTLci3[ci0-1]
        ci2=n
        if(ci0<=length(QTLci3)){ci2=QTLci3[ci0]-1}
        sdi=((ci2-ci1)/(2*1.96))
        k1=pnorm((1:n)+0.5,pii,sdi)-pnorm((1:n)-0.5,pii,sdi)
        k0=k0+k1
      }
      ka=c(ka,k0)
    }
    EQF[i,]=ka
  }
  
  EQF.all=apply(EQF,2,sum)
  
  plot(EQF.all,type = "h",ylab="EQF",xlab="chromosome",main=paste("LOD thresholds =",thre),xaxt="n",ylim=c(0,max(EQF.all)*1.2))
  axis(1,1:nc,at=ncr-(lcr)/2)
  abline(v=c(0,ncr),col="green")

  return(list(detect.QTL.number=qtldetect,QTL.matrix=det,EQF.matrix=EQF,linkage.QTL.number=link,threshole=thre,chro=chro))
}

EQF.permu=function(LOD.QTLdetect.result,ptime=1000,alpha=0.05,Q=TRUE,plot.all=TRUE,plot.cr=TRUE){
  dat=LOD.QTLdetect.result
  detect=dat$QTL.matrix
  EQF=dat$EQF.matrix
  nt=nrow(detect)
  ns=ncol(detect)
  thre=dat$threshole
  chro=dat$chro
  nc=nrow(chro)
  lcr=chro[,2]
  ncr=c()
  for(i in 1:nc){
    ncr[i]=sum(chro[1:i,2])
  }
  cr0=c()
  for(i in 1:nc){
    cr0=c(cr0,rep(i,chro[i,2]))
  }
  
  clu.det=detect
  clu.det[is.na(clu.det)]=0
  rownames(clu.det)=1:nt
  det.all=apply(clu.det,1,sum)
  clu.det=clu.det[det.all>0,]
  
  clu.eqf=EQF
  clu.eqf[is.na(clu.det)]=0
  rownames(clu.eqf)=1:nt
  clu.eqf=clu.eqf[det.all>0,]
  
  qtlind=as.numeric(rownames(clu.det))
  cluster=list()
  
  maxeqf=1
  add.group=c()
  cat("cluster group", "\n")
  while(class(clu.eqf)[[1]]=="matrix"){
    peqf=c()
    clu.eqf1=apply(clu.eqf,2,sum)
    peqf=c(peqf,which.max(clu.eqf1))
    group=which(clu.det[,peqf]==1)
    if(length(group)==0){
      clu.eqf[,peqf]=0
      next
    }
    add=1
    while(add>0){
      group1=group
      group.det=clu.det[group1,]
      if(length(group1)>1){
        clu.det1=apply(group.det,2,sum)
      } else {clu.det1=group.det}
      group.pos=which(clu.det1>0)
      group.pos=group.pos[!group.pos%in%peqf]
      peqf=c(peqf,group.pos)
      if(length(group.pos)>0){
        for(i in 1:length(group.pos)){
          group2=which(clu.det[,group.pos[i]]==1)
          group1=c(group1,group2[!group2%in%group1])
        }
      }
      add=length(group1)-length(group)
      group=group1
    }
    ind=sort(as.numeric(rownames(clu.det)[group]))
    cluster[[maxeqf]]=ind
    add.group=c(add.group,ind)
    clu.det=clu.det[-group,]
    clu.eqf=clu.eqf[-group,]
    if(maxeqf%%10==0){
      cat(maxeqf, "\n")
    }
    maxeqf=maxeqf+1
  }
  cat(maxeqf, "\n")
  
  cluster[[maxeqf]]=qtlind[!qtlind%in%add.group]
  Lagend=c()
  for(i in 1:maxeqf){
    Lagend[i+2]=length(cluster[[i]])
  }
  Lagend[1]=maxeqf
  Lagend[2]=sum(Lagend[3:(maxeqf+2)])
  names(Lagend)=c("ngroup","nQTL",(paste("ng",1:maxeqf,sep = "")))
  
  clumatrix=matrix(0,maxeqf,ns)
  for(i in 1:maxeqf){
    submatrix=EQF[cluster[[i]],]
    if(class(submatrix)[[1]]=="matrix"){
      clumatrix[i,]=apply(submatrix,2,sum)
    } else {clumatrix[i,]=submatrix}
  }
  
  eqf.premutation=function(eqfmatrix,ptime=ptime){
    cat("permutation time", "\n")
    n=ncol(eqfmatrix)
    g=nrow(eqfmatrix)
    out=matrix(0,ptime,n)
    for(i in 1:ptime){
      p.matrix=matrix(0,g,n)
      for(j in 1:g){
        eqfpos=eqfmatrix[j,eqfmatrix[j,]>0]
        k=sample(1:n,length(eqfpos))
        p.matrix[j,k]=eqfpos
      }
      out[i,]=sort(apply(p.matrix,2,sum),decreasing = T)
      if(i%%50==0 | i==ptime){
        cat(paste(i, ptime, sep = "/"), "\n", sep = "\t")
      }
    }
    return(out)
  }
  
  cat("cluster group permutation", "\n")
  permu.clu=eqf.premutation(clumatrix,ptime)

  sortEQF=function(x,a=ceiling(ptime*alpha)){
    return(sort(x,decreasing = T)[a])
  }
  
  thre.clu=apply(permu.clu,2,sortEQF)

  a1=thre.clu[1]
  a2=thre.clu[2]
  a5=thre.clu[5]
  a20=thre.clu[20]
  eqfthre=c(a20,a5,a2,a1)
  names(eqfthre)=c(20,5,2,1)

  permu.q=NULL
  if(Q){
    cat("Q-method permutation", "\n")
    q0=EQF[apply(EQF,1,sum)>0,]
    permu.q=eqf.premutation(q0,ptime)
    thre.q=apply(permu.q,2,sortEQF)
    b1=thre.q[1]
    beta=max(which(b1<thre.clu))
    eqfthre=c(b1,a20,a5,a2,a1)
    names(eqfthre)=c(paste("β:",beta),20,5,2,1)
  }
  
  eqfthre=round(eqfthre,1)
  
  eqf.all=apply(EQF,2,sum)
  nseqf=c()
  for(k in 1:length(eqfthre)){
    nseqf[k]=length(which(eqf.all>eqfthre[k]))
  }
  
  if(plot.all){
    par(mai=c(1,1,1,2))
    plot(eqf.all,type = "h",ylab="EQF",xlab="chromosome",main=paste("LOD thresholds =", thre,"  # of group =",nrow(clumatrix)),
         xaxt="n",ylim=c(0,max(eqf.all)*1.2),yaxt="n",cex.main=1.5)
    axis(1,1:nc,at=ncr-(lcr)/2)
    for(j in 1:length(eqfthre)){
      axis(2,eqfthre[j],las=2)
      axis(4,paste(names(eqfthre[j])," (",nseqf[j],")",sep=""),at=eqfthre[j],las=2)
    }
    abline(v=c(0,ncr),col="green")
    abline(h=eqfthre,col="red")
  }
  
  if(plot.cr){
    par(mfrow=c(2,2))
    for(i in 1:16){
      par(mai=c(0.5,0.5,0.5,1))
      eqf=eqf.all[cr0==i]
      plot(eqf,type = "h",ylab="EQF",main=paste("cr",i,"   LOD thresholds = ",thre),
           ylim=c(0,max(eqf.all)*1.2),yaxt="n",cex.main=1.2)
      for(k in 1:length(eqfthre)){
        axis(2,eqfthre[k],las=2)
        axis(4,paste(names(eqfthre)[k]," (",nseqf[k],")",sep=""),at=as.numeric(eqfthre[k]),las=2)
      }
      abline(h=eqfthre,col="red")
    }
    par(mfrow=c(1,1))
  }

  return(list(cluster.number=Lagend,cluster.id=cluster,cluster.matrix=clumatrix,permu.matrix.cluster=permu.clu,permu.matrix.Q=permu.q))
}