library(boot)
library(MASS)
library(base)
library(stats)


##################################################################                                                                                                                                                                                    
##     TOWA is an association test to test association        ##
##     between a phenotype and an optimally weighted          ##
##     combination of variants (can be rare and common)       ##                                              
##     of a genomic region in admixed populations.            ##                                   
##     y represents disease status (1=case;0=control).        ##   
##     a1 represents the allele matrix for the first allele   ##
##     of SNPs. a2 represents the allele matrix for the       ##
##     second allele of SNPs. p1 represents the ancestry      ##
##     matrix of ancestral population 1 for the first allele  ##
##     of SNPs. p2 represents the ancestry matrix of          ##
##     ancestral population 1 for the second allele of SNPs.  ## 
##     Each row of a1, a2, p1, and p2 represents an           ##
##     individual and each column represents a SNP. numper    ## 
##     represents the number of permutation times to          ## 
##     obtain a p value of the test TOWA.                     ##                                                                                                    
##                                                            ##
################################################################

TOWA=function(y,a1,a2,p1,p2,numper)
{
     ee=0.000001
     m=ncol(a1)
     n=nrow(a1)
     i=matrix(rep(1,m*n), nrow = n, ncol = m, byrow = TRUE)
     p11=i-p1
     p22=i-p2
     x=a1*p1+a2*p2
     z=a1*p11+a2*p22	
     ############################
     v1=(n-1)*apply(x,2,var)+ee 
     v2=(n-1)*apply(z,2,var)+ee 
     u1=(y-mean(y))%*%x
     u2=(y-mean(y))%*%z
     u1v1=(u1^2)/v1
     u2v2=(u2^2)/v2
     uv12=u1v1+u2v2
	T1=sum(uv12)
	#############################
	T1per=rep(0,numper)
	for(i in 1:numper)
		{
		  yper=sample(y)
                  u1=(yper-mean(yper))%*%x
                  u2=(yper-mean(yper))%*%z
                  u1v1=(u1^2)/v1
                  u2v2=(u2^2)/v2
                  uv12=u1v1+u2v2
	          T1per[i]=sum(uv12)
		}
	##############################	
pv=sum(T1per[1:numper]>T1)/numper+sum(T1per[1:numper]==T1)/numper/2
     pvalue=pv
}


##################################################################                                                            
##                                                            ##                                                             
##     TOWAs is TOWA used to test association in each         ##
##     ancestral population.                                  ##                                                                                                                                                 
##     y represents disease status (1=case;0=control).        ## 
##     a1 represents the allele matrix for the first allele   ##
##     of SNPs. a2 represents the allele matrix for the       ##
##     second allele of SNPs. p1 represents the ancestry      ##
##     matrix of ancestral population 1 for the first allele  ##
##     of SNPs. p2 represents the ancestry matrix of ancestral##
##     population 1 for the second allele of SNPs. Each row   ##
##     of a1, a2, p1, and p2 represents an individual and     ##
##     each column represents a SNP. numper represents the    ##
##     number of permutation times to obtain a p value of     ##
##     the test TOWAs.                                        ##                                                                                                                                                                                                                                                                                                 
################################################################



TOWAs=function(y,a1,a2,p1,p2,numper)
{
	ee=0.000001
	m=ncol(a1)
	n=nrow(a1)
     i=matrix(rep(1,m*n), nrow = n, ncol = m, byrow = TRUE)
     p11=i-p1
     p22=i-p2
     x=a1*p1+a2*p2
     z=a1*p11+a2*p22	
	############################
     v1=(n-1)*apply(x,2,var)+ee 
     v2=(n-1)*apply(z,2,var)+ee 
     u1=(y-mean(y))%*%x
     u2=(y-mean(y))%*%z
     u1v1=(u1^2)/v1
     u2v2=(u2^2)/v2
     T1=sum(u1v1)
     T2=sum(u2v2)
     #############################
     T1per=rep(0,numper)
     T2per=rep(0,numper)
	for(i in 1:numper)
		{
		  yper=sample(y)
                  u1=(yper-mean(yper))%*%x
                  u2=(yper-mean(yper))%*%z
                  u1v1=(u1^2)/v1
                  u2v2=(u2^2)/v2
                  T1per[i]=sum(u1v1)
                  T2per[i]=sum(u2v2)
		}
	##############################		
pv1=sum(T1per[1:numper]>T1)/numper+sum(T1per[1:numper]==T1)/numper/2
pv2=sum(T2per[1:numper]>T2)/numper+sum(T2per[1:numper]==T2)/numper/2
     pvalue=c(pv1,pv2)
}