# What is the best way to compare differences in DEG expression between two populations?

I am interested in differential response to a drug in the "RED" group versus the "GRN" group. I mean, I need to compare differences in DEG expression between two populations (RED vs GRN)

I have 3 pairs of control vs treated samples. 3 pairs for each population.

Here is an example of results for a gene named "gen1".

Gname   RED_ctl1  RED_drg1  RED_ctl2  RED_drg2  RED_ctl3  RED_drg3
gen1           4        10         5        10         6         9
Gname   GRN_ctl1  GRN_drg1  GRN_ctl2  GRN_drg2  GRN_ctl3  GRN_drg3
gen1           4         9         5        10         4         9


In the end, I expect to plot a boxplot that would show me (or not) the difference between RED and GRN groups in terms of DEG expression. I expect to find a group (either RED or GRN) with DEGs that respond to drugs much more (or less) in terms of expression.

Since this is my first time doing this analysis, I would like to know what is the classical way of analyzing that kind of data? What are the steps?

A related question was presented here on Bioinformatics SE

A t-test or eyeballing a box plot might be okay for something very informal, but that's not the statistically right way to do it.

If this is based on RNASeq data, the answer is clear. You need DESeq2 or EdgeR. Read their tutorials, and pick one to use.

If this is microarray, limma is probably best.

If it's something else, you'll need to read some papers to see what people use.

It doesn't matter what data you are dealing with classical statistics can handle this data and non-parametric are very versatile, but there is one notable limitation.

If this is strictly paired data (which what I suspect):

GRN 1 | GRN 2 | GRN 3
RED 1 | RED 2 | RED 3


This would be a two-way ANOVA and is a very powerful test it is fulfils the criteria (below). Interpreting the output needs a bit of background reading.

However, if it is simply GRN vs RED then it is one-way ANOVA (or T-test) and Kruskal-Wallis test (non-parametric). In addition, if you are comparing

RED1
RED2
RED3
GRN1
GRN2
GRN3


Thus 6 data sets a T-test doesn't work and its one-way ANOVA or KW test.

ANOVA The Fmax test is needed for ANOVA analysis to assess the homogeneity of variance.

• The test works best when the number of samples drawn from each population is the same (check in this case). This is not an essential requirement but if the sampling between groups is heavily that's not good.
• The underlying populations are normally distributed (same is true of the T-test). If Fmax says 'yes' skip this bit.

I honestly don't know if Fmax is available in Python, it is definitely not numpy.fmax. It is available in R:

hartleyTest(x, ...)


where x is the vector. You can do it manually.

For 1-way ANOVA

from scipy.stats import f_oneway

f_oneway(GRN123, RED123)
f_oneway(GRN1,GRN2,GRN3,RED1,RED2,RED3) # This makes it like the second example


The data is within a list e.g. GRN1 = [1,2,3,1,1,14,1]

Easy.

Two-way ANOVA is complicated

import statsmodels.api as sm
from statsmodels.formula.api import old

model = ols('treatment ~ C(RED1) + C(GRN1) ... : C(GRN1)', data=df).fit() # see note
sm.stats.anova_lm(model, typ=2)


You should be able to circumvent the ols - because its complicated. It is however useful, particularly is Fmax fails.

Here for details.

You might substitute your one-way ANOVA for a T-test in this instance, but it depends how you structure you comparison.

If Fmax fails and you were wanting to perform a one-way ANOVA you shift to the non-parametric Kruskal-Wallis test. In python,

from scipy import stats

stats.kruskal(GRN123, RED123)


Work around if its paired data ...

stats.kruskal (GRN1, RED1)
stats.kruskal (GRN2, RED2)
....


Its a bit cheeky but you can cover yourself with Bonferoni correction and strictly Wilcoxan is better here.

So its easy.

HOWEVER if your data is paired and fails an Fmax test, strictly you go to the Friedman test, this is the non-parametric version of the 2-way ANOVA (by Milton Friedman). Now for some reason no-one ever does this, but I don't remember what the reason why. Thus your options are:

• in that situation you could just use the ols regression model fitted into the ANOVA - but you need to think hard about his model because it needs to present nicely paired data to ANOVA.
• You could of course do the Friedman test and its up to the reviewer to understand the theoretical criticism of this (I don't remember it and generally the quantitive standards in molecular peer-review isn't great)
• You could use multiple KW

Alternative the 'correct' way to get ANOVA working when Fmax says 'no' is transformation. This gets complicated however and is one step before a custom model.

Final points

1. Exception If this is microarray data - there is an issue about standardisation and that is complicated (comparing different microarray slides against experimental variation).
2. What algorithms such as DSeq2 do will likely do is apply a fixed transformation to the data, but I don't know about standardisation. Non-parametric statistics will circumvent the requirement for transformation ... their whole point is you don't transform and must not transform - must be raw data.
3. Again non-parametrics become difficult in replicating 2-way ANOVA ... there are several approaches here and of course there is nothing to stop you using it and leave to someone else explain the limitation.

Note the Wilcoxon is the non-parametric equivalent of the T-test, which I think we covered in a previous post (these are mice right?).

What would I do? Personally, if a non-parametric fits I use it every time ... its simple, everyone understand it and a significant result is not in dispute. The reality is you are probably into a skewed distribution. Furthermore, no-one gets worried about a few KW tests (type 2 error) and you can always deploy bonferoni's correction.

I should understand the limitation of Friedman test and I might be working under a false assumption. If a simple transformation results in good Fmax I'd use ANOVA (if raw data isn't homogeneous). Beyond its statistical modelling and that's complicated.