**110**wrote:

Dear Community,

I stole the codes from Stephen Turner to analyze RNAseq data using DESeq2 (out put of featureCount).

Everything worked well, except that I kept getting error when plotting heatmap:

Error in plot.new() : figure margins too large Error in par(op) : invalid value specified for graphical parameter "pin"

Below is the heatmap code:

sampleDists <- as.matrix(dist(t(assay(rld)))) library(gplots) png("qc-heatmap-samples.png", w=1500, h=2500, pointsize=1500) heatmap.2(as.matrix(sampleDists), key=F, trace="none", col=colorpanel(100, "black", "white"), ColSideColors=mycols[condition], RowSideColors=mycols[condition], margin=c(10, 10), main="Sample Distance Matrix") dev.off()

I played with the margin but did not make any difference...

Any idea guys? Any input will be appreciated!

Best,

Wenhan

Below is the full version of code for reference:

countdata <- read.table("B6_vs_PPARA_count.txt", header=TRUE, row.names=1) countdata <- countdata[ ,6:ncol(countdata)] colnames(countdata) <- gsub("\.[sb]am$", "", colnames(countdata)) countdata <- as.matrix(countdata) head(countdata) (condition <- factor(c(rep("ctl", 3), rep("exp", 3)))) library(DESeq2) (coldata <- data.frame(row.names=colnames(countdata), condition)) dds <- DESeqDataSetFromMatrix(countData=countdata, colData=coldata, design=~condition) dds dds <- DESeq(dds)

# Plot

dispersions png("qc-dispersions.png", 2000, 2000, pointsize=20) plotDispEsts(dds, main="Dispersion plot") dev.off()

# Regularized log transformation for clustering/heatmaps, etc

rld <- rlogTransformation(dds) head(assay(rld)) hist(assay(rld))

# Colors for plots below

## Ugly:

## (mycols <- 1:length(unique(condition)))

## Use RColorBrewer, better

library(RColorBrewer) (mycols <- brewer.pal(8, "Dark2")[1:length(unique(condition))])

# Sample distance heatmap

sampleDists <- as.matrix(dist(t(assay(rld)))) library(gplots) png("qc-heatmap-samples.png", w=1500, h=2500, pointsize=1500) heatmap.2(as.matrix(sampleDists), key=F, trace="none", col=colorpanel(100, "black", "white"), ColSideColors=mycols[condition], RowSideColors=mycols[condition], margin=c(10, 10), main="Sample Distance Matrix") dev.off()

# Principal components analysis

## Could do with built-in DESeq2 function:

## DESeq2::plotPCA(rld, intgroup="condition")

## I like mine better:

rld_pca <- function (rld, intgroup = "condition", ntop = 500, colors=NULL, legendpos="bottomleft", main="PCA Biplot", textcx=1, ...) {
require(genefilter)
require(calibrate)
require(RColorBrewer)
rv = rowVars(assay(rld))
select = order(rv, decreasing = TRUE)[seq_len(min(ntop, length(rv)))]
pca = prcomp(t(assay(rld)[select, ]))
fac = factor(apply(as.data.frame(colData(rld)[, intgroup, drop = FALSE]), 1, paste, collapse = " : "))
if (is.null(colors)) {
if (nlevels(fac) >= 3) {
colors = brewer.pal(nlevels(fac), "Paired")
} else {
colors = c("black", "red")
}
}
pc1var <- round(summary(pca)$importance[2,1]*100, digits=1)
pc2var <- round(summary(pca)$importance[2,2]*100, digits=1)
pc1lab <- paste0("PC1 (",as.character(pc1var),"%)")
pc2lab <- paste0("PC1 (",as.character(pc2var),"%)")
plot(PC2~PC1, data=as.data.frame(pca$x), bg=colors[fac], pch=21, xlab=pc1lab, ylab=pc2lab, main=main, ...)
with(as.data.frame(pca$x), textxy(PC1, PC2, labs=rownames(as.data.frame(pca$x)), cex=textcx))
legend(legendpos, legend=levels(fac), col=colors, pch=20)
# rldyplot(PC2 ~ PC1, groups = fac, data = as.data.frame(pca$rld),
# pch = 16, cerld = 2, aspect = "iso", col = colours, main = draw.key(key = list(rect = list(col = colours),
# terldt = list(levels(fac)), rep = FALSE)))
}
png("qc-pca.png", 1500, 1500, pointsize=25)
rld_pca(rld, colors=mycols, intgroup="condition", xlim=c(-75, 35))
dev.off()

# Get differential expression results

res <- results(dds) table(res$padj<0.05)

## Order by adjusted p-value

res <- res[order(res$padj), ]

## Merge with normalized count data

resdata <- merge(as.data.frame(res), as.data.frame(counts(dds, normalized=TRUE)), by="row.names", sort=FALSE) names(resdata)[1] <- "Gene" head(resdata)

## Write results

write.csv(resdata, file="diffexpr-results.csv")

## Examine plot of p-values

hist(res$pvalue, breaks=50, col="grey")

## Examine independent filtering

attr(res, "filterThreshold") plot(attr(res,"filterNumRej"), type="b", xlab="quantiles of baseMean", ylab="number of rejections")

## MA plot

## Could do with built-in DESeq2 function:

## DESeq2::plotMA(dds, ylim=c(-1,1), cex=1)

## I like mine better:

maplot <- function (res, thresh=0.05, labelsig=TRUE, textcx=1, ...) { with(res, plot(baseMean, log2FoldChange, pch=20, cex=.5, log="x", ...)) with(subset(res, padj<thresh), points(basemean,="" log2foldchange,="" col="red" ,="" pch="20," cex="1.5))" if="" (labelsig)="" {="" require(calibrate)="" with(subset(res,="" padj<thresh),="" textxy(basemean,="" log2foldchange,="" labs="Gene," cex="textcx," col="2))" }="" }="" png("diffexpr-maplot.png",="" 1500,="" 1500,="" pointsize="25)" maplot(resdata,="" main="MA Plot" )="" dev.off()<="" p="">

## Volcano plot with "significant" genes labeled

volcanoplot <- function (res, lfcthresh=2, sigthresh=0.05, main="Volcano Plot", legendpos="bottomright", labelsig=TRUE, textcx=1, ...) { with(res, plot(log2FoldChange, -log10(pvalue), pch=20, main=main, ...)) with(subset(res, padj<sigthresh ),="" points(log2foldchange,="" -log10(pvalue),="" pch="20," col="red" ,="" ...))="" with(subset(res,="" abs(log2foldchange)>lfcthresh),="" points(log2foldchange,="" -log10(pvalue),="" pch="20," col="orange" ,="" ...))="" with(subset(res,="" padj<sigthresh="" &="" abs(log2foldchange)>lfcthresh),="" points(log2foldchange,="" -log10(pvalue),="" pch="20," col="green" ,="" ...))="" if="" (labelsig)="" {="" require(calibrate)="" with(subset(res,="" padj<sigthresh="" &="" abs(log2foldchange)>lfcthresh),="" textxy(log2foldchange,="" -log10(pvalue),="" labs="Gene," cex="textcx," ...))="" }="" legend(legendpos,="" xjust="1," yjust="1," legend="c(paste("FDR<",sigthresh,sep="")," paste("|logfc|>",lfcthresh,sep="" ),="" "both"),="" pch="20," col="c("red","orange","green"))" }="" png("diffexpr-volcanoplot.png",="" 1500,="" 1500,="" pointsize="25)" volcanoplot(resdata,="" lfcthresh="1," sigthresh="0.05," textcx=".8," xlim="c(-5," 5))="" dev.off()<="" p="">

make the figure smaller and try again

4.3kThanks TriS! I tried and it only works when I made the figure ridiculously large: png("qc-heatmap-samples.png", w=25000, h=27500, pointsize=1500) but it looked super ugly and I can`t see any genes labeled.

Any ideas?

110The error message is indeed related to the dimensions of the image.

Firstly, I would switch to using the ComplexHeatmap package - since I have 'mastered' it, I rarely go back to using heatmap, heatmap.2, or wrappers of these anymore. With ComplexHeatmap, you have greater flexibility in how you arrange pretty much everything on the plot, including labeling.

Generally, I would also switch to using PDF and not png, tiff or jpeg. A PDF is created using vectors and not pixels, meaning better quality and important for when you want to publish. With the pdf() function, you also only need to specify inches of the canvas, i.e., width=5, height=5 is a 5x5" plotting canvas and is sufficient for most graphics.

Kevin

69k