From f1c36fc29b88105969a5ff9209993acc7d28c309 Mon Sep 17 00:00:00 2001 From: Daniel Gatti Date: Sat, 22 Jul 2017 09:31:38 -0400 Subject: [PATCH 1/5] Update qtl_heatmap.R --- qtl_heatmap.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/qtl_heatmap.R b/qtl_heatmap.R index 71abaa7..c1e7bae 100644 --- a/qtl_heatmap.R +++ b/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) From 07470a794a993f0e96d7e7d0f562022eeec28609 Mon Sep 17 00:00:00 2001 From: Daniel Gatti Date: Sat, 22 Jul 2017 09:38:55 -0400 Subject: [PATCH 2/5] Changed to ggplot --- qtl_histogram.R | 39 ++++++++++++++++----------------------- 1 file changed, 16 insertions(+), 23 deletions(-) diff --git a/qtl_histogram.R b/qtl_histogram.R index 26cf1af..96c3316 100644 --- a/qtl_histogram.R +++ b/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 # dan.gatti@jax.org # July 20, 2017 @@ -8,7 +8,7 @@ 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) @@ -21,39 +21,32 @@ 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() - - From 4342aa9fcafbe59f6ca6d163776823dfcd3e7c04 Mon Sep 17 00:00:00 2001 From: Daniel Gatti Date: Sat, 22 Jul 2017 09:41:50 -0400 Subject: [PATCH 3/5] Removed LOD threshold. --- qtl_histogram.R | 3 --- 1 file changed, 3 deletions(-) diff --git a/qtl_histogram.R b/qtl_histogram.R index 96c3316..7d39f36 100644 --- a/qtl_histogram.R +++ b/qtl_histogram.R @@ -10,15 +10,12 @@ library(tidyverse) # Arguments: # 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 = read.csv(input.file) From 90205376d7484a2cbe082cbcffc93bef40f7f343 Mon Sep 17 00:00:00 2001 From: Daniel Gatti Date: Sat, 22 Jul 2017 09:42:29 -0400 Subject: [PATCH 4/5] Removed threshold. --- qtl_histogram.sh | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/qtl_histogram.sh b/qtl_histogram.sh index 03ae0de..7f97c93 100644 --- a/qtl_histogram.sh +++ b/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 From fa984474d31e8e20acc6c4637b2c7c1b007d71b7 Mon Sep 17 00:00:00 2001 From: Daniel Gatti Date: Sat, 22 Jul 2017 09:53:09 -0400 Subject: [PATCH 5/5] Update README.md --- README.md | 19 +++++++++++++++++-- 1 file changed, 17 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index efa72f6..663f198 100644 --- a/README.md +++ b/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. + +