# Map QTLs for shape phenotypes measured in CFW outbred mice using a linear
# mixed model as implemented in GEMMA. It corrects for relatedness between 
# samples using a kinship matrix.
# The method "leave one chromosome out". The kinship matrix doesn't include 
# the SNPs of the chromosome being mapped. 


# load packages and source code for extra functions
library(abind)
library(qtl)
library(data.table)
source("run.gemma.R")
source ("other.code.R") 
source("analyses.PCs.R")

# define the path where 1) the files required to run gemma are gonna be created,
# 2) the results are gonna be stored, 3) gemma.exe is stored.
gemmadir       <- "/home/pallares/Chicago"  
resultdir      <- "/home/pallares/Chicago"   
gemma.exe      <- "/home/pallares/bin/GEMMA/gemma"  

# Load the phenotype data
cat("Loading phenotype data.\n")
pheno <- read.pheno.pl("../data/pheno.mandskull.csv") 
pheno <- prepare.pheno(pheno) 

# Choose the "phenovector" depending on the trait to analyse, mandible or skull.

	# phenovector for mandible traits
phenovector <- c("mCS","mPC1", "mPC2", "mPC3", "mPC4","mPC5", "mPC6","mPC7","mPC8","mPC9","mPC10","mPC11","mPC12","mPC13","mPC14","mPC15","mPC16","mPC17","mPC18","mPC19","mPC20", "mPC21")

	# phenovector for skull traits
#phenovector <- c("sCS","sPC1","sPC2","sPC3","sPC4","sPC5", "sPC6","sPC7","sPC8","sPC9","sPC10","sPC11","sPC12","sPC13","sPC14","sPC15","sPC16","sPC17","sPC18","sPC19","sPC20", "sPC21", "sPC22")

# This loop will run gemma for each PC within the phenovector
# The phenovector for mandible has 22 items (1:22) & for skull 23 items (1:23)
for (i in 1:2){
  trait <- phenovector[i]

# SCRIPT PARAMETERS

which.analysis <- trait 

# Get the phenotype and covariates used in the QTL mapping, and the
# pathname of the file for saving the results.
analysis   <- analyses[[which.analysis]]
phenotype  <- analysis$pheno
covariates <- analysis$cov
outliers   <- analysis$outliers

if (!is.null(outliers))
  pheno <- remove.outliers(pheno,phenotype,covariates,outliers)

# Only analyze samples for which the phenotype and all the covariates are
# observed.
pheno <- pheno[which(none.missing.row(pheno[c(phenotype,covariates)])),]


# LOAD GENOTYPE DATA

# Load the "mean genotypes", or the mean alternative allele counts.
cat("Loading genotype data.\n")
load("../data/geno.dosage.RData")  
# Discard genotype samples from mislabeled flowcell samples.
X <- X[which(discard == "no"),]

# Align the phenotypes and genotypes
ids   <- intersect(pheno$id,rownames(X))
pheno <- pheno[match(ids,pheno$id),]
X     <- X[match(ids,rownames(X)),]

# Discard SNPs with low "imputation quality" assessed by inspecting
# the genotype probabilities. Retain SNPs for which: (1) at least 95%
# of the samples have a maximum probability genotype greater than than
# 0.5; (2) the minor allele frequency is greater than 2%.
f       <- apply(X,2,computemaf)  # computemaf is in misc.r
markers <- which(quality > 0.95 & f > 0.02)
map     <- map[markers,]
X       <- X[,markers]

# MAP QTLs 

# Calculate p-values using GEMMA.
gwscan.gemma <- run.gemma(phenotype,covariates,pheno,X,map,
                          gemmadir,gemma.exe)

#Save results to file.
cat("Saving results to file.\n")
save(list = c("analysis","gwscan.gemma"),
     file = paste0(resultdir,"/",analysis$file))
	 
rm(X)
}