#####################################################################
#
# mqmplots.R
#
# Copyright (c) 2009-2020, Danny Arends
# Copyright polyplot routine (c) 2009, Rutger Brouwer
#
# Modified by Pjotr Prins and Karl Broman
#
#
# first written Februari 2009
# last modified Feb 2020
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License,
# version 3, as published by the Free Software Foundation.
#
# This program is distributed in the hope that it will be useful,
# but without any warranty; without even the implied warranty of
# merchantability or fitness for a particular purpose. See the GNU
# General Public License, version 3, for more details.
#
# A copy of the GNU General Public License, version 3, is available
# at http://www.r-project.org/Licenses/GPL-3
#
# Part of the R/qtl package
# Contains: mqmplot.directedqtl
# mqmplot.cistrans
# addloctocross
# polyplot
# getThird
# getChr
# mqmplot.multitrait
# mqmplot.permutations
# mqmplot.singletrait
#
#
#
#####################################################################
mqmplot.directedqtl <- function(cross, mqmresult, pheno.col=1, draw = TRUE)
{
if(is.null(cross)){
stop("No cross object. Please supply a valid cross object.")
}
if(!is.null(cross$mqm$Nind)){
stop("Augmented crossobject. Please supply the original unaugmented dataset.")
}
if(is.null(mqmresult)){
stop("No mqmresult object. Please supply a valid scanone object.")
}
if(!inherits(mqmresult, "scanone")){
stop("No mqmresult object. Please supply a valid scanone object.")
}
onlymarkers <- mqmextractmarkers(mqmresult)
eff <- effectscan(sim.geno(cross),pheno.col=pheno.col,draw=FALSE)
if(any(eff[,1]=="X")){
eff <- eff[-which(eff[,1]=="X"),]
}
onlymarkers[,3] <- onlymarkers[,3]*(eff[,3]/abs(eff[,3]))
if(draw) plot(ylim=c((min(onlymarkers[,3])*1.1),(max(onlymarkers[,3])*1.1)),onlymarkers)
class(onlymarkers) <- c("scanone",class(onlymarkers))
if(!is.null(attr(mqmresult,"mqmmodel"))) attr(onlymarkers,"mqmmodel") <- attr(mqmresult,"mqmmodel")
invisible(onlymarkers)
}
mqmplot.heatmap <- function(cross, result, directed=TRUE, legend=FALSE, breaks = c(-100,-10,-3,0,3,10,100), col = c("darkblue","blue","lightblue","yellow","orange","red"), ...)
{
if(is.null(cross)){
stop("No cross object. Please supply a valid cross object.")
}
if(directed && !is.null(cross$mqm$Nind)){
stop("Augmented crossobject. Please supply the original unaugmented dataset.")
}
if(is.null(result)){
stop("No result object. Please supply a valid scanone object.")
}
if(!inherits(result,"mqmmulti")){
stop("Not a mqmmulti object. Please supply a valid scanone object.")
}
cross <- sim.geno(cross)
names <- NULL
for(x in 1:nphe(cross)){
result[[x]] <- mqmextractpseudomarkers(result[[x]])
if(directed){
effect <- effectscan(sim.geno(cross,step=stepsize(result[[x]])), pheno.col=x, draw=FALSE)
for(y in 1:nrow(result[[x]])){
effectid <- which(rownames(effect)==rownames(result[[x]])[y])
if(!is.na(effectid&&1)){
result[[x]][y,3] <- result[[x]][y,3] *(effect[effectid,3]/abs(effect[effectid,3]))
}
}
}
names <- c(names,substring(colnames(result[[x]])[3],5))
}
chrs <- unique(lapply(result,getChr))
data <- NULL
for(x in 1:length(result)){
data <- rbind(data,result[[x]][,3])
}
rownames(data) <- names
if(nphe(cross) < 100){
image(seq(0,nrow(result[[1]])),seq(0,nphe(cross)),t(data),xlab="Markers",ylab="Traits",breaks=breaks,col=col,yaxt="n",...)
axis(2,at=seq(1,nphe(cross)),labels=colnames(pull.pheno(cross)),las=1)
}else{
image(seq(0,nrow(result[[1]])),seq(0,nphe(cross)),t(data),xlab="Markers",ylab="Traits",breaks=breaks,col=col,...)
}
abline(v=0)
for(x in unique(chrs[[1]])){
abline(v=sum(as.numeric(chrs[[1]])<=x))
}
for(x in 1:nphe(cross)){
abline(h=x)
}
if(legend){
leg <- NULL
for(x in 2:length(breaks)){
leg <- c(leg,paste("LOD",breaks[x-1],"to",breaks[x]))
}
legend("bottom",leg,col=col,lty=1,bg="white")
}
invisible(data)
}
mqmplot.clusteredheatmap <- function(cross, mqmresult, directed=TRUE, legend=FALSE, Colv=NA, scale="none", verbose=FALSE, breaks = c(-100,-10,-3,0,3,10,100), col = c("darkblue","blue","lightblue","yellow","orange","red"), ...)
{
if(is.null(cross)){
stop("No cross object. Please supply a valid cross object.")
}
if(directed && !is.null(cross$mqm$Nind)){
stop("Augmented crossobject. Please supply the original unaugmented dataset.")
}
if(is.null(mqmresult)){
stop("No mqmresult object. Please supply a valid mqmmulti object.")
}
if(!inherits(mqmresult,"mqmmulti")){
stop("Not a mqmmulti object. Please supply a valid mqmmulti object.")
}
cross <- sim.geno(cross)
names <- NULL
for(x in 1:nphe(cross)){
mqmresult[[x]] <- mqmextractpseudomarkers(mqmresult[[x]])
if(directed){
effect <- effectscan(sim.geno(cross,step=stepsize(mqmresult[[x]])), pheno.col=x, draw=FALSE)
if(verbose) cat(".")
for(y in 1:nrow(mqmresult[[x]])){
effectid <- which(rownames(effect)==rownames(mqmresult[[x]])[y])
if(!is.na(effectid&&1)){
mqmresult[[x]][y,3] <- mqmresult[[x]][y,3] *(effect[effectid,3]/abs(effect[effectid,3]))
}
}
}
names <- c(names,substring(colnames(mqmresult[[x]])[3],5))
}
if(verbose && directed) cat("\n")
chrs <- unique(lapply(mqmresult,getChr))
data <- NULL
for(x in 1:length(mqmresult)){
data <- rbind(data,mqmresult[[x]][,3])
}
colnames(data) <- rownames(mqmresult[[1]])
rownames(data) <- names
if(length(names) < 100){
retresults <- heatmap(data,Colv=Colv,scale=scale, xlab="Markers",main="Clustered heatmap",breaks=breaks,col=col,keep.dendro = TRUE, ...)
}else{
retresults <- heatmap(data,Colv=Colv,scale=scale, xlab="Markers",main="Clustered heatmap",breaks=breaks,col=col,keep.dendro = TRUE,labRow=1:length(names), ...)
}
if(legend){
leg <- NULL
for(x in 2:length(breaks)){
leg <- c(leg,paste("LOD",breaks[x-1],"to",breaks[x]))
}
legend("bottom",leg,col=col,lty=1,bg="white")
}
invisible(retresults)
}
mqmplot.cistrans <- function(result,cross,threshold=5,onlyPEAK=TRUE,highPEAK=FALSE,cisarea=10,pch=22,cex=0.5, verbose=FALSE, ...)
{
if(is.null(cross)){
stop("No cross object. Please supply a valid cross object.")
}
if(is.null(cross$locations)){
stop("Please add trait locations to the cross file\n")
}
if(inherits(result, "mqmmulti")){
sum.map <- 0
chr.breaks <- NULL
for(j in 1:nchr(cross)){
l.chr <- max(result[[1]][result[[1]][,1]==j,2])
chr.breaks <- c(chr.breaks,sum.map)
sum.map <- sum.map+l.chr
}
sum.map <- ceiling(sum.map)
if(verbose) cat("Total maplength:",sum.map," cM in ",nchr(cross),"Chromosomes\nThe lengths are:",chr.breaks,"\n")
locations <- do.call(rbind,cross$locations)
QTLs <- do.call(rbind,lapply(FUN=getThird,result))
colnames(QTLs) <- rownames(result[[1]])
bmatrix <- QTLs>threshold
pmatrix <- NULL
for(j in 1:nrow(QTLs)){
if(verbose && (j%%1000 == 0)) cat("QTL row:",j,"\n")
temp <- as.vector(bmatrix[j,])
tempv <- QTLs[j,]
value = 0
index = -1
for(l in 1:length(temp)){
if(temp[l]){
if( tempv[l] > value){
#New highest marker set the old index to false
if(index != -1){
temp[index] <- FALSE
}
#and store the new one
value = tempv[l]
index <- l
}else{
#Lower marker
temp[l] <- FALSE
}
}else{
index = -1
value = 0
}
}
if(onlyPEAK){
bmatrix[j,] <- temp
}
if(highPEAK){
pmatrix <- rbind(pmatrix,temp)
}
}
locz <- NULL
for(marker in 1:ncol(QTLs)){
if(verbose && (j%%1000 == 0)) cat("QTL col:",marker,"\n")
pos <- find.markerpos(cross, colnames(QTLs)[marker])
if(!is.na(pos[1,1])){
locz <- c(locz,round(chr.breaks[as.numeric(pos[[1]])] + as.numeric(pos[[2]])))
}else{
mark <- colnames(QTLs)[marker]
mchr <- substr(mark,sum(regexpr("c",mark)+attr(regexpr("c",mark),"match.length")),regexpr(".loc",mark)-1)
mpos <- as.numeric(substr(mark,sum(regexpr("loc",mark)+attr(regexpr("loc",mark),"match.length")),nchar(mark)))
locz <- c(locz,round(chr.breaks[as.numeric(mchr)] + as.numeric(mpos)))
}
}
axi <- 1:sum.map
plot(x=axi,y=axi,type="n",main=paste("Cis/Trans QTL plot at LOD",threshold),xlab="Markers (in cM)",ylab="Location of traits (in cM)",xaxt="n",yaxt="n")
trait.locz <- NULL
for(j in 1:nrow(QTLs)){
if(verbose && (j%%10 == 0)) cat("QTL row:",j,"\n")
values <- rep(NA,sum.map)
aa <- locz[bmatrix[j,]]
trait.locz <- c(trait.locz,chr.breaks[locations[j,1]] + locations[j,2])
values[aa] = chr.breaks[locations[j,1]] + locations[j,2]
if(!highPEAK){
points(values,pch=pch,cex=cex)
}else{
points(values,pch=24,cex=1.25*cex,col="black",bg="red")
}
}
points(axi,type="l")
points(axi+(cisarea/2),type="l",col="green")
points(axi-(cisarea/2),type="l",col="green")
chr.breaks <- c(chr.breaks,sum.map)
axis(1,at=chr.breaks,labels=FALSE)
axis(2,at=chr.breaks,labels=FALSE)
axis(1,at=locz,line=1,pch=24)
axis(2,at=seq(0,sum.map,25),line=1)
}else{
stop("invalid object supplied\n")
}
}
addloctocross <- function(cross,locations=NULL,locfile="locations.txt", verbose=FALSE)
{
if(is.null(cross)){
stop("No cross object. Please supply a valid cross object.")
}
if(is.null(locations)){
locations <- read.table(locfile,row.names=1,header=TRUE, stringsAsFactors=TRUE)
}
if(verbose) {
cat("Phenotypes in cross:",nphe(cross),"\n")
cat("Phenotypes in file:",nrow(locations),"\n")
}
if(max(as.numeric(rownames(locations))) != nphe(cross)){
stop("ID's of traits in file are larger than # of traits in crossfile.")
}
if(nphe(cross)==nrow(locations)){
locs <- vector(mode = "list", length = nphe(cross))
for(x in as.numeric(rownames(locations))){
if(names(cross$pheno)[x] == locations[x,1]){
locs[[x]] <- locations[x,2:3]
rownames(locs[[x]]) <- locations[x,1]
}else{
warning("Mismatch between name of trait in cross & file.\n")
}
}
}else{
stop("Number of traits in cross & file don't match.")
}
cross$locations <- locs
cross
}
polyplot <- function( x, type='b', legend=TRUE,legendloc=0, labels=NULL, cex = par("cex"), pch = 19, gpch = 21, bg = par("bg"), color = par("fg"), col=NULL, ylim=range(x[is.finite(x)]), xlim = NULL,
main = NULL, xlab = NULL, ylab = NULL, add=FALSE, ... )
{
#Addition by Danny Arends
if(legend){
if(legendloc){
op <- par(mfrow = c(1,2))
}else{
op <- par(mfrow = c(1,1))
}
}else{
op <- par(mfrow = c(1,1))
}
#End of addition
if ( is.vector(x) ) {
x <- t( as.matrix(x) )
if (is.null(labels) ) {
rownames(x) <- c("unspecified gene")
} else {
rownames(x) <- labels
}
} else {
x <- as.matrix(x)
}
if (is.null(labels) ) labels = rownames( x )
if (is.null(col) ) col = rainbow( nrow(x),alpha=0.35 )
if (is.null(xlab) ) xlab="Markers"
if (is.null(ylab) ) ylab="QTL (LOD)"
timepoints <- as.numeric( colnames(x) )
tps <- sort( unique( timepoints ) )
if(is.null(xlim)) xlim = c(min(tps),max(tps))
plot.new() # make a new plot
plot.window(xlim=xlim, ylim=ylim, log="") # add the plot windows size
#grid()
for( k in 1:nrow( x ) ) {
max.p <- NULL # the expression of the maximum
min.p <- NULL #
med.p <- NULL #
for( i in 1:length(tps)) { #
idx <- ( 1:length( timepoints ) )[timepoints==tps[i] ] # get the indeces of the
work <- x[k, idx]
pmax <- max(work[is.finite(work)], na.rm=TRUE)
pmin <- min(work[is.finite(work)], na.rm=TRUE)
pmed <- median(work[is.finite(work)], na.rm=TRUE)
max.p <- append(max.p, pmax) #
min.p <- append(min.p, pmin) #
med.p <- append(med.p, pmed) #
}
lines( x=tps, y=max.p, type='l', col=col[k] ) # add the lines if requested
lines( x=tps, y=min.p, type='l', col=col[k] ) # add the lines if requested
xp <- append(tps, rev(tps))
yp <- append(max.p, rev(min.p) )
polygon(xp, y=yp, col=col[k], border=FALSE)
lines( x=tps, y=med.p, type='l', col=col[k] ) # add the lines if requested
}
axis(1) # add the x axis
axis(2) # add the y axis
title(main=main, xlab=xlab, ylab=ylab, ...) # add the title and axis labels
#Addition by Danny Arends
if (legend){
if(legendloc){
plot.new()
}
legend("topright", labels, col=col, pch=pch) # add a legend if requested
}
#End of addition
op <- par(mfrow = c(1,1))
invisible() # return the plot
}
getThird <- function(x){
x[,3]
}
getChr <- function(x){
x[,1]
}
mqmplot.multitrait <- function(result, type=c("lines","image","contour","3Dplot"), group=NULL, meanprofile=c("none","mean","median"), theta=30, phi=15, ...)
{
#Helperfunction to plot mqmmulti objects made by doing multiple mqmscan runs (in a LIST)
type <- match.arg(type)
meanprofile <- match.arg(meanprofile)
if(!inherits(result, "mqmmulti")){
stop("Wrong type of result file, please supply a valid mqmmulti object.")
}
n.pheno <- length(result)
temp <- lapply(result,getThird)
chrs <- unique(lapply(result,getChr))
qtldata <- do.call("rbind",temp)
if(!is.null(group)){
qtldata <- qtldata[group,]
colors <- rep("blue",n.pheno)
}else{
group <- 1:n.pheno
colors <- rainbow(n.pheno)
}
qtldata <- t(qtldata)
if(type=="contour"){
#Countour plot
contour(x=seq(1,dim(qtldata)[1]),
y=seq(1,dim(qtldata)[2]),
qtldata,
xlab="Markers",ylab="Trait", ...)
for(x in unique(chrs[[1]])){
abline(v=sum(as.numeric(chrs[[1]])<=x))
}
}
if(type=="image"){
#Image plot
image(x=1:dim(qtldata)[1],y=1:dim(qtldata)[2],qtldata,xlab="Markers",ylab="Trait",...)
for(x in unique(chrs[[1]])){
abline(v=sum(as.numeric(chrs[[1]])<=x))
}
}
if(type=="3Dplot"){
#3D perspective plot
persp(x=1:dim(qtldata)[1],y=1:dim(qtldata)[2],qtldata,
theta = theta, phi = phi, expand = 1,
col="gray", xlab = "Markers", ylab = "Traits", zlab = "LOD score")
}
if(type=="lines"){
#Standard plotting option, Lineplot
first <- TRUE
for(i in group) {
if(first){
plot(result[[i]],ylim=c(0,max(qtldata)),col=colors[i],lwd=1,ylab="LOD score",xlab="Markers",main="Multiple profiles", ...)
first <- FALSE
}else{
plot(result[[i]],add=TRUE,col=colors[i],lwd=1,...)
}
}
if(meanprofile != "none"){
temp <- result[[1]]
if(meanprofile=="median"){
temp[,3] <- apply(qtldata,1,median)
legend("topright",c("QTL profiles","Median profile"),col=c("blue","black"),lwd=c(1,3))
}
if(meanprofile=="mean"){
temp[,3] <- rowMeans(qtldata)
legend("topright",c("QTL profiles","Mean profile"),col=c("blue","black"),lwd=c(1,3))
}
plot(temp,add=TRUE,col="black",lwd=3,...)
}
}
}
mqmplot.permutations <- function(permutationresult, ...)
{
#Helperfunction to show mqmmulti objects made by doing multiple mqmscan runs (in a LIST)
#This function should only be used for bootstrapped data
matrix <- NULL
row1 <- NULL
row2 <- NULL
i <- 1
if(!inherits(permutationresult, "mqmmulti"))
ourstop("Wrong type of result file, please supply a valid mqmmulti object.")
for( j in 1:length( permutationresult[[i]][,3] ) ) {
row1 <- NULL
row2 <- NULL
for(i in 1:length( permutationresult ) ) {
if(i==1){
row1 <- c(row1,rep(permutationresult[[i]][,3][j],(length( permutationresult )-1)))
names(row1) <- rep(j,(length( permutationresult )-1))
}else{
row2 <- c(row2,permutationresult[[i]][,3][j])
}
}
names(row2) <- rep(j,(length( permutationresult )-1))
matrix <- cbind(matrix,rbind(row1,row2))
}
rownames(matrix) <- c("QTL trait",paste("# of bootstraps:",length(permutationresult)-1))
#Because bootstrap only has 2 rows of data we can use black n blue
polyplot(matrix,col=c(rgb(0,0,0,1),rgb(0,0,1,0.35)),...)
#PLot some lines so we know what is significant
perm.temp <- mqmprocesspermutation(permutationresult) #Create a permutation object
numresults <- dim(permutationresult[[1]])[1]
lines(x=1:numresults,y=rep(summary(perm.temp)[1,1],numresults),col="green",lwd=2,lty=2)
lines(x=1:numresults,y=rep(summary(perm.temp)[2,1],numresults),col="blue",lwd=2,lty=2)
chrs <- unique(lapply(permutationresult,getChr))
for(x in unique(chrs[[1]])){
abline(v=sum(as.numeric(chrs[[1]])<=x),lty="dashed",col="gray",lwd=1)
}
}
mqmplot.singletrait <- function(result, extended=0,...)
{
#Helperfunction to show scanone objects made by doing mqmscan runs
if(!inherits(result, "scanone")){
stop("Wrong type of result file, please supply a valid scanone (from MQM) object.")
}
if(is.null(result$"info")){
stop("Wrong type of result file, please supply a valid scanone (from MQM) object.")
}
if(is.null(attr(result,"mqmmodel"))){
op <- par(mfrow = c(1,1))
}else{
op <- par(mfrow = c(2,1))
}
info.l <- result
info.l[,3] <- result[,4]
if(extended){
if(!is.null(attr(result,"mqmmodel"))){
plot(attr(result,"mqmmodel"))
}
plot(result,lwd=1,col=c("black"),ylab="QTL (LOD)",...)
par(new=TRUE)
plot(info.l,lwd=1,col=c("red"),ylab="QTL (LOD)",yaxt="n",lty=1,...)
grid(length(result$chr),5)
labels <- c(colnames(result)[3],"Information Content")
mtext("Information Content",side=4,col="red",line=4)
axis(4, ylim=c(0,1), col="red",col.axis="red",las=1)
legend("right", labels,col=c("black","red"),lty=c(1,1,1))
}else{
if(!is.null(attr(result,"mqmmodel"))){
plot(attr(result,"mqmmodel"))
}
plot(result,lwd=1,ylab="QTL (LOD)",...)
grid(length(result$chr),5)
labels <- c(colnames(result)[3])
legend("right", labels,col=c("black","blue"),lty=c(1,1))
}
op <- par(mfrow = c(1,1))
}
# end of mqmplots.R
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.