DESeq2 Heatmap issue
Entering edit mode
3.6 years ago
dazhudou1122 ▴ 110

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 : 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")

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

Any idea guys? Any input will be appreciated!



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)


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

Regularized log transformation for clustering/heatmaps, etc

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

Colors for plots below


(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")

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([, 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,$x), bg=colors[fac], pch=21, xlab=pc1lab, ylab=pc2lab, main=main, ...) with($x), textxy(PC1, PC2, labs=rownames($x)), cex=textcx)) legend(legendpos, legend=levels(fac), col=colors, pch=20) # rldyplot(PC2 ~ PC1, groups = fac, data =$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))

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(,, 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" )=""<="" 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)&gt;lfcthresh),="" points(log2foldchange,="" -log10(pvalue),="" pch="20," col="orange" ,="" ...))="" with(subset(res,="" padj<sigthresh="" &amp;="" abs(log2foldchange)&gt;lfcthresh),="" points(log2foldchange,="" -log10(pvalue),="" pch="20," col="green" ,="" ...))="" if="" (labelsig)="" {="" require(calibrate)="" with(subset(res,="" padj<sigthresh="" &amp;="" abs(log2foldchange)&gt;lfcthresh),="" textxy(log2foldchange,="" -log10(pvalue),="" labs="Gene," cex="textcx," ...))="" }="" legend(legendpos,="" xjust="1," yjust="1," legend="c(paste("FDR&lt;",sigthresh,sep="")," paste("|logfc|&gt;",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))=""<="" p="">

RNA-Seq heatmap DEseq2 • 3.6k views
Entering edit mode
figure margins too large

make the figure smaller and try again

png("qc-heatmap-samples.png", w=500, h=750)
Entering edit mode

Thanks 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.

looks ugly

Any ideas?

Entering edit mode

The 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.



Login before adding your answer.

Traffic: 1688 users visited in the last hour
Help About
Access RSS

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6