Six flies: three normal and three without the Adar gene. For each one we know how many reads landed on each of 16,642 genes. The task is to work out what changed in the fly. Everything is computed on this page from the raw numbers: nothing is prepared in advance, and the result changes when you move the sliders.
A row is a gene, a column is one fly. The number is how many pieces of that gene's RNA reached the machine. Not molecules — reads. A zero means not «there is no RNA at all», but «none of the millions of reads landed on it».
Find three genes: Adar — the one that was switched off; DptA — an antimicrobial peptide; Act5C — actin, an ordinary housekeeping gene that should not change.
counts <- read.csv("muhi.csv", row.names = 1)
dim(counts)
counts[c("Adar", "DptA", "Act5C"), ]They cannot. One fly gave 6.9 million reads, another 21.8 — more than a threefold difference, and it is about the machine, not about biology.
How to think about it. Two classes sat the same test. Class A solved 18 problems, class B solved 30. Class B looks stronger. But class A has six pupils and class B fifteen: per pupil that is 3 against 2, and class A is actually stronger. The number of reads in a sample is «how many pupils are in the class». So every number is divided by its column total and multiplied by a million. That gives CPM — «how many reads this gene would get per million».
The rhodopsin ninaE looks as if it had doubled. After normalisation you can see it did not change: what grew was the file size, not the gene.
colSums(counts) # depth of each sample
cpm <- t(t(counts) / colSums(counts)) * 1e6For a gene with three reads, «doubled» means it now has six. That is noise. Such genes are thrown out before any statistics — otherwise they fill the list with random jumps.
keep <- rowSums(counts) >= 10
counts <- counts[keep, ]
l <- log2(cpm[keep, ] + 1) # see step 4 on the logarithmWhy a logarithm. The «times» scale is lopsided: twice as much is 2,
half as much is 0.5. Two is one unit away from one, but 0.5 is only half a unit away,
although the change means the same thing. log2 counts not «how many times
more» but «how many doublings»: +1 is doubled, −1 is halved, 0 is unchanged. Now a rise
and a fall weigh the same.
At the bottom of the «went down» list you should find Adar itself. That is the best check that everything was done right: the gene that was switched off in the fly must switch off in our numbers too.
m1 <- rowMeans(l[, 1:3]) # normal flies m2 <- rowMeans(l[, 4:6]) # flies without Adar lfc <- m2 - m1 head(sort(lfc, decreasing = TRUE), 10)
There are only six flies. Even for a gene that did not change, three numbers rarely match three others exactly. A real difference has to be told apart from scatter.
What p is. Suppose this gene did not change at all, and the whole difference between the six numbers is natural scatter. Then it does not matter which three flies we call mutants: any three of the six will do. p is the share of such reshuffles in which the difference comes out at least as large as the one we saw. A small p means «it is hard to get this by chance».
For Def the gene went up several fold, but the points inside each group are more scattered than the groups differ from each other — so p is large. A big change and a reliable change are not the same thing.
# row-wise, without apply — otherwise 16 seconds instead of 40 milliseconds
s2 <- (rowSums((l[, 1:3] - m1)^2) + rowSums((l[, 4:6] - m2)^2)) / 4
t <- lfc / sqrt(s2 * (1/3 + 1/3))
p <- 2 * pt(-abs(t), df = 4)Three numbers give a poor estimate of a gene's scatter: the estimate itself is noisy. A gene whose three numbers happened to land close together will get a tiny scatter and slip into the list for no good reason.
The way out, used by every real program: borrow the missing precision from the other genes. Most genes have roughly the same scatter, so one takes the average of the gene's own scatter and the typical scatter across all eleven thousand.
At the left end of the slider it is the ordinary test, which finds almost nothing on three replicates. At the right end only the common scatter is used, which is too crude. Real packages sit in the middle, and so do we.
s2 <- (s2 + median(s2)) / 2 # half is borrowed from the neighbours t <- lfc / sqrt(s2 * (1/3 + 1/3)) p <- 2 * pt(-abs(t), df = 8) # borrowing added degrees of freedom
Why a correction is needed. A threshold of p < 0.05 means «I accept being wrong 5 % of the time». But the question was not asked once — it was asked for every gene. 5 % of 11,855 is 593 genes that will slip into the list even if nothing at all changed in the flies. The Benjamini—Hochberg correction rescales the thresholds so that the share of junk among the selected genes stays below what you allow. The rescaled number is called q.
q <- p.adjust(p, method = "BH") sig <- names(which(q < 0.05 & abs(lfc) > 1)) length(sig)
A list of several hundred names says nothing on its own. There is a database called Gene Ontology: for every gene it records which processes that gene takes part in. We look for the process that shows up among the selected genes more often than it would in a random selection of the same size.
The «expected» column is how many genes of that process would land in the list by chance. If we found 22 and expected 1.5, the process is genuinely affected.
go <- read.csv("go.csv") # gene — process pairs
tab <- table(go$process[go$gene %in% sig])
# Fisher's exact test for one process:
fisher.test(matrix(c(22, 27, 329, 11477), 2), alternative = "greater")$p.valueNow the same thing as a picture. A dot is a gene. A line between two genes means they take part in shared processes. Genes that work together gather into a clump on their own.
Hover over a dot to see the gene name and its processes.
Everything above was computed by JavaScript right in the browser. But the same analysis can be run in real R — it also works inside a page. The first load is about 13 MB and takes a few seconds; after that R comes from the cache.
The data and the script can be taken away and run on your own machine: muhi.csv, go.csv, urok.R.
A fly has no interferon, and losing Adar does not kill it. A mouse does have interferon, and a mouse without editing dies as an embryo. The accepted explanation is that what kills is not the missing edits themselves but the sensor MDA5, which mistakes unedited RNA for viral RNA and raises the alarm.
This can be checked: take a mouse whose enzyme is present but dead, and compare two cases — the sensor present, and the sensor removed as well. Run the same pipeline with both buttons and compare the numbers.
The answers are computed from the same numbers as everything else on the page, so the «right answer» cannot disagree with what you saw.