mageck2 mle design matrix meaning of multiple baselines

[carlosj-rr]

Hello,
I want to run the mle command but I am getting confused with the documentation [edit: of the design matrix]. I can understand that for an experiment column, setting more than one row to “1” will return the beta scores for the combined samples with 1s in relation to the baseline sample (the one with a 1 on the baseline column and zeroes everywhere else). However, on the example documentation for the mle subcommand, there are two baselines and it is not clear how MAGeCK2 treats them - does it know they pertain to different cell lines? How? Or does it join the two .initial samples and uses that as a baseline?

Many thanks both for developing this tool and for your time and attention. Please do let me know if you need more information on my question.

Edit: I’ve also been through the Advanced Design Matrix Tutorial but all the examples there use a single baseline.

[davidliwei]

Good question — your second guess is right: MAGeCK2 pools the two initial samples into one shared baseline. It has no concept of cell line; the HL60/KBM7 names are just labels it never uses for grouping.

The rule is purely mechanical: any sample whose design-matrix row sums to 1 (a 1 in baseline, 0 elsewhere) is treated as a baseline sample. Both HL60.initial and KBM7.initial qualify, so the baseline reference is estimated as the average of the two. Then:

• HL60|beta = change in HL60.final vs. that shared baseline,
• KBM7|beta = change in KBM7.final vs. the same shared baseline.

This is intended for the demo's setup, where both lines come from the same library plasmid pool, so their day-0 representations should match and pooling just stabilizes the reference.

If your two cell lines instead have genuinely different baselines (different libraries, or you don't want to assume a shared start), don't pool them — run mle separately for each line, each with its own initial as the baseline. Tell me which case you're in and I'll confirm the right matrix.

[carlosj-rr]

Hi @davidliwei, thanks a lot. That clarifies it all for me, thank you! I have a combination of the two cases you mention. I am testing 4 genotypes, same library, across 3 time points. I’m interested in each genotype’s beta values (like the last case you explain).

However, my collaborator is also asking to compare different cases against one baseline. Briefly, our lines have mutations on either of two genes, gene1 or gene2 (the 4 lines being WT, gene1-mut, gene2-mut, and double mut), and we wish to identify which gene dropouts that happen in a double mutant do not happen in a single mutant of one of the genes. Thinking about this, I assumed a matrix in which the WT was the baseline, then the samples to be analysed would be the the mutant of the gene of interest and the double mutant?

Please let me know if I need to be more clear.

[davidliwei]

Great — that's a clean 2x2 factorial, and it's worth separating your two goals, because they want different design matrices.

Goal 1: each genotype's own beta (vs. day 0)

--day0-label builds this for you, no design matrix needed:

mageck2 mle -k counts.txt --day0-label T0 -n perGenotype

One caveat: it emits one beta per sample, not per genotype — replicates and time points are not pooled. If you want a single beta per genotype, write the matrix by hand and give the samples that share a genotype the same row.

Goal 2: dropouts in the double mutant that aren't in either single

Your instinct (WT as baseline, singles + double as columns) gets you differential essentiality vs. WT, but it can't quite answer the question your collaborator is asking. It will flag any gene that drops out more in the double mutant than in WT — including genes where the double mutant is simply the sum of the two single-mutant effects. What you want is the subset where the double mutant deviates from that additive expectation, which is a genetic-interaction term.

Encode it as an explicit interaction column:

Samples   baseline  WT  g1_vs_WT  g2_vs_WT  interaction
T0            1      0      0         0          0
WT            1      1      0         0          0
gene1_mut     1      1      1         0          0
gene2_mut     1      1      0         1          0
double_mut    1      1      1         1          1

Repeat a row verbatim for each replicate of that genotype. Then:

• WT|beta — essentiality in WT, vs. T0 (your common essentials)
• g1_vs_WT|beta — extra effect of the gene1 mutation
• g2_vs_WT|beta — extra effect of the gene2 mutation
• interaction|beta — the column you want. How far the double mutant departs from what the two singles predict. Strongly negative = drops out in the double mutant beyond the additive expectation.

The advantage over comparing betas from a --day0-label run is that interaction is a fitted coefficient, so MAGeCK gives it its own z, p-value and FDR directly — there's no valid significance test for eyeballing the difference between two separately-fitted betas.

Note betas are on the natural-log scale, so a 4-fold dropout is about -1.39.

Your three time points

The simplest starting point is to run the matrix above once per time point (T0 plus the four genotypes at that time), giving you three sets of betas so you can watch the interaction emerge over time.

You can also fold all three time points into a single fit. The matrix gets more complicated, but it's a supported design — see section 9 of the tutorial, which covers exactly this kind of time-series encoding. Following "Option C" there (time as a linear covariate), your factorial version would be:

Samples   baseline  WT  g1_vs_WT  g2_vs_WT  interaction
T0            1      0      0         0          0
WT_t1         1      1      0         0          0
WT_t2         1      2      0         0          0
WT_t3         1      3      0         0          0
g1_t1         1      1      1         0          0
g1_t2         1      2      2         0          0
g1_t3         1      3      3         0          0
g2_t1         1      1      0         1          0
g2_t2         1      2      0         2          0
g2_t3         1      3      0         3          0
dm_t1         1      1      1         1          1
dm_t2         1      2      2         2          2
dm_t3         1      3      3         3          3

Each beta is now a per-time-point rate rather than a total effect, estimated from all 12 samples at once, which buys you power. The assumption is that the log fold change scales linearly with time point — worth confirming with the per-time-point runs first before committing to it. If your time points aren't evenly spaced, use the actual doubling times or days in place of 1/2/3.