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

#Run the Ordinal 1 Factor Model GWAS

snpCovs(FacModelData = "Ord1fac.txt", vars = c("i1","i2","i3","i4","i5"), controls = c("c1","c2","c3"), SNPdata = "SNP.txt", output = "Ord1fac")
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")
end <- Sys.time()


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


# Make a Manhattan Plot for the 1 Factor SEM
results <- read.table("Ord1Facout.txt")
plot(-log10(results$p.value), ylim = c(0, 9), ylab = "-log10 of the 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)



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


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