Atypical volcano plot RNA-seq

Question

[Teresa]

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

[andrew.j.skelton73]

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?

[Teresa]

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)

Sorry for my english. And thank you so much.