## Demonstration
require(OpenMx)
source("GW-SEM.R")

#Run the Ordinal 1 Factor Residuals GWAS

snpCovs(FacModelData = "Ord1fac.txt", vars = c("i1","i2","i3","i4","i5"), controls = c("c1","c2","c3"), SNPdata = "SNP.txt", output = "Ord1fac")
resDWLS(itemData = "Ord1fac.txt", snpCov = "Ord1faccov.txt", snpWei = "Ord1facwei.txt", VarNames =  c("i1","i2","i3","i4","i5"), Controls = c("c1","c2","c3"), resids = c("i1","i2","i3","i4","i5"), output = "OrdResOut.txt")
gwasDWLS(itemData= "Ord1fac.txt", snpCov = "Ord1faccov.txt", snpWei = "Ord1facwei.txt", VarNames = c("i1","i2","i3","i4","i5"), Controls = c("c1","c2","c3"), output = "Ord1Facout.txt")

# Read the output into R and look at the first few SNPs
resultsFac <- read.table("Ord1Facout.txt")
resultsRes <- read.table("OrdResOut.txt")
head(resultsFac)
head(resultsRes)


# Make a Manhattan Plot for the 1 Factor SEM
resultsFac <- read.table("Ord1Facout.txt")
par(mar = c(4, 5, 2, 2) + 0.1)
plot(-log10(resultsFac$p.value), ylim = c(0, 9), ylab = "One Factor GW-SEM  \n -Log10 (P value)", xlab = "SNP", 
     main = "Manhattan Plot of a One Factor GW-SEM", pch = 16)
abline(h = 5, col = "blue", lwd = 2)
abline(h = -log10(5 * 10^-8), col = "red", lwd = 2)


# Make a Manhattan Plot for the 1 Factor Residuals GWAS
resultsRes <- read.table("OrdResOut.txt")

par(mfrow = c(5,1), mar = c(0, 5, 0, 0) + 0.1)
plot(-log10(resultsRes$p_i1), ylim = c(0, 9), ylab = "i1  \n -Log10 (P value)", 
     main = " ", xlab = " ", xaxt = "n", pch = 16)
abline(h = 5, col = "blue", lwd = 2)
abline(h = -log10(5 * 10^-8), col = "red", lwd = 2)

plot(-log10(resultsRes$p_i2), ylim = c(0, 9), ylab = "i2  \n -Log10 (P value)", 
     main = " ", xlab = " ", xaxt = "n", pch = 16)
abline(h = 5, col = "blue", lwd = 2)
abline(h = -log10(5 * 10^-8), col = "red", lwd = 2)

plot(-log10(resultsRes$p_i3), ylim = c(0, 9), ylab = "i3  \n -Log10 (P value)", 
     main = " ", xlab = " ", xaxt = "n", pch = 16)
abline(h = 5, col = "blue", lwd = 2)
abline(h = -log10(5 * 10^-8), col = "red", lwd = 2)

plot(-log10(resultsRes$p_i4), ylim = c(0, 9), ylab = "i4  \n -Log10 (P value)", 
     main = " ", xlab = " ", xaxt = "n", pch = 16)
abline(h = 5, col = "blue", lwd = 2)
abline(h = -log10(5 * 10^-8), col = "red", lwd = 2)

plot(-log10(resultsRes$p_i5), ylim = c(0, 9), ylab = "i5  \n -Log10 (P value)", 
     main = " ", xlab = " ", xaxt = "n", pch = 16)
abline(h = 5, col = "blue", lwd = 2)
abline(h = -log10(5 * 10^-8), col = "red", lwd = 2)


# Select the SNPs with the lowest p-values
resultsFac <- read.table("Ord1Facout.txt")
head(resultsFac[order(resultsFac$p.value),], 10)

resultsRes <- read.table("OrdResOut.txt")
resultsRes["rs5218",]


# Examine the Estimates from the one factor GW-SEM and the residuals GW-SEM
oneFacEst  <- read.table("estimates.txt")
resEst     <- read.table("resEstimates.txt")

oneFacEst["rs5218",]
resEst["rs5218",]



# Make a QQ plot for the 1 Factor GW-SEM

require(car)
resultsFac <- read.table("Ord1Facout.txt")
qqPlot(resultsFac$t.statistic, grid = F, main = "QQ-Plot for the One Factor GW-SEM", 
       ylab = "Sample Quantiles", xlab = "Theoretical Quantiles")

# Make a QQ plot for the Residuals GW-SEM
resultsRes <- read.table("OrdResOut.txt")
par(mfrow= c(2,3))
qqPlot(resultsRes$t_i1, grid = F, main = "QQ-Plot for Item 1", 
       ylab = "Sample Quantiles", xlab = "Theoretical Quantiles")
qqPlot(resultsRes$t_i2, grid = F, main = "QQ-Plot for Item 2", 
       ylab = "Sample Quantiles", xlab = "Theoretical Quantiles")
qqPlot(resultsRes$t_i3, grid = F, main = "QQ-Plot for Item 3", 
       ylab = "Sample Quantiles", xlab = "Theoretical Quantiles")
qqPlot(resultsRes$t_i4, grid = F, main = "QQ-Plot for Item 4", 
       ylab = "Sample Quantiles", xlab = "Theoretical Quantiles")
qqPlot(resultsRes$t_i5, grid = F, main = "QQ-Plot for Item 5", 
       ylab = "Sample Quantiles", xlab = "Theoretical Quantiles")