##################################################################
#############            Read DNA Sequences                ########
##################################################################

DNAread<-function(DNAdata)
{
	
	file1=scan(file=DNAdata,what="",sep="@",quiet=TRUE)
        file2=scan(file=DNAdata,what="",sep="",quiet=TRUE)
	R.no=as.numeric(file2[[1]])
	C.no=as.numeric(file2[[2]])
	Dname=matrix(0,R.no,1)
	Data1=matrix(0,R.no,C.no)
	for(i in 2:(R.no+1)){
	aa=c(unlist(strsplit(file1[[i]],NULL)))
	Dname[i-1,]=paste(aa[1],aa[2],aa[3],aa[4],aa[5],aa[6],aa[7],aa[8],aa[9],aa[10],sep="")
	Data1[i-1,]=aa[-(1:10)]
		}
	Data2=matrix(Data1,R.no,C.no,dimnames = list(c(Dname), 1:(C.no)))
	Data2
	
dd=NULL
for (i in 1:R.no)
{
	for (j in 1:C.no)

{
	if(c(Data2[[i,j]]!="A"&&Data2[[i,j]]!="C"&&Data2[[i,j]]!="G"&&Data2[[i,j]]!="T"))
		{dd=cbind(dd,j)}
}}
if(length(dd!=0)){dd=sort(dd)}
dd3=NULL
L1=length(dd)
if(L1 > 1) {
	
	
	for (i in 1:(L1-1)) {
		if((dd[i ] - dd[(i+1)]) != 0)
			{dd3 = cbind(dd3,dd[i])}
	}
	
	dd3=cbind(dd3,dd[L1])
	
}
else {dd3 = dd}
if(length(dd3!=0)){dd3=-1*dd3}
if((length(dd3)!=0)){Data3=Data2[,dd3]}
	else{Data3=Data2}
Data3=as.matrix(Data3)
if(dim(Data3)[[2]]==0){print("All the sites have at least one missing nucleotide")}else {
C1.no=dim(Data3)[[2]]	
Data3=matrix(Data3,R.no,dimnames = list(c(Dname), 1:(C1.no)))
Data3}

}
##################################################################
#############  Test for Symmetry of Matched DNA Sequences   ######
##################################################################

Testpairs <- function(Fmatrix)
{
	f=Fmatrix
	n <- dim(f)[[1]]
	B1 <- NULL
	S1 <- NULL
	I1 <- NULL
	for(i in 1:n) {
		B <- matrix(0, n, 1)
		S <- matrix(0, n, 1)
		I <- matrix(0, n, 1)
		for(j in 1:n) {
			if(j != i) {
				k <- table(f[i,  ], f[j,  ])
				testr <- TEST(k)$pvals
			}
			else {
				testr <- matrix(0, 3, 1)
			}
			B[j] <- testr[1]
			S[j] <- testr[2]
			I[j] <- testr[3]
		}
		
		
		B1 <- cbind(B1, B)
		S1 <- cbind(S1, S)
		I1 <- cbind(I1, I)
	}
    colnames(I1)=colnames(S1)=colnames(B1)=rownames(f)
	list(Bowker.test = as.dist(round(B1, 4)), Stuart.test = as.dist(round(S1, 4)), Internal.test = as.dist(round(I1, 4)))
}

TEST=function(Nt){
	QB <- 0
	QS <- 0
	QR <- 0
	PB <- 0
	PR <- 0
	PS <- 0
	ZB <- 0
	ZS <- 0
	ZR <- 0
	V <- matrix(0, 4, 4)
	d <- matrix(0, 4, 1)
	zcount=0
	zzcount=0
	for(i in 1:4) {
		for(j in 1:4) {
			if((i < j))
			if(Nt[i, j] + Nt[j, i] == 0) {
				zcount=zcount+1
			}
			else {
				QB<-QB+(((Nt[i, j]-Nt[j, i])^2)/(Nt[i, j] + Nt[j, i]))
			}
		}
	}
	for(i in 1:4) {
		d[i] <- (sum(Nt[i,  ]) - sum(Nt[, i]))
		for(j in 1:4) {
			if(i == j) V[i, j]<-(sum(Nt[i,  ]) + sum(Nt[, i]) - 2 * sum(Nt[i, i])) else V[i, j]<--1*(Nt[i, j] + Nt[j, i])
		}
	}
        r=qr(V)$rank
        if (r==0) QS=NA else{if (r==1){
             L<-1/eigen(V)$values[1]
	   U<-eigen(V)$vectors[, 1]
           QS<-c(t(d)%*%U%*%L%*%t(U)%*%d)} else{
           L<-diag(1/eigen(V)$values[1:r])
	   U<-eigen(V)$vectors[, 1:r]
           QS<-c(t(d)%*%U%*%L%*%t(U)%*%d)}}	
        QR <- QB - QS
	if (6-zcount==0) PB=NA else PB <- 1 - pchisq(QB, 6-zcount)
	if (r==0) PS=NA else PS <- 1 - pchisq(QS, r)
        if (6-zcount-r==0) PR=NA else PR <- 1 - pchisq(QR, 6-zcount-r)
        list(pvals=matrix(c(PB, PS, PR),3,1),df=matrix(c(6-zcount,r,6-zcount-r),3,1))
}




##################################################################
#############      Overall Test for Marginal Symmetry     ########
##################################################################
Testoverall<-function(Fm)
{
    v<-NULL
    m <- dim(Fm)[1]
    r <- 4
    n <- NULL
    one <- rep(1, r * (m - 1))
    
    L1 <- t(matrix(t(rep(diag(c(1, 1, 1, 1)), (m - 1))), r,
        (r * (m - 1))))
    L <- cbind(L1,  - (diag(one)))
    for(i in 1:m) {
        vc1 <- NULL
        for(k in 1:m) {
            vc <- NULL
            if(i == k) {
                fi <- table(Fm[i,])
                fi= fi/sum(fi)
                vc <- diag(fi)
                n1 <- matrix(fi, r, 1)
                n <- rbind(n, n1)
            }
            if(i < k) {
                fij <- table(Fm[i,],Fm[k,])
                fij = fij/sum(fij)
                fr <- apply(fij, 1, sum)
                fc <- apply(fij, 2, sum)
                vc <- fij
            }
            if(i > k) {
                fij <- table(Fm[i,],Fm[k,])
                fij = fij/sum(fij)
                fr <- apply(fij, 1, sum)
                fc <- apply(fij, 2, sum)
                vc <- fij
            }
            vc1 <- cbind(vc1, vc)
        }
        v <- rbind(v, vc1)
    }
    rk=qr(L%*%v%*%t(L))$rank
    ev=eigen(L%*%v%*%t(L))
    
    Ts <- (t(L %*% n) %*% ev$vectors[,1:rk]%*%diag(1/ev$values[1:rk])%*%t(ev$vectors[,1:rk]) %*% (
        L %*% n)) * dim(Fm)[2]
    pval<-1-pchisq(c(Ts),rk)
    list(chisqstat=c(Ts),df=rk,pvalue=pval)
}