Язык R

и его применение в биоинформатике

Лекция 9

Анастасия Жарикова

8 ноября 2024

Случайная величина

Величина, принимающая в ходе эксперимента то или иное значение

Численные исходы некого случайного эксперимента

Точно сказать, какое это будет значение, мы не можем

Случайные величины

  • Число, выпавшее на шестигранном кубике

  • Пол ребенка

  • Время прибытия автобуса

  • Время работы лампочки

  • Количество клиентов в день

  • Число бракованных товаров в партии

Бывают дискретными и непрырывными

Дискретные случайные величины

Принимает не более, чем счетное множество значений

Распределение можно задавать при помощи таблицы: каждому значению сопоставлена его вероятность.

В сумме вероятности равны единице

Задача

Есть 2 случайные величины

Насколько согласовано их изменение?

Ковариация

Мера совместной изменчивости двух случайных величин

Большие значения одной переменной в основном соответствуют большим значениям в другой переменной, аналогично для меньших значений

Ковариация

library(tidyverse)
x = rnorm(10)
y = rnorm(10)
ggplot(mapping = aes(x = x, y = y)) + geom_point(size = 3) + theme_bw() + theme(text = element_text(size = 20))
cov(x, y)
[1] 0.4784603

Ковариация

x = 1:10
y = 4:13
ggplot(mapping = aes(x = x, y = y)) + geom_point(size = 3) + theme_bw() + theme(text = element_text(size = 20))
cov(x, y)
[1] 9.166667

Ковариация

x = 1:10
y = 13:4
ggplot(mapping = aes(x = x, y = y)) + geom_point(size = 3) + theme_bw() + theme(text = element_text(size = 20))
cov(x, y)
[1] -9.166667

Ковариация

Не нормированная величина

Зависит от величин переменных

Сложно интерпретировать

Ковариация

x = 1:10
y = x*0.1
ggplot(mapping = aes(x = x, y = y)) + geom_point(size = 3) + theme_bw() + theme(text = element_text(size = 20))
cov(x, y)
[1] 0.9166667

Ковариация

x = 1:10
y = x*10
ggplot(mapping = aes(x = x, y = y)) + geom_point(size = 3) + theme_bw() + theme(text = element_text(size = 20))
cov(x, y)
[1] 91.66667

Корреляция

Можно отнормировать

Получим корреляцию Пирсона

Безразмерная величина

Изменяется от -1 до 1

Можно применять не только к нормально распределенным данным

Корреляция = 0

Из этого не следует независимость переменных

Корреляция != 0

Нельзя рассуждать о причинно-следственной связи

Корреляция 0,66%

ref

Корреляция != 0

Корреляция 0,99%

Визуализируй!

set.seed(123)
x = rnorm(1000)
y = rnorm(1000)
ggplot(mapping = aes(x = x, y = y)) + geom_point(size = 3) + theme_bw() + theme(text = element_text(size = 20))
cor(x, y)
[1] 0.08647944

Визуализируй!

set.seed(123)
x = c(rnorm(100), rpois(5, 5))
y = c(rnorm(100), rpois(5, 6))
ggplot(mapping = aes(x = x, y = y)) + geom_point(size = 3) + theme_bw() + theme(text = element_text(size = 20))
cor(x, y)
[1] 0.5275666

Визуализируй!

set.seed(123)
x = c(rnorm(100), rnorm(100, 6))
y = c(rnorm(100), rnorm(100, 6))
ggplot(mapping = aes(x = x, y = y)) + geom_point(size = 3) + theme_bw() + theme(text = element_text(size = 20))
cor(x, y)
[1] 0.8977794

Корреляция

Если зависимость не линейная, но монотонная?

Корреляция Спирмена

Считаем корреляцию между рангами наблюдений

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

mt_filt <- mtcars[, c(1,3,4,5,6,7)]

cor(mt_filt$mpg, mt_filt$disp)
[1] -0.8475514
cor(mt_filt$mpg, mt_filt$disp, method = 'spearman')
[1] -0.9088824

Визуализация корреляционной матрицы

library(ggcorrplot)
correlation_matrix <- round(cor(mt_filt),1)
head(correlation_matrix[, 1:4])
      mpg disp   hp drat
mpg   1.0 -0.8 -0.8  0.7
disp -0.8  1.0  0.8 -0.7
hp   -0.8  0.8  1.0 -0.4
drat  0.7 -0.7 -0.4  1.0
wt   -0.9  0.9  0.7 -0.7
qsec  0.4 -0.4 -0.7  0.1
ggcorrplot(correlation_matrix, method ="square")

Визуализация корреляционной матрицы

ggcorrplot(correlation_matrix, method ="circle")

Визуализация корреляционной матрицы

ggcorrplot(correlation_matrix, lab =TRUE)

Визуализация корреляционной матрицы

ggcorrplot(correlation_matrix, type ="lower", lab =TRUE)

Визуализация корреляционной матрицы

ggcorrplot(correlation_matrix, type ="lower", lab =TRUE, hc.order = T)

Визуализация корреляционной матрицы

library(corrplot)
corr_mat <- cor(mt_filt)
corrplot(corr_mat)

Визуализация корреляционной матрицы

corrplot(corr_mat, col=colorRampPalette(c("blue","white","red"))(200))

Визуализация корреляционной матрицы

corrplot(corr_mat, method = 'number')

Визуализация корреляционной матрицы

corrplot(corr_mat, method = 'ellipse')

Задание 5

Пусть виды ирисов будут “до” обработки супер-средством и “после”.

Влияет ли супер-средство на длину лепестков ирисов?

Решение 5.1

df <- iris %>%
  filter(Species %in% c('versicolor', 'virginica'))
t.test(Petal.Length ~ Species, df, paired = T)

    Paired t-test

data:  Petal.Length by Species
t = -12.091, df = 49, p-value = 2.562e-16
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
 -1.506744 -1.077256
sample estimates:
mean difference 
         -1.292

Решение 5.2

Визуализируем

library(ggpubr)
ggpaired(df, x = "Species", y = "Petal.Length", 
          fill = ("Species"), 
          add = c('mean', 'jitter'),
          palette = c("#9e2a2b", "#62929e"),
         line.color = "gray", line.size = 0.4) +
  stat_compare_means(method = "t.test")

t.test(Petal.Length ~ Species, df, paired = T)$p.value -> a

ggpaired(df, x = "Species", y = "Petal.Length", 
          fill = ("Species"), 
          add = c('mean', 'jitter'),
          palette = c("#9e2a2b", "#62929e"),
         line.color = "gray", line.size = 0.4) + 
annotate("text", x=1.3, y=7.5, 
         label = paste('p.value (t.test): ', a, sep = ''), size = 4)

immer

Урожайность ячменя в 1931(Y1) и 1932(Y2) годах.

Задание 6

Различается ли урожайность?

Смотрим

library(MASS)
head(immer)
  Loc Var    Y1    Y2
1  UF   M  81.0  80.7
2  UF   S 105.4  82.3
3  UF   V 119.7  80.4
4  UF   T 109.7  87.2
5  UF   P  98.3  84.2
6   W   M 146.6 100.4

Поправим данные

im <- immer %>%
  pivot_longer(Y1:Y2, names_to = 'year', values_to = 'value')

Рисуем

im %>% 
  ggplot() + geom_boxplot(aes(value, fill = year)) + theme_bw()

Еще рисуем

im %>%
  ggplot(aes(x = year, y = value)) + 
  geom_bar(position='dodge', stat='summary', fun='mean', 
                      fill = 'lightblue', color = 'brown') + 
  geom_jitter(width = 0.25) + theme_bw()

Решение 6.1

wilcox.test(value ~ year, im, paired = F)

    Wilcoxon rank sum test with continuity correction

data:  value by year
W = 589, p-value = 0.04059
alternative hypothesis: true location shift is not equal to 0

Задание 7

Есть ли зависимость между длиной лепестка и длиной чашелистика у ирисов?

Рисуем

Рисуем

ggplot(data = iris) + 
  geom_point(aes(x = Sepal.Length, y = Petal.Length)) + theme_bw()
cor(iris$Sepal.Length, iris$Petal.Length)
[1] 0.8717538

Пирсона или Спирмена?

cor(iris$Sepal.Length, iris$Petal.Length)
[1] 0.8717538
cor(iris$Sepal.Length, iris$Petal.Length, method = 'spearman')
[1] 0.8818981

А значимость?

cor.test(iris$Sepal.Length, iris$Petal.Length)

    Pearson's product-moment correlation

data:  iris$Sepal.Length and iris$Petal.Length
t = 21.646, df = 148, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.8270363 0.9055080
sample estimates:
      cor 
0.8717538 
cor.test(iris$Sepal.Length, iris$Petal.Length, method = 'spearman')

    Spearman's rank correlation rho

data:  iris$Sepal.Length and iris$Petal.Length
S = 66429, p-value < 2.2e-16
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.8818981

Матрица корреляций

cor(mtcars)
            mpg        cyl       disp         hp        drat         wt
mpg   1.0000000 -0.8521620 -0.8475514 -0.7761684  0.68117191 -0.8676594
cyl  -0.8521620  1.0000000  0.9020329  0.8324475 -0.69993811  0.7824958
disp -0.8475514  0.9020329  1.0000000  0.7909486 -0.71021393  0.8879799
hp   -0.7761684  0.8324475  0.7909486  1.0000000 -0.44875912  0.6587479
drat  0.6811719 -0.6999381 -0.7102139 -0.4487591  1.00000000 -0.7124406
wt   -0.8676594  0.7824958  0.8879799  0.6587479 -0.71244065  1.0000000
qsec  0.4186840 -0.5912421 -0.4336979 -0.7082234  0.09120476 -0.1747159
vs    0.6640389 -0.8108118 -0.7104159 -0.7230967  0.44027846 -0.5549157
am    0.5998324 -0.5226070 -0.5912270 -0.2432043  0.71271113 -0.6924953
gear  0.4802848 -0.4926866 -0.5555692 -0.1257043  0.69961013 -0.5832870
carb -0.5509251  0.5269883  0.3949769  0.7498125 -0.09078980  0.4276059
            qsec         vs          am       gear        carb
mpg   0.41868403  0.6640389  0.59983243  0.4802848 -0.55092507
cyl  -0.59124207 -0.8108118 -0.52260705 -0.4926866  0.52698829
disp -0.43369788 -0.7104159 -0.59122704 -0.5555692  0.39497686
hp   -0.70822339 -0.7230967 -0.24320426 -0.1257043  0.74981247
drat  0.09120476  0.4402785  0.71271113  0.6996101 -0.09078980
wt   -0.17471588 -0.5549157 -0.69249526 -0.5832870  0.42760594
qsec  1.00000000  0.7445354 -0.22986086 -0.2126822 -0.65624923
vs    0.74453544  1.0000000  0.16834512  0.2060233 -0.56960714
am   -0.22986086  0.1683451  1.00000000  0.7940588  0.05753435
gear -0.21268223  0.2060233  0.79405876  1.0000000  0.27407284
carb -0.65624923 -0.5696071  0.05753435  0.2740728  1.00000000

Рисуем

corrplot(cor(mtcars), col=colorRampPalette(c("#1b4965","white","#c1121f"))(200))

libraries

  • library(tidyverse)

  • library(corrplot)

  • library(psych)

  • library(ggcats)

  • library(ggdogs)

  • library(CatterPlots)

  • library(ggimage)

  • library(rsvg)

ggcats

install.packages("remotes")
remotes::install_github("R-CoderDotCom/ggcats@main")
library(ggcats)
grid <- expand.grid(1:5, 3:1)

ggcats

df <- data.frame(x = grid[, 1],
                 y = grid[, 2],
                 image = c("nyancat", "bongo", "colonel", "grumpy", "hipster", "lil_bub",
                           "maru", "mouth", "pop", "pop_close",  "pusheen", "pusheen_pc",
                           "toast", "venus", "shironeko"))

ggplot(df) +
   geom_cat(aes(x, y, cat = image), size = 5) +
   geom_text(aes(x, y - 0.5, label = image), size = 2.5) +
   xlim(c(0.25, 5.5)) + ylim(c(0.25, 3.5)) + theme_bw()

ggcats

iris$cat <- factor(iris$Species,
                   labels = c("pusheen", "toast", "venus"))
ggplot(iris, aes(Petal.Length, Petal.Width)) +
 geom_cat(aes(cat = cat), size = 4) + theme_bw()

ggdogs

install.packages("remotes")
remotes::install_github("R-CoderDotCom/ggdogs@main")
library(ggdogs)

ggdogs

df <- data.frame(x = grid[, 1],
                 y = grid[, 2],
                 image = c("doge", "doge_strong", "chihuahua", "eyes", "gabe", "glasses",
                           "tail", "surprised", "thisisfine", "hearing", "pug", "ears",
                           "husky", "husky_2", "chilaquil"))
                           
ggplot(df) +
 geom_dog(aes(x, y, dog = image), size = 5) +
 geom_text(aes(x, y - 0.5, label = image), size = 2.5) +
 xlim(c(0.25, 5.5)) + ylim(c(0.25, 3.5)) + theme_bw()

ggdogs

set.seed(1)
df <- data.frame(x = 1:10, y = rnorm(10))

ggplot(df, aes(x = x, y = y)) +
  geom_point(size = 3, color = 4) +
  geom_dog(aes(x = 7, y = -0.5), dog = "thisisfine", size = 5) +
  geom_label(aes(x = 7.75, y = -0.1, label = "This is fine")) + theme_bw()

CatterPlot

install.packages("remotes")
remotes::install_github("Gibbsdavidl/CatterPlots")
library(CatterPlots)

CatterPlot

x <- seq(0, 15, 0.5)
y <- sin(x)

obj <- catplot(x, y, cat = 2, catcolor = 2, type = "line")

cats(obj, 4.5, 0.8, cat = 3, catcolor = 3)

Custom dots

set.seed(123)

d <- data.frame(x = rnorm(10), y = rnorm(10), z=1:10,
  image = 'https://www.svgrepo.com/show/69771/dna.svg')

ggplot(d, aes(x, y)) + 
  geom_image(aes(image=image, color=z), size = 0.1) +
  scale_color_gradient(low='blue', high='red') +
  theme_bw()

Немного биологии

У человека ~20000 белок-кодирующих генов

Для каждого гена знаем уровень его экспрессии в группе больных и здоровых

Есть ли различия?

Сделаем t-test для каждого гена

Получили ~2000 дифференциально экспрессирующихся генов

Вопрос

Сколько же генов ассоциированы с заболеванием?

У нас же есть уже p-value, сравненные с уровнем значимости…

Другая ситуация

Есть 100 чашек Петри, в которых растут клетки при одинаковых условиях

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

Получаем 100 реплик

Делим в произвольном порядке 100 реплик на 2 группы

Проверяем каждый из 1000 генов на изменение уровня экспрессии

Ожидаем, что гены не должны менять свою экспрессию

Ожидаем, что \(H_{0}\) верна для всех генов

p-value и ошибка первого рода

Если \(p-value \le \alpha\), отвергаем нулевую гипотезу

Поддерживаем ошибку первого рода на уровне значимости \(\alpha\)

Если вернуться к чашкам, что получается, что мы ожидаем:

1000 * 0.05 = 50 дифференциально экспрессирующихся генов

Мы говорим, что есть статистически значимый результат, но на самом деле его нет

Вернемся к больным

Было 20000 генов

Сделали 20000 статистических тестов

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

Если данные долго анализировать…

… в конце концов они в чем-то сознаются

Scanning Dead Salmon in fMRI Machine Highlights Risk of Red Herrings

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

Проводим много тестов

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

Получаем ложноположительные значения

Хотелось бы как-то контролировать уровень ошибки первого рода сразу для всех проведенных тестов

Поправки на множественное тестирование

Их очень много

Выделяют два основных семейства:

  • FWER - family-wise error rate

  • FDR - false discovery rate

FWER

Вероятность совершить хотя бы одну ошибку первого рода

Контроль FWER на уровне значимости \(\alpha\): \(FWER \le \alpha\)

Достаточно жесткие поправки

FWER

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

  • Поправка Холм-Бонферрони: выше мощность, меньше ложноотрицательных результатов

Неважна независимость

Применять, если важно не получить ни одного ложноположительного результата

Помним, что потеряем много истинноположительных

FDR

  • Поправка Беньямини — Хохберга

Готовы к “мусорным” данным в результате

Порог на FDR зависит от амбиций исследователя ## Дома

Повторяем ANOVA

Конец!