search

btmonier / vidger

Make rapid visualizations of RNA-seq data in R
GNU General Public License v3.0
18 stars 6 forks source link

vsVolcano not displaying DESeq2 lfcShrink values? #5

Open diegobah opened 5 years ago

diegobah commented 5 years ago

Can vsVolcano display shrunken lfc values from DESeq2?

Currently supplying such an object (e.g. dds <- DESeq(ddsCounts,betaPrior=TRUE)) to vsVolcano still produces output from lfcMLE of an unshrunk DESeq2 data object.

vsMAPlot displays shrunken LFC values correctly. Does vsVolcano generate a DESEq2 results object as part of its operation and if so how does it generate unshrunken lfc values?

diegobah commented 5 years ago

Looking into the code I noticed the LogFC value is generated rather than obtained from the DESeq object and thus the problem with representing pre-generated shrunken LFC values. The data aquisition for MAplot can be seen to pull the values directly.

getDeseqVolcano <- function(x, y, data, d.factor) { if(missing(d.factor)) { stop( 'This appears to be a DESeq object. Please state d.factor variable.' ) } dat1 <- as.data.frame(colData(data)) dat2 <- fpm(data) dat3 <- results(data, contrast = c(d.factor, y, x)) nam_x <- row.names(dat1[which(dat1[d.factor] == x),]) nam_y <- row.names(dat1[which(dat1[d.factor] == y),]) x <- rowMeans(dat2[, nam_x]) y <- rowMeans(dat2[, nam_y]) id <- row.names(data) dat4 <- data.frame(id, x, y) dat4$logFC <- log2(dat4$y / dat4$x) dat4$pval <- dat3$pvalue dat4$padj <- dat3$padj dat4 <- dat4[complete.cases(dat4),] return(dat4) }

.getDeseqMA <- function(x, y, data, d.factor) { if(missing(d.factor)) { stop( paste( "This appears to be a DESeq object.", "Please state \"d.factor\" variable." ) ) } dat1 <- as.data.frame(colData(data)) dat2 <- fpm(data) dat3 <- results(data, contrast = c(d.factor, y, x)) nam_x <- row.names(dat1[which(dat1[d.factor] == x),]) nam_y <- row.names(dat1[which(dat1[d.factor] == y),]) x <- log10(rowMeans(dat2[, nam_x]) + 1) y <- log10(rowMeans(dat2[, nam_y]) + 1) id <- row.names(data) dat4 <- data.frame(id, x, y) dat4$A <- log10(dat3$baseMean) #dat4$A <- 0.5 log2(dat4$x dat4$y) dat4$M <- dat3$log2FoldChange dat4$pval <- dat3$pvalue dat4$padj <- dat3$padj dat4 <- dat4[complete.cases(dat4), ] return(dat4)
}