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
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.
Any ideas?
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