Шесть мух: три обычных и три без гена Adar. Для каждой посчитано, сколько прочтений попало на каждый из 16 642 генов. Задача — выяснить, что в мухе изменилось. Всё считается прямо на этой странице, из сырых чисел: ничего не заготовлено заранее, и результат меняется, когда вы двигаете ползунки.
Строка — ген, столбец — одна муха. Число — сколько кусочков РНК этого гена попало в прибор. Не молекул, а именно прочтений: ноль означает не «РНК нет совсем», а «ни одно из миллионов прочтений на него не попало».
Найдите три гена: Adar — тот, который выключили; DptA — антимикробный пептид; Act5C — актин, обычный рабочий ген, который меняться не должен.
counts <- read.csv("muhi.csv", row.names = 1)
dim(counts)
counts[c("Adar", "DptA", "Act5C"), ]Нельзя. У одной мухи прочтений собрали 6.9 миллиона, у другой 21.8 — разница более чем втрое, и она про работу прибора, а не про биологию.
Как это понять. Два класса писали контрольную. 9А решил 18 задач, 9Б — 30. Кажется, 9Б сильнее. Но в 9А шесть человек, а в 9Б пятнадцать: на человека выходит 3 против 2, и сильнее на самом деле 9А. Число прочтений в образце — это «сколько человек в классе». Поэтому каждое число делят на итог по столбцу и умножают на миллион. Получается CPM — «сколько прочтений этого гена пришлось бы на миллион».
У родопсина ninaE сырые числа как будто выросли вдвое. После нормировки видно, что он не изменился: выросла не работа гена, а размер файла.
colSums(counts) # глубина каждого образца
cpm <- t(t(counts) / colSums(counts)) * 1e6У гена, на который попало три прочтения, «вырос в два раза» означает, что прочтений стало шесть. Это шум. Такие гены выкидывают до всякой статистики — иначе они засоряют список случайными скачками.
keep <- rowSums(counts) >= 10
counts <- counts[keep, ]
l <- log2(cpm[keep, ] + 1) # см. шаг 4 про логарифмПочему логарифм. Шкала «во сколько раз» кривая: вдвое больше — это
2, вдвое меньше — 0.5. От единицы двойка отстоит на 1, а 0.5 всего на 0.5, хотя по смыслу
изменение одинаковое. log2 считает не «во сколько раз», а
«сколько раз удвоили»: +1 — удвоили, −1 — уполовинили, 0 — не изменилось. Теперь рост и
падение весят одинаково.
Внизу списка упавших должен стоять сам Adar. Это лучшая проверка, что мы всё сделали правильно: ген, который в мухе выключили, обязан выключиться и в наших числах.
m1 <- rowMeans(l[, 1:3]) # обычные мухи m2 <- rowMeans(l[, 4:6]) # мухи без Adar lfc <- m2 - m1 head(sort(lfc, decreasing = TRUE), 10)
Мух всего шесть. Даже если ген вообще не менялся, три числа редко совпадают с тремя другими. Надо отличить настоящую разницу от разброса.
Что такое p. Представим, что этот ген на самом деле не менялся, и вся разница между шестью числами — природный разброс. Тогда неважно, каких трёх мух назвать мутантами: годится любая тройка из шести. p — доля таких перетасовок, при которых разница выходит не меньше наблюдаемой. Маленькое p значит «случайно так разложиться трудно».
У Def ген вырос в несколько раз, но точки внутри групп разбросаны сильнее, чем группы отличаются друг от друга — и p большое. Большое изменение и надёжное изменение это не одно и то же.
# по строкам, без apply — иначе 16 секунд вместо 40 миллисекунд
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)По трём числам разброс гена оценивается плохо: сама оценка получается шумной. У гена, которому случайно повезло и три числа сошлись близко, разброс выйдет крошечным, и он пролезет в список без всяких оснований.
Выход, которым пользуются все настоящие программы: занять недостающую точность у остальных генов. Разброс у большинства генов примерно одинаковый, так что берут среднее между разбросом самого гена и типичным разбросом по всем 11 тысячам.
Слева ползунка — обычный тест, который на трёх повторах почти ничего не находит. Справа — только общий разброс, это уже слишком грубо. Настоящие пакеты стоят посередине, и мы тоже.
s2 <- (s2 + median(s2)) / 2 # половину занимаем у соседей t <- lfc / sqrt(s2 * (1/3 + 1/3)) p <- 2 * pt(-abs(t), df = 8) # занятый разброс добавил степеней свободы
Почему нужна поправка. Порог p < 0.05 означает «согласен ошибаться в 5 % случаев». Но вопрос задан не один раз, а для каждого гена. 5 % от 11 855 — это 593 генов, которые пролезут в список, даже если в мухах не изменилось ровным счётом ничего. Поправка Бенджамини—Хохберга пересчитывает пороги так, чтобы доля мусора среди отобранных не превышала заданной. Пересчитанное число называют q.
q <- p.adjust(p, method = "BH") sig <- names(which(q < 0.05 & abs(lfc) > 1)) length(sig)
Список из сотен названий сам по себе ничего не говорит. Есть база Gene Ontology: для каждого гена там записано, в каких процессах он участвует. Смотрим, какой процесс встречается среди отобранных генов чаще, чем вышло бы при случайном выборе такого же размера.
Колонка «ждали» — сколько генов этого процесса попало бы в список случайно. Если нашли 22, а ждали 1.5 — процесс задет по-настоящему.
go <- read.csv("go.csv") # пары ген — процесс
tab <- table(go$process[go$gene %in% sig])
# точный тест Фишера для одного процесса:
fisher.test(matrix(c(22, 27, 329, 11477), 2), alternative = "greater")$p.valueТеперь то же самое картинкой. Точка — ген. Линия между двумя генами — они участвуют в общих процессах. Гены, которые работают вместе, сами собираются в комок.
Наведите на точку — увидите название гена и его процессы.
Всё, что вы прошли выше, посчитал JavaScript прямо в браузере. Но тот же разбор можно выполнить в настоящем R — он тоже умеет работать внутри страницы. Первая загрузка около 13 МБ и занимает несколько секунд, дальше R берётся из кэша.
Данные и скрипт можно забрать и запустить у себя: muhi.csv, go.csv, urok.R.
У мухи нет интерферона, и потеря Adar её не убивает. У мыши интерферон есть, и мышь без редактирования гибнет ещё эмбрионом. Считается, что убивает не отсутствие правок само по себе, а сенсор MDA5, который принимает неотредактированную РНК за вирусную и поднимает тревогу.
Проверить это можно так: взять мышь, у которой фермент цел, но не работает, и сравнить два случая — когда сенсор на месте и когда его тоже убрали. Прогоните тот же пайплайн обеими кнопками и сравните числа.
Ответы считаются из тех же чисел, что и всё на странице, поэтому «правильный ответ» не может разойтись с тем, что вы видели.