#other useful links:
#https://github.com/cran/qtlbim.git
#http://www.ssg.uab.edu/qtlbim/
#http://www.rqtl.org/rqtltour2.pdf

#sets working directory to read files
setwd("/Users/guillaume/Desktop/Rqtl_DataSets/rqtl_Carotenoids")
#identify and load packages needed
install.packages("qtl")
install.packages("devtools")
devtools::install_github("cran/qtlbim")
library(qtl)
library(qtlbim)

#######################################################################################
##########################***1-IMPORTING DATA ***######################################
#######################################################################################

#to get help for the module
help("qtl-package")

#Maize data read input file
pop11 = read.cross(format='csv', file='ASCarotFamily091014.csv')
pop11 = read.cross(format='csv', file='WACarotFamily091014.csv')

#Tomato data read input file
tomate <- read.cross("mm", ".", "tomateCL_raw.txt", , "tomateCL_map.txt")

#since some markers are at the same position based on composite maps adjust
pop11 <- jittermap(pop11)

#######################################################################################
##########################***2-DATA CHECKING ***#######################################
#######################################################################################

#By default the "class" of the cross is "f2" or "bc", these are F6 or S5 lines
#so this is not correct for this material
#the support for RILs is strange
#the best method is to load into bc version,
#so AA=A, AB=NA (hets are converted to missing data), and BB=H
#class(pop11)[1]="riself"

#sim.geno takes care of the missing data - hets and real missing
# = missing genotypes are simulated given observed marker data
pop11 <- sim.geno(pop11)

#Data summary 
pop11
summary(pop11) 
nind(pop11) 
nphe(pop11) 
nchr(pop11) 
nmar(pop11) 
totmar(pop11)


######## Plot Linkage map, missing values and phenotypes on the same output ###########

par(mfrow=c(2,2))
#to get pairwise recombination fraction LOD scores for tests of linkage between pairs of markers
plot.rf(pop11)
#to get various plots of the data
plot(pop11)
# plot markers physical pos only
plot.map(pop11)

############################ Check genotypic data #####################################

#Uses the hidden Markov model technology to calculate the probabilities of the true underlying 
#genotypes given the observed multipoint marker data, with possible allowance for genotyping errors. 
pop11 <- calc.genoprob(pop11, step=2.5, error.prob=0.01)

#Identify likely genotyping errors
pop11 <- calc.errorlod(pop11,error.prob=0.01)

## plot observed genotypes, flagging likely errors
plot.geno(pop11,chr=c("2"),ind=c(40:80),cutoff=4,min.sep=2,cex=0.8)
#Likely errors are indicated by red squares. In a backcross, genotypes AA and AB are indicated by white and black circles, respectively.
#In an intercross, genotypes AA, AB and BB are indicated by white, gray, and black circles, respectively, 
#and the partially missing genotypes "not BB" (D in mapmaker) and "not AA" (C in mapmaker) are indicated by green and orange circles, respectively.

############################ plot observed phenotypes #################################

#plot one phenotype
plot.pheno(pop11, pheno.col=1)

#plot multiple phenotypes
par(mfrow=c(2,3))
for (i in 1:6)
{
  plot.pheno(pop11, pheno.col=i)
}

############################ correlation between traits ##############################

phe.pop11<-pull.pheno(pop11,c(1:6))
pairs(phe.pop11)


# Estimate genetic maps
newmap<-est.map(pop11,error.prob=0.01)
# where: 
# steps =Maximum distance (in cM) between positions at which the genotype probabilities are calculated, 
# error prob: Assumed genotyping error rate used in the calculation of the penetrance 

#######################################################################################
##########################***3-SINGLE QTL ANALYSIS ***#################################
#######################################################################################

##arguments
#cross: cross object
#pheno.col: phenotypic trait column to analyse
#model: (normal model, a model for binary traits, a two-part model or non-parametric analysis)
#method: (EM algorithm, imputation, Haley-Knott regression, the extended Haley-Knott method, or marker regression)
#addcovar:	(Additive covariates; allowed only for the normal and binary models.)

#Single marker using EM algorithm (Lander and Botstein 1989) on trait #2
#output format: chr, pos, lod
smout.em2 <- scanone(pop11, pheno.col=2, method="em")
#Single marker using Haley-Knott regression (Haley and Knott 1992) on trait #2
smout.hk2 <- scanone(pop11, pheno.col=2, method="hk")

#Single marker using EM algorithm (Lander and Botstein 1989)  on trait #6
smout.em6 <- scanone(pop11, pheno.col=6, method="em")
#Single marker using Haley-Knott regression (Haley and Knott 1992)  on trait #6
smout.hk6 <- scanone(pop11, pheno.col=6, method="hk")

#plot QTL results according to different models on trait #2
plot(smout.em2, smout.hk2, lty=1, col=c("blue","red"))
#plot differences between the two models
plot(smout.hk2 - smout.em2, ylim=c(-0.3, 0.3), ylab="LOD(HK)-LOD(EM)")

#plot QTL results according to different models on trait #6
plot(smout.em6, smout.hk6, lty=1, col=c("blue","red"))

#plot QTL results according for trait #1, trait #2 and trait #6
plot(smout.hk1, smout.hk2, smout.hk6,lty=1, col=c("blue","red", "green")) 

#Summary provides the top QTL peak on each chromosome for trait #2 
#this can be a real problem where there are multiple QTL per chromosome
summary(smout.em6)
summary(smout.hk6)
#add LOD threshold
summary(smout.hk2,5)

########################### Plot QTL effects across the genome ########################
##infos:
# For a backcross, only one line displayed: additive effect. 
# For an intercross, two lines displayed: the additive and dominance effects.
# additive: contributions to the final phenotype from more than 1 gene, or from alleles of a single gene 
# dominance: effect on phenotype of one allele masks the contribution of a second allele at the same locus
addeff=effectscan(pop11)
addeff

#Do permutations to get a genome-wide significance threshold (=genome-scan-adjusted p-values)
smout.hk.perm <- scanone(pop11, method="hk", n.perm=500)
mean(smout.hk.perm)
hkthreshold=summary(smout.hk.perm, alpha=0.05, pvalues=TRUE)
plot(smout.hk2, lty=1, col=c("red"))
add.threshold(smout.hk2, perms=smout.hk.perm, alpha=0.05) #OR abline(h=hkthreshold)

#To obtain the 1.5-LOD support interval and 95% Bayes interval for the QTL on chromosome 9, type:
lodint(out.hk, chr=9)
bayesint(out.hk, chr=9)
#The first and last rows define the ends of the intervals; the middle row is the estimated QTL location.
#It is sometimes useful to identify the closest flanking markers; use expandtomarkers=TRUE:
lodint(out.hk, chr=9, expandtomarkers=TRUE)
bayesint(out.hk, chr=9, expandtomarkers=TRUE)
#We can calculate the 2-LOD support interval and the 99% Bayes interval as follows.
lodint(out.hk, chr=9, drop=2)
bayesint(out.hk, chr=9, prob=0.99)

#############################################################################################
##########################***4- COMPOSITE INTERVAL MAPPING ***###############################
#############################################################################################

out.cim5 <- cim(pop11, n.marcovar=5)
summary(out.cim5)
#Permutations in CIM (takes some time)
out.cim5.perm <- cim(pop11, n.marcovar=5, n.perm=10)
mean(out.cim5.perm)
#define LOD threshold
cimthreshold=summary(out.cim5.perm, alpha=0.05)
#plot composite (green) and single locus analysis (red) analysis
plot(out.cim5, smout.hk2, lty=1, col=c("darkgreen","red"))
# add threshold from single locus mapping
add.threshold(smout.hk2, perms=smout.hk.perm, alpha=0.05)
# add threshold from CIM
add.threshold(out.cim5, perms=out.cim5.perm, alpha=0.05)

#Extract out the large effect QTL and create a model
bigQTL=summary(smout.hk2, perms=smout.hk.perm, alpha=0.2)
bigQTL
#pull genotype probabilities/imputed genotypes at the nearest pseudomarkers for bigQTL
qtl <- makeqtl(pop11, bigQTL$chr, bigQTL$pos, what="prob")
#Fit a user-specified multiple-QTL model.
#fitqtl(pop11, pheno.col=1, qtl, formula=y~Q1+Q2+Q3+Q4, get.ests=TRUE, method="hk")
fitqtl(pop11, pheno.col=1, qtl, formula=y~Q1+Q2, get.ests=TRUE, method="hk")

#r QTL also has a model building approach that includes marginal 
#and epistatic effects - I am not sure I agree that it is working
#had to change back to BC to make it work
class(pop11)[1]="bc"
swQTLmodel <- stepwiseqtl(pop11, max.qtl=5, method="hk", keeptrace=TRUE)
swQTLmodel 
#fitqtl(pop11, pheno.col=1, swQTLmodel, formula=y~Q1+Q2+Q3+Q4+Q5, get.ests=TRUE, method="hk")
fitqtl(pop11, pheno.col=1, swQTLmodel, formula=y~Q1+Q2+Q3+Q4, get.ests=TRUE, method="hk")
plotModel(swQTLmodel)

#return it back to riself for any other playing
class(pop11)[1]="riself"

# QTL Cartographer format
write.cross(pop11, "qtlcart", "pop11qtlcart")
#try this out at home on QTLCartographer - there is a windows version with a
#graphical interface and also command line versions

#######################################################################################
##########################***5-NON NORMAL PHENOTYPES ***###############################
#######################################################################################

#adapt your model to trait distribution
# binaires : model="binary"
# normally distributed :continuous trait   : model="normal"
# not normally distributed: needs nonparametric model   : model="np"

# Shapiro test: normal distribution if W close to 1, p-value is non significant 
#get pheno 
phen.tomate <- pull.pheno(tomate,c(1:8))
phen.pop11 <- pull.pheno(pop11,c(1:6))

# exemple: 
shapiro.test(phen.pop11[,1])

# Boucle sur chaque trait:

par(mfrow=c(3,2))
for (i in c(1:6))
  
{
  hist(phen.pop11[,i],main=dimnames(phen.pop11)[[2]][i])
  mtext(round(shapiro.test(phen.pop11[,i])$p.value,6))
  print(dimnames(phen.pop11)[[2]][i])
  print(shapiro.test(phen.pop11[,i])$statistic)
  print(shapiro.test(phen.pop11[,i])$p.value)
}

####################################################################################
##########################***6-SINGLE QTL EFFECTS ***###############################
####################################################################################

#We may obtain plots indicating the estimated effects of the QTL via plot PXG , which creates a dot plot, or effectplot,which plots the average phenotype for each genotype group.
#For plotPXG , we must first identify the marker closest to the QTL peak. Use find.marker.
max(out.hk)
mar <- find.marker(pop11, chr=2, pos=68)
plotPXG(pop11, marker=mar)
#Note that red dots correspond to inferred genotypes (based on a single imputation).
#The function effectplot uses the multiple imputation results from sim.geno effectplot(sug, mname1=mar)
#We may use effectplot at a position on the “grid” between markers, using "7@47.7" to indicate the position at 47.7 cM on chr 7.
effectplot(pop11, mname1="2@68")
#Similar plots may be obtained for the locus on chr 9.
max(out.hk, chr=9)
mar2 <- find.marker(pop11, chr=9, pos=58)
plotPXG(pop11, marker=mar2)
effectplot(pop11, mname1="9@58")
#We may plot the joint effects of the two loci via plotPXG as follows:
plotPXG(pop11, marker=c(mar, mar2))
plotPXG(pop11, marker=c(mar2, mar))
#The function effectplot gives more readable figures in this case; it’s often useful to look at it in both ways.
effectplot(pop11, mname1="2@68", mname2="9@58")
effectplot(pop11, mname2="9@58", mname1="2@68")
#The two loci do not appear to interact.