Comparing standardised values of microbial colony perimeters

Question

I’m having a statistical problem (a rather major one) and I was wondering if you could help. I’m researching microbial chemotaxis and analysing colony perimeters by scanning their fluorescence. Chemotaxis is the ability of a microbial cell to swim towards an attractant, in this case succinate. I’ve a set of ~10 scans (A,B,C..) of these colonies. with some having errors.

I’m comparing the perimeter (cm) of:

Alpha chemotactic mono-culture. Beta non-chemotactic mono-culture. Gamma co-culture. I have two other strain combinations (Delta, Epsilon), however these are the ones I’m focussing on now. In three different treatment conditions; 0 mM sugar, 1 mM sugar and 10 mM sugar (succinate).

Statistically comparing the perimeters of three different strain mixtures;

Alpha vs. Beta Alpha vs. Gamma Beta vs. Gamma So far I’ve tried comparing the raw values (cm) which are generally non-normally distributed; upon inspection of density plots and using a Shapiro-Wilks test. I’ve then used a Wilcox rank test to compare them (Fig. 1), I’ve done this with and without outliers removed after assessing their IQR (Fig. 2, Unpair.xlsx).

There is quite a lot of variability, so I decided to standardise the perimeter values of Beta, Gamma (Delta, Epsilon) against the Alpha chemotactic control strain within each scan/batch (A,B,C...). This requires pairs within each scan and successful Alpha scans, if there was an issue with Alpha strain growth; such as contamination, successful scans of Beta for example within the same batch were deleted (Fig. 3). I then removed outliers based on their IQR (Fig. 4) and conducted standardisation: eg. Example: (Alpha perimeter: 7cm / Alpha perimeter: 7cm = 1 and Beta perimeter: 2 cm/ Alpha perimeter: 7 cm = 0.286) (Fig. 5, Pair.xlsx).

Unpair Excel columns: UniqueID:, Bacteria, Treatment, Date, Perimeter (cm), Scan, Biological samples (this is the exact culture number I used). Pair Excel columns: UniqueID:, Bacteria, Treatment, Date, Perimeter (cm), Standardised (cm), Scan, Biological samples (this is the exact culture number I used).

Fig. 1-4.: a), b), c) Treatment comparison. Fig. 1-5.: d), e), f) Strain comparison. Link: https://drive.google.com/file/d/113-u7nVtAuCyVpBUU-cyW_mwvJn4mHAI/view?usp=sharing

As you can see from the figures the results are fairly inconsistent and I’m unsure if the standardisation method I applied works, so I was wondering your take on the most robust statistical test I can use to analyse the strain and treatments.

For example, if I conduct a Wilcox rank test of the Alpha vs. Beta there are all ties, because every Beta standardised is being compared against 1 (Alpha). Is there any kind of paired batch effect standardisation across scans, biological samples or comparison to the mean or average values you’d suggest, or any smoothing techniques I could apply to the perimeter (cm) values. Or any kind of way I could say model the likelihood of one perimeter being greater than the other ?

Answer 1

So far I’ve tried comparing the raw values (cm) which are generally non-normally distributed...

The raw values don't have to be normally distributed. A strict requirement (typically too strict) is that the residuals between observations and modeled values should be normally distributed.

What's important for statistical interpretation of regression coefficients is that the distributions of the coefficient estimates are normally distributed. A normal distribution of error terms is sufficient but not necessary for that. Having no associations between residual distributions and modeled outcome values (homoscedasticity) in a large enough study is often good enough. See this page for more details.

Often in biology and biochemistry the magnitudes of residuals tend to increase as a function of modeled outcomes. That can happen if error magnitudes are proportional to observed values instead of constant. A log transformation of the outcome values can sometimes solve that problem.

Is there any kind of paired batch effect standardisation across scans, biological samples or comparison to the mean or average values you’d suggest, or any smoothing techniques I could apply to the perimeter (cm) values. Or any kind of way I could say model the likelihood of one perimeter being greater than the other ?

A mixed model, in this case with biological sample as a random effect, is one good way to deal with systematic differences among samples that might need "standardization." Random intercepts allow for variation among biological samples in terms of the estimated baseline outcomes (here, at 0 mM sugar). They also can deal with missing data for particular combinations of samples and treatments. That's better than removing all observations from a biological sample just because of contamination in one trial involving it.

These considerations seem to solve your problems. I show a start on your data below. You should work this through on your own, make sure that you understand what each step involves, and incorporate your understanding of the subject matter if there's something else that needs to be addressed. There also are more extensive tests of mixed-model quality, for example in the R DHARMa package, than what I describe.

Start at modeling

I took your data and changed some of the column names to fit better into R data frames. With only 3 treatment levels, it's best to model treatment with a categorical factor.

colonyData <- read.delim("colonyData.txt")
## set "Treatment" to factor called "Sugar" with values "0", "1", "10"
colonyData$Sugar <- factor(colonyData$Sugar)
colonyData$BioSample <- factor(colonyData$BioSample)
colonyData$Bacterium <- factor(colonyData$Bacterium)

I tried a simple mixed model without transforming outcomes. The interaction * allows the response to sugar to differ among bacterial strains.

library(lme4)
lme1 <- lmer(Perimeter~Sugar*Bacterium + (1|BioSample),data=colonyData)
plot(lme1) ## not shown; suggests increasing residual magnitude with modeled values

That plot indicated that the magnitudes of residuals tended to increase with modeled values. Working with log-transformed perimeter values worked better.

lme2 <- lmer(log(Perimeter)~Sugar*Bacterium + (1|BioSample),data=colonyData)
plot(lme2) ## not shown; much better

When there are multiple levels of categorical predictors then the usual model summary (not shown here) can be difficult to interpret. For a categorical predictor, it displays coefficients for the difference between each of the individual factor levels and the reference level. Thus the apparent "significance" of one level can depend on the choice of the reference level.

Use post-modeling tools to estimate the combined significance of all levels of a categorical predictor. The standard R anova() function doesn't handle unbalanced data well. The Anova() function in the R car package is one good alternative.

car::Anova(lme2)
# Analysis of Deviance Table (Type II Wald chisquare tests)
# 
# Response: log(Perimeter)
#                    Chisq Df Pr(>Chisq)    
# Sugar            59.9987  2  9.364e-14 
# Bacterium       133.7860  4  < 2.2e-16 
# Sugar:Bacterium   7.8283  8     0.4504    

This indicates that there are differences among levels of Sugar (treatment) and among bacterial strains. The overall Sugar:Bacterium interaction isn't "significant" but that doesn't mean that it's necessarily unimportant. That's illustrated by detailed analysis of the model predictions.

The emmeans package can provide reports of detailed model predictions. Its "revparwise" comparison method, in this case, evaluates all 3 differences among Sugar levels for each of your bacterial strains. The type="response" specification lets these differences be expressed in terms of perimeter ratios. That makes sense for a model based on log-transformed perimeter values, as a difference in logs is the log of a corresponding ratio.

emm2pairwise <- emmeans(lme2,revpairwise~Sugar|Bacterium, type="response")
emm2pairwise$contrasts
# Bacterium = Alpha:
#  contrast         ratio    SE   df null t.ratio p.value
#  Sugar1 / Sugar0   1.87 0.306 86.9    1   3.819  0.0007
#  Sugar10 / Sugar0  2.19 0.359 86.9    1   4.800  <.0001
#  Sugar10 / Sugar1  1.17 0.198 88.8    1   0.951  0.6095
# 
# Bacterium = Beta:
#  contrast         ratio    SE   df null t.ratio p.value
#  Sugar1 / Sugar0   1.32 0.228 87.1    1   1.625  0.2407
#  Sugar10 / Sugar0  1.62 0.270 85.3    1   2.912  0.0126
#  Sugar10 / Sugar1  1.23 0.211 87.1    1   1.188  0.4634
# 
# Bacterium = Delta:
#  contrast         ratio    SE   df null t.ratio p.value
#  Sugar1 / Sugar0   1.45 0.250 87.1    1   2.175  0.0812
#  Sugar10 / Sugar0  1.49 0.266 89.3    1   2.248  0.0688
#  Sugar10 / Sugar1  1.03 0.177 87.1    1   0.151  0.9874
# 
# Bacterium = Epsilon:
#  contrast         ratio    SE   df null t.ratio p.value
#  Sugar1 / Sugar0   1.35 0.233 87.1    1   1.758  0.1899
#  Sugar10 / Sugar0  1.37 0.237 87.1    1   1.850  0.1597
#  Sugar10 / Sugar1  1.02 0.169 85.3    1   0.095  0.9950
# 
# Bacterium = Gamma:
#  contrast         ratio    SE   df null t.ratio p.value
#  Sugar1 / Sugar0   1.92 0.342 89.3    1   3.661  0.0012
#  Sugar10 / Sugar0  2.19 0.389 88.2    1   4.414  0.0001
#  Sugar10 / Sugar1  1.14 0.203 88.2    1   0.742  0.7393
# 
# Degrees-of-freedom method: kenward-roger 
# P value adjustment: tukey method for comparing a family of 3 estimates 
# Tests are performed on the log scale 

The report incorporates an appropriate correction for multiple comparisons within each bacterial strain.

No strain showed a "statistically significant" difference at p < 0.05 between 1 mM and 10 mM sugar. Only 2 strains showed such a difference between 1 and 0 mM sugar, but a third showed such a difference between 10 and 0 mM. With this size of study, it looks like you need a perimeter ratio of about 1.5 to meet that standard (if arbitrary) criterion of "significance."

The above doesn't deal with differences among assay dates. In principle, you could include them as random effects, also.