Merge branch 'master' of ssh://github.com/dmgatti/AttieMetabolomics

README.md

@@ -30,7 +30,22 @@ Plasma lipids: 1767 phenotypes.
    + The output is a \*.rds file containing LOD scores for all analytes (num_markers X num_analytes).
  4. Collect the LOD scores for all chunks into one file.
    + Modify file 'gather_qtl_output.R'.
-   + This file will gather the chunked QTL files an write them out to a single file.
+   + This file will gather the chunked QTL files an write them out to a single \*.rds file.
    + It will also create a PDF containing all QTL plots.
  5. Harvest the QTL peaks above your desired threshold.
- 6. Calculate founder allele effects and association mapping LODs at selected peaks.
+   + Modify 'harvest_thr_qtl.R'.
+   + The input file is the \*.rds file containing all LOD scores created in step 4.
+   + Set a LOD threshold.
+   + The output will be a \*.csv file containing the highest peak on each chromosome above the LOD threshold.
+ 6. Create QTL heatmap and histogram.
+   + Modify 'qtl_heatmap.R' and 'qtl_histogram.R'.
+   + The input for the QTL heatmap is the \*.rds file created in step 4 containing all LOD scores at all markers.
+   + The output is a large PNG.
+   + The input for the QTL histogram is the \*.csv file created in step 5 containing the harvested QTL peaks.
+   + The output is a PNG that plots the chromosomes that have peaks within 2 Mb.
+ 7. Calculate founder allele effects and association mapping LODs at selected peaks.
+  + Modify 'qtl2_coef_assoc_engine.R'.
+  + The input is the \*.rds file created in step 4 and the harvested QTL file created in step 5.
+  + The output will be one founder allele effects plots and one association mapping plot for each peak in the harvested QTL file.
+
+

qtl_heatmap.R

@@ -10,7 +10,7 @@ library(RColorBrewer)
 
 # Arguments:
 # input.file: The *.rds file containing the LOD curves.
-# output.file: full path to the figure file to output.
+# output.file: full path to the figure file to output (PNG).
 # lod.thr: LOD threshold above which LODs will be truncated to prevent large
 #          peaks from dominating the coloring of the heatmap.
 args = commandArgs(trailingOnly = TRUE)

qtl_histogram.R

@@ -1,5 +1,5 @@
 ################################################################################
-# Make a histogram of the top QTL.
+# Given the harvested QTL, make a histogram of the top QTL.
 # Daniel Gatti
 # [email protected]
 # July 20, 2017
@@ -8,52 +8,42 @@ options(stringsAsFactors = F)
 library(tidyverse)
 
 # Arguments:
-# input.file: full path to the *.rds qtl summary file.
+# input.file: full path to the *.csv qtl summary file.
 # output.file: full path to the output figure file as a PNG.
-# thr: LOD threshold to use when selecting QTL peaks.
 args = commandArgs(trailingOnly = TRUE)
 input.file  = args[1]
 output.file = args[2]
-thr = as.numeric(args[3])
 
 print(paste("INPUT.FILE=", input.file))
 print(paste("OUPUT.FILE=", output.file))
-print(paste("THR=", thr))
 
 # Load in the data.
-data = readRDS(input.file)
+data = read.csv(input.file)
+data$chr = factor(data$chr, levels = c(1:19, "X"))
 
+# Load in markers.
 load(url("ftp://ftp.jax.org/MUGA/GM_snps.Rdata"))
-markers = GM_snps[rownames(data), 1:4]
+markers = GM_snps[, 1:4]
 markers[,2] = factor(markers[,2], levels = c(1:19, "X"))
 
 map = split(markers[,3], markers[,2])
 chrlen = sapply(map, length)
-chrlen = chrlen[order(as.numeric(names(chrlen)))]
+chrlen.mb = sapply(map, max)
 chrsum = cumsum(chrlen)
+chrsum.mb = cumsum(c(0, chrlen.mb[-length(chrlen.mb)]))
+names(chrsum.mb) = names(chrlen)
 chrmid = c(1, chrsum[-length(chrsum)]) + chrlen / 2
 names(chrmid) = names(chrlen)
 col = rep(c("grey50", "black"), 10)
 
-png(output.file, width = 1000, height = 800, res = 128)
-
-rs = rowSums(data > thr)
-ylim = range(rs)
-
-print(paste("YLIM = ", ylim))
+# Add genome Mb to data.
+data = data.frame(data, gmb = data$pos + chrsum.mb[data$chr])
 
-rs.df = data.frame(pos = 1:length(rs), rs = rs)
-rs.df = split(rs.df, markers[,2])
+png(output.file, width = 1000, height = 800, res = 128)
 
-par(plt = c(0.08, 0.99, 0.12, 0.95))
-plot(-1, -1, type = "l", las = 1, col = col, xlim = c(0, nrow(markers)), 
-     xaxt = "n", xlab = "", xaxs = "i", ylim = c(0, ylim[2]),
-     ylab = "Number of Analytes")
-for(i in 1:length(rs.df)) {
-  lines(rs.df[[i]][,1], rs.df[[i]][,2], col = col[i])
-}
-mtext(side = 1, line = 1, at = chrmid, text = names(chrmid), cex = 1.5)
+ggplot(data, aes(x = pos)) +
+  geom_histogram(binwidth = 2) +
+  facet_grid(~chr) + 
+  theme(panel.spacing.x = unit(0.1, "lines"))
 
 dev.off()
-
-

qtl_histogram.sh

@@ -9,26 +9,24 @@ module load R/3.3.2
 cd /hpcdata/gac/projects/Attie_DO_Liver_Metabolomics/scripts
 BASEDIR=/hpcdata/gac/projects/Attie_DO_Liver_Metabolomics/
 
-THR=6
-
 ##########
 # Liver lipids: JAX: sex, gen & batch
 INFILE=${BASEDIR}QTL/Liver/lipids_norm_jax/liver_lipids_jax_norm_all_qtl.rds
 OUTFILE=${BASEDIR}figures/QTL/liver_lipids_norm_jax/liver_lipids_jax_norm_qtl_histogram.png
 
-R --no-save --args ${INFILE} ${OUTFILE} ${THR} < qtl_histogram.R > qtl_histogram_liver_lipids.Rout
+R --no-save --args ${INFILE} ${OUTFILE} < qtl_histogram.R > qtl_histogram_liver_lipids.Rout
 
 ##########
 # Liver metabolites: JAX: sex, gen & batch
 INFILE=${BASEDIR}QTL/Liver/metabolites_norm_jax/liver_metabolites_jax_norm_all_qtl.rds
 OUTFILE=${BASEDIR}figures/QTL/liver_metabolites_norm_jax/liver_metabolites_jax_norm_qtl_histogram.png
 
-R --no-save --args ${INFILE} ${OUTFILE} ${THR} < qtl_histogram.R > qtl_histogram_liver_metabolites.Rout
+R --no-save --args ${INFILE} ${OUTFILE} < qtl_histogram.R > qtl_histogram_liver_metabolites.Rout
 
 ##########
 # Plasma metabolites: JAX: sex, gen & batc
 INFILE=${BASEDIR}QTL/Plasma/lipids_norm_jax/plasma_lipids_jax_norm_all_qtl.rds
 OUTFILE=${BASEDIR}figures/QTL/plasma_lipids_norm_jax/plasma_lipids_jax_norm_qtl_histogram.png
 
-R --no-save --args ${INFILE} ${OUTFILE} ${THR} < qtl_histogram.R > qtl_histogram_plasma_lipids.Rout
+R --no-save --args ${INFILE} ${OUTFILE} < qtl_histogram.R > qtl_histogram_plasma_lipids.Rout