Двенадцать шагов от таблицы с числами до ответа на вопрос «что случилось с мухой». Перед каждым шагом — картинка, объясняющая идею; потом несколько строк настоящего R, который выполняется прямо в браузере. Код можно менять: любая правка меняет результат, и это единственный способ понять, что каждая строчка делает.
Опыт, который мы разбираем. У мухи выключили ген Adar — фермент, который правит собственную РНК клетки. Шесть мух: три обычных и три без этого гена. Вопрос: что в мухе изменилось? Ответ придётся достать из чисел самим.
R надо загрузить один раз — около 13 МБ и несколько секунд. Это происходит при первом запуске любого шага. Дальше всё считается мгновенно и без интернета.
Шесть мух: три обычных и три без гена Adar. Из каждой достали всю РНК, прибор прочёл миллионы коротких обрывков, и каждый обрывок опознали — от какого он гена. Получилась таблица: строка — ген, столбец — муха, число — сколько обрывков насчитали.
Ноль в такой таблице означает не «гена нет», а «ни один из миллионов обрывков на него не попал».
Что видно. Adar — тот ген, который сломали: 2860 прочтений у обычных мух и 5 у мутантов. DptA — противомикробный белок: у обычных ноль, у мутантов сотни. Вот с этого «ноль против сотен» и начнётся вся история.
У одной мухи прибор прочёл 6.9 миллиона обрывков, у другой — 21.8. Это про то, сколько материала попало в прибор, а не про биологию. Если сравнить столбцы как есть, мы измерим работу прибора.
Что видно. Синие — обычные мухи, красные — мутанты. Обратите внимание: и высокие, и низкие столбики есть в обеих группах. Значит дело точно не в биологии.
Лечится это делением. Каждое число делим на итог по своему столбцу и умножаем на миллион — получается CPM, «сколько прочтений этого гена пришлось бы на миллион». Теперь столбцы сравнимы.
t() переворачивает таблицу. Он написан дважды не для красоты: R делит матрицу по строкам, а нам нужно по столбцам — поэтому таблицу переворачивают, делят и переворачивают обратно.
Что видно. Было от 4478 до 12740 — почти втрое. Стало от 519 до 882. Разброс не исчез совсем, но актин перестал выглядеть так, будто у половины мух он работает втрое сильнее.
Больше трети генов в голове мухи почти не работают: ноль, два, ноль. Про такой ген нельзя сказать, вырос он или упал. А в общей куче он вредит: каждый лишний ген — ещё одна возможность случайно «попасть».
Правило: ген остаётся, если хотя бы у трёх мух его CPM не ниже единицы. Три — потому что в каждой группе три мухи.
Что видно. Поставьте >= 5 вместо >= 1 и пройдите страницу до конца. Список значимых генов изменится. Порог — решение человека, а не свойство данных, и в статьях его обязаны указывать.
Считаем среднее по трём обычным мухам, среднее по трём мутантам и вычитаем. Разница логарифмов — это отношение: 1 означает «вдвое», 3 — «в восемь раз», -1 — «вдвое меньше».
Что видно. В верхушке списка — Dro, Mtk, AttC, DptB, DptA. Это всё противомикробные пептиды: белки, которыми муха убивает бактерий. А Adar упал в 217 раз — так и должно быть, его же выключили.
Представьте пушку. Даже если её не трогать, снаряды ложатся не в одну точку — есть разброс. Теперь ствол чуть повернули, и залп лёг выше. Как понять, что ствол правда повернули, а не снаряды сами случайно легли выше?
Сравнить сдвиг с разбросом. Это и есть t: во сколько раз сдвиг больше разброса. А p отвечает на второй вопрос — как часто неподвижная пушка дала бы такой же сдвиг сама по себе. Три мухи в группе — это и есть один залп.
На схеме та же мысль записана формулой. Сверху — сдвиг: среднее у мутантов минус среднее у обычных мух. Снизу — разброс, пересчитанный на этот сдвиг: чем больше мух в группе, тем меньше сдвиг гуляет сам по себе, поэтому там и появляется деление на три. t.test считает ровно эту дробь.
Функция печатает не одно число, а целый абзац — R вообще любит рассказывать всё, что посчитал. Полезных строк тут шесть, и вот что означает каждая.
Что видно. Замените в коде DptA на Act5C — обычный рабочий ген. t упадёт до −0.25, а p-value вырастет до 0.82: такой сдвиг случайность выдаёт в восьми случаях из десяти. Так выглядит ген, который ничего не почувствовал.
Вызывать t.test десять тысяч раз медленно, поэтому ту же формулу пишут матрицами. Первая строка кода — это разброс: насколько три мухи внутри группы расходятся между собой.
Зачем нужна вторая строка — «заём у соседей». Мух всего три. Иногда три числа случайно совпадают почти точно, разброс выходит крошечным, и по формуле «сдвиг ÷ разброс» ерундовое изменение вылезает в чемпионы — просто потому, что мы делим на почти ноль.
Лечится это так: берут половину собственного разброса гена и половину типичного разброса по всем десяти тысячам генов. Гену с невероятно маленьким разбросом добавляют реализма, гену с огромным — наоборот, сбавляют.
Что видно. Уберите вторую строку совсем и запустите заново — «значимых» станет 962 вместо 1163. Разница в двести генов и есть цена трёх мух на группу. Сравните два числа в выводе: 1163 и 499. Почти половина «значимых» — просто везение.
Поправка не вычитает эти пятьсот из списка — она не знает, какие именно гены случайные. Вместо этого она требует от каждого гена доказательства посильнее: ген на месте k проходит, только если его p не больше, чем k × 0.05 / n.
Поэтому для 367-го гена порог оказывается не 0.05, а 0.0018 — почти в тридцать раз строже. Всё это делает одна функция p.adjust.
Что видно. Три числа: 1163 — сколько казалось значимым; 367 — сколько выжило после поправки; 353 — сколько осталось, когда мы вдобавок потребовали изменения хотя бы вдвое. Последнее требование — уже не статистика, а биология: изменение на 3 % может быть безупречно достоверным и при этом никому не нужным.
Классическая картинка всей области. По горизонтали — во сколько раз изменился ген, по вертикали — насколько мы в этом уверены. Красное — выросло, синее — упало, серое — не прошло порог.
Что видно. 278 генов вверх и 75 вниз. Act5C и ninaE сидят внизу посередине — там, где и положено генам, с которыми ничего не случилось.
Триста имён вроде AttC и CG3270 сами по себе ничего не значат. Gene Ontology — огромный список пар «ген — чем он занимается», собранный вручную по тысячам статей.
Вопрос к нему простой: если бы мы выбрали 300 генов наугад, сколько из них случайно оказались бы про иммунитет? На это отвечает phyper — это задача про шары в урне, никакой биологии.
Что видно. Первая строка — humoral immune response, «гуморальный иммунный ответ»: 22 наших гена из 42, что есть в этом процессе. Вероятность такого совпадения — 9 на 1024. Мухе никто не вводил бактерий: она подняла защиту сама, потому что своя неотредактированная РНК стала похожа на чужую.
Таблицу с шестью строками читать легко, а с двадцатью — уже нет. Поэтому обогащение почти всегда рисуют одной картинкой, где сразу три величины: по горизонтали — какая доля наших генов попала в процесс, размер точки — сколько их там, цвет — насколько уверенно.
В настоящих работах это делает dotplot() из пакета clusterProfiler. Здесь то же самое собрано из plot и axis — чтобы было видно, что ничего волшебного внутри нет.
Что видно. Три верхние строки — humoral immune response, antimicrobial humoral response, antibacterial humoral response — это почти один и тот же набор генов, вложенный сам в себя: Gene Ontology устроена как дерево, и один ген числится сразу и в узком термине, и во всех, что над ним. Поэтому топ такой картинки часто занят синонимами, и читать её надо не как список из двенадцати открытий, а как один вывод: иммунитет. У нас почти-дубли заранее схлопнуты, иначе их было бы вдвое больше.
Последний шаг — увидеть не список, а картину. Два гена соединяем линией, если у них много общих процессов. Иммунные гены собираются в плотный ком, всё остальное остаётся по краям.
M %*% t(M) — умножение таблицы на саму себя. В каждой клетке получается число процессов, общих у пары генов. Одна строчка вместо двойного цикла.
Что видно. Поставьте >= 3 вместо >= 5 — линий станет намного больше и картинка превратится в кашу. Поставьте >= 8 — останется только самое плотное ядро.
Никто не считает RNA-seq так, как считали вы. В настоящей работе всё, что было выше, делает один пакет — 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 складывает вместе счётчики, таблицу образцов и формулу дизайна ~ group — «объясняй различия группой». Если бы мухи были ещё и разного пола, написали бы ~ sex + group, и различия между полами перестали бы мешать. Вот этого мы руками не умеем.
DESeq() — это три шага, и все три вы делали: нормировка (шаг 3), разброс с займом у соседей (шаг 7) и проверка (шаг 6). results() достаёт таблицу и делает поправку BH (шаг 8).
Запустить этот блок в браузере нельзя: DESeq2 написан на C++ и для браузера его никто не собирал. Поэтому мы прогнали его заранее на обычном компьютере, а сюда положили таблицу с его ответом — она читается в следующем блоке. Скрипт целиком лежит внизу страницы, дома он запустится как есть.
Главный вопрос к любому упрощённому методу: он ошибается или просто видит меньше? Проверяется это в две строчки.
Что видно. «Только у нас» — ноль. Все 353 наших гена есть и у DESeq2, а он нашёл ещё 325, которые нам оказались не по зубам. Простой метод ничего не выдумал: он пропускает настоящее, но не выдаёт ложное. Из двух видов ошибок это та, с которой можно жить.
Все переменные остались на месте: counts, cpm, lg,
fc, p, q, sig, go,
D. Спрашивайте что угодно.
Все ответы даёт сам R — запустите нужный шаг и посмотрите в вывод. Числа верны, пока вы не меняли код (кнопка «вернуть код» возвращает исходный).
Всё, что вы написали, — это базовый R: ни одного установленного пакета. Настоящие работы считают тем же самым через DESeq2 или edgeR — аккуратнее, но по той же логике: нормировка, фильтр, разница, разброс, поправка, наборы генов.
Ту же историю про мышь — почему у неё потеря редактирования смертельна, а муха выживает — смотрите на соседней вкладке.
Забрать файлы и запустить дома в RStudio: muhi.csv, go.csv, deseq_muha.csv, deseq.R.