Best statistical test when plotting normalised counts

Question

[Mozart]

Hi there, I know it can be considered naive but I am trying to plot normalised counts on selected genes, outputted from my deseq2, that I have further confirmed with rt-qpcr experiments (the idea is just to confirm that these genes follow the same, expected, trend).

since these genes are supposed to be significantly regulated upon my treatment, I was wondering if a statistical test should be performed on those or this should be considered as 'wrong step' to do.

Answer 1

[jared.andrews07]

How are you plotting the counts? Side by side with treated vs untreated? If so, why not just use the p-value from DESeq2? If not, you need to clarify what/how exactly you're plotting.

[Kevin Blighe]

Yes, the p-values generated by DESeq2 are derived from tests performed on the normalised counts. So, please use those, Mozart. It is a p-value from the Wald test.

[Mozart]

Thanks both for very useful comments, I will be using adjusted p values so in that case I should write in the methods section: p value from Wald test adjusted with benjiamini hochberg correction? (if DESeq2 is used, of course)

is that correct?

[Kevin Blighe]

Sounds good to me / Sembra buono, assuming that you used results() in the standard way (?), and assuming that you are comparing, for example, Condition A vs Condition B and also showing Condition A and Condition B in the plot (?)

[Mozart]

Thanks a lot for this information, Kevin! I am also in doubt with another problem related to this topic. I was wondering if one could plot (e.g. PCA plot etc.) and publish results generated with surrogate variable analysis (SVA package) while just using the batch effect DESeq2 built-in compensation strategy (i.e. ~intercept+batch) for differentially expressed genes. at the end of the day, so far, the only way to visualise batch effect-corrected plots is by borrowing limma package..so that's why I am asking this. If this sounds reasonable to you, guys.

[Kevin Blighe]

Bongiorno, you mean that you used limma::removeBatchEffect()? At which point in the pipeline did you use it? - on the r-log expression levels?

[Mozart]

So, after having created a tximport object - which imports counts from kallisto - I then create the sample table needed for the dds object. A new column is then added that contains the days in which the experiment was performed (i.e. day1, day2 etc). then I create the dds object

dds <- DESeqDataSetFromTximport(txi, coldata, ~batch+condition)
results=results(dds)
rld <- rlog(dds)

and then

assay(rld) <- limma::removeBatchEffect(assay(rld), rld$batch)

[jared.andrews07]

I generally dislike limma's removeBatchEffect function, as it directly affects your base counts. Accounting for this in your design formula is typically enough for DGE, though I suppose for visualization, it may prove useful to use the limma function at times. I don't see anything particularly wrong with your methodology. Kevin may have additional insight here, though I feel using the variance stabilized counts from DESeq2 should typically be good enough for most visualizations.

[Mozart]

Waiting for Kevin's reply, I have another doubt (so hopefully you and Kevin could help me on both topics!). If I wanted to use combat or sva, should I implement this function just for visualisation purposes or should I use dds corrected object for generating the list of differentially expressed genes? sorry for my English, hope this makes some sort of sense to you guys!

[jared.andrews07]

DESeq2 should handle any batch effects in the model itself for differential gene expression, so there is no need to use combat or sva for getting your list of DEGs. Use the dds object. I would try the DESeq2 output for visualization first too - there isn't much point trying to deal with batch effects if they aren't present.