install.packages('qtl')

install.packages('ggcorrplot')

library(qtl)
library(ggplot2)
library(ggcorrplot)
source("myFunctions.R")

# Read genetic data from the data set
cross1 <- read.cross("csv", file = "Lyons1.csv", 
                     genotypes = c("D", "H", "C"), 
                     na.strings = "-", alleles = c("D", "C"))

# Plot a jittermap of the genetic data - jitters the data so no two markers are on the same spot
jittermap(cross1)

# Summarize the genetic data
summary(cross1)

# Extract the names of phenotypes
names(cross1$pheno)

# Estimate recombination frequencies - Helps us to understand how clean the data is
cross1 <- est.rf(cross1)

plotRF(cross1) #Plots the pairwise recombination data to understand how clean the data is.
#Clean if there are no large red spots

plot.map(cross1) #Plots the positions of the genetic markers on the chromosomes

plotMissing(cross1) #Plots the locations of missing genotypes on the chromosome map, we don't want may missing spots

gallscore <- cross1$pheno$gallstone_score
gallcount <- cross1$pheno$gallstone_count
gallweight <- cross1$pheno$gallstone_weight
chmcagg <- cross1$pheno$ChMC_agg

#Looking to see if the distribution of the data is approximately normal. 
par(mfrow = c(2, 2))
qqnorm(gallscore)
qqline(gallscore)

qqnorm(gallcount)
qqline(gallcount)

qqnorm(gallweight)
qqline(gallweight)

qqnorm(chmcagg)
qqline(chmcagg)
par(mfrow = c(1,1))

#Setting up the basic settings for the main scan for all the variables. 
#We are using a step size of 2 and the Haldane method to map the data. 
cross1 <- calc.genoprob(cross1, step=2.0, off.end = 0.0, error.prob = 1.0e-4, 
                        map.function = "haldane", stepwidth = "fixed")
cross1 <- sim.geno(cross1, step=2.0, off.end = 0.0, error.prob = 1.0e-4, 
                        map.function = "haldane", stepwidth = "fixed")

#Creating a mainscan analysis for the gallscore data using 100 permutations
cross1.scan.gallscore <- scanone(cross1, pheno.col = gallscore, model = "normal", method = "em")
cross1.scan.gallscore.perm <- scanone(cross1, pheno.col = gallscore, model = "normal", method = "em", n.perm = 100) 

#Plotting the gallscore main scan data with threshold lines
plot(cross1.scan.gallscore, main = "Main Scan Plot of Gallscore")
thresh <- summary(cross1.scan.gallscore.perm, alpha = c(0.37, 0.10, 0.05))
abline(h=thresh[1], col = "blue")
abline(h=thresh[2], col = "green")
abline(h=thresh[3], col = "red")

#Creating a summary of the main scan
summary(cross1.scan.gallscore, perm = cross1.scan.gallscore.perm, alpha = 0.05, pvalues = TRUE)
summary(cross1.scan.gallscore, perm = cross1.scan.gallscore.perm, alpha = 0.10, pvalues = TRUE)
summary(cross1.scan.gallscore, perm = cross1.scan.gallscore.perm, alpha = 0.37, pvalues = TRUE)

score <- find.marker(cross1, 2, 53.8)

#Creating an effect plot for the gallscore mainscan run
effectplot(cross1, pheno.col=gallscore, mname1=score)

#Creating a mainscan analysis for the gall count data using 100 permutations
cross1.scan.gallcount <- scanone(cross1, pheno.col = gallcount, model = "normal", method = "hk") #question about pheno.col #

cross1.scan.gallcount.perm <- scanone(cross1, pheno.col = gallcount, model = "normal", method = "hk", n.perm = 100) #question about pheno.col #

#Plotting the gallscore main scan data with threshold lines
plot(cross1.scan.gallcount, main = "Main Scan Plot of Gall Count")
thresh <- summary(cross1.scan.gallcount.perm, alpha = c(0.37, 0.10, 0.05))
abline(h=thresh[1], col = "blue")
abline(h=thresh[2], col = "green")
abline(h=thresh[3], col = "red")

#Creating a summary of the main scan
summary(cross1.scan.gallcount)

#Creating a summary of the main scan with 95% confidence
summary(cross1.scan.gallcount, perm = cross1.scan.gallcount.perm, alpha = 0.05, pvalues = T)

#Creating an effect plot for the gall count mainscan run
par(mfrow = c(1,2))
count1 <- find.marker(cross1, 6, 7.81)
effectplot(cross1, pheno.col=gallcount, mname1=count1)
count2<- find.marker(cross1, 8, 60.14)
effectplot(cross1, pheno.col=gallcount, mname1=count2)
par(mfrow = c(1,1))

cross1.scan.gallweight <- scanone(cross1, pheno.col = gallweight, model = "normal", method = "em") #question about pheno.col #
cross1.scan.gallweight.perm <- scanone(cross1, pheno.col = gallweight, model = "normal", method = "em", n.perm = 100) #question about pheno.col #

plot(cross1.scan.gallweight, main = "Main Scan Plot of Gall Weight")
thresh <- summary(cross1.scan.gallweight.perm, alpha = c(0.37, 0.10, 0.05))
abline(h=thresh[1], col = "blue")
abline(h=thresh[2], col = "green")
abline(h=thresh[3], col = "red")

summary(cross1.scan.gallweight)
summary(cross1.scan.gallweight, perm = cross1.scan.gallweight.perm, alpha = 0.05)
summary(cross1.scan.gallweight, perm = cross1.scan.gallweight.perm, alpha = 0.10)
summary(cross1.scan.gallweight, perm = cross1.scan.gallweight.perm, alpha = 0.37, pvalue = TRUE)

weight <- find.marker(cross1, 8, 61.7)
effectplot(cross1, pheno.col=gallweight, mname1=weight)

cross1.scan.chmcagg <- scanone(cross1, pheno.col = chmcagg, model = "normal", method = "em") #question about pheno.col #
cross1.scan.chmcagg.perm <- scanone(cross1, pheno.col = chmcagg, model = "normal", method = "em", n.perm = 100) #question about pheno.col #

plot(cross1.scan.chmcagg, main = "Main Scan Plot of Severity of Cholesterol Monohydrate Crystals, aggregates")
thresh <- summary(cross1.scan.chmcagg.perm, alpha = c(0.37, 0.10, 0.05))
abline(h=thresh[1], col = "blue")
abline(h=thresh[2], col = "green")
abline(h=thresh[3], col = "red")

summary(cross1.scan.chmcagg)
summary(cross1.scan.chmcagg, perm = cross1.scan.chmcagg.perm, alpha = 0.05, pvalues = T)

chol <- find.marker(cross1, 6, 57.8)
effectplot(cross1, pheno.col=chmcagg, mname1=chol)

augmented <-  mqmaugment(cross1, minprob=0.1, verbose = TRUE)

par(mfrow=c(1,2))
geno.image(cross1)
geno.image(augmented)

summary(cross1)
summary(augmented)

ex <- mqmscan(augmented)
ex1 <- scanone(cross1)

max(ex)
max(ex1)

real_markers <- mqmextractmarkers(ex)
find.marker(augmented, chr=18, pos=5)

markers2 <-mqmextractmarkers(ex1)
find.marker(cross1, chr=18, pos=8.46)


multitoset <- find.markerindex(augmented, "D18Mit64")
setcofactors <- mqmsetcofactors(augmented, cofactors=multitoset)
mqm_co1 <- mqmscan(augmented, setcofactors)

heatmap(cor(cross1$pheno[,4:15], use="complete.obs"))

par(mfrow = c(1,1))
plot(ex, col = "red")
plot(ex1, add=TRUE, col = "green")
plot(mqm_co1, add= TRUE, col = "blue")

data <- read.csv("LyonsNoQTL.csv", header = TRUE)
summary(data)
data[] <- lapply(data, as.numeric)
summary(data)
data <- na.omit(data)

corr <- round(cor(data), 1)
corrchart <- 
  ggcorrplot(corr, hc.order = TRUE,
             type = "lower",
             lab = TRUE,
             lab_size = 3,
             method="circle",
             colors = c("tomato2", "white", "springgreen3"),
             title="Correlogram of Factors for Gallstone Suscebtibility",
             ggtheme=theme_bw)
corrchart