R/pLARmEB.R

In mrMLM: Multi-Locus Random-SNP-Effect Mixed Linear Model Tools for GWAS

pLARmEB<-function(gen,phe,outATCG,genRaw,kk,psmatrix,CriLOD,lars1,Genformat,Bootstrap,CLO){
  
    lodvalue<-CriLOD
    gene.data<-gen
    
    rm(gen)
    gc()
    
    inputform<-Genformat
    
    if(is.null(psmatrix)){
      flagps<-1
    }else{
      flagps<-0
    }
    
    if(is.null(lodvalue)==TRUE||is.null(lars1)==TRUE){
      warning("Please set parameter!")
    }
    if(lodvalue<0)
    {
      warning("Please input critical LOD score: > 0 !")
    }
    if(lars1<0||lars1>=nrow(phe))
    {
      warning("Please input the number of most relevant variables select by LARS: >0 and less than numbers of sample!")
    }
    if(is.null(gene.data)==TRUE)
    {
      warning("Please input correct genotypic data !")
      
    }
    if(is.null(phe)==TRUE)
    {
      warning("Please input correct phenotypic data !")
    }
    if((is.null(gene.data)==FALSE)&&(is.null(phe)==FALSE)&&(ncol(gene.data)!=(nrow(phe)+2)))
    {
      warning("Sample size in genotypic dataset doesn't equal to the sample size in phenotypic dataset !")
    }
    
    if((is.null(gene.data)==FALSE)&&(is.null(phe)==FALSE)&&((ncol(gene.data)==(nrow(phe)+2)))&&(lodvalue>=0)&&(lars1>0))
    {
      
      wan<-NULL
      result<-NULL
      
      multinormal<-function(y,mean,sigma)
      {
        pdf_value<-(1/sqrt(2*3.14159265358979323846*sigma))*exp(-(y-mean)*(y-mean)/(2*sigma));
        return (pdf_value)
      }
      
      ebayes_EM<-function(x,z,y)
      {
        n<-nrow(z);k<-ncol(z)
        
        if(abs(min(eigen(crossprod(x,x))$values))<1e-6){
          b<-solve(crossprod(x,x)+diag(ncol(x))*1e-8)%*%crossprod(x,y)
        }else{
          b<-solve(crossprod(x,x))%*%(crossprod(x,y))
        }
        
        v0<-as.numeric(crossprod((y-x%*%b),(y-x%*%b))/n)
        u<-matrix(rep(0,k),k,1)
        v<-matrix(rep(0,k),k,1)
        s<-matrix(rep(0,k),k,1)
        for(i in 1:k)
        {
          zz<-z[,i]
          s[i]<-((crossprod(zz,zz)+1e-100)^(-1))*v0
          u[i]<-s[i]*crossprod(zz,(y-x%*%b))/v0
          v[i]<-u[i]^2+s[i]
        }
        
        vv<-matrix(rep(0,n*n),n,n);
        for(i in 1:k)
        {
          zz<-z[,i]
          vv=vv+tcrossprod(zz,zz)*v[i]
        }
        vv<-vv+diag(n)*v0
        
        iter<-0;err<-1000;iter_max<-500;err_max<-1e-8
        tau<-0;omega<-0
        while((iter<iter_max)&&(err>err_max))
        {
          iter<-iter+1
          v01<-v0
          v1<-v
          b1<-b
          vi<-solve(vv)
          xtv<-crossprod(x,vi)
          
          if(ncol(x)==1)
          {
            b<-((xtv%*%x)^(-1))*(xtv%*%y)
          }else{
            if(abs(min(eigen(xtv%*%x)$values))<1e-6){
              b<-solve((xtv%*%x)+diag(ncol(x))*1e-8)%*%(xtv%*%y)
            }else{
              b<-solve(xtv%*%x)%*%(xtv%*%y)
            }
          }
          r<-y-x%*%b
          ss<-matrix(rep(0,n),n,1)
          for(i in 1:k)
          {
            zz<-z[,i]
            zztvi<-crossprod(zz,vi)
            u[i]<-v[i]*zztvi%*%r
            s[i]<-v[i]*(1-zztvi%*%zz*v[i])
            v[i]<-(u[i]^2+s[i]+omega)/(tau+3)
            ss<-ss+zz*u[i]
          }
          v0<-as.numeric(crossprod(r,(r-ss))/n)
          
          vv<-matrix(rep(0,n*n),n,n)
          for(i in 1:k)
          {
            zz<-z[,i]
            vv<-vv+tcrossprod(zz,zz)*v[i]
          }
          vv<-vv+diag(n)*v0
          
          err<-(crossprod((b1-b),(b1-b))+(v01-v0)^2+crossprod((v1-v),(v1-v)))/(2+k)
          beta<-t(b)
          sigma2<-v0
        }
        
        wang<-matrix(rep(0,k),k,1)
        for (i in 1:k){
          stderr<-sqrt(s[i]+1e-20)
          t<-abs(u[i])/stderr
          f<-t*t
          p<-pchisq(f,1,lower.tail = F)
          wang[i]<-p
        }
        
        return(list(u=u,sigma2=sigma2,wang=wang))
      }
      
      likelihood<-function(xxn,xxx,yn,bbo)
      {
        nq<-ncol(xxx)
        ns<-nrow(yn)
        at1<-0
        
        if(is.null(bbo)==TRUE){
          ww1<-1:ncol(xxx)
          ww1<-as.matrix(ww1)
        }else{
          ww1<-as.matrix(which(abs(bbo)>1e-5))
        }
        at1<-dim(ww1)[1]
        lod<-matrix(rep(0,nq),nq,1)
        if(at1>0.5)
          ad<-cbind(xxn,xxx[,ww1])
        else
          ad<-xxn
        if(abs(min(eigen(crossprod(ad,ad))$values))<1e-6)
          bb<-solve(crossprod(ad,ad)+diag(ncol(ad))*0.01)%*%crossprod(ad,yn)
        else
          bb<-solve(crossprod(ad,ad))%*%crossprod(ad,yn)
        vv1<-as.numeric(crossprod((yn-ad%*%bb),(yn-ad%*%bb))/ns);
        ll1<-sum(log(abs(multinormal(yn,ad%*%bb,vv1))))
        sub<-1:ncol(ad);
        if(at1>0.5)
        {
          for(i in 1:at1)
          {
            ij<-which(sub!=sub[i+ncol(xxn)])
            ad1<-ad[,ij]
            if(abs(min(eigen(crossprod(ad1,ad1))$values))<1e-6)
              bb1<-solve(crossprod(ad1,ad1)+diag(ncol(ad1))*0.01)%*%crossprod(ad1,yn)
            else
              bb1<-solve(crossprod(ad1,ad1))%*%crossprod(ad1,yn) 
            vv0<-as.numeric(crossprod((yn-ad1%*%bb1),(yn-ad1%*%bb1))/ns);
            ll0<-sum(log(abs(multinormal(yn,ad1%*%bb1,vv0))))
            lod[ww1[i]]<--2.0*(ll0-ll1)/(2.0*log(10))
          }
        }
        return (lod)
      }
      
      
      emma.eigen.L <- function(Z,K,complete=TRUE) {
        if ( is.null(Z) ) {
          return(emma.eigen.L.wo.Z(K))
        }
        else {
          return(emma.eigen.L.w.Z(Z,K,complete))
        }
      }
      #likelihood
      emma.eigen.L.wo.Z <- function(K) {
        eig <- eigen(K,symmetric=TRUE)
        return(list(values=eig$values,vectors=eig$vectors))
      }
      #likelihood
      emma.eigen.L.w.Z <- function(Z,K,complete=TRUE) {
        if ( complete == FALSE ) {
          vids <- colSums(Z)>0
          Z <- Z[,vids]
          K <- K[vids,vids]
        }
        eig <- eigen(K%*%crossprod(Z,Z),symmetric=FALSE,EISPACK=TRUE)
        return(list(values=eig$values,vectors=qr.Q(qr(Z%*%eig$vectors),complete=TRUE)))
      }
      
      #restricted likelihood
      emma.eigen.R <- function(Z,K,X,complete=TRUE) {
        if ( ncol(X) == 0 ) {
          return(emma.eigen.L(Z,K))
        }
        else if ( is.null(Z) ) {
          return(emma.eigen.R.wo.Z(K,X))
        }
        else {
          return(emma.eigen.R.w.Z(Z,K,X,complete))
        }
      }
      #restricted likelihood
      emma.eigen.R.wo.Z <- function(K, X) {
        n <- nrow(X)
        q <- ncol(X)
        S <- diag(n)-X%*%solve(crossprod(X,X))%*%t(X)
        eig <- eigen(S%*%(K+diag(1,n))%*%S,symmetric=TRUE)
        stopifnot(!is.complex(eig$values))
        return(list(values=eig$values[1:(n-q)]-1,vectors=eig$vectors[,1:(n-q)]))
      }
      
      emma.eigen.R.w.Z <- function(Z, K, X, complete = TRUE) {
        if ( complete == FALSE ) {
          vids <-  colSums(Z) > 0
          Z <- Z[,vids]
          K <- K[vids,vids]
        }
        n <- nrow(Z)
        t <- ncol(Z)
        q <- ncol(X)
        
        SZ <- Z - X%*%solve(crossprod(X,X))%*%crossprod(X,Z)
        eig <- eigen(K%*%crossprod(Z,SZ),symmetric=FALSE)
        if ( is.complex(eig$values) ) {
          eig$values <- Re(eig$values)
          eig$vectors <- Re(eig$vectors)    
        }
        qr.X <- qr.Q(qr(X))
        return(list(values=eig$values[1:(t-q)],
                    vectors=qr.Q(qr(cbind(SZ%*%eig$vectors[,1:(t-q)],qr.X)),
                                 complete=TRUE)[,c(1:(t-q),(t+1):n)]))   
      }
      
      emma.delta.ML.LL.wo.Z <- function(logdelta, lambda, etas, xi) {
        n <- length(xi)
        delta <- exp(logdelta)
        return( 0.5*(n*(log(n/(2*pi))-1-log(sum((etas*etas)/(delta*lambda+1))))-sum(log(delta*xi+1))) )  
      }
      
      emma.delta.ML.LL.w.Z <- function(logdelta, lambda, etas.1, xi.1, n, etas.2.sq ) {
        delta <- exp(logdelta)
        return( 0.5*(n*(log(n/(2*pi))-1-log(sum(etas.1*etas.1/(delta*lambda+1))+etas.2.sq))-sum(log(delta*xi.1+1)) ))
        
      }
      
      emma.delta.ML.dLL.wo.Z <- function(logdelta, lambda, etas, xi) {
        n <- length(xi)
        delta <- exp(logdelta)
        etasq <- etas*etas
        ldelta <- delta*lambda+1
        return( 0.5*(n*sum(etasq*lambda/(ldelta*ldelta))/sum(etasq/ldelta)-sum(xi/(delta*xi+1))) )
      }
      
      emma.delta.ML.dLL.w.Z <- function(logdelta, lambda, etas.1, xi.1, n, etas.2.sq ) {
        delta <- exp(logdelta)
        etasq <- etas.1*etas.1
        ldelta <- delta*lambda+1
        return( 0.5*(n*sum(etasq*lambda/(ldelta*ldelta))/(sum(etasq/ldelta)+etas.2.sq)-sum(xi.1/(delta*xi.1+1))) )
      }
      
      emma.delta.REML.LL.wo.Z <- function(logdelta, lambda, etas) {
        nq <- length(etas)
        delta <-  exp(logdelta)
        return( 0.5*(nq*(log(nq/(2*pi))-1-log(sum(etas*etas/(delta*lambda+1))))-sum(log(delta*lambda+1))) )
      }
      
      emma.delta.REML.LL.w.Z <- function(logdelta, lambda, etas.1, n, t, etas.2.sq ) {
        tq <- length(etas.1)
        nq <- n - t + tq
        delta <-  exp(logdelta)
        return( 0.5*(nq*(log(nq/(2*pi))-1-log(sum(etas.1*etas.1/(delta*lambda+1))+etas.2.sq))-sum(log(delta*lambda+1))) ) 
      }
      
      emma.delta.REML.dLL.wo.Z <- function(logdelta, lambda, etas) {
        nq <- length(etas)
        delta <- exp(logdelta)
        etasq <- etas*etas
        ldelta <- delta*lambda+1
        return( 0.5*(nq*sum(etasq*lambda/(ldelta*ldelta))/sum(etasq/ldelta)-sum(lambda/ldelta)) )
      }
      
      emma.delta.REML.dLL.w.Z <- function(logdelta, lambda, etas.1, n, t1, etas.2.sq ) {
        t <- t1
        tq <- length(etas.1)
        nq <- n - t + tq
        delta <- exp(logdelta)
        etasq <- etas.1*etas.1
        ldelta <- delta*lambda+1
        return( 0.5*(nq*sum(etasq*lambda/(ldelta*ldelta))/(sum(etasq/ldelta)+etas.2.sq)-sum(lambda/ldelta) ))
      }
      
      
      emma.MLE <- function(y, X, K, Z=NULL, ngrids=100, llim=-10, ulim=10,
                           esp=1e-10, eig.L = NULL, eig.R = NULL)
      {
        n <- length(y)
        t <- nrow(K)
        q <- ncol(X)
        stopifnot(ncol(K) == t)
        stopifnot(nrow(X) == n)
        if ( det(crossprod(X,X)) == 0 ) {
          warning("X is singular")
          return (list(ML=0,delta=0,ve=0,vg=0))
        }
        
        if ( is.null(Z) ) {
          if ( is.null(eig.L) ) {
            eig.L <- emma.eigen.L.wo.Z(K)
          }
          if ( is.null(eig.R) ) {
            eig.R <- emma.eigen.R.wo.Z(K,X)
          }
          etas <- crossprod(eig.R$vectors,y)
          logdelta <- (0:ngrids)/ngrids*(ulim-llim)+llim
          m <- length(logdelta)
          delta <- exp(logdelta)
          Lambdas.1<-matrix(eig.R$values,n-q,m)    
          Lambdas <- Lambdas.1 * matrix(delta,n-q,m,byrow=TRUE)+1
          Xis.1<-matrix(eig.L$values,n,m)
          Xis <- Xis.1* matrix(delta,n,m,byrow=TRUE)+1
          Etasq <- matrix(etas*etas,n-q,m)
          dLL <- 0.5*delta*(n*colSums(Etasq*Lambdas.1/(Lambdas*Lambdas))/colSums(Etasq/Lambdas)-colSums(Xis.1/Xis))
          optlogdelta <- vector(length=0)
          optLL <- vector(length=0)
          if ( dLL[1] < esp ) {
            optlogdelta <- append(optlogdelta, llim)
            optLL <- append(optLL, emma.delta.ML.LL.wo.Z(llim,eig.R$values,etas,eig.L$values))
          }
          if ( dLL[m-1] > 0-esp ) {
            optlogdelta <- append(optlogdelta, ulim)
            optLL <- append(optLL, emma.delta.ML.LL.wo.Z(ulim,eig.R$values,etas,eig.L$values))
          }
          
          for( i in 1:(m-1) )
          {
            if ( ( dLL[i]*dLL[i+1] < 0-esp*esp ) && ( dLL[i] > 0 ) && ( dLL[i+1] < 0 ) ) 
            {
              r <- uniroot(emma.delta.ML.dLL.wo.Z, lower=logdelta[i], upper=logdelta[i+1], lambda=eig.R$values, etas=etas, xi=eig.L$values)
              optlogdelta <- append(optlogdelta, r$root)
              optLL <- append(optLL, emma.delta.ML.LL.wo.Z(r$root,eig.R$values, etas, eig.L$values))
            }
          }
        }
        else {
          if ( is.null(eig.L) ) {
            eig.L <- emma.eigen.L.w.Z(Z,K)
          }
          if ( is.null(eig.R) ) {
            eig.R <- emma.eigen.R.w.Z(Z,K,X)
          }
          etas <- crossprod(eig.R$vectors,y)
          etas.1 <- etas[1:(t-q)]
          etas.2 <- etas[(t-q+1):(n-q)]
          etas.2.sq <- sum(etas.2*etas.2)
          logdelta <- (0:ngrids)/ngrids*(ulim-llim)+llim
          m <- length(logdelta)
          delta <- exp(logdelta)
          Lambdas.1<-matrix(eig.R$values,t-q,m)
          Lambdas <- Lambdas.1 * matrix(delta,t-q,m,byrow=TRUE) + 1
          Xis.1<-matrix(eig.L$values,t,m)
          Xis <- Xis.1 * matrix(delta,t,m,byrow=TRUE) + 1
          Etasq <- matrix(etas.1*etas.1,t-q,m)
          dLL <- 0.5*delta*(n*colSums(Etasq*Lambdas.1/(Lambdas*Lambdas))/(colSums(Etasq/Lambdas)+etas.2.sq)-colSums(Xis.1/Xis))
          optlogdelta <- vector(length=0)
          optLL <- vector(length=0)
          if ( dLL[1] < esp ) {
            optlogdelta <- append(optlogdelta, llim)
            optLL <- append(optLL, emma.delta.ML.LL.w.Z(llim,eig.R$values,etas.1,eig.L$values,n,etas.2.sq))
          }
          if ( dLL[m-1] > 0-esp ) {
            optlogdelta <- append(optlogdelta, ulim)
            optLL <- append(optLL, emma.delta.ML.LL.w.Z(ulim,eig.R$values,etas.1,eig.L$values,n,etas.2.sq))
          }
          
          for( i in 1:(m-1) )
          {
            if ( ( dLL[i]*dLL[i+1] < 0-esp*esp ) && ( dLL[i] > 0 ) && ( dLL[i+1] < 0 ) ) 
            {
              r <- uniroot(emma.delta.ML.dLL.w.Z, lower=logdelta[i], upper=logdelta[i+1], lambda=eig.R$values, etas.1=etas.1, xi.1=eig.L$values, n=n, etas.2.sq = etas.2.sq )
              optlogdelta <- append(optlogdelta, r$root)
              optLL <- append(optLL, emma.delta.ML.LL.w.Z(r$root,eig.R$values, etas.1, eig.L$values, n, etas.2.sq ))
            }
          }
        }
        maxdelta <- exp(optlogdelta[which.max(optLL)])
        optLL=replaceNaN(optLL)
        maxLL <- max(optLL)
        if ( is.null(Z) ) {
          maxve <- sum(etas*etas/(maxdelta*eig.R$values+1))/n    
        }
        else {
          maxve <- (sum(etas.1*etas.1/(maxdelta*eig.R$values+1))+etas.2.sq)/n
        }
        maxvg <- maxve*maxdelta
        
        return (list(ML=maxLL,delta=maxdelta,ve=maxve,vg=maxvg))
      }
      
      emma.REMLE <- function(y, X, K, Z=NULL, ngrids=100, llim=-10, ulim=10,
                             esp=1e-10, eig.L = NULL, eig.R = NULL) {
        n <- length(y)
        t <- nrow(K)
        q <- ncol(X)
        stopifnot(ncol(K) == t)
        stopifnot(nrow(X) == n)
        if ( det(crossprod(X,X)) == 0 ) {
          warning("X is singular")
          return (list(REML=0,delta=0,ve=0,vg=0))
        }
        
        if ( is.null(Z) ) {
          if ( is.null(eig.R) ) {
            eig.R <- emma.eigen.R.wo.Z(K,X)
          }
          etas <- crossprod(eig.R$vectors,y)
          logdelta <- (0:ngrids)/ngrids*(ulim-llim)+llim
          m <- length(logdelta)
          delta <- exp(logdelta)
          Lambdas.1<-matrix(eig.R$values,n-q,m)
          Lambdas <- Lambdas.1 * matrix(delta,n-q,m,byrow=TRUE) + 1
          Etasq <- matrix(etas*etas,n-q,m)
          dLL <- 0.5*delta*((n-q)*colSums(Etasq*Lambdas.1/(Lambdas*Lambdas))/colSums(Etasq/Lambdas)-colSums(Lambdas.1/Lambdas))
          optlogdelta <- vector(length=0)
          optLL <- vector(length=0)
          if ( dLL[1] < esp ) {
            optlogdelta <- append(optlogdelta, llim)
            optLL <- append(optLL, emma.delta.REML.LL.wo.Z(llim,eig.R$values,etas))
          }
          if ( dLL[m-1] > 0-esp ) {
            optlogdelta <- append(optlogdelta, ulim)
            optLL <- append(optLL, emma.delta.REML.LL.wo.Z(ulim,eig.R$values,etas))
          }
          
          for( i in 1:(m-1) )
          {
            if ( ( dLL[i]*dLL[i+1] < 0-esp*esp ) && ( dLL[i] > 0 ) && ( dLL[i+1] < 0 ) ) 
            {
              r <- uniroot(emma.delta.REML.dLL.wo.Z, lower=logdelta[i], upper=logdelta[i+1], lambda=eig.R$values, etas=etas)
              optlogdelta <- append(optlogdelta, r$root)
              optLL <- append(optLL, emma.delta.REML.LL.wo.Z(r$root,eig.R$values, etas))
            }
          }
        }
        else {
          if ( is.null(eig.R) ) {
            eig.R <- emma.eigen.R.w.Z(Z,K,X)
          }
          etas <- crossprod(eig.R$vectors,y)
          etas.1 <- etas[1:(t-q)]
          etas.2 <- etas[(t-q+1):(n-q)]
          etas.2.sq <- sum(etas.2*etas.2)
          logdelta <- (0:ngrids)/ngrids*(ulim-llim)+llim
          m <- length(logdelta)
          delta <- exp(logdelta)
          Lambdas.1 <- matrix(eig.R$values,t-q,m) 
          Lambdas <- Lambdas.1 * matrix(delta,t-q,m,byrow=TRUE) + 1
          Etasq <- matrix(etas.1*etas.1,t-q,m)
          dLL <- 0.5*delta*((n-q)*colSums(Etasq*Lambdas.1/(Lambdas*Lambdas))/(colSums(Etasq/Lambdas)+etas.2.sq)-colSums(Lambdas.1/Lambdas))
          optlogdelta <- vector(length=0)
          optLL <- vector(length=0)
          if ( dLL[1] < esp ) {
            optlogdelta <- append(optlogdelta, llim)
            optLL <- append(optLL, emma.delta.REML.LL.w.Z(llim,eig.R$values,etas.1,n,t,etas.2.sq))
          }
          if ( dLL[m-1] > 0-esp ) {
            optlogdelta <- append(optlogdelta, ulim)
            optLL <- append(optLL, emma.delta.REML.LL.w.Z(ulim,eig.R$values,etas.1,n,t,etas.2.sq))
          }
          
          for( i in 1:(m-1) )
          {
            if ( ( dLL[i]*dLL[i+1] < 0-esp*esp ) && ( dLL[i] > 0 ) && ( dLL[i+1] < 0 ) ) 
            {
              r <- uniroot(emma.delta.REML.dLL.w.Z, lower=logdelta[i], upper=logdelta[i+1], lambda=eig.R$values, etas.1=etas.1, n=n, t1=t, etas.2.sq = etas.2.sq )
              optlogdelta <- append(optlogdelta, r$root)
              optLL <- append(optLL, emma.delta.REML.LL.w.Z(r$root,eig.R$values, etas.1, n, t, etas.2.sq ))
            }
          }
        }  
        maxdelta <- exp(optlogdelta[which.max(optLL)])
        optLL=replaceNaN(optLL)
        maxLL <- max(optLL)
        if ( is.null(Z) ) {
          maxve <- sum(etas*etas/(maxdelta*eig.R$values+1))/(n-q)    
        }
        else {
          maxve <- (sum(etas.1*etas.1/(maxdelta*eig.R$values+1))+etas.2.sq)/(n-q)
        }
        maxvg <- maxve*maxdelta
        return (list(REML=maxLL,delta=maxdelta,ve=maxve,vg=maxvg))
      }
      
      
      emma.maineffects.B<-function(Z=NULL,K,deltahat.g,complete=TRUE){
        if( is.null(Z) ){
          return(emma.maineffects.B.Zo(K,deltahat.g))
        }
        else{
          return(emma.maineffects.B.Z(Z,K,deltahat.g,complete))
        }
      }
      
      
      emma.maineffects.B.Zo <-function(K,deltahat.g){
        t <- nrow(K)
        stopifnot(ncol(K) == t)
        B<-deltahat.g*K+diag(1,t)
        eig<-eigen(B,symmetric=TRUE)
        qr.B<-qr(B)
        q<-qr.B$rank
        stopifnot(!is.complex(eig$values))
        A<-diag(1/sqrt(eig$values[1:q]))
        Q<-eig$vectors[,1:q]
        C<-Q%*%A%*%t(Q)
        return(list(mC=C,Q=Q,A=A))
      }
      
      emma.maineffects.B.Z <- function(Z,K,deltahat.g,complete=TRUE){
        if ( complete == FALSE ) {
          vids <- colSums(Z)>0
          Z <- Z[,vids]
          K <- K[vids,vids]
        }
        n <- nrow(Z)  
        B <- deltahat.g*Z%*%K%*%t(Z)+diag(1,n)
        eig <- eigen(B,symmetric=TRUE,EISPACK=TRUE)
        qr.B<-qr(B)
        q<-qr.B$rank
        stopifnot(!is.complex(eig$values))
        A<-diag(1/sqrt(eig$values[1:q]))
        Q<-eig$vectors[,1:q]
        C<-Q%*%A%*%t(Q)
        return(list(mC=C,Q=Q,A=A,complete=TRUE))
      }
      emma.MLE0.c <- function(Y_c,W_c){
        n <- length(Y_c)
        stopifnot(nrow(W_c)==n)
        M_c<-diag(1,n)-W_c%*%solve(crossprod(W_c,W_c))%*%t(W_c)
        etas<-crossprod(M_c,Y_c)
        LL <- 0.5*n*(log(n/(2*pi))-1-log(sum(etas*etas)))
        return(list(ML=LL))
      }
      
      emma.REMLE0.c <- function(Y_c,W_c){
        n <- length(Y_c)
        stopifnot(nrow(W_c)==n)
        M_c <-diag(1,n)-W_c%*%solve(crossprod(W_c,W_c))%*%t(W_c)
        eig <-eigen(M_c)
        t <-qr(W_c)$rank
        v <-n-t
        U_R <-eig$vector[,1:v]
        etas<-crossprod(U_R,Y_c)
        LL <- 0.5*v*(log(v/(2*pi))-1-log(sum(etas*etas)))
        return(list(REML=LL))
      }
      
      replaceNaN<-  function(LL) {
        index=(LL=="NaN")
        if(length(index)>0) theMin=min(LL[!index])
        if(length(index)<1) theMin="NaN"
        LL[index]=theMin
        return(LL)    
      }
      
      lars <-  function(x, y, type = c("lasso", "lar", "forward.stagewise","stepwise"), trace = FALSE,
                        normalize=TRUE, intercept=TRUE, Gram, 
                        eps = .Machine$double.eps,  max.steps, use.Gram = TRUE)
      {
        
        call <- match.call()
        type <- match.arg(type)
        TYPE <- switch(type,
                       lasso = "LASSO",
                       lar = "LAR",
                       forward.stagewise = "Forward Stagewise",
                       stepwise = "Forward Stepwise")
        if(trace)
          cat(paste(TYPE, "sequence\n"))
        
        nm <- dim(x)
        n <- nm[1]
        m <- nm[2]
        im <- inactive <- seq(m)
        one <- rep(1, n)
        vn <- dimnames(x)[[2]]  
        ### Center x and y, and scale x, and save the means and sds
        if(intercept){
          meanx <- drop(one %*% x)/n
          x <- scale(x, meanx, FALSE)  # centers x
          mu <- mean(y)
          y <- drop(y - mu)
        }
        else {
          meanx <- rep(0,m)
          mu <- 0
          y <- drop(y)
        }
        if(normalize){
          normx <- sqrt(drop(one %*% (x^2)))
          nosignal<-normx/sqrt(n) < eps
          if(any(nosignal))# ignore variables with too small a variance
          {
            ignores<-im[nosignal]
            inactive<-im[-ignores]
            normx[nosignal]<-eps*sqrt(n)
            if(trace)
              cat("LARS Step 0 :\t", sum(nosignal), "Variables with Variance < eps; dropped for good\n")  #
          }
          else ignores <- NULL #singularities; augmented later as well
          names(normx) <- NULL
          x <- scale(x, FALSE, normx)	# scales x
        }
        else {
          normx <- rep(1,m)
          ignores <- NULL
        }
        if(use.Gram & missing(Gram)) {
          if(m > 500 && n < m)
            cat("There are more than 500 variables and n<m;\nYou may wish to restart and set use.Gram=FALSE\n"
            )
          if(trace)
            cat("Computing X'X .....\n")
          Gram <- t(x) %*% x	#Time saving
        }
        Cvec <- drop(t(y) %*% x)
        ssy <- sum(y^2)	### Some initializations
        residuals <- y
        if(missing(max.steps))
          max.steps <- 8*min(m, n-intercept)
        beta <- matrix(0, max.steps + 1, m)	# beta starts at 0
        lambda=double(max.steps)
        Gamrat <- NULL
        arc.length <- NULL
        R2 <- 1
        RSS <- ssy
        first.in <- integer(m)
        active <- NULL	# maintains active set
        actions <- as.list(seq(max.steps))	
        
        drops <- FALSE
        Sign <- NULL
        R <- NULL	###
        
        k <- 0
        while((k < max.steps) & (length(active) < min(m - length(ignores),n-intercept)) )
        {
          action <- NULL
          C <- Cvec[inactive]	#
          
          Cmax <- max(abs(C))
          if(Cmax<eps*100){ 
            if(trace)cat("Max |corr| = 0; exiting...\n")
            break
          }
          k <- k + 1
          lambda[k]=Cmax
          
          if(!any(drops)) {
            new <- abs(C) >= Cmax - eps
            C <- C[!new]	# for later
            new <- inactive[new]	# Get index numbers
            
            for(inew in new) {
              if(use.Gram) {
                R <- updateR(Gram[inew, inew], R, drop(Gram[
                  inew, active]), Gram = TRUE,eps=eps)
              }
              else {
                R <- updateR(x[, inew], R, x[, active], Gram
                             = FALSE,eps=eps)
              }
              if(attr(R, "rank") == length(active)) {
                
                nR <- seq(length(active))
                R <- R[nR, nR, drop = FALSE]
                attr(R, "rank") <- length(active)
                ignores <- c(ignores, inew)
                action <- c(action,  - inew)
                if(trace)
                  cat("LARS Step", k, ":\t Variable", inew, 
                      "\tcollinear; dropped for good\n")	#
              }
              else {
                if(first.in[inew] == 0)
                  first.in[inew] <- k
                active <- c(active, inew)
                Sign <- c(Sign, sign(Cvec[inew]))
                action <- c(action, inew)
                if(trace)
                  cat("LARS Step", k, ":\t Variable", inew, 
                      "\tadded\n")
              }
            }
          }
          else action <-  - dropid
          Gi1 <- backsolve(R, backsolvet(R, Sign))	
          
          dropouts<-NULL
          if(type == "forward.stagewise") {
            directions <- Gi1 * Sign
            if(!all(directions > 0)) {
              if(use.Gram) {
                nnls.object <- nnls.lars(active, Sign, R, 
                                         directions, Gram[active, active], trace = 
                                           trace, use.Gram = TRUE,eps=eps)
              }
              else {
                nnls.object <- nnls.lars(active, Sign, R, 
                                         directions, x[, active], trace = trace, 
                                         use.Gram = FALSE,eps=eps)
              }
              positive <- nnls.object$positive
              dropouts <-active[-positive]
              action <- c(action, -dropouts)
              active <- nnls.object$active
              Sign <- Sign[positive]
              Gi1 <- nnls.object$beta[positive] * Sign
              R <- nnls.object$R
              C <- Cvec[ - c(active, ignores)]
            }
          }
          A <- 1/sqrt(sum(Gi1 * Sign))
          w <- A * Gi1	# note that w has the right signs
          if(!use.Gram) u <- drop(x[, active, drop = FALSE] %*% w)	###
          
          if( (length(active) >=  min(n-intercept, m - length(ignores) ) )|type=="stepwise") {
            gamhat <- Cmax/A
          }
          else {
            if(use.Gram) {
              a <- drop(w %*% Gram[active,  - c(active,ignores), drop = FALSE])
            }
            else {
              a <- drop(u %*% x[,  - c(active, ignores), drop=FALSE])
            }
            gam <- c((Cmax - C)/(A - a), (Cmax + C)/(A + a))	
            
            gamhat <- min(gam[gam > eps], Cmax/A)	
          }
          if(type == "lasso") {
            dropid <- NULL
            b1 <- beta[k, active]	# beta starts at 0
            z1 <-  - b1/w
            zmin <- min(z1[z1 > eps], gamhat)
            if(zmin < gamhat) {
              gamhat <- zmin
              drops <- z1 == zmin
            }
            else drops <- FALSE
          }
          beta[k + 1,  ] <- beta[k,  ]
          beta[k + 1, active] <- beta[k + 1, active] + gamhat * w
          if(use.Gram) {
            Cvec <- Cvec - gamhat * Gram[, active, drop = FALSE] %*% w
          }
          else {
            residuals <- residuals - gamhat * u
            Cvec <- drop(t(residuals) %*% x)
          }
          Gamrat <- c(Gamrat, gamhat/(Cmax/A))
          arc.length <- c(arc.length, gamhat)	
          if(type == "lasso" && any(drops)) {
            dropid <- seq(drops)[drops]	
            for(id in rev(dropid)) {
              if(trace)
                cat("Lasso Step", k+1, ":\t Variable", active[
                  id], "\tdropped\n")
              R <- downdateR(R, id)
            }
            dropid <- active[drops]
            beta[k+1,dropid]<-0  
            active <- active[!drops]
            Sign <- Sign[!drops]
          }
          if(!is.null(vn))
            names(action) <- vn[abs(action)]
          actions[[k]] <- action
          inactive <- im[ - c(active, ignores)]
          if(type=="stepwise")Sign=Sign*0
        }
        beta <- beta[seq(k + 1), ,drop=FALSE ]	
        lambda=lambda[seq(k)]
        dimnames(beta) <- list(paste(0:k), vn)
        if(trace)
          cat("Computing residuals, RSS etc .....\n")
        residuals <- y - x %*% t(beta)
        beta <- scale(beta, FALSE, normx)
        RSS <- apply(residuals^2, 2, sum)
        R2 <- 1 - RSS/RSS[1]
        actions=actions[seq(k)]
        netdf=sapply(actions,function(x)sum(sign(x)))
        df=cumsum(netdf)### This takes into account drops
        if(intercept)df=c(Intercept=1,df+1)
        else df=c(Null=0,df)
        rss.big=rev(RSS)[1]
        df.big=n-rev(df)[1]
        if(rss.big<eps|df.big<eps)sigma2=NaN
        else
          sigma2=rss.big/df.big
        Cp <- RSS/sigma2 - n + 2 * df
        attr(Cp,"sigma2")=sigma2
        attr(Cp,"n")=n
        object <- list(call = call, type = TYPE, df=df, lambda=lambda,R2 = R2, RSS = RSS, Cp = Cp, 
                       actions = actions[seq(k)], entry = first.in, Gamrat = Gamrat, 
                       arc.length = arc.length, Gram = if(use.Gram) Gram else NULL, 
                       beta = beta, mu = mu, normx = normx, meanx = meanx)
        class(object) <- "lars"
        object
      }
      
      Y.data<-as.matrix(phe)
      if(is.null(psmatrix)==FALSE){
        psmatrix<-as.matrix(psmatrix)
      }
      nsam <-ncol(gene.data)-2
      chrnum<-length(unique(gene.data[,1]))
      
      W.orig<-matrix(1,nsam,1)
      if(is.null(psmatrix)==FALSE){
        W1 <-cbind(W.orig,psmatrix)
      }else{
        W1<-W.orig
      }
      
      kk<-list(NULL)
      cc<-list(NULL)
      kktotal<-matrix(0,nsam,nsam)
      
      for(i in 1:chrnum){
        xxot <- as.matrix(gene.data[gene.data[,1]==i,3:ncol(gene.data)])
        xot <-t(xxot)
        nmarkot <-ncol(xot)
        # kk[[i]]<-matrix(0,nsam,nsam)
        # for(k in 1:nmarkot)
        # {
        #   z<-as.matrix(xot[,k])
        #   kk[[i]]<-kk[[i]]+z%*%t(z)
        # }
        kk[[i]]<-mrMLM::multiplication_speed(xot,t(xot))
        cc[[i]]<-mean(diag(kk[[i]]))
        kktotal<-kktotal+kk[[i]]
      }
      
      rm(xot)
      gc()
      
      
      cl.cores <- detectCores()
      if((cl.cores<=2)||(is.null(CLO)==FALSE)){
        cl.cores<-1
      }else if(cl.cores>2){
        if(cl.cores>10){
          cl.cores<-10
        }else {  
          cl.cores <- detectCores()-1
        }
      }   
      
      cl <- makeCluster(cl.cores)
      registerDoParallel(cl)
      
      
      larsres<-foreach(i=1:chrnum,.multicombine=TRUE,.combine='rbind')%dopar%{
        
        requireNamespace("foreach")
        requireNamespace("lars")
        requireNamespace("sampling") 
        
        xxot <- as.matrix(gene.data[(gene.data[,1])!=i,3:ncol(gene.data)])
        
        xx1 <- as.matrix(gene.data[gene.data[,1]==i,3:ncol(gene.data)])
        YY1 <- matrix(Y.data,,1)
        
        K1 <- (kktotal-kk[[i]])/(sum(unlist(cc))-as.numeric(cc[i]))
        
        repl<-numeric()
        if(Bootstrap==TRUE){
          
          res1<-foreach(repl=1:5,.multicombine=TRUE,.combine='cbind')%do%{
            
            if(repl==1){
              YY<-YY1
              xx<-xx1
              K<-K1
              W<-W1
            }else{
              s<-srswr(nrow(YY1),nrow(YY1))
              ind<-(1:nrow(YY1))[s!=0]
              n<-s[s!=0]
              ind<-rep(ind,times=n)
              YY<-as.matrix(YY1[ind,])
              xx<-xx1[,ind]
              
              xxot2<-xxot[,ind]
              xot2 <-t(xxot2)
              nmarkot2 <-ncol(xot2)
              kk2<-matrix(0,nsam,nsam)
              for(k in 1:nmarkot2)
              {
                z2<-as.matrix(xot2[,k])
                kk2<-kk2+z2%*%t(z2)
              }
              cc2<-mean(diag(kk2))
              K <- numeric()
              K <- kk2/cc2
              W<-as.matrix(W1[ind,])
            }
            
            remle2<-emma.REMLE(YY, W, K, Z=NULL, ngrids=100, llim=-10, ulim=10,esp=1e-10, eig.L = NULL, eig.R = NULL)
            remle1.B1<-emma.maineffects.B(Z=NULL,K,remle2$delta)
            rm(K)
            gc()
            C2<-remle1.B1$mC
            Y_c <- C2%*%YY
            W_c <- C2%*%W
            G_c <- C2%*%t(xx)
            GGG <- t(G_c)
            rm(G_c)
            gc()
            ylars <- as.matrix(Y_c)
            xlars <- cbind(W_c,t(GGG))
            rm(GGG)
            gc()
            LAR <- lars(xlars,ylars,type="lar",use.Gram=F,max.steps=lars1)
            rm(xlars)
            gc()
            LAR$beta[nrow(LAR$beta),]
          }
          
        }else if(Bootstrap==FALSE){
          res1 <- numeric()
          remle2<-emma.REMLE(YY1, W1, K1, Z=NULL, ngrids=100, llim=-10, ulim=10,esp=1e-10, eig.L = NULL, eig.R = NULL)
          remle1.B1<-emma.maineffects.B(Z=NULL,K1,remle2$delta)
          rm(K1)
          gc()
          C2<-remle1.B1$mC
          Y_c <- C2%*%YY1
          W_c <- C2%*%W1
          G_c <- C2%*%t(xx1)
          rm(xx1)
          gc()
          GGG <- t(G_c)
          rm(G_c)
          gc()
          ylars <- as.matrix(Y_c)
          xlars <- cbind(W_c,t(GGG))
          rm(GGG)
          gc()
          LAR <- lars(xlars,ylars,type="lar",use.Gram=F,max.steps=lars1)
          rm(xlars)
          gc()
          res1<-cbind(res1,LAR$beta[nrow(LAR$beta),])
        }  
        
        if(is.null(psmatrix)==FALSE){
          rr <- as.matrix(res1[-c(1:(ncol(psmatrix)+1)),])
        }else{
          rr <- as.matrix(res1[-1,])
        }
      }
      
      stopCluster(cl)
      
      rm(kk,kktotal)
      gc()
      
     if(Bootstrap==TRUE){
        count <- matrix(rep(0,nrow(larsres)),nrow(larsres),1)
        
        ttt <- numeric()
        for(ii in 1:nrow(larsres))
        {
          tt <- 0
          for(jj in 1:ncol(larsres))
          {
            if ((abs(larsres[ii,jj]))>0){tt <- tt+1}
          }
          count[ii] <-tt
        }
        larsres <-cbind(larsres,count)
        
        for(ii in 1:nrow(larsres))
        {
          if(larsres[ii,ncol(larsres)]>=3){ttt <- cbind(ttt,ii)}
        }
        
        countnum <- ttt
        
      }else{
        
        countnum <- numeric()
        for(ii in 1:nrow(larsres))
        {
          if ((abs(larsres[ii]))>0){countnum <- cbind(countnum,ii)}
        }
      }
      if(ncol(countnum)>nrow(phe)){
        
        if(length(countnum)==1){
          xx2 <- matrix(gene.data[c(countnum),3:ncol(gene.data)],1,)
          
        }else{
          xx2 <- as.matrix(gene.data[c(countnum),3:ncol(gene.data)])
        }
        YY2 <- matrix(Y.data,,1)
        
        ylars <- as.matrix(YY2)
        xlars <- cbind(W1,t(xx2))
        LAR <- lars(xlars,ylars,type="lar",use.Gram=F)
        
        res1<-as.matrix(LAR$beta[nrow(LAR$beta),])
        
        rm(xlars,xx2)
        gc()
        
        if(is.null(psmatrix)==FALSE){
          rr <- as.matrix(res1[-c(1:(ncol(psmatrix)+1)),])
        }else{
          rr <- as.matrix(res1[-1,])
        }
        
        ct <- numeric()
        for(ii in 1:nrow(rr))
        {
          if ((abs(rr[ii]))>0){ct <- cbind(ct,ii)}
        }
        
        inct<-c(ct)
        countnum<-countnum[,inct]
        
      }
      
      if(length(countnum)==1){
        xeb <- matrix(gene.data[c(countnum),],1,)
        ebrow <-matrix(xeb[,1:2],,2)
        xeb1<-matrix(xeb[,3:ncol(xeb)],1,)
        xxeb <- as.matrix(t(xeb1))
        nmak <- ncol(xxeb)
        
      }else{
        xeb <- as.matrix(gene.data[c(countnum),])
        ebrow <-as.matrix(xeb[,1:2])
        xeb1<-as.matrix(xeb[,3:ncol(xeb)])
        xxeb <- as.matrix(t(xeb1))
        nmak <- ncol(xxeb)
      }
      bayeslodres <- numeric()
      
      genname<-gene.data[,1:2] 
      
      rm(xeb,gene.data)
      gc()
      
      yeb <- as.matrix(phe)
      
      if(is.null(psmatrix)==FALSE){
        u1<-ebayes_EM(cbind(matrix(1,nrow(xxeb),1),psmatrix),xxeb,yeb)
        xb<-u1$u
      }else{
        u1<-ebayes_EM(matrix(1,nrow(xxeb),1),xxeb,yeb)
        xb<-u1$u
      }
      xb<-as.matrix(xb)
      if(is.null(psmatrix)==FALSE){
        temp<-cbind(matrix(1,nrow(xxeb),1),psmatrix)
      }else{
        temp<-matrix(1,nrow(xxeb),1)
      }
      
      lodres<-likelihood(temp,xxeb,yeb,xb)
      lodres<-as.matrix(lodres)
      #### compute heredity#######
      ch_er <- as.numeric()
      ch_x <- cbind(matrix(1,nrow(xxeb),1),xxeb)
      
      ch_bb <- rbind(mean(yeb),as.matrix(xb))
      
      rm(xxeb)
      gc()
      
      
      for(i in 1:(ncol(ch_x)-1))
      {
        ch_xi <- ch_x[,(1+i)]
        as1 <- length(which(ch_xi==1))/nrow(ch_x)
        as2 <- 1-as1
        ch_er <- rbind(ch_er,(1-(as1-as2)*(as1-as2))*ch_bb[i+1]*ch_bb[i+1])
      }
      ch_v0 <- (1/(nrow(ch_x)-1))*(t(yeb-ch_x%*%ch_bb)%*%(yeb-ch_x%*%ch_bb))
      
      rm(ch_x)
      gc()
      
      
      if(var(yeb)>=sum(ch_er)+ch_v0){
        hered <- (ch_er/as.vector(var(yeb)))*100 
      }else{
        hered <- (ch_er/as.numeric(sum(ch_er)+ch_v0))*100
      }
      
      bayeslodres<-cbind(ebrow,xb,lodres,hered)
      
      
      lodid<-which(bayeslodres[,4]>lodvalue)
      if(length(lodid)!=0){
        
        if(length(lodid)==1){
          lastres<-matrix(bayeslodres[lodid,],1,)
          xeb2<-matrix(xeb1[lodid,],1,)
        }else{
          lastres<-bayeslodres[lodid,]
          xeb2<-as.matrix(xeb1[lodid,])
        }
        
        rm(xeb1)
        gc()
        
        xxmaf<- xeb2
        leng.maf<-dim(xxmaf)[2]
        
        maf.fun<-function(snp){
          leng<-length(snp)
          snp1<-length(which(snp==1))
          snp11<-length(which(snp==-1))
          snp0<-length(which(snp==0))
          ma1<-(2*snp1+snp0)/(2*leng)
          ma2<-(2*snp11+snp0)/(2*leng)
          maf<-min(ma1,ma2)
          return(maf)
        }
        
        maf<-apply(xxmaf,1,maf.fun)
        maf<-as.matrix(round(maf,4))
        vee<-round(u1$sigma2,4)
        pee<-round(var(yeb),4)
        
        vees<-matrix("",nrow = nrow(lastres),1)
        pees<-matrix("",nrow = nrow(lastres),1)
        pees[1,1]<-pee
        vees[1,1]<-vee
        result<-lastres
        result<-result
        
        if(nrow(result)>1){
          temp<-as.matrix(result[,3:5])
          temp[which(abs(temp)>=1e-4)]<-round(temp[abs(temp)>=1e-4],4)
          temp[which(abs(temp)<1e-4)]<-as.numeric(sprintf("%.4e",temp[abs(temp)<1e-4]))
          wan<-cbind(result[,1:2],temp)
        }else{
          temp<-t(as.matrix(result[,3:5]))
          temp[which(abs(temp)>=1e-4)]<-round(temp[abs(temp)>=1e-4],4)
          temp[which(abs(temp)<1e-4)]<-as.numeric(sprintf("%.4e",temp[abs(temp)<1e-4]))
          wan<-cbind(t(as.matrix(result[,1:2])),temp)  
        }
        
        if(inputform==1){
          genRaw<-as.data.frame(genRaw)
          genraw<-genRaw[-1,1:4]
          
          wan_len<-dim(wan)[1]
          marker<-character()
          snp<-character()
          
          for(i in 1:wan_len){
            chr_pos<-which(genraw[,2]==wan[i,1])
            new_matrix<-genraw[chr_pos,]
            posi_pos<-which(new_matrix[,3]==wan[i,2])
            mark<-matrix(new_matrix[posi_pos,1],1,)
            marker<-rbind(marker,mark)
            sn<-matrix(new_matrix[posi_pos,4],1,)
            snp<-rbind(snp,sn)
          }
        }
        if(inputform==2){
          
          genRaw<-as.data.frame(genRaw)
          genraw<-genRaw[-1,1:4]
          
          wan_len<-dim(wan)[1]
          marker<-character()
          snp<-character()
          for(i in 1:wan_len){
            chr_pos<-which(genraw[,2]==wan[i,1])
            new_matrix<-genraw[chr_pos,]
            posi_pos<-which(new_matrix[,3]==wan[i,2])
            mark<-matrix(new_matrix[posi_pos,1],1,)
            marker<-rbind(marker,mark)
            sn<-matrix(new_matrix[posi_pos,4],1,)
            snp<-rbind(snp,sn)
          }
        }
        if(inputform==3){
          genRaw<-as.data.frame(genRaw)
          genraw<-genRaw[-1,c(1,3,4,12)]
          
          wan_len<-dim(wan)[1]
          marker<-character()
          snp<-character()
          for(i in 1:wan_len){
            chr_pos<-which(genraw[,2]==wan[i,1])
            new_matrix<-genraw[chr_pos,]
            posi_pos<-which(new_matrix[,3]==wan[i,2])
            mark<-matrix(new_matrix[posi_pos,1],1,)
            marker<-rbind(marker,mark)
            sn<-matrix(new_matrix[posi_pos,4],1,)
            snp<-rbind(snp,sn)
          }
        }
        
        wan<-cbind(marker,wan,maf,snp,vees,pees)
        tempwan <- wan
        lodscore1 <- as.numeric(tempwan[,5])
        log10P <- as.matrix(round(-log10(pchisq(lodscore1*4.605,1,lower.tail = F)),4))
        if(nrow(tempwan)>1){
          tempwan1 <- cbind(tempwan[,1:5],log10P,tempwan[,6:10])
        }else{
          tempwan1 <- cbind(t(as.matrix(tempwan[,1:5])),log10P,t(as.matrix(tempwan[,6:10])))  
        }
        wan <- tempwan1
        colnames(wan)<-c("RS#","Chromosome","Marker position (bp)","QTN effect","LOD score","'-log10(P)'","r2 (%)","MAF","Genotype for code 1","Var_error","Var_phen (total)")
        wan<-as.data.frame(wan)
      }
      output<-list(result=wan)
    } 
    return(output) 
  }