Hi,

I have run an RNAseq analysis and am stuck on what would be the most appropriate thing to report in results.

I first decided to report the results of IHW tested results from DESeq2, from reading the vignettes this output is suggested to be a better method than the default results output by DESeq2. Is this correct?

Secondly, I applied apeglm shrinkage to my dds object for plotting a MA plot, because from reading the vignettes this output is suggested to be a better shrinkage than the default shrinkage output by DESeq2.

However I am getting confused by which results to finally report.

I am hoping to sort my results by the most robust changes (significnace value) in my dataset as many of my genes have a LFC of >2 in all results.

Should I output apeglm results and report the associated s values? or should I use IHW results and report the associated adj.pvalues?

OR should I output apeglm results and report the associated adj.pvalues?

Which is the 'gold standard' or best course of action?

Fot the latter, svalues = FALSE, doesn't seem to be working for me to produce padj values when I run my dds object through lfcshrink()

> res.sLFC <- lfcShrink(dds.s, coef="comparison_disease_vs_Control", svalue = FALSE, type="apeglm", lfcThreshold=1)
using 'apeglm' for LFC shrinkage. If used in published research, please cite:
    Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
    sequence count data: removing the noise and preserving large differences.
    Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895
computing FSOS 'false sign or small' s-values (T=1)
> res.sLFC
log2 fold change (MAP): comparison PPCD vs Control 

DataFrame with 31172 rows and 4 columns
                   baseMean log2FoldChange     lfcSE      svalue
                  <numeric>      <numeric> <numeric>   <numeric>
ENSG00000000003 30899.27765      -2.980110  0.313910 1.14037e-11
ENSG00000000005     2.92703      -3.052505  1.084273 4.83066e-03
ENSG00000000419  2457.80966       0.288835  0.164714 7.25880e-01
ENSG00000000457  1310.70709      -0.930851  0.182720 2.58388e-01
ENSG00000000460   417.26650       0.919899  0.334406 2.23694e-01
...                     ...            ...       ...         ...
ENSG00000288584    2.418664    -0.21016742  0.713188    0.536877
ENSG00000288585   10.382099     0.00861317  0.381028    0.678976
ENSG00000288586   13.772947    -0.34706693  0.427470    0.613719
ENSG00000288587    8.671689     0.01334065  0.859633    0.553893
ENSG00000288588    0.715984    -0.34188076  0.952566    0.358753

Many thanks for any assistance,

Nathan

"better method than the default results output by DESeq2. Is this correct?"

Yes, IHW is optimal adjustment compared to the procedure built into results.

"I am hoping to sort my results by the most robust changes (significance value)"

If you want to focus on significance, then go with the IHW adjusted p-value.

"Should I output apeglm results and report the associated s values?"

You could report the MLE LFC, the IHW adjusted p-values, and the posterior LFC as an effect size suitable for ranking genes by largest effects (whereas the MLE is not suitable for ranking).

svalue=FALSE has to be swapped because you asked lfcShrink to compute s-value on posterior LFC when setting an lfcThreshold.

If you want p-values applying lfcThreshold use results. Then you can use lfcShrink to produce the posterior LFC and you can add this column to the results object from results.