Question: Atypical volcano plot RNA-seq
1
gravatar for Teresa
5 months ago by
Teresa20
Teresa20 wrote:

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

ADD COMMENTlink modified 5 months ago by WouterDeCoster38k • written 5 months ago by Teresa20
1

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 REPLYlink written 5 months ago by andrew.j.skelton735.6k
1

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 REPLYlink modified 5 months ago • written 5 months ago by Teresa20
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: 986 users visited in the last hour