Лекция #3

Анализ bulk RNA-Seq

Серёжа Исаев

​

аспирант MedUni Vienna

«Анализ данных NGS»

Дорожная карта анализа RNA-Seq

2

QC прочтений

Выравнивание

Подсчёт�экспрессий

Псевдовыравнивание

Нормализация

Дифференциальная экспрессия

GO, GSEA, ssGSEA и прочее

Распределение каунтов генов

Экспрессии генов TP53 и EGFR в образцах рака лёгкого

Какое это распределение?

3

Распределение Пуассона

Распределение Пуассона отражает число событий, произошедших за фиксированное время, при условии, что данные события происходят с некоторой фиксированной средней интенсивностью и независимо друг от друга

4

Распределение Пуассона

Представим, что у нас есть бесконечно большая шляпа, в которой есть несколько типов шариков — красные, синие, зелёные, … Сфокусируемся на красном шарике, доля красных шариков 0.01 (то есть вероятность вытащить красный шарик — 1 из 100).

Мы забираем из шляпы 300 шариков, то есть в среднем мы увидим красный шарик 3 раза

Какое будет распределение вероятности различного количества красных шариков, которые мы увидим? Это как раз Пуассон

• Шарики = прочтения
• Цвет шарика = ген

5

Среднее и дисперсия распределения Пуассона

В распределении Пуассона среднее равно дисперсии, а потому достаточно легко понять, если несколько случайных величин распределены по Пуассону

6

Отрицательное биномиальное распределение

Отрицательное биномиальное распределение определяется как количество произошедших неудач в последовательности испытаний Бернулли с вероятностью успеха p, проводимой до r-го успеха.

7

Отрицательное биномиальное распределение

Несложно заметить, что можно таким же образом подсчитать число удач до n-ой неудачи, только теперь в вероятность мы подставим не p, а 1 — p

• Допустим, я беру по одному прочтению из образца X
• Если прочтение будет из гена g, то это успех (число удач = число каунтов гена)
• Если нет, то неудача (число неудач = глубина секвенирования)
• p — вероятность успеха (= экспрессия гена)

8

Отрицательное биномиальное распределение

• Допустим ген имеет не очень высокую экспрессию, например, p = 10e-6, а мы секвенируем прочтения по штучке за раз
• Сколько прочтений из этого гена я получу пока не отсеквенирую r = 1e7 прочтений не из этого гена?

Для этого воспользуемся формулой NB(r, 1 — p), которое будет показывать число удач до r-ой неудачи

9

Отрицательное биномиальное распределение

• Допустим ген имеет не очень высокую экспрессию, например, p = 10e-6, а мы секвенируем прочтения по штучке за раз
• Сколько прочтений из этого гена я получу пока не отсеквенирую r = 1e6 прочтений не из этого гена?

Для этого воспользуемся формулой NB(r, 1 — p), которое будет показывать число удач до r-ой неудачи

10

Отрицательное биномиальное распределение

• Допустим ген имеет не очень высокую экспрессию, например, p = 10e-6, а мы секвенируем прочтения по штучке за раз
• Сколько прочтений из этого гена я получу пока не отсеквенирую r = 1e5 прочтений не из этого гена?

Для этого воспользуемся формулой NB(r, 1 — p), которое будет показывать число удач до r-ой неудачи

11

Среднее и дисперсия NB-распределения

Среднее и дисперсия отрицательного биномиального распределения связаны, благодаря чему мы можем инспектировать наши распределения даже без каких-либо тестов на Goodness of Fit

Это свойство называют овердисперсией

12

Среднее и дисперсия NB-распределения

Среднее и дисперсия отрицательного биномиального распределения связаны, благодаря чему мы можем инспектировать наши распределения даже без каких-либо тестов на Goodness of Fit

Это свойство называют овердисперсией

13

Как понять распределение наших данных?

1. Допустим, мы считаем, что наши значения описываются некоторым распределением X(a, b)
2. При помощи MLE мы можем оценить наиболее правдоподобные значения параметров этого распределения a и b
3. После этого мы можем посчитать правдоподобие того, что наши данные порождены данной моделью
4. В итоге, используя информацию о правдоподобии данных в контексте данного распределения и числе параметров распределения, мы можем сравнить Goodness of Fitness наших данных различными распределениями

14

Нормализации

Количество каунтов гена, которые мы видим, зависит от нескольких параметров:

• от длины гена,
• от глубины библиотеки,
• от экспрессии гена,
• от дополнительных факторов, которые сложно оценить.

Для того, чтобы убрать влияние глубины секвенирования и длины (а в особенности чтобы суммировать информацию по экспрессии транскриптов в экспрессию гена, отнормировав на длину каждого из транскриптов), придумали ряд метрик

15

RPKM и TPM

16

В чём разница?

Связь TPM и RPKM

17

Распределение CPM / TPM

18

CPM

TPM

Проблемы TPM и RPKM

Нормализация на глубину библиотеки предполагает, что суммарное “истинное” количество РНК в клетке константно

Это не работает в случае, когда, например, экспрессия одного набора генов увеличилась, а других — не поменялась

19

Evans et al., Brief Bioinform, 2017

Корректная нормализация

При корректной нормализации (которую, например, выполняет DESeq2 или edgeR) мы принимаем во внимание, что большая часть генов не меняет свою экспрессию между образцами

20

Нормализация в DESeq2 (RLE)

21

Нормализация в DESeq2 (RLE)

22

Нормализация в DESeq2 (RLE)

23

Нормализация в DESeq2 (RLE)

24

Итого по нормализациям

• CPM — простое сравнение одинаковых генов каунтов между разными образцами, грубая нормировка только на глубину библиотеки. Не для DE
• RPKM — сравнение генов внутри одного образца (например, для ранговых методов, о которых поговорим дальше). Не для DE
• TMP — сравнение генов как внутри одного образца (для ранговых методов), так и грубого — между образцами (но не для DE!)
• RLE и TMM — сравнение генов между разными образцами (в том числе и для DE), но не внутри одного образца (отсутствует нормировка на длину)

25

Дорожная карта анализа RNA-Seq

26

QC прочтений

Выравнивание

Подсчёт�экспрессий

Псевдовыравнивание

Нормализация

Дифференциальная экспрессия

GO, GSEA, ssGSEA и прочее

Суть задачи

Нам необходимо статистически сравнить среднее экспрессий между двумя выборками образцов

Что бы мы сделали в классическом случае?

1. Тест Манна-Уитни,
2. t-test

Проблема в том, что тест Манна-Уитни будет слишком слабый, так как чаще всего у нас мало точек в каждой из выборок, а t-test просто не подойдёт потому, что наши данные распределены не нормально

Что делать?

27

Причём тут регрессия?

С одной стороны, регрессионные модели могут позволить нам оценить статистическую достоверность разниц в средних

С другой стороны, GLM позволяют обобщить регрессию на ненормальные распределения

28

Причём тут регрессия?

Статистический вопрос, который мы будем извлекать из регрессии, — значимо ли различаются параметры β1 и β2?

Это можно сказать, сравнив правдоподобия моделей или при помощи других подходов (будет оговорено дальше)

29

Причём тут регрессия?

Линейную модель можно обобщить и добавить более двух уровней фактора, чтобы сравнивать сразу несколько категорий

30

Intercept

Вместо того, чтобы сравнивать значимость разницы между β1 и β2, обычно используют модель со свободным членом β0 и после этого вычисляют значимость β1

Свободный член в данном случае называют словом intercept

31

Intercept

Эту же логику можно обобщить и на модели с несколькими категориями в таргетной переменной

32

Линейные модели

y ~ 0 + feature1 + feature2 + …�без intercept

y ~ 1 + feature1 + feature 2 + …�с intercept

33

Какие переменные включают в модель?

Таргет:

• экспериментальные условия,

сопутствующие факторы:

• пациент,
• пол,
• возраст,
• … (всё, что может иметь влияние на экспрессию)

Что не включают:

• техническую повторность

34

Обобщённые линейные модели (GLM)

В обобщённой линейной модели нет требования к нормальности и гомоскедамтичности остатков

Коэффициенты определяются при помощи MLE

35

Модель DESeq2

Модель, которая вшита в DESeq2, может описываться следующим образом:

36

Последовательность действий DESeq2

1. Сначала происходит оценка size factor’а,
2. потом происходит оценка дисперсии и затем
3. происходит оценка параметров β модели при помощи GLM

37

Подрезание дисперсии

При малых размерах выборки оценка дисперсии становится достаточно неточной, поэтому используют процедуру подрезание дисперсии

38

Взаимодействие переменных

Удобным способом понимания и отображения того, что с чем сравнивается в дизайне экспериментов по секвенированию РНК могут служить модельные матрицы

Модельные матрицы содержат 0 или 1 для каждого из элементов линейной модели� model.matrix(~1+condition+time+condition:time, samples)

Рассмотрим примеры модельных матриц для разных дизайнов (по материалам Hugo Tavares)

39

Один фактор, два уровня

40

Один фактор, два уровня

41

Один фактор, два уровня

42

Один фактор, два уровня

43

Один фактор, три уровня

44

Два фактора и взаимодействие

45

Три фактора с вложенностью

46

Три фактора с вложенностью

47

P-value

48

Способы определения достоверности коэффициентов линейной модели

Likelihood-Ratio Test (LRT)

​

Рассматривает отношение правдоподобий H₀ и Hₐ, логарифм их отношения распределён как χ²

Тест Вальда

​

Похож на LRT, но в явном виде сравнивает не правдоподобия моделей, а коэффициенты

p-value = NA?

Если в строке все значения = 0, что изменение экспрессии и дисперсию не посчитать

Если в строке есть очень большой выброс, то p-value назначается NA

Строка не прошла фильтрацию по средней экспрессии

49

Проблема множественного сравнения

50

Принципы принятия решений

Некоторые обобщения ошибки первого рода:

• FWER — family-wise error rate, групповая вероятность ошибки первого рода. Используется при поправке методом Бонферрони
• FDR — false discovery rate, средняя доля ложных отклонений гипотез (среди всех отклонений). Используется при поправке методом Бенджамини — Хохберга

51

Поправка Бонферрони

52

Поправка Бенджамини-Хохберга

53

Volcano plot

54

От генов к транскриптам: tximport

Как мы уже говорили ранее, самой правильной стратегией будет проводить анализ дифференциальной экспрессии на уровне транскриптов, а потом уже агрегировать информацию до уровня генов

55