#' Plot founder haplotype dosages/probabilities, ordered by phenotype
#'
#' This function produces a probability heatmap plot, ordered by the phenotype. This gives an idea of what the regression
#' procedure is actually being handed as inputs.
#'
#' @param marker A marker that is contained in the genome cache. In general this should be a marker of interest, such as one
#' beneath a putative QTL peak.
#' @param p.value DEFAULT: NULL. Includes the observed p-value in the plot title.
#' @param genomecache The path to the genome cache that contains founder haplotype information.
#' @param model DEFAULT: "additive". If "additive", dosages are plotted. If "full", probabilities are plotted.
#' @param phenotype The name of the phenotype column in the data set.
#' @param phenotype.data A data.frame object that contains the phenotype information. Should also have a column that matches
#' genomes in the genome cache.
#' @param merge.by DEFAULT: "SUBJECT.NAME". Specifies the columns to merge phenotype and haplotype data.
#' @param founder.labels DEFAULT: "NULL". If NULL, will default to the labels in the genome cache.
#' @param founder.cex DEFAULT: 1. Defines the text size of the founder labels.
#' @param include.ramp DEFAULT: TRUE. If TRUE, spectrum ramp for dosages or probabilities is included.
#' @param ramp.label.cex DEFAULT: 0.7. Specifies the size of the dosage/probability spectrum label.
#' @param ramp.label.line DEFAULT: 0.5. Specifies the number of margin lines that separates the spectrum label from the bottom axis.
#' @param prob.axis.cex DEFAULT: 1. Specifies the size of the dosage/probability tick labels.
#' @param include.marker DEFAULT: TRUE. If TRUE, the marker name is included in the title.
#' @param alternative.phenotype.label DEFAULT: NULL. Allows for an alternative label for the phenotype.
#' @param alternative.marker.label DEFAULT: NULL. Allows for an alternative label for the marker.
#' @export
#' @examples prob.heatmap()
prob.heatmap = function(marker, p.value=NULL, genomecache, model="additive",
phenotype, phenotype.data, merge.by="SUBJECT.NAME",
founder.labels=NULL, founder.cex=1, founder.line=1, founder.col=rep("black", 8),
include.phenotype.axis=TRUE, phenotype.lab.cex=1, phenotype.num.cex=1, phenotype.num.padj=NA,
phenotype.line=NA, phenotype.num.line=NA,
include.ramp=TRUE, ramp.label.cex=0.7, ramp.label.line=0.5, prob.axis.cex=1,
include.marker=TRUE, marker.line=2.5,
alternative.phenotype.label=NULL, alternative.marker.label=NULL){
h <- DiploprobReader$new(genomecache)
X <- h$getLocusMatrix(locus=marker, model=model)
subjects <- h$getSubjects()
rownames(X) <- subjects
if(!include.marker){ marker <- NULL }
if(!is.null(alternative.marker.label)){ marker <- alternative.marker.label }
prob.heatmap.from.matrix(geno.matrix=X, marker=marker, p.value=p.value, model=model, phenotype=phenotype,
phenotype.data, merge.by=merge.by,
founder.labels=founder.labels, founder.cex=founder.cex, founder.line=founder.line, founder.col=founder.col,
include.phenotype.axis=include.phenotype.axis, phenotype.lab.cex=phenotype.lab.cex, phenotype.num.cex=phenotype.num.cex, phenotype.num.padj=phenotype.num.padj,
phenotype.line=phenotype.line, phenotype.num.line=phenotype.num.line,
include.ramp=include.ramp, ramp.label.cex=ramp.label.cex, ramp.label.line=ramp.label.line, prob.axis.cex=prob.axis.cex,
marker.line=marker.line,
alternative.phenotype.label=alternative.phenotype.label)
}
#' @export
prob.heatmap.from.matrix = function(geno.matrix, marker=NULL, marker.line=2.5,
p.value=NULL, model="additive",
phenotype, phenotype.data,
merge.by="SUBJECT.NAME",
founder.labels=NULL, founder.cex=1, founder.line=0.5, founder.col=rep("black", 8),
include.phenotype.axis=TRUE, phenotype.lab.cex=1, phenotype.num.cex=1, phenotype.num.padj=NA,
phenotype.line=NA, phenotype.num.line=NA,
include.ramp=TRUE, ramp.label.cex=0.7, ramp.label.line=0.5, prob.axis.cex=1,
alternative.phenotype.label=NULL){
if(!is.null(p.value)){
p.value <- round(-log10(p.value), 4)
}
X <- geno.matrix
X.data <- data.frame(rownames(X), X)
names(X.data)[1] <- merge.by
if(is.null(founder.labels)){
founder.labels <- colnames(X)
}
num.col <- length(founder.labels)
# Allow function of phenotype
phenotype.data <- model.frame(formula(paste(phenotype, "~ 1 +", merge.by)), data=phenotype.data)
names(phenotype.data)[1] <- "this.y"
final.data <- merge(x=phenotype.data, y=X.data, by=merge.by, all=FALSE)
final.data <- final.data[order(final.data$this.y),] # sort by phenotypic value
probs <- as.matrix(final.data[,-(1:2)]) # Just keep prob of 8 strains for a certain marker
probs <- probs[,rev(1:ncol(probs))]
s <- summary(as.vector(final.data[,"this.y"]))
s1 <- as.character(round(s[1], 2))
s2 <- as.character(round(s[2], 2))
s3 <- as.character(round(s[3], 2))
s5 <- as.character(round(s[5], 2))
s6 <- as.character(round(s[6], 2))
# plot heatmap
## save original par settings
cols <- rev(gray(10000:1/10000))
if (include.ramp) {
op <- par(no.readonly=TRUE)
oplt <- par()$plt
par(fig=c(0.1, 0.85, 0.05, 0.95),
mai=c(0.4, 0.05, 0.7, 0.05)) ##set the margin
}
if (model == "additive") { z.val <- 2 - probs; z.lim <- c(0, 2) }
else { z.val <- 1 - probs; z.lim <- c(0, 1) }
image(z=z.val, axes=FALSE, col=cols, zlim=z.lim)
box()
# axis(2, at=seq(0, num.col, 1+1/num.col)/num.col,
# labels=rev(founder.labels), cex.axis=founder.cex, line=founder.line,
# lty=0, srt=90, las=2, col.axis=founder.col) # add txt on the strain
axis(2, at=seq(0, num.col, 1+1/num.col)/num.col,
labels=FALSE,
lty=0, srt=90)
mtext(text=rev(founder.labels), side=2, at=seq(0, num.col, 1+1/num.col)/num.col, col=founder.col,
las=2, line=founder.line, cex=founder.cex)
phenotype <- ifelse(is.null(alternative.phenotype.label), phenotype, alternative.phenotype.label)
axis(1, at=0.5, labels=parse(text=paste('"-" %<-%', gsub(x=phenotype, pattern=" ", replacement="~"), '%->% "+"')),
tick=FALSE, cex.axis=phenotype.lab.cex, line=phenotype.line)
if (include.phenotype.axis) {
axis(3, at=c(0, 0.25, 0.5, 0.75, 1), labels=c(s1, s2, s3, s5, s6), cex.axis=phenotype.num.cex, line=phenotype.num.line, padj=phenotype.num.padj)
}
if (!is.null(marker)) {
if (is.null(p.value)) {
this.title <- marker
}
else {
this.title <- bquote(.(paste0(marker, ":")) ~ -log[10]*P ~ .(paste("=", p.value)))
}
title(this.title, line=marker.line)
}
if (include.ramp) {
ramp.label <- ifelse(model == "additive", "Haplotype Dose", "Diplotype Prob")
par(fig=c(0.89, 0.92, 0.33, 0.66),
mai=c(0.1, 0.05, 0.5, 0.05),
new=TRUE)
if (model == "additive") { image(y=seq(from=0, to=2, length.out=length(cols)), z=matrix(seq(from=0, to=2, length.out=length(cols)), nrow=1),
zlim=c(0, 2), ylim=c(0, 2), axes=FALSE, col=rev(cols), main="", cex.main=0.77) } #for the legend
if (model == "full") { image(y=seq(from=0, to=1, length.out=length(cols)), z=matrix(seq(from=0, to=1, length.out=length(cols)), nrow=1),
zlim=c(0, 1), ylim=c(0, 1), axes=FALSE, col=rev(cols), main="", cex.main=0.77) }
box()
axis(4, las=1, cex.axis=prob.axis.cex)
mtext(text=ramp.label, side=1, line=ramp.label.line, cex=ramp.label.cex)
par(op)
}
#else{ par(plt <- oplt) }
}
#' @export
prob.image = function(marker.data, marker=NULL, p.value=NULL,
phenotype, phenotype.data, column.labels=NULL){
if(!is.null(p.value)){ p.value <- round(-log10(p.value), 4) }
if(is.null(column.labels)){ column.labels <- colnames(marker.data)[-1] }
if(is.null(marker)){ marker <- "locus" }
num.col <- length(column.labels)
if(sum(marker.data[1,-1], na.rm=TRUE) != 1){
marker.data[,-1] <- marker.data[,-1]/2
}
# allow function of phenotype
phenotype.data <- model.frame(formula(paste0(phenotype, " ~ 1 + SUBJECT.NAME")), data=phenotype.data)
names(phenotype.data)[1] <- "y"
final.data <- merge(x=phenotype.data, y=marker.data, by="SUBJECT.NAME", all=FALSE)
final.data <- final.data[order(final.data$y),] # sort by phenotypic value
probs <- as.matrix(final.data[,-(1:2)]) # Just keep prob of 8 strains for a certain marker
s <- summary(final.data[,"y"])
s1 <- as.character(round(s[1], 2))
s2 <- as.character(round(s[2], 2))
s3 <- as.character(round(s[3], 2))
s5 <- as.character(round(s[5], 2))
s6 <- as.character(round(s[6], 2))
# plot image
## save original par settings
oplt <- par()$plt
cols <- rev(gray(10000:1/10000))
par(plt=c(0.1,.75,0.1,.8)) ##set the margin
image(z=1-probs, axes=F, col=cols, zlim=c(0, 1))
box()
axis(2, at=seq(0, num.col, 1+1/num.col)/num.col,
labels=column.labels,
lty=0, srt=90, las=2) # add txt on the strain
axis(1, at=0.5, labels=phenotype, tick=FALSE)
axis(3, at=c(0,0.25,0.5,0.75,1), labels=c(s1,s2,s3,s5,s6))
this.title <- ifelse(is.null(p.value), marker, paste0(marker, ": -log10P=", p.value))
title(this.title, line=2.5)
par(plt <- oplt)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.