Atypical volcano plot RNA-seq
0
1
Entering edit mode
5.5 years ago
Teresa ▴ 20

Hi all,

I am doing two RNA-seq analysis with unequal sample size both. In the first one, I use 10 healthy controls and 52 patients and I obtain good results when I use logFC=|2| y p-adj 0.01. Nevertheless, when I analyze the population of 7 healthy controls and 58 patients I obtain this volcano plot, which looks really strange. I have never had problems plotting it, so I think that the code I am using is correct. What might be the problem?

alpha <- 0.01
cols <- densCols(res$log2FoldChange, -log10(res$pvalue))
plot(res$log2FoldChange, -log10(res$padj), col=cols, panel.first=grid(),
     main="Volcano plot", xlab="Effect size: log2(fold-change)", ylab="-log10(adjusted p-value)",pch=20, cex=0.6)
abline(v=0)
abline(v=c(-1,1), col="brown")
abline(h=-log10(alpha), col="brown")

gn.selected <- abs(res$log2FoldChange) > 2 & res$padj < alpha 
text(res$log2FoldChange[gn.selected],
     -log10(res$padj)[gn.selected],
     lab=rownames(res)[gn.selected ], cex=0.4)

Thanks so much,
Teresa

enter image description here

volcano-plot RNA-seq deseq2 • 4.4k views
ADD COMMENT
1
Entering edit mode

It would be helpful if you included code for your analysis, did you use lfcShrink? Are your samples all from the same experiment? What did you PCA show?

ADD REPLY
1
Entering edit mode

Thank you for your fast answer. My experiment consists in samples from two centers what I call GOE and BOL. The samples from both centers have been proccess together in the same machine, the only difference between them is that they come from different centers. As the samples from these two centers are of different condition (healthy controls and patients), I want to perform the DE analysis between GOE and UNIBO.

As you can see, the PCA plot shows that the samples are clustering by gender. I have discussed this issue previously with some experts and they have told me that this might be because the condition of interest does not represent a clear difference between these two groups. For that reason, I have used for the design ~ gender + condition.

I have just modified the code as you have mentioned lfcShrink, obtaining the following:

library(DESeq2)
directory = "/Users/Desktop/ESTUDIOS/RNA_SEQ/COUNTS_HTSEQ/"
setwd(directory)

# Load data directly from HTSeq Count
sampleFiles=grep('COUNT',list.files(directory),value=TRUE)
sampleFiles
sampleNames=sub('.txt','',sampleFiles)
sampleNames
sampleCondition=c("BOL","BOL","BOL","BOL","BOL","BOL","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE","GOE")
sampleGender=c("M","M","F","F","M","F","M","F","M","M","M","M","F","M","M","F","F","M","M","F","F","M","F","M","F","F","F","F","F","F","M","F","M","F","M","M","M","M","F","M","M","M","F","F","F","M","M","M","M","M","M","M","F","M","F","F","F","M","M","M","M","M")

sampleTable=data.frame(sampleName=sampleNames, fileName=sampleFiles, condition=sampleCondition, gender=sampleGender)
sampleTable
ddsHTSeq = DESeqDataSetFromHTSeqCount(sampleTable=sampleTable, directory=directory, design= ~ gender + condition)
ddsHTSeq<-estimateSizeFactors(ddsHTSeq)

ddsHTSeq<-estimateSizeFactors(ddsHTSeq)

sizeFactors(ddsHTSeq)

vsd <- varianceStabilizingTransformation(ddsHTSeq, blind=FALSE)

plotPCA(vsd, intgroup=c("gender","condition"))

DESeq.ds <- DESeq(ddsHTSeq, betaPrior = FALSE)

res1 <- results(DESeq.ds, contrast = c("condition", "BOL", "GOE"))
res1 <- lfcShrink(DESeq.ds, contrast = c("condition", "BOL", "GOE"), res = res1)

rownames(res1)

library(EnhancedVolcano)

library(magrittr)

EnhancedVolcano(res1,

                lab = rownames(res1),

                x = "log2FoldChange",

                y = "padj",

                pCutoff = 0.0001,

                transcriptPointSize = 1.5,

                transcriptLabSize = 3.0)

enter image description here enter image description here

Sorry for my english. And thank you so much.

ADD REPLY

Login before adding your answer.

Traffic: 1943 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

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

Powered by the version 2.3.6