Question: DESeq2 Heatmap issue
0
gravatar for dazhudou1122
3 months ago by
dazhudou112240
dazhudou112240 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)&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))="" dev.off()<="" p="">

heatmap rna-seq deseq2 • 327 views
ADD COMMENTlink modified 3 months ago • written 3 months ago by dazhudou112240
figure margins too large

make the figure smaller and try again

png("qc-heatmap-samples.png", w=500, h=750)
ADD REPLYlink written 3 months ago by TriS3.2k

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?

ADD REPLYlink modified 3 months ago • written 3 months ago by dazhudou112240

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.

Kevin

ADD REPLYlink written 3 months ago by Kevin Blighe9.0k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1336 users visited in the last hour