Differential expression with limma on small microarray dataset: design, contrasts, and lack of significant genes

Hi everyone,

I’m here again with some questions regarding differential expression analysis (DEG), contrasts, and limma.

I’m working with the dataset GSE118337, which contains human proximal tubular cells (HK-2 and RPTEC/TERT1) under different conditions: control, TGF-beta, empagliflozin (EMPA), and canagliflozin (CANA), each with ~2 replicates. The main goal of my study is to understand the difference in action between empagliflozin and canagliflozin.

First, when I perform PCA, I observe a clear outlier (HK2_TGFB). Since I am working with a very small number of samples, does it still make sense to remove this outlier? (see PCA here: https://imgur.com/a/P9GK6hY )

Also, from the PCA, I cannot clearly determine whether there is any replicate/batch effect, or if what I am seeing is mainly driven by differences between the two cell types. Is there a recommended way to formally assess this?

For the DEG analysis using limma, I tried two different approaches:

Using a combined group variable (e.g., RPTEC.EMPA, RPTEC.TGFB) and performing contrasts within each cell type (e.g., RPTEC_EMPA - RPTEC_TGFB).

This approach gives me very few or no genes with FDR < 0.05. Using an additive model like ~0 + Condition + Cell (I’m not sure whether I should also include replicate here). With this approach, I obtain many more significant genes.

This makes me unsure about which approach is more appropriate biologically.

Another issue is that for some contrasts, I obtain reasonable p-values, but after multiple testing correction, all adjusted p-values are ~1. I assume this is due to the small sample size. In this scenario, does it still make sense to rely on limma results? Or would it be more appropriate to use other methods?

Overall, I’m struggling to understand what kind of analysis makes the most sense given such a small dataset, and whether limma is still the right tool here. In the end, what I am most interested in are pathway-level changes — are approaches like GSVA reliable with such small sample sizes?

I would really appreciate any guidance. Sorry if some of these questions sound basic — I currently have limited supervision, and this has been quite frustrating as there seem to be many different ways to approach the same problem.

Thanks in advance!

pca differential-expression rna limma

• 809 views • link

updated 2 months ago by Gordon Smyth

• written 3 months ago by fnepo •

The PCA plot shows an absolutely huge outlier that will obviously prevent any meaningful statistical results from the dataset. Have you done anything about the outlier, either removing it, or downweighting using limma's methods for handling outliers? Or have you just continued the analysis without taking any action? Getting non-significant results with such an huge outlier is absolutely predictable, and indeed giving non-significant results for such data is what the statistical methods are supposed to do.

• 0 views • link

written 3 months ago by Gordon Smyth

Hi Gordon, thank you very much for taking the time to reply, I really appreciate it.

Yes, I removed the outlier after confirming it with PCA, heatmap and MA plots.

The main issue is that I only have 2 biological replicates per cell × condition. For HK2_TGFB this means that after removing the outlier I am left with only one sample. In a first approach I removed the entire HK2_TGFB group and analyzed each cell line independently, but then I can only do treatment vs TGFB contrasts in RPTEC.

An alternative I was considering is to use a model like ~ 0 + Condition + Cell, remove only the outlier, and keep the remaining HK2_TGFB sample. This would allow me to still do contrasts vs TGF-beta. Do you think this makes sense, or would it be inappropriate?

Also, would you say there is any replicate effect? Is there a good way to assess this? (edit: then I did individual PCA for each cell and I can clearly see replicate effect) When I include replicate in the model (~ 0 + Group + Replicate) I get around 120–200 DE genes, but without it I get around ~40. I’m concerned that including replicate might be not ok because of the small sample size.

Thank you again for your help!

• 0 views • link

written 3 months ago by fnepo •

I think it is a bit unrealistic to do standard pathway analyses for a small dataset like this with few DE genes. GSVA won't give results, partly for the same reasons as limma, but also because it is designed for simpler two-group comparisons with moderate to large samples sizes.

My feeling it that the best hope for pathway analysis of such a small data would be to use self-contained gene set tests like limma::fry. The original paper for this dataset used DAVID, but the DE results that the DAVID analysis is based on look pretty dodgy, and the test by which they were obtained is not explained in the paper.

• 0 views • link

written 3 months ago by Gordon Smyth

My main interest is comparing canagliflozin vs empagliflozin, so I was thinking that it might make more sense to focus only on the RPTEC cell line and perform those contrasts there.I thought GSVA and GSEA could be more robust here. Even though the statistical power is limited, I would treat the results more as exploratory.

I also tried to look into the limma::fry function you mentioned, but I couldn’t find a specific article describing it in detail. Do you recommend any references or papers that explain its use?

More generally, are there any key or “must-read” papers you would suggest regarding: limma and small sample sizes or pathway analysis approaches in low-powered datasets?

I had a look at some of your work and other limma-related papers, but I’m trying to better understand what is considered good practice in this type of scenario.

Thank you again for your time, I really appreciate it.

• 0 views • link

written 3 months ago by fnepo •

1 answer

The limma documentation and workflows give many example analyses, and they focus on small datasets very much like yours. The limma User's guide gives detailed advice. It includes a number of cases studies, including one with Agilent microarrays, just like your dataset, that you could follow closely. The example case study has rat Agilent arrrays instead of human, but the code is almost the same.

Just out of curiosity, I analysed the GSE118337 dataset myself. You can see my analysis at this link: https://smythlab.github.io/limma/GSE118337.html. I found lots of DE probes (between 400 and 2000) for every treatment vs control comparison.

• 0 views • link

written 3 months ago by Gordon Smyth

Thank you so much, Gordon. I really appreciate the time you took to help me.

I also made a short report to quickly show my analysis. https://gse118337.netlify.app/ When analysing the two cell types separately, I obtained results similar to yours, although my main interest is the contrast against TGFB, as the disease-induced condition. (Edit: After re-reading the article, it seems that they treated control cells separately with TGFB, CANA, and EMPA, rather than applying CANA/EMPA on top of TGFB. If that is the case,I guess that comparing directly against TGFB might not be the most appropriate approach.

I have just a few questions:

Initially, I applied avereps() before the differential expression analysis, because some probes are duplicated. Is it recommended to run the differential expression analysis using all probes, including duplicated ones, and only summarise/filter afterwards? If so, what would be the correct step for doing this?

Related to this, when the same gene has different probes and consequently different logFC and p-values, what is the recommended way to deal with it? For example, should I keep the one with the largest absolute logFC, or is there a better/recommended approach?

Another question is about the design. Why is it preferable here to analyse each cell type separately? What is the main difference between this approach and using a combined Group factor, as I tried in the report, although my results there did not seem reasonable?

Finally, I noticed something strange when checking the differentially expressed results: some GeneName values appear as negative controls. I assume these should be removed when filtering out control probes.

Thank you again for all your help and for your time. Best regards.

• 0 views • link

written 3 months ago by fnepo •

I applied avereps() before the differential expression analysis, because some probes are duplicated. Is it recommended to run the differential expression analysis using all probes, including duplicated ones, and only summarise/filter afterwards? If so, what would be the correct step for doing this?

I don't see any problem with retaining multiple probes for the same gene. To me, it gives more confidence in the results. Reducing to a unique probe for each gene is only needed for specific applications, for example when you need to match results across datasets. If you really do need to uniquify, then I recommend keeping the probe with the highest average expression, see:

I do not recommend averaging probes, because that risks diluting informative probes with uninformative probes. Averaging also introduces complications for the empirical Bayes estimation, because the probe-averages need to be treated differently to probes that are not averaged. I recommend consolidating genewise results at the end, after the linear modeling and eBayes, rather than at the beginning.

Related to this, when the same gene has different probes and consequently different logFC and p-values, what is the recommended way to deal with it? For example, should I keep the one with the largest absolute logFC, or is there a better/recommended approach?

I generally consider a gene to be DE if any of the probes are DE, and DE is judged by p-value and FDR rather than by logFC.

In principle, you could used fry() to consolidate probe results for each gene in a more optimized fashion, but I generally view that as overkill.

Why is it preferable here to analyse each cell type separately? What is the main difference between this approach and using a combined Group factor, as I tried in the report, although my results there did not seem reasonable?

If you are only interested in results are that are consistent between the cell lines, then you could use a combined model with design <- model.matrix(~Treat+CellLine.Replicate). Here, CellLine.Replicate would have four levels, one for each cell line experiment.

We do see that the HK2 and RPTEC results are correlated:

> cor(fit.HK2$coef[,2:4], fit.RPTEC$coef[,2:4])
          TreatTBFB TreatEmpa TreatCana
TreatTBFB   0.34272  -0.03858  -0.06735
TreatEmpa  -0.08143   0.17730   0.16444
TreatCana  -0.06692   0.17518   0.15385

so a combined analysis could work.

The choice depends on your scientific aims and on the underlying biology rather than on any question of statistical "best practice". I tend to view cell lines as independent examples rather than as biological replicates, but it is entirely up to you.

I noticed something strange when checking the differentially expressed results: some GeneName values appear as negative controls. I assume these should be removed when filtering out control probes.

You're right, I should have removed both negative and positive controls by:

NotControl <- y$genes$ControlType==0L
IsExpr <- rowSums(y$other$gIsWellAboveBG > 0) >= 2
yfilt <- y[NotControl & IsExpr, ]

• 0 views • link

written 3 months ago by Gordon Smyth