Nothing
mrMLMFun<-function(gen,phe,outATCG,genRaw,kk,psmatrix,svpal,svrad,svmlod,Genformat,CLO){
inputform<-Genformat
if(is.null(kk)){
if(is.null(gen)==TRUE)
{
showModal(modalDialog(title = "Warning", strong("Please input correct genotypic dataset !"), easyClose = TRUE))
}else{
envgenq<-deepcopy(gen,3:ncol(gen))
envgenq2<-t(envgenq[,])
# envgen<-big.matrix(nrow(envgenq2),ncol(envgenq2),type='double',shared = FALSE)
# envgen[,]<-envgenq2[,]
rm(envgenq)
gc()
# m<-ncol(envgen)
# n<-nrow(envgen)
#kk1<-matrix(0,n,n)
# for(k in 1:m){
# z<-as.matrix(envgen[,k])
# kk1<-kk1+z%*%t(z)
# }
kk1<-mrMLM.GUI::multiplication_speed(envgenq2,t(envgenq2))
cc<-mean(diag(kk1))
kk1<-kk1/cc
kk<-as.matrix(kk1)
}
rm(kk1)
gc()
}
if(is.null(psmatrix)){
flagps<-1
}else{
flagps<-0
}
if(is.null(svrad)==TRUE||is.null(svmlod)==TRUE){
showModal(modalDialog(title = "Warning", strong("Please set parameters!"), easyClose = TRUE))
}
if(svrad<0)
{
showModal(modalDialog(title = "Warning", strong("Please input search radius (kb) of candidate gene: > 0 !"), easyClose = TRUE))
}
if(svmlod<0)
{
showModal(modalDialog(title = "Warning", strong("Please input critical LOD score: > 0 !"), easyClose = TRUE))
}
if(exists("gen")==FALSE)
{
showModal(modalDialog(title = "Warning", strong("Please input correct genotypic dataset!"), easyClose = TRUE))
}
if(exists("phe")==FALSE)
{
showModal(modalDialog(title = "Warning", strong("Please input correct phenotypic dataset !"), easyClose = TRUE))
}
if(exists("kk")==FALSE)
{
showModal(modalDialog(title = "Warning", strong("Please input correct kinship (K) dataset !"), easyClose = TRUE))
}
if((exists("gen")==TRUE)&&(exists("phe")==TRUE)&&(ncol(gen)!=(nrow(phe)+2)))
{
showModal(modalDialog(title = "Warning", strong("Sample size in genotypic dataset doesn't equal to the sample size in phenotypic dataset !"), easyClose = TRUE))
}
if((exists("gen")==TRUE)&&(exists("phe")==TRUE)&&(exists("kk")==TRUE)&&((ncol(gen)==(nrow(phe)+2)))&&(svrad>0)&&(svmlod>=0))
{
parmsShow<-NULL
wan<-NULL
parms<-NULL
parms.pchange<-NULL
mannewp<-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)
}
mixed<-function(x,y,kk){
loglike<-function(theta){
lambda<-exp(theta)
logdt<-sum(log(lambda*delta+1))
h<-1/(lambda*delta+1)
yy<-sum(yu*h*yu)
yx<-matrix(0,q,1)
xx<-matrix(0,q,q)
for(i in 1:q){
yx[i]<-sum(yu*h*xu[,i])
for(j in 1:q){
xx[i,j]<-sum(xu[,i]*h*xu[,j])
}
}
loglike<- -0.5*logdt-0.5*(n-q)*log(yy-t(yx)%*%solve(xx)%*%yx)-0.5*log(det(xx))
return(-loglike)
}
fixed<-function(lambda){
h<-1/(lambda*delta+1)
yy<-sum(yu*h*yu)
yx<-matrix(0,q,1)
xx<-matrix(0,q,q)
for(i in 1:q){
yx[i]<-sum(yu*h*xu[,i])
for(j in 1:q){
xx[i,j]<-sum(xu[,i]*h*xu[,j])
}
}
beta<-solve(xx,yx)
sigma2<-(yy-t(yx)%*%solve(xx)%*%yx)/(n-q)
sigma2<-drop(sigma2)
var<-diag(solve(xx)*sigma2)
stderr<-sqrt(var)
return(c(beta,stderr,sigma2))
}
qq<-eigen(kk)
delta<-qq[[1]]
uu<-qq[[2]]
q<-ncol(x)
n<-ncol(kk)
vp<-var(y)
yu<-t(uu)%*%y
xu<-t(uu)%*%x
theta<-0
parm<-optim(par=theta,fn=loglike,hessian = TRUE,method="L-BFGS-B",lower=-50,upper=10)
lambda<-exp(parm$par)
conv<-parm$convergence
fn1<-parm$value
fn0<-loglike(-Inf)
lrt<-2*(fn0-fn1)
hess<-parm$hessian
parmfix<-fixed(lambda)
beta<-parmfix[1:q]
stderr<-parmfix[(q+1):(2*q)]
sigma2<-parmfix[2*q+1]
lod<-lrt/4.61
p_value<-pchisq(lrt,1,lower.tail = F)
sigma2g<-lambda*sigma2
goodness<-(vp-sigma2)/vp
par<-data.frame(lrt,beta,stderr,sigma2,lambda,sigma2g,lod,p_value)
return(par)
}
loglike<-function(theta){
xi<-exp(theta)
tmp0<-zz*xi+1
tmp<-xi*solve(tmp0)
yHy<-yy-t(zy)%*%tmp%*%zy
yHx<-yx-zx%*%tmp%*%zy
xHx<-xx-zx%*%tmp%*%t(zx)
logdt2<-log(det(tmp0))
loglike<- -0.5*logdt2-0.5*(n-s)*log(yHy-t(yHx)%*%solve(xHx)%*%yHx)-0.5*log(det(xHx))
return(-loglike)
}
fixed<-function(xi){
tmp0<-zz*xi+diag(1)
tmp<-xi*solve(tmp0)
yHy<-yy-t(zy)%*%tmp%*%zy
yHx<-yx-zx%*%tmp%*%zy
xHx<-xx-zx%*%tmp%*%t(zx)
zHy<-zy-zz%*%tmp%*%zy
zHx<-zx-zx%*%tmp%*%zz
zHz<-zz-zz%*%tmp%*%zz
beta<-solve(xHx,yHx)
tmp2<-solve(xHx)
sigma2<-(yHy-t(yHx)%*%tmp2%*%yHx)/(n-s)
gamma<-xi*zHy-xi*t(zHx)%*%tmp2%*%yHx
var<-abs((xi*diag(1)-xi*zHz*xi)*as.numeric(sigma2))
stderr<-sqrt(diag(var))
result<-list(gamma,stderr,beta,sigma2)
return(result)
}
name<-gen[,1:2]
genq<-deepcopy(gen,3:ncol(gen))
genq2<-t(genq[,])
gen<-big.matrix(nrow(genq2),ncol(genq2),type='double',shared = FALSE)
gen[,]<-genq2[,]
rm(genq2)
gc()
n<-nrow(gen)
m<-ncol(gen)
if((flagps==1)||(exists("psmatrix")==FALSE))
{
x<-matrix(1,n,1)
}else if(flagps==0)
{
x<-cbind(matrix(1,n,1),psmatrix)
}
ll<-numeric()
s<-ncol(x)
kk<-as.matrix(kk)
qq<-eigen(kk)
delta<-qq[[1]]
uu<-qq[[2]]
xu<-t(uu)%*%x
rm(qq)
gc()
yy<-phe[,1]
y<-as.matrix(yy)
parm<-mixed(x=x,y=y,kk=kk)
lambda<-parm$lambda[1]
h<-1/(delta*lambda+1)
yu<-t(uu)%*%y
xx<-matrix(0,s,s)
for(i in 1:s){
for(j in 1:s){
xx[i,j]<-sum(xu[,i]*h*xu[,j])
}
}
yy<-sum(yu*h*yu)
yx<-matrix(0,s,1)
for(i in 1:s){
yx[i]<-sum(yu*h*xu[,i])
}
genf<-gen[,]
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)
if((flagps==1)||(is.null("psmatrix")))
{
k<-numeric()
ff=foreach(k=1:m, .multicombine=TRUE, .combine = 'rbind')%dopar%
{
#browser()
z<-as.matrix(genf[,k])
zu<-t(uu)%*%z
zy<-as.matrix(sum(yu*h*zu))
zz<-as.matrix(sum(zu*h*zu))
zx<-matrix(0,s,1)
for(i in 1:s){
zx[i]<-sum(xu[,i]*h*zu)
}
theta<-c(0)
par<-optim(par=theta,fn=loglike,hessian = TRUE,method="L-BFGS-B",lower=-10,upper=10)
xi<-exp(par$par)
conv<-par$convergence
fn1<-par$value
hess<-par$hessian
parmfix<-fixed(xi)
gamma<-parmfix[[1]]
stderr<-parmfix[[2]]
beta<-parmfix[[3]]
sigma2<-parmfix[[4]]
lambda<-xi
sigma2g<-lambda*sigma2
fn0<-loglike(-Inf)
lrt<-2*(fn0-fn1)
p_lrt<-pchisq(lrt,1,lower.tail = F)
wald<-(gamma/stderr)^2
p_wald<-pchisq(wald,1,lower.tail = F)
parm0<-c(1,name[k,1],name[k,2],beta,sigma2,sigma2g,gamma,stderr,wald,p_wald)
}
stopCluster(cl)
ll<-rbind(ll,ff)
}else if(flagps==0){
k<-numeric()
ff=foreach(k=1:m, .multicombine=TRUE, .combine = 'rbind')%dopar%
{
#browser()
z<-as.matrix(genf[,k])
zu<-t(uu)%*%z
zy<-as.matrix(sum(yu*h*zu))
zz<-as.matrix(sum(zu*h*zu))
zx<-matrix(0,s,1)
for(i in 1:s){
zx[i]<-sum(xu[,i]*h*zu)
}
theta<-c(0)
par<-optim(par=theta,fn=loglike,hessian = TRUE,method="L-BFGS-B",lower=-10,upper=10)
xi<-exp(par$par)
conv<-par$convergence
fn1<-par$value
hess<-par$hessian
parmfix<-fixed(xi)
gamma<-parmfix[[1]]
stderr<-parmfix[[2]]
beta<-parmfix[[3]][1]
sigma2<-parmfix[[4]]
lambda<-xi
sigma2g<-lambda*sigma2
fn0<-loglike(-Inf)
lrt<-2*(fn0-fn1)
p_lrt<-pchisq(lrt,1,lower.tail = F)
wald<-(gamma/stderr)^2
p_wald<-pchisq(wald,1,lower.tail = F)
parm0<-c(1,name[k,1],name[k,2],beta,sigma2,sigma2g,gamma,stderr,wald,p_wald)
}
stopCluster(cl)
ll<-rbind(ll,ff)
}
rm(uu,kk,genf)
gc()
parms<-ll
parms<-matrix(parms,,10)
chr_pos<-parms[,2:3]
pfit<-which(parms[,10]<=(svpal))
pfit<-as.matrix(pfit)
pfitrow<-nrow(pfit)
no_p<-cbind((1:(nrow(parms))),parms[,10])
no_porder<-order(no_p[,2])
no_p<-no_p[no_porder,]
choose_orderp<-no_p[1:pfitrow,]
orderno<-no_p[1:pfitrow,1]
orderno<-as.matrix(orderno)
sigma2g_SNPerr<-cbind(parms[,6],parms[,8])
correct_each<-matrix(1,(nrow(sigma2g_SNPerr)),1)-sigma2g_SNPerr[,2]*sigma2g_SNPerr[,2]/sigma2g_SNPerr[,1]
k0<-which(correct_each<0)
k0<-as.matrix(k0)
if(nrow(k0)>0){
correct_each[k0,1]<-matrix(0,(nrow(k0)),1)
}
correct_sum<-sum(correct_each)
newp<-0.05/correct_sum
mannewp<-newp
manstandchoice<-1
no_porder<-which(no_p[,2]<=newp)
no_porder<-as.matrix(no_porder)
no_porderrow<-nrow(no_porder)
gg<-orderno
if(nrow(orderno)>1){
for (ii in 1:(nrow(orderno)-1)){
for (jj in (ii+1):(nrow(orderno))){
ci<- chr_pos[orderno[ii],1]
cj<- chr_pos[orderno[jj],1]
if (ci==cj){
ye<-abs(chr_pos[orderno[ii],2]-chr_pos[orderno[jj],2])
if (ye<=((svrad)*1000)){
gg[jj,1]<-0
}
}
}
}
}
parms.pchange<-parms
parmsp<-as.matrix(parms.pchange[,10])
locsub<-which(parmsp==0)
if(length(locsub)!=0){
pmin<-min(parmsp[parmsp!=0])
subvalue<-10^(1.1*log10(pmin))
parms.pchange[locsub,10]<-subvalue
}else{
parms.pchange<-parms
}
if(inputform==1){
#output result1 using mrMLM numeric format
parmsShow<-parms[,-1]
meadd<-matrix(1,nrow(parms),1)
meadd[which(parms[,10]<newp),1]<-sprintf("%.4e",newp)
meadd[which(parms[,10]>=newp),1]<-" "
tempparms<-parms[,4:10]
tempparms[,7]<--log10(tempparms[,7])
tempparms[which(abs(tempparms)>=1e-4)]<-round(tempparms[which(abs(tempparms)>=1e-4)],4)
tempparms[which(abs(tempparms)<1e-4)]<-as.numeric(sprintf("%.4e",tempparms[which(abs(tempparms)<1e-4)]))
parmsShow<-cbind(genRaw[-1,1],parms[,2:3],tempparms,genRaw[-1,4],meadd)
colnames(parmsShow)<-c("RS#","Chromosome","Marker position (bp)","Mean","Sigma2","Sigma2_k","SNP effect (mrMLM)","Sigma2_k_posteriori","Wald","'-log10(P) (mrMLM)'","Genotype for code 1","Significance")
}
if(inputform==2){
#output result1 using mrMLM character format
parmsShow<-parms[,-1]
outATCG<-matrix(outATCG,,1)
meadd<-matrix(1,nrow(parms),1)
meadd[which(parms[,10]<newp),1]<-sprintf("%.4e",newp)
meadd[which(parms[,10]>=newp),1]<-" "
tempparms<-parms[,4:10]
tempparms[,7]<--log10(tempparms[,7])
tempparms[which(abs(tempparms)>=1e-4)]<-round(tempparms[which(abs(tempparms)>=1e-4)],4)
tempparms[which(abs(tempparms)<1e-4)]<-as.numeric(sprintf("%.4e",tempparms[which(abs(tempparms)<1e-4)]))
parmsShow<-cbind(genRaw[-1,1],parms[,2:3],tempparms,outATCG,meadd)
colnames(parmsShow)<-c("RS#","Chromosome","Marker position (bp)","Mean","Sigma2","Sigma2_k","SNP effect (mrMLM)","Sigma2_k_posteriori","Wald","'-log10(P) (mrMLM)'","Genotype for code 1","Significance")
}
if(inputform==3){
#output result1 using TASSEL format
parmsShow<-parms[,-1]
outATCG<-matrix(outATCG,,1)
#outATCG<-unlist(strsplit(outATCG,""))
#outATCG<-matrix(outATCG[c(TRUE,FALSE)],,1)
meadd<-matrix(1,nrow(parms),1)
meadd[which(parms[,10]<newp),1]<-sprintf("%.4e",newp)
meadd[which(parms[,10]>=newp),1]<-" "
tempparms<-parms[,4:10]
tempparms[,7]<--log10(tempparms[,7])
tempparms[which(abs(tempparms)>=1e-4)]<-round(tempparms[which(abs(tempparms)>=1e-4)],4)
tempparms[which(abs(tempparms)<1e-4)]<-as.numeric(sprintf("%.4e",tempparms[which(abs(tempparms)<1e-4)]))
parmsShow<-cbind(genRaw[-1,1],parms[,2:3],tempparms,outATCG,meadd)
colnames(parmsShow)<-c("RS#","Chromosome","Marker position (bp)","Mean","Sigma2","Sigma2_k","SNP effect (mrMLM)","Sigma2_k_posteriori","Wald","'-log10(P) (mrMLM)'","Genotype for code 1","Significance")
}
rm(genRaw)
gc()
gg<-as.matrix(gg)
misfit<-numeric()
kk<- numeric()
kk0<- numeric()
l0<- numeric()
bong<-no_porderrow
if (bong>0){
g0<-gg[1:no_porderrow,1]
g0<-as.matrix(g0)
kk0<-no_porderrow
no_porderrow<-which(g0>0)
no_porderrow<-as.matrix(no_porderrow)
g0<-g0[no_porderrow,1]
g0<-as.matrix(g0)
xxx0<-genq[c(g0),]
if(dim(g0)[1]==1){
xxx0<-as.matrix(xxx0)
}
if(dim(g0)[1]>1)
{
xxx0<-as.matrix(xxx0)
xxx0<-t(xxx0)
}
phe<-as.matrix(phe)
if((flagps==1)||(exists("psmatrix")==FALSE))
{
par<-likelihood(matrix(1,(nrow(xxx0)),1),xxx0,phe,bbo=NULL)
lod<-par
}else if(flagps==0)
{
temp<-cbind(matrix(1,(nrow(xxx0)),1),psmatrix)
par<-likelihood(temp,xxx0,phe,bbo=NULL)
lod<-par
}
kk<-which(lod>=1.5)
kk<-as.matrix(kk)
kk1<-which(lod<1.5)
kk1<-as.matrix(kk1)
if ((nrow(kk1))>0){
misfit<-g0[kk1,1]
misfit<-as.matrix(misfit)
}
if ((nrow(kk))>0){
g0<-as.matrix(g0)
g0<-g0[kk,1]
xx0<-xxx0[,kk]
lo<-lod[kk,1]
}
if ((nrow(kk))==0){kk<-0}
}
if (bong==0){
kk0<-0
kk<-0
}
nleft<-as.matrix(gg[(kk0+1):(nrow(gg)),1])
if ((length(misfit))>0){gg<-rbind(nleft,misfit)}
if ((length(misfit))==0){gg<-nleft}
a1<-which(gg>0)
a1<-as.matrix(a1)
a2<-gg[a1,1]
a2<-as.matrix(a2)
if(nrow(a2)>1){
xx<-t(genq[c(a2),])
}else{
xx<-genq[c(a2),]
}
xx<-as.matrix(xx)
if((flagps==1)||(exists("psmatrix")==FALSE))
{
if (length(kk)>1){xin<-cbind(matrix(1,(nrow(xx)),1),xx0)}
if (length(kk)==1){
if(kk==0){
xin<- matrix(1,(nrow(xx)),1)
}
if(kk>0){
xin<-cbind(matrix(1,(nrow(xx)),1),xx0)
}
}
}else if(flagps==0)
{
temp<-cbind(matrix(1,(nrow(xx)),1),psmatrix)
if (length(kk)>1){xin<-cbind(temp,xx0)}
if (length(kk)==1){
if(kk==0){
xin<-temp
}
if(kk>0){
xin<-cbind(temp,xx0)
}
}
}
xin<-as.matrix(xin)
par1<-ebayes_EM(xin,xx,phe)
par<-par1$wang
w2<-which(par[,1]<=0.01)
if(length(w2)!=0){
w2<-as.matrix(w2)
ww<- numeric()
if ((nrow(w2))>0){
orderno<-a2[w2,1]
orderno<-as.matrix(orderno)
x3<-cbind(xin,xx[,w2])
x3<-as.matrix(x3)
lodfix<-matrix(x3[,1],nrow(x3),)
lodrand<-matrix(x3[,2:(ncol(x3))],nrow(x3),)
if((flagps==1)||(exists("psmatrix")==FALSE))
{
lod<-likelihood(lodfix,lodrand,phe,bbo=NULL)
}else if(flagps==0)
{
temp<-cbind(psmatrix,lodfix)
lod<-likelihood(temp,lodrand,phe,bbo=NULL)
}
w3<-which(lod[,1]>=(svmlod))
w3<-as.matrix(w3)
if ((kk[1])>0){
g0<-as.matrix(g0)
orderno<-rbind(g0,orderno)
orderno<-as.matrix(orderno)
}
#if ((nrow(w3))==0){ww<-0}change20190125
if ((nrow(w3)!=0)&&(w3[1]>0)){
if((flagps==1)||(exists("psmatrix")==FALSE))
{
lo<-lod[w3,1]
ww<-orderno[w3,]
}else if(flagps==0)
{
lo<-lod[w3,1]
no_loc<-w3-ncol(psmatrix)
ww<-orderno[no_loc,]
}
}
}
if ((nrow(w2))==0){
g0<-as.matrix(g0)
lo<-as.matrix(lo)
yang<-which(lo>=(svmlod))
yang<-as.matrix(yang)
if ((nrow(yang))>0){
ww<-g0[yang,1]
lo<-lo[yang,1]
}
#if ((nrow(yang))==0){ww<-0}
}
#ww<-as.matrix(ww)
needww<-ww
if (length(ww)>=1){
#ww<-as.matrix(ww)
if (length(ww)>1){
ww<-as.matrix(ww)
if((flagps==1)||(exists("psmatrix")==FALSE))
{
ex<-cbind(matrix(1,(nrow(xx)),1),t(genq[c(ww),]))
}else if(flagps==0)
{
ex<-cbind(cbind(matrix(1,(nrow(xx)),1),psmatrix),t(genq[c(ww),]))
}
}else{
if((flagps==1)||(exists("psmatrix")==FALSE))
{
ex<-cbind(matrix(1,(nrow(xx)),1),as.matrix(genq[c(ww),]))
}else if(flagps==0)
{
ex<-cbind(cbind(matrix(1,(nrow(xx)),1),psmatrix),as.matrix(genq[c(ww),]))
}
}
rm(genq)
gc()
ex<-as.matrix(ex)
cui<-det(t(ex)%*%ex)
p1<-rep(1,ncol(ex))
p2<-diag(p1)
if (cui<1e-6){bbbb<-solve(t(ex)%*%ex+p2*0.01)%*%t(ex)%*%phe}
if (cui>=1e-6){ bbbb<-solve(t(ex)%*%ex)%*%t(ex)%*%phe }
if((flagps==1)||(exists("psmatrix")==FALSE))
{
eeff<-bbbb[2:(nrow(bbbb)),1]
}else if(flagps==0)
{
eeff<-bbbb[(2+ncol(psmatrix)):(nrow(bbbb)),1]
}
eeff<-as.matrix(eeff)
er<-as.numeric()
her<-as.numeric()
if((flagps==1)||(exists("psmatrix")==FALSE))
{
excol<-ncol(ex)
for(i in 1:(excol-1))
{
em<-ex[,(1+i)]
as1<-length(which(em==1))/nrow(ex)
as2<-1-as1
er<-rbind(er,(1-(as1-as2)*(as1-as2))*eeff[i]*eeff[i])
}
v0<-(1/(nrow(ex)-1))*(t(phe-ex%*%bbbb)%*%(phe-ex%*%bbbb))
if(var(phe)>=(sum(er)+v0)){
her<-(er/as.vector(var(phe)))*100
}else{
her<-(er/as.numeric(sum(er)+v0))*100
}
}else if(flagps==0)
{
excol<-ncol(ex)
for(i in 1:(excol-1-ncol(psmatrix)))
{
em<-ex[,(1+ncol(psmatrix)+i)]
as1<-length(which(em==1))/nrow(ex)
as2<-1-as1
er<-rbind(er,(1-(as1-as2)*(as1-as2))*eeff[i]*eeff[i])
}
v0<-(1/(nrow(ex)-1))*(t(phe-ex%*%bbbb)%*%(phe-ex%*%bbbb))
if(var(phe)>=(sum(er)+v0)){
her<-(er/as.vector(var(phe)))*100
}else{
her<-(er/as.numeric(sum(er)+v0))*100
}
}
vee<-round(v0,4)
pee<-round(var(y),4)
if(nrow(her)>1){
vees<-matrix("",nrow = nrow(her),1)
pees<-matrix("",nrow = nrow(her),1)
pees[1,1]<-pee
vees[1,1]<-vee
}else{
pees<-as.matrix(pee)
vees<-as.matrix(vee)
}
x<-deepcopy(gen,,3:nrow(gen))
xxxx<-as.matrix(x[,ww])
rm(x)
gc()
xxmaf<-t(xxxx)
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))
eeff[which(abs(eeff)>=1e-4)] <- round(eeff[which(abs(eeff)>=1e-4)],4)
eeff[which(abs(eeff)<1e-4)] <- as.numeric(sprintf("%.4e",eeff[which(abs(eeff)<1e-4)]))
lo[which(abs(lo)>=1e-4)] <- round(lo[which(abs(lo)>=1e-4)],4)
lo[which(abs(lo)<1e-4)] <- as.numeric(sprintf("%.4e",lo[which(abs(lo)<1e-4)]))
her[which(abs(her)>=1e-4)] <- round(her[which(abs(her)>=1e-4)],4)
her[which(abs(her)<1e-4)] <- as.numeric(sprintf("%.4e",her[which(abs(her)<1e-4)]))
log10P <- as.matrix(-log10(pchisq(lo*4.605,1,lower.tail = F)))
log10P[which(abs(log10P)>=1e-4)] <- round(log10P[which(abs(log10P)>=1e-4)],4)
log10P[which(abs(log10P)<1e-4)] <- as.numeric(sprintf("%.4e",her[which(abs(log10P)<1e-4)]))
if (length(ww)>1){
wan<-data.frame(parmsShow[needww,1],chr_pos[ww,],eeff,lo,log10P,her,maf,parmsShow[needww,11])
wan<-wan[order(wan[,2]),]
wan<-data.frame(wan,vees,pees)
}else{
wan<-data.frame(parmsShow[needww,1],t(as.matrix(chr_pos[ww,])),eeff,lo,log10P,her,maf,parmsShow[needww,11],vees,pees)
}
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)")
}
}
if(is.null(parmsShow)==FALSE){
parmsShow<-parmsShow[,-c(4,5,6,8,9,12)]
}
output<-list(result1=parmsShow,result2=wan)
return(output)
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.