# This is written in R # The input file "dat" is a matrix # The first column of dat are the phenotype values # Each of the other columns are the coded genotypes at a SNP # The number of SNPs is supposed to be greater than 3 # The output is the p-value calculated by 10,000 permutations # Before running the program, the following should be included # library(mvtnorm) # library(wavethresh) # library(EbayesThresh) function(dat){ np=4 x=dat[,-1] y=dat[,1] n=dim(x)[1] m=dim(x)[2] if(m<4){ print("The number of SNPs must be greater than 3") } mind=(log(m)/log(2)) if(floor(log(m)/log(2))==mind){ Ans=matrix(0,n,m) for (i in 1:n){ x1=as.vector(x[i,1:m]) x2=wd(x1,family="DaubLeAsymm",filter.number=8) if(sum(x1)==0) x3=x2 else{ x3=ebayesthresh.wavelet(x2,vscale="level") } x4=wr(x3) Ans[i,]=x4 } Xt=Ans U=numeric(m) V=numeric(m) for ( j in 1:m){ xmean=mean(Xt[,j]) U=c(U,sum(y*(Xt[,j]-xmean))) V=c(V,((n/(4*(n-1)))*sum((Xt[,j]-xmean)^2))) } sumv=sum(V) if(sumv>0) t=((sum(U))^2)/(sumv) else t=0 tmax=t nsim=(10^np) tx=numeric(nsim) ic=0 for( i in i:nsim){ U=numeric(m) V=numeric(m) y1=sample(y) for ( j in 1:m){ xmean=mean(Xt[,j]) U=c(U,sum(y1*(Xt[,j]-xmean))) V=c(V,(n/(4*(n-1))*sum((Xt[,j]-xmean)^2))) } sumv=sum(V) if(sumv>0) tx[i]=((sum(U))^2)/sumv else tx[i]=0 if(tx[i]>tmax) ic=(ic+1) } pval=(ic/nsim) return(pval) } else{ m1=floor(log(m)/log(2)) m1=2^m1 Xt=matrix(0,n,m) Ans=matrix(0,n,m1) for (i in 1:n){ x1=as.vector(x[i,1:m1]) x2=wd(x1,family="DaubLeAsymm",filter.number=8) if(sum(x1)==0) x3=x2 else{ x3=ebayesthresh.wavelet(x2,vscale="level") } x4=wr(x3) Ans[i,]=x4 } Xt[,(1:m1)]=Ans Xt[,((m1+1):m)]=x[,((m1+1):m)] U=numeric(m) V=numeric(m) for ( j in 1:m){ xmean=mean(Xt[,j]) U=c(U,sum(y*(Xt[,j]-xmean))) V=c(V,((n/(4*(n-1)))*sum((Xt[,j]-xmean)^2))) } sumv=sum(V) if(sumv>0) t=((sum(U))^2)/(sumv) else t=0 tmax=t nsim=(10^np) tx=numeric(nsim) ic=0 for( i in i:nsim){ U=numeric(m) V=numeric(m) y1=sample(y) for ( j in 1:m){ xmean=mean(Xt[,j]) U=c(U,sum(y1*(Xt[,j]-xmean))) V=c(V,(n/(4*(n-1))*sum((Xt[,j]-xmean)^2))) } sumv=sum(V) if(sumv>0) tx[i]=((sum(U))^2)/sumv else tx[i]=0 if(tx[i]>tmax) ic=(ic+1) } pval=(ic/nsim) return(pval) } }