# ========================================================================
# Example 
# ========================================================================
library(matrixStats)
library(mvtnorm)
library(MASS)
source("summaryaspu.r")


## Examples
# -- Simulating data
# -- n.snp: number of SNPs
# -- n.trait: number of traits
# -- n.subject: number of subjects
set.seed(136)
n.snp <- 1000
n.traits <- 10
n.subjects <- 500
traits <- matrix(rnorm(n.subjects*n.traits), n.subjects, n.traits)
v <- cov(traits)
allZ <- rmvnorm(n.snp, sigma=v)
colnames(allZ) <- paste("trait", 1:n.traits, sep="")
rownames(allZ) <- paste("snp", 1:n.snp, sep="")
head(allZ)









# =========================================================================
# estimating covariance across multiple traits based on the NULL SNPs
# This v is applied to minP, chi-squared, aSPU test.
v <- estcov(allZ)






# =========================================================================
# applying methods
# =========================================================================
# minP test
# =========================================================================

n.snp <- dim(allZ)[1]
minP_pval <- numeric(n.snp)
for (i in 1:n.snp){
	minP_pval[i] = minP(Zi = allZ[i,], v = v)
}

head(minP_pval)


# =========================================================================
# chi-squared test
# =========================================================================

n.snp <- dim(allZ)[1]
chi_pval <- numeric(n.snp)
for (i in 1:n.snp){
	u0 = allZ[i,]; 
	Ts = sum(u0 * ginv(v) %*% u0)
	chi_pval[i] = pchisq(Ts, df = length(u0),lower = FALSE)
}

head(chi_pval)



# ========================================================================
# aSPU 
# gradually increase permutation number B=1e3, 1e4, upto 1e8
# apply aSPU with larger B to SNPs only when the p-value less than 5/B
# Directly applying B=10^8 requires large memory for 2.6m SNPs
# If have parallel cores, it would be better to split data in parallels
# ========================================================================
n.snp <- dim(allZ)[1]
mtaspuS_pval <- numeric(n.snp)


##  Screening steps
##  In each step, summary Z statistics for single SNP is applied to the 
##  funtion MTaSPUs()


##  B = 1e3
n.snp <- dim(allZ)[1]
B <- 1e3
aspuB1e3 <- numeric(n.snp)
for (i in 1:n.snp){
	P <-MTaSPUs(allZ[i,], v = v, B = B, pow = c(1:8, Inf), transform = FALSE)
	aspuB1e3[i] <-P[,"MTaSPUs"]
}

##  update aspu pvalue with B=1e3
mtaspuS_pval <- aspuB1e3





##  B = 1e4
subset <- matrix(allZ[mtaspuS_pval < 5e-3,], ncol = ncol(allZ))
n.snp2 <- dim(subset)[1]
B <- 1e4
aspuB1e4 <- numeric(n.snp2)
for (i in n.snp2){
	P <- MTaSPUs(subset[i,], v = v, B = B, pow = c(1:8, Inf), transform = FALSE)
	aspuB1e4[i] <- P[,"MTaSPUs"]
}
##  update aspu pvalue with B=1e4
mtaspuS_pval[mtaspuS_pval < 5e-3] <- aspuB1e4





##  B = 1e5
subset <- matrix(allZ[mtaspuS_pval < 5e-4,], ncol=ncol(allZ))
n.snp3 <- dim(subset)[1]
B <- 1e5
aspuB1e5 <- numeric(n.snp3)
for (i in n.snp3){
	P <- MTaSPUs(subset[i,], v = v, B = B, pow = c(1:8, Inf), transform = FALSE)
	aspuB1e5[i] <- P[,"MTaSPUs"]
}
##  update aspu pvalue with B=1e5
mtaspuS_pval[mtaspuS_pval < 5e-4] <- aspuB1e5





##  B = 1e6
subset <- matrix(allZ[mtaspuS_pval < 5e-5,], ncol=ncol(allZ))
n.snp4 <- dim(subset)[1]
B <- 1e6
aspuB1e6 <- numeric(n.snp4)
for (i in n.snp4){
  P <- MTaSPUs(subset[i,], v = v, B = B, pow = c(1:8, Inf), transform = FALSE)
  aspuB1e6[i] <- P[,"MTaSPUs"]
}
##  update aspu pvalue with B=1e6
mtaspuS_pval[mtaspuS_pval < 5e-5] <- aspuB1e6





##  Keep increasing B
##  Apply chank of summary Z statistics for multiple SNPs is applied to the
##  where each columns are individual trait; rows for single SNP
##  funtion MTaSPUs()
##  P are results
##  each row represents each snp, and each column represents 
##  the pvalues obtained from the SPU and aSPU tests.


##  B = 1e7
subset <- matrix(allZ[mtaspuS_pval < 5e-6,], ncol=ncol(allZ))
B <- 1e7
P <- MTaSPUs(subset, v = v, B = B, pow = c(1:8, Inf), transform = FALSE)
##  update aspu pvalue with B=1e6
mtaspuS_pval[mtaspuS_pval < 5e-6] <- P[,"MTaSPUs"]





# -- In this example no SNP has p value of < 5e-7
# no run 
# subset <- matrix(allZ[mtaspuS_pval < 5e-7, ], ncol=ncol(allZ))
# P <-  MTaSPUs(subset, v = v, B = B, pow = c(1:8, Inf), transform = FALSE)
# mtaspuS_pval[mtaspuS_pval < 5e-7] <- P[,"MTaSPUs"]



##  Summarize results
names(chi_pval) <- rownames(allZ)
names(minP_pval) <- rownames(allZ)
names(mtaspuS_pval) <- rownames(allZ)

head(mtaspuS_pval)