Twelve steps from a table of numbers to an answer to the question “what happened to the fly”. Each step opens with a picture explaining the idea, then a few lines of real R that run inside the browser. The code can be edited: every change changes the result, and that is the only way to learn what each line does.
The experiment. A fly had its Adar gene switched off — the enzyme that edits the cell’s own RNA. Six flies: three normal and three without that gene. The question: what changed in the fly? The answer has to be dug out of the numbers.
R has to be downloaded once — about 13 MB and a few seconds. It happens when you run the first step. After that everything runs instantly and offline.
Six flies: three normal and three without the Adar gene. All the RNA was extracted from each, the machine read millions of short fragments, and every fragment was matched to a gene. The result is a table: a row is a gene, a column is a fly, the number is how many fragments were counted.
A zero here does not mean “the gene is absent”, but “none of the millions of fragments landed on it”.
What to see. Adar is the gene that was broken: 2860 reads in normal flies and 5 in the mutants. DptA is an antimicrobial protein: zero in normal flies, hundreds in the mutants. That “zero against hundreds” is where the whole story begins.
For one fly the machine read 6.9 million fragments, for another 21.8. That is about how much material got into the machine, not about biology. Comparing the columns as they are would measure the machine.
What to see. Blue are the normal flies, red the mutants. Notice that tall and short bars appear in both groups — so this is definitely not biology.
The cure is division. Divide every number by its own column total and multiply by a million — that gives CPM, “how many reads of this gene there would be per million”. Now the columns are comparable.
t() flips the table. It appears twice not for decoration: R divides a matrix row-wise, and we need it column-wise — so the table is flipped, divided and flipped back.
What to see. It was 4478 to 12740 — almost threefold. Now it is 519 to 882. The scatter has not vanished entirely, but actin no longer looks as if half the flies run it three times harder.
More than a third of the genes in a fly head barely work: zero, two, zero. Nothing can be said about such a gene, and in the pile it does harm: every extra gene is one more chance to hit something by accident.
The rule: a gene stays if its CPM is at least one in at least three flies. Three, because each group has three flies.
What to see. Try >= 5 instead of >= 1 and run through to the end. The list of significant genes will change. The threshold is a human decision, not a property of the data, and papers are required to state it.
Take the mean over the three normal flies, the mean over the three mutants, and subtract. A difference of logarithms is a ratio: 1 means twofold, 3 means eightfold, -1 means half.
What to see. At the top of the list: Dro, Mtk, AttC, DptB, DptA. These are all antimicrobial peptides — proteins a fly uses to kill bacteria. And Adar dropped 217-fold, exactly as it should: it was the gene that was disabled.
Picture a cannon. Even if nobody touches it, the shots never land on one point — there is scatter. Now the barrel is turned a little and the volley lands higher. How do you tell that the barrel really moved, rather than the shots landing higher by chance?
Compare the shift with the scatter. That is t: how many times bigger the shift is than the scatter. And p answers the second question — how often a motionless cannon would produce the same shift by itself. Three flies in a group are one volley.
The schema says the same thing as a formula. On top is the shift: the mutant mean minus the normal mean. Below is the scatter converted to that shift: the more flies in a group, the less the shift wanders on its own, which is where the division by three comes from. t.test computes exactly this fraction.
The function prints a whole paragraph rather than a single number — R likes to report everything it computed. Six lines matter here, and this is what each of them says.
What to see. Replace DptA with Act5C in the code — an ordinary working gene. t falls to −0.25 and p-value rises to 0.82: chance alone produces a shift like that eight times out of ten. This is what a gene which felt nothing looks like.
Calling t.test ten thousand times is slow, so the same formula is written with matrices. The first line is the scatter: how much the three flies within a group disagree.
Why the second line — “borrowing from the neighbours”. There are only three flies. Sometimes three numbers happen to agree almost exactly, the scatter comes out tiny, and in the “shift ÷ scatter” formula a trivial change is crowned champion — simply because we divided by almost nothing.
The cure: take half of the gene’s own scatter and half of the typical scatter across all ten thousand genes. A gene with an implausibly small scatter gets a dose of realism; a gene with a huge one is toned down.
What to see. Delete the second line and rerun — “significant” drops to 962 from 1163. Those two hundred genes are the price of having only three flies per group. And compare the two numbers in the output: 1163 against 499. Almost half of the “significant” genes are pure luck.
The correction does not subtract those five hundred — it does not know which genes are the accidental ones. Instead it demands stronger evidence from each: the gene in position k passes only if its p is no more than k × 0.05 / n.
So for the 367th gene the threshold is not 0.05 but 0.0018 — nearly thirty times stricter. All of it is done by one function, p.adjust.
What to see. Three numbers: 1163 looked significant; 367 survived the correction; 353 remained once we also demanded at least a twofold change. That last requirement is not statistics but biology: a 3 % change can be impeccably certain and still useless.
The classic picture of the field. Horizontally, how much a gene changed; vertically, how sure we are. Red went up, blue went down, grey did not pass.
What to see. 278 genes up, 75 down. Act5C and ninaE sit at the bottom in the middle — exactly where genes that nothing happened to belong.
Three hundred names like AttC and CG3270 mean nothing on their own. Gene Ontology is a huge list of “gene — what it does” pairs, curated by hand from thousands of papers.
The question is simple: if we picked 300 genes at random, how many would be about immunity by accident? phyper answers that — it is balls in an urn, no biology involved.
What to see. The first line is humoral immune response: 22 of our genes out of the 42 in that process. The chance of that is 9 in 1024. Nobody injected the fly with bacteria — it raised its defences on its own, because its own unedited RNA started to look foreign.
A six-row table is easy to read; a twenty-row one is not. So enrichment is almost always drawn as a single figure carrying three quantities at once: horizontally, what share of our genes landed in the process; dot size, how many of them; colour, how certain it is.
In real work this is dotplot() from the clusterProfiler package. Here the same thing is built out of plot and axis — to show there is no magic inside.
What to see. The top three rows — humoral immune response, antimicrobial humoral response, antibacterial humoral response — are nearly the same set of genes nested inside itself: Gene Ontology is a tree, and a gene belongs both to the narrow term and to everything above it. So the top of such a figure is often filled with synonyms, and it should be read not as twelve discoveries but as one conclusion: immunity. Here the near-duplicates were collapsed in advance, otherwise there would be twice as many.
The last step is to see a picture rather than a list. Two genes get a line if they share many processes. The immune genes gather into a dense clot, everything else stays at the edges.
M %*% t(M) multiplies the table by itself. Each cell ends up holding the number of processes a pair of genes shares. One line instead of a double loop.
What to see. Try >= 3 instead of >= 5 — many more lines and the picture turns to mush. Try >= 8 and only the densest core survives.
Nobody analyses RNA-seq the way you just did. In real work everything above is done by a single package — DESeq2:
library(DESeq2)
counts <- as.matrix(read.csv("muhi.csv", row.names = 1))
info <- data.frame(row.names = colnames(counts),
group = factor(c(rep("norm", 3), rep("mut", 3)),
levels = c("norm", "mut")))
dds <- DESeqDataSetFromMatrix(counts, info, design = ~ group)
dds <- DESeq(dds)
res <- results(dds, contrast = c("group", "mut", "norm"))
DESeqDataSetFromMatrix brings together the counts, the sample table and the design formula ~ group — “explain the differences by group”. If the flies also differed in sex, one would write ~ sex + group and the sex difference would stop interfering. That is the part we cannot do by hand.
DESeq() is three steps, and you did all three: normalisation (step 3), scatter with borrowing (step 7) and the test (step 6). results() pulls out the table and applies the BH correction (step 8).
This block cannot be run in the browser: DESeq2 is written in C++ and nobody has built it for the browser. So we ran it beforehand on an ordinary computer and put its answer here as a table, read in the next block. The full script is at the bottom of the page; at home it runs as is.
The main question about any simplified method: does it get things wrong, or does it just see less? Two lines settle it.
What to see. “Ours only” is zero. All 353 of our genes are in the DESeq2 list, and it found 325 more that were beyond us. The simple method invented nothing: it misses real things but does not produce false ones. Of the two kinds of error, that is the one you can live with.
Every variable is still there: counts, cpm, lg,
fc, p, q, sig, go,
D. Ask anything you like.
R gives every answer — run the step in question and read the output. The numbers hold as long as you have not edited the code (“reset code” puts it back).
Everything you wrote is base R: not a single installed package. Real analyses do the same through DESeq2 or edgeR — more carefully, but by the same logic: normalise, filter, difference, scatter, correction, gene sets.
The same story about the mouse — why losing editing kills it while the fly survives — is on the neighbouring tab.
Take the files and run them at home in RStudio: muhi.csv, go.csv, deseq_muha.csv, deseq.R.