## ---- echo = FALSE, message=FALSE---------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
library("PhenotypeSimulator")
library("ggplot2")

## ---- eval=FALSE--------------------------------------------------------------
#  # load kinship data, add small value to diagonal for numerical stability and
#  # do an eigendecomposition
#  indir <- "~/data/hmeyer/Supplementary/CEU.0908.impute.files"
#  
#  kinship <- read.table(paste(indir,
#                              "/genotypes_genome_hapgen.controls.grm.rel",sep=""),
#                        sep="\t", header=FALSE)
#  
#  kinship <- as.matrix(kinship)
#  diag(kinship) <- diag(kinship) + 1e-4
#  
#  kinship_decomposed <- eigen(kinship)
#  write.table(kinship,
#  			paste(indir, "/genotypes_genome_hapgen.controls.grm.rel",
#  				  sep=""),
#  			sep="\t", col.names=FALSE, row.names=FALSE)
#  
#  write.table(kinship_decomposed$values,
#  			paste(indir, "/genotypes_eigenvalues",
#  				  sep=""),
#  			sep="\t", col.names=FALSE, row.names=FALSE)
#  
#  write.table(kinship_decomposed$vectors,
#  			paste(indir, "/genotypes_eigenvectors",
#  				  sep=""),
#  			sep="\t", col.names=FALSE, row.names=FALSE)

## ---- eval=FALSE--------------------------------------------------------------
#  indir <- "/homes/hannah/data/hmeyer/Supplementary/CEU.0908.impute.files"
#  
#  # specify directory to save data; if it doesn't exist yet, create i.e. needs
#  # writing permissions
#  datadir <- '/homes/hannah/tmp/phenotypeSimulator'
#  if (!dir.exists(datadir)) dir.create(datadir, recursive=TRUE)
#  
#  # specify filenames and parameters
#  totalGeneticVar <- 0.4
#  totalSNPeffect <- 0.1
#  h2s <- totalSNPeffect/totalGeneticVar
#  kinshipfile <- paste(indir, "/genotypes_genome_hapgen.controls.grm.rel",
#                                  sep="")
#  
#  genoFilePrefix <- paste(indir, "/genotypes_", sep="")
#  genoFileSuffix <- "_hapgen.controls.gen"
#  
#  # simulate phenotype with three phenotype components
#  simulation <- runSimulation(N = 1000, P = 3, cNrSNP=10, seed=43,
#                             kinshipfile = kinshipfile,
#                             format = "oxgen",
#                             genoFilePrefix = genoFilePrefix,
#                             genoFileSuffix = genoFileSuffix,
#                             chr = c(5,7,11),
#                             mBetaGenetic = 0, sdBetaGenetic = 0.2,
#                             theta=1,
#                             NrSNPsOnChromosome=c(533966, 503129, 428863),
#                             genVar = totalGeneticVar, h2s = h2s,
#                             phi = 0.6, delta = 0.2, rho=0.2,
#                             NrFixedEffects = 2, NrConfounders = c(2, 2),
#                             distConfounders = c("bin", "norm"),
#                             probConfounders = 0.2,
#                             genoDelimiter=" ",
#                             kinshipDelimiter="\t",
#                             kinshipHeader=FALSE,
#                             verbose = TRUE )

## ---- eval=FALSE, echo=FALSE--------------------------------------------------
#  savedir <- paste(system.file("extdata", package="PhenotypeSimulator"),
#                              "/resultsSimulationAndLinearModel", sep="")
#  saveRDS(simulation$phenoComponentsFinal,
#          paste(savedir, "/simulation_phenoComponentsFinal.rds", sep=""))

## ---- echo=FALSE--------------------------------------------------------------
savedir <- paste(system.file("extdata", package="PhenotypeSimulator"), 
                            "/resultsSimulationAndLinearModel", sep="")
simulation <- list(
    phenoComponentsFinal=readRDS(paste(savedir, 
                                       "/simulation_phenoComponentsFinal.rds", 
                            sep="")))

## ----heatmaps, fig.keep='all', fig.height=3, fig.width=7, eval=TRUE, fig.cap="\\label{fig:heatmaps}Heatmaps of the trait-by-trait correlation (Pearson correlation) of the simulated phenotype and its five phenotype components."----
cor_phenotype <- reshape2::melt(cor(simulation$phenoComponentsFinal$Y))
cor_phenotype$type <- "Y" 

cor_genBg <- reshape2::melt(cor(simulation$phenoComponentsFinal$Y_genBg))
cor_genBg$type <- "U" 

cor_genFixed <- reshape2::melt(cor(simulation$phenoComponentsFinal$Y_genFixed))          
cor_genFixed$type <- "XB" 

cor_noiseFixed <- reshape2::melt(cor(simulation$phenoComponentsFinal$Y_noiseFixed))          
cor_noiseFixed$type <- "WA" 

cor_noiseBg <- reshape2::melt(cor(simulation$phenoComponentsFinal$Y_noiseBg))          
cor_noiseBg$type <- "Psi" 

cor_correlatedBg <- reshape2::melt(cor(simulation$phenoComponentsFinal$Y_correlatedBg))          
cor_correlatedBg$type <- "T" 

cor_components <- rbind(cor_phenotype, cor_genBg, cor_genFixed, cor_noiseBg, 
                        cor_noiseFixed, cor_correlatedBg)
cor_components$type <- factor(cor_components$type, levels=c("Y", 
                                                            "WA",
                                                            "XB", 
                                                            "U", 
                                                            "T",
															"Psi"),
													labels=c("bold(Y)", 
                                                            "bold(WA)",
                                                            "bold(XB)", 
                                                            "bold(U)", 
                                                            "bold(T)",
															"bold(Psi)"))
colnames(cor_components) <- c("TraitA", "TraitB", "Correlation", "type")

# For prettier plotting
cor_components$TraitA <- gsub("_", " ", cor_components$TraitA)
cor_components$TraitB <- gsub("_", " ", cor_components$TraitB)

p_corr <- ggplot(data=cor_components, aes(x=TraitA, y=TraitB, fill=Correlation))  
p_corr <- p_corr + geom_tile() +
    facet_grid(~ type, labeller = label_parsed) +
    scale_x_discrete(expand = c(0,0)) +
    scale_y_discrete(expand = c(0,0)) +
    scale_fill_gradient2(low="#F2AD00", high="#5BBCD6" , mid="white",
                    limits = c(-1,1), 
                    breaks = c(-1, -0.5, 0, 0.5, 1),
                    labels = c(-1, -0.5, 0, 0.5, 1),
                    guide=guide_colorbar(direction = "horizontal",
                                         title.position="top",
                                         title.hjust =0.5)) + 
    xlab(" ") +
    ylab(" ") +
    theme_bw() +
    coord_fixed() +
    theme(axis.text.x = element_text(angle=90, vjust=0.5),
			legend.position = "bottom",
			strip.background = element_rect(fill="white", color="white"),
            axis.text  =element_text(size=8),
            legend.title  =element_text(size=8),
            legend.text  =element_text(size=8),
            strip.text  =element_text(size=8))
print(p_corr)

## ---- echo=FALSE, eval=FALSE--------------------------------------------------
#  ggsave(plot=p_corr,
#          filename=paste(savedir, "/correlation_phenotype.pdf", sep=""),
#          units="mm", width=86, height=60)

## ---- eval=FALSE--------------------------------------------------------------
#  # Save phenotypes, kinship and non-genetic covariates in csv and
#  # GEMMA-specific format
#  outdirectory <- savePheno(simulation, directory = datadir,
#                              intercept_gemma=TRUE, format=c("csv", "gemma"),
#                              saveIntermediate=TRUE)

## ---- eval=FALSE--------------------------------------------------------------
#  datadir="/homes/hannah/tmp/phenotypeSimulator"
#  
#  causalSNPs <- fread(paste(datadir, "/SNP_NrSNP10.csv", sep=""),
#  					header=TRUE, sep=",", stringsAsFactors=FALSE,
#  					data.table=FALSE)
#  causalSNPnames <- colnames(causalSNPs)
#  
#  LMM_mt <- fread(paste(datadir, "/GWAS_gemma_LMM_mt_pvalues.txt", sep=""),
#  				data.table=FALSE, header=FALSE, stringsAsFactors=FALSE)
#  LMM_st_trait1 <- fread(paste(datadir, "/GWAS_gemma_LMM_st_trait1_pvalues.txt",
#                              sep=""),
#  				data.table=FALSE, header=FALSE, stringsAsFactors=FALSE)
#  LMM_st_trait2 <- fread(paste(datadir, "/GWAS_gemma_LMM_st_trait2_pvalues.txt",
#                              sep=""),
#  				data.table=FALSE, header=FALSE, stringsAsFactors=FALSE)
#  LMM_st_trait3 <- fread(paste(datadir, "/GWAS_gemma_LMM_st_trait3_pvalues.txt",
#                              sep=""),
#  				data.table=FALSE, header=FALSE, stringsAsFactors=FALSE)
#  
#  LMM_mt <- LMM_mt[order(LMM_mt[,2]),]
#  LMM_mt$expected <- stats::ppoints(nrow(LMM_mt))
#  LMM_mt$model <- "mvLMM"
#  LMM_mt$trait <- "all"
#  colnames(LMM_mt)[1:2] <- c("rsID", "observed")
#  
#  LMM_st_trait1 <- LMM_st_trait1[order(LMM_st_trait1[,2]),]
#  LMM_st_trait1$expected <- LMM_mt$expected
#  LMM_st_trait1$model <- "uvLMM"
#  LMM_st_trait1$trait <- "trait1"
#  colnames(LMM_st_trait1)[1:2] <- c("rsID", "observed")
#  
#  LMM_st_trait2 <- LMM_st_trait2[order(LMM_st_trait2[,2]),]
#  LMM_st_trait2$expected <- LMM_mt$expected
#  LMM_st_trait2$model <- "uvLMM"
#  LMM_st_trait2$trait <- "trait2"
#  colnames(LMM_st_trait2)[1:2] <- c("rsID", "observed")
#  
#  LMM_st_trait3 <- LMM_st_trait3[order(LMM_st_trait3[,2]),]
#  LMM_st_trait3$expected <- LMM_mt$expected
#  LMM_st_trait3$model <- "uvLMM"
#  LMM_st_trait3$trait <- "trait3"
#  colnames(LMM_st_trait3)[1:2] <- c("rsID", "observed")
#  
#  pvalues <- rbind(LMM_mt, LMM_st_trait1, LMM_st_trait2, LMM_st_trait3)
#  pvalues$causal <- "all"
#  pvalues$causal[which(pvalues$rsID %in% causalSNPnames)] <- "simulated causal"
#  pvalues$causal <- factor(pvalues$causal, levels=c("all", "simulated causal"))
#  
#  p_qq <- ggplot(data=pvalues,
#  				aes(x=-log10(expected), y=-log10(observed)))
#  p_qq <- p_qq + geom_point(aes(color=causal, shape=trait)) +
#      scale_color_manual(values=c('darkgrey', '#1b9e77'),
#                         name="SNPs",
#  					guide=guide_legend(nrow = 2, title.position="top",
#                                         title.hjust =0.5,
#                                         override.aes = list(shape=21))) +
#      scale_shape_manual(values=c(18, 15:17),
#                         name="Traits",
#  					guide=guide_legend(nrow = 2,
#                                         override.aes = list(colour='darkgrey'),
#                                         title.position="top",
#                                         title.hjust =0.5)) +
#      geom_abline(intercept=0, slope=1, col="black") +
#  	geom_point(data=dplyr::filter(pvalues, causal == "simulated causal"),
#  				color='#1b9e77', aes(shape=trait)) +
#      xlab(expression(Expected~~-log[10](italic(p))) ) +
#      ylab(expression(Observed~~-log[10](italic(p)))) +
#      facet_grid(~ model) +
#      theme_bw() +
#  	theme(strip.background = element_rect(fill="white"),
#  			axis.title = element_text(size=12),
#  			axis.text  =element_text(size=10),
#  			legend.title  =element_text(size=12),
#  			legend.text  =element_text(size=10),
#  			strip.text  =element_text(size=10),
#  			legend.position = "bottom",
#              legend.direction = "horizontal")

## ----echo=FALSE, eval=FALSE, out.width='100%'---------------------------------
#  ggsave(plot=p_qq,
#  		filename=paste(savedir, "/gwas_results_gemma_qqplot.png", sep=""),
#  		units="mm", width=86, height=100)

## ----qqplot, echo=FALSE, fig.cap="\\label{fig:qqplot}Quantile-quantile plots of p-values observed from the multivariate linear mixed model (mvLMM, traits:all) and the univariate linear mixed models (uvLMM, traits: trait1/trait2/trait3) fitted to each of about eight million genome-wide SNPs (grey), including the ten SNPs for which a phenotype effect was modelled (green)."----
#knitr::include_graphics("../docs/articles/Simulation-and-LinearModel_files/figure-html/qqplot-1.png")
knitr::include_graphics("gwas_results_gemma_qqplot.png")

## ---- eval=FALSE--------------------------------------------------------------
#  sigSNPs <- dplyr::filter(pvalues, causal == "simulated causal",
#                           observed < 5*10^-8)
#  sigSNPs$observed_adjusted <- sigSNPs$observed
#  sigSNPs$observed_adjusted[which(sigSNPs$trait != "all")] <-
#      3*sigSNPs$observed[which(sigSNPs$trait != "all")]
#  sigSNPs <- dplyr::filter(sigSNPs, observed_adjusted < 5*10^-8)

## ---- echo=FALSE, eval=FALSE--------------------------------------------------
#  write.table(sigSNPs,
#              paste(savedir, "/sigSNPs.csv", sep=""),
#              col.names=TRUE, row.names=FALSE, quote=FALSE, sep=",")

## ----sigSNPs, echo=FALSE------------------------------------------------------
                        
sigSNPs <- read.table(paste(savedir, "/sigSNPs.csv", sep=""), header=TRUE, 
                      sep=",")
knitr::kable(sigSNPs[, c(1,4,5,7)], 
             caption="\\label{tab:sigSNPs}P-values from multi-trait GWAS (model == mvLMM) and 
                        single-trait GWAS (model == uvLMM)",
             table.attr = "id=\"sigSNPs\"",
             digits=40)

## ---- eval=FALSE--------------------------------------------------------------
#  # Read the SNPs simulated to have an effect on the phenotype and their SNP
#  # effects
#  SNPs <- read.csv(paste(datadir,"/SNP_NrSNP10.csv", sep=""),
#                   row.names=1)
#  SNPeffect <- read.csv(paste(datadir,"/SNP_effects_NrSNP10.csv", sep=""),
#                        row.names=1)
#  
#  
#  # HapGen simulated SNPs might be encoded as 7-25481146 , which R converts to
#  # X7.25481146 via make.names in read.table; substitute X and . to get original
#  # names
#  SNPnames <- gsub("\\.", "-", gsub("X", "", colnames(SNPs)))
#  SNPnames <- SNPnames[order(SNPnames)]
#  simulated_withEffect <- dplyr::filter(pvalues, rsID %in% SNPnames)

## ---- echo=FALSE, eval=FALSE--------------------------------------------------
#  write.table(simulated_withEffect,
#              paste(savedir,"/simulated_withEffect.csv", sep=""), col.names=TRUE,
#              quote=FALSE,row.names=FALSE)

## ---- echo=FALSE, eval=TRUE---------------------------------------------------
SNPeffect <- read.csv(paste(savedir,"/SNP_effects_NrSNP10.csv", sep=""),
                      row.names=1)
SNPs <- read.csv(paste(savedir,"/SNP_NrSNP10.csv", sep=""),
                 row.names=1)
SNPnames <- gsub("\\.", "-", gsub("X", "", colnames(SNPs)))
SNPnames <- SNPnames[order(SNPnames)]
simulated_withEffect <- read.csv(paste(savedir,"/simulated_withEffect.csv", sep=""))

## ----freq-beta, fig.cap="\\label{fig:freq-beta}P-values, allele frequencies and simulated effect sizes of the ten SNPs simulated to have an effect on phenotype."----
# get the allele frequencies of the SNPs simulated to have an effect on the 
# phenotype
freq <- apply(SNPs, 2, getAlleleFrequencies)
minor <- apply(freq, 2, min)
minor <- minor[order(gsub("X", "", names(minor)))]

# format the effect size table and combine with p-value and frequency
# information
effects <- reshape2::melt(SNPeffect, value.name="beta", variable.name="name")
effects$type <- gsub('\\..*', '', effects$name)
effects$rsID <- gsub("\\.", "-", gsub('.*_', '', effects$name))
effects$trait <- rep(paste("trait", 1:3, sep=""), 10)
effects <- effects[order(effects$rsID),]
effects <- dplyr::filter(effects, rsID %in% simulated_withEffect$rsID)

simulated_withEffect <- simulated_withEffect[order(simulated_withEffect[, 1]),]
simulated_withEffect$beta <- NA
simulated_withEffect$beta[simulated_withEffect$trait != "all"] <- effects$beta
simulated_withEffect$freq <- rep(minor[which(SNPnames %in% effects$rsID)], each=4)

# plot association p-values in relation to the absolute value of the simulated 
# effect sizes and the SNPs allele frequencies
p <- ggplot(dplyr::filter(simulated_withEffect, trait != "all"))
p <- p +
    geom_point(aes(x=freq, y=-log10(observed), color=abs(beta), shape=trait)) +
    geom_hline(yintercept=-log10(5*10^(-8)), color="darkgrey") +
    scale_color_gradientn(colours=c('#f0f9e8','#ccebc5','#a8ddb5','#7bccc4',
                                    '#4eb3d3','#2b8cbe','#08589e'),
                          name=expression("| simulated" ~beta~ "|")) +
    scale_shape_discrete(name="Phenotype in GWAS") +
    xlab("Allele frequencies") +
    ylab(expression(-log[10](p-value))) +
    theme_bw()
print(p)

## ----echo=FALSE, eval=FALSE, out.width='100%'---------------------------------
#  ggsave(plot=p,
#  		filename=paste(savedir, "/effectsizes-freq-pvalues.pdf", sep=""),
#  		units="mm", width=150, height=120)

