Язык R

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

Лекция 5

Анна Валяева

4 октября 2024

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

Данные из пакета palmerpenguins.

penguins
# A tibble: 344 × 8
   species island    bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
   <chr>   <chr>              <dbl>         <dbl>             <dbl>       <dbl>
 1 Adelie  Torgersen           39.1          18.7               181        3750
 2 Adelie  Torgersen           39.5          17.4               186        3800
 3 Adelie  Torgersen           40.3          18                 195        3250
 4 Adelie  Torgersen           NA            NA                  NA          NA
 5 Adelie  Torgersen           36.7          19.3               193        3450
 6 Adelie  Torgersen           39.3          20.6               190        3650
 7 Adelie  Torgersen           38.9          17.8               181        3625
 8 Adelie  Torgersen           39.2          19.6               195        4675
 9 Adelie  Torgersen           34.1          18.1               193        3475
10 Adelie  Torgersen           42            20.2               190        4250
# ℹ 334 more rows
# ℹ 2 more variables: sex <chr>, year <dbl>

Столбчатая диаграмма

Столбчатая диаграмма - geom_bar

Сколько было пингвинов каждого вида?

ggplot(penguins, aes(x = species)) +
  geom_bar()

Цвета по столбцам

ggplot(penguins, 
  aes(
    x = species, 
    fill = species)) +
  geom_bar()

Убрать легенду

ggplot(penguins, 
  aes(
    x = species, 
    fill = species)) +
  geom_bar() +
  theme(legend.position = "none")

Что делает geom_bar()?

penguins %>% 
  group_by(species) %>% 
  summarise(n = n())
# A tibble: 3 × 2
  species       n
  <chr>     <int>
1 Adelie      152
2 Chinstrap    68
3 Gentoo      124

summarise()

penguins %>% 
  summarise(
    n = n(),
    avg_flipper = mean(flipper_length_mm, na.rm = TRUE)
  )
# A tibble: 1 × 2
      n avg_flipper
  <int>       <dbl>
1   344        201.

group_by() + summarise()

penguins %>% 
  group_by(species) %>% 
  summarise(
    n = n(),
    avg_flipper = mean(flipper_length_mm, na.rm = TRUE)
  )
# A tibble: 3 × 3
  species       n avg_flipper
  <chr>     <int>       <dbl>
1 Adelie      152        190.
2 Chinstrap    68        196.
3 Gentoo      124        217.

group_by() + mutate()

penguins %>% 
  group_by(species) %>% 
  mutate(
    n = n(),
    avg_flipper = mean(flipper_length_mm, na.rm = TRUE),
    .before = 1)
# A tibble: 344 × 10
# Groups:   species [3]
       n avg_flipper species island    bill_length_mm bill_depth_mm
   <int>       <dbl> <chr>   <chr>              <dbl>         <dbl>
 1   152        190. Adelie  Torgersen           39.1          18.7
 2   152        190. Adelie  Torgersen           39.5          17.4
 3   152        190. Adelie  Torgersen           40.3          18  
 4   152        190. Adelie  Torgersen           NA            NA  
 5   152        190. Adelie  Torgersen           36.7          19.3
 6   152        190. Adelie  Torgersen           39.3          20.6
 7   152        190. Adelie  Torgersen           38.9          17.8
 8   152        190. Adelie  Torgersen           39.2          19.6
 9   152        190. Adelie  Torgersen           34.1          18.1
10   152        190. Adelie  Torgersen           42            20.2
# ℹ 334 more rows
# ℹ 4 more variables: flipper_length_mm <dbl>, body_mass_g <dbl>, sex <chr>,
#   year <dbl>

summarise(.by = ...) 🏠

Аналогично mutate(.by = ...).

penguins %>% 
  summarise(
    n = n(),
    avg_flipper = mean(flipper_length_mm, na.rm = TRUE),
    .by = species)
# A tibble: 3 × 3
  species       n avg_flipper
  <chr>     <int>       <dbl>
1 Adelie      152        190.
2 Gentoo      124        217.
3 Chinstrap    68        196.

count()

Краткая замена сочетанию функций group_by() и summarise().

penguins %>% 
  count(species)
# A tibble: 3 × 2
  species       n
  <chr>     <int>
1 Adelie      152
2 Chinstrap    68
3 Gentoo      124
penguins %>% 
  group_by(species) %>% 
  summarise(n = n())
# A tibble: 3 × 2
  species       n
  <chr>     <int>
1 Adelie      152
2 Chinstrap    68
3 Gentoo      124

Визуализировать значения из таблицы

Если у нас есть таблица, в которой уже посчитаны пингвины разных видов:

penguin_species <- penguins %>% 
  count(species)

penguin_species
# A tibble: 3 × 2
  species       n
  <chr>     <int>
1 Adelie      152
2 Chinstrap    68
3 Gentoo      124
ggplot(penguin_species, 
  aes(
    x = species, 
    y = n,
    fill = species)) +
  geom_bar(stat = "identity") +
  theme(legend.position = "none")

Визуализировать значения из таблицы

Можно использовать geom_bar(stat = "identity") или geom_col():

ggplot(penguin_species, 
  aes(
    x = species, 
    y = n,
    fill = species)) +
  geom_col() +
  theme(legend.position = "none")

Высота столбца как среднее по группе

Хотим визуализировать среднюю длину крыла у пингвинов разных видов.

ggplot(penguins, 
  aes(
    x = species, 
    y = flipper_length_mm,
    fill = species)) +
  geom_bar(
    stat = "summary", 
    fun = "mean") +
  theme(legend.position = "none")

Высота столбца как среднее по группе с разбросом

Хотим визуализировать среднюю длину крыла у пингвинов разных видов и показать разброс в виде “среднее ± стандартное отклонение”.

Сначала нужно подготовить таблицу со значениями средних и SD по группам:

penguin_flippers <- penguins %>% 
  group_by(species) %>% 
  summarise(
    avg_flipper_length = mean(flipper_length_mm, na.rm = TRUE),
    sd_flipper_length = sd(flipper_length_mm, na.rm = TRUE))

penguin_flippers
# A tibble: 3 × 3
  species   avg_flipper_length sd_flipper_length
  <chr>                  <dbl>             <dbl>
1 Adelie                  190.              6.54
2 Chinstrap               196.              7.13
3 Gentoo                  217.              6.48

geom_bar() + geom_errorbar()

ggplot(penguin_flippers, 
  aes(
    x = species, 
    y = avg_flipper_length,
    ymin = avg_flipper_length - sd_flipper_length,
    ymax = avg_flipper_length + sd_flipper_length,
    fill = species)) +
  geom_bar(stat = "identity") +
  geom_errorbar(width = 0.5) +
  theme(legend.position = "none")

stat_summary()

ggplot(penguins, 
       aes(x = species, y = flipper_length_mm, fill = species)) +
  geom_bar(stat = "summary", fun = "mean") +
  stat_summary(fun = mean, 
               fun.max = function(x) mean(x) + sd(x),
               fun.min = function(x) mean(x) - sd(x),
               geom = "errorbar", width = 0.5)

ggbarplot() из пакета ggpubr

Пакет ggpubr не из коллекции tidyverse, его нужно установить:

install.packages("ggpubr")


library(ggpubr)

ggbarplot(
  penguins, 
  x = "species", 
  y = "flipper_length_mm", 
  fill = "species", 
  add = "mean_sd")

Какой разброс показывать?

  • SD - стандартное отклонение
  • SE, или SEM - стандартная ошибка среднего - \(SE = \frac{SD}{\sqrt{n}}\)
  • CI - доверительный интервал
  • IQR - интерквартильный размах
  • …
?ggpubr::ggbarplot

Альтернативы столбчатой диаграммы для количественного признака

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

  • малые выборки (\(n<20\)): на столбчатую диаграмму добавляйте точки - индивидуальные наблюдения

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

Разделение по двум категориальным переменным

По умолчанию рисуется стековая столбчатая диаграмма - используется параметр position = "stack".

ggplot(penguins, 
  aes(
    x = species, 
    fill = sex)) +
  geom_bar()

Группированная диаграмма

Чтобы получить группированную столбчатую диаграмма, нужно изменить параметр на position = "dodge".

ggplot(penguins, 
  aes(
    x = species, 
    fill = sex)) +
  geom_bar(
    position = "dodge")

Столбчатая диаграмма с процентами

Чтобы вместо числа наблюдений визуализировать процент, нужно изменить параметр position = "fill".

ggplot(penguins, 
  aes(
    x = species, 
    fill = sex)) +
  geom_bar(
    position = "fill")

Изменить порядок столбцов

По умолчанию столбцы располагаются в алфавитном порядке.

ggplot(penguins, 
  aes(
    x = species, 
    fill = species)) +
  geom_bar() +
  theme(legend.position = "none")

Факторы

Категориальные переменные как факторы 🏠

  • Фактор - fct - особый тип данных в R
  • Уровни фактора [levels] - ограниченное число известных значений категориальной переменной, отсортированные в определенном порядке
set.seed(123)

ten_penguins <- penguins %>% 
  slice_sample(n = 10) # случайные 10 строк из датафрейма

ten_penguins$species
 [1] "Gentoo"    "Adelie"    "Gentoo"    "Chinstrap" "Adelie"    "Chinstrap"
 [7] "Gentoo"    "Gentoo"    "Chinstrap" "Gentoo"

factor() 🏠

Преобразовать какой-то вектор в вектор с факторами:

factor(ten_penguins$species)
 [1] Gentoo    Adelie    Gentoo    Chinstrap Adelie    Chinstrap Gentoo   
 [8] Gentoo    Chinstrap Gentoo   
Levels: Adelie Chinstrap Gentoo


Можем сами задать уровни фактора:

factor(ten_penguins$species, levels = c("Chinstrap", "Adelie", "Gentoo"))
 [1] Gentoo    Adelie    Gentoo    Chinstrap Adelie    Chinstrap Gentoo   
 [8] Gentoo    Chinstrap Gentoo   
Levels: Chinstrap Adelie Gentoo

Задать уровни фактора при отрисовке графика 🏠

ggplot(penguins, 
  aes(
    x = factor(species, levels = c("Chinstrap", "Adelie", "Gentoo")), 
    y = flipper_length_mm,
    fill = species)) +
  geom_boxplot() +
  theme(legend.position = "none")

Создать столбец с факторами 🏠

Обычно оптимальнее заранее создать столбец с факторами или преобразовать столбец с категориями в факторы:

penguins_fct <- penguins %>% 
  mutate(
    species = factor(species, levels = c("Chinstrap", "Adelie", "Gentoo")))

penguins_fct
# A tibble: 344 × 8
   species island    bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
   <fct>   <chr>              <dbl>         <dbl>             <dbl>       <dbl>
 1 Adelie  Torgersen           39.1          18.7               181        3750
 2 Adelie  Torgersen           39.5          17.4               186        3800
 3 Adelie  Torgersen           40.3          18                 195        3250
 4 Adelie  Torgersen           NA            NA                  NA          NA
 5 Adelie  Torgersen           36.7          19.3               193        3450
 6 Adelie  Torgersen           39.3          20.6               190        3650
 7 Adelie  Torgersen           38.9          17.8               181        3625
 8 Adelie  Torgersen           39.2          19.6               195        4675
 9 Adelie  Torgersen           34.1          18.1               193        3475
10 Adelie  Torgersen           42            20.2               190        4250
# ℹ 334 more rows
# ℹ 2 more variables: sex <chr>, year <dbl>

График по датафрейму с факторами 🏠

ggplot(penguins_fct, 
  aes(
    x = species, 
    y = flipper_length_mm,
    fill = species)) +
  geom_boxplot() +
  theme(legend.position = "none")

Пакет forcats 🏠

  • Для работы с факторами есть пакет forcats в составе tidyverse
library(tidyverse)
# или
library(forcats)

fcts <- factor(ten_penguins$species)
fcts
 [1] Gentoo    Adelie    Gentoo    Chinstrap Adelie    Chinstrap Gentoo   
 [8] Gentoo    Chinstrap Gentoo   
Levels: Adelie Chinstrap Gentoo
# В обратном порядке
fct_rev(fcts) %>% levels()
[1] "Gentoo"    "Chinstrap" "Adelie"   
# По частоте встречаемости
fct_infreq(fcts) %>% levels()
[1] "Gentoo"    "Chinstrap" "Adelie"   
# Случайно перемешать
fct_shuffle(fcts) %>% levels()
[1] "Chinstrap" "Adelie"    "Gentoo"

base R vs forcats 🏠

  • as.factor() сортирует уровни фактора по алфавиту
as.factor(ten_penguins$species)
 [1] Gentoo    Adelie    Gentoo    Chinstrap Adelie    Chinstrap Gentoo   
 [8] Gentoo    Chinstrap Gentoo   
Levels: Adelie Chinstrap Gentoo
  • forcats::as_factor() сортирует уровни фактора по порядку их встречи
as_factor(ten_penguins$species)
 [1] Gentoo    Adelie    Gentoo    Chinstrap Adelie    Chinstrap Gentoo   
 [8] Gentoo    Chinstrap Gentoo   
Levels: Gentoo Adelie Chinstrap

Что такое фактор? 🏠

Строки:

typeof(ten_penguins$species)
[1] "character"
class(ten_penguins$species)
[1] "character"
as.integer(ten_penguins$species)
 [1] NA NA NA NA NA NA NA NA NA NA
sort(ten_penguins$species[c(6:10,1:5)])
 [1] "Adelie"    "Adelie"    "Chinstrap" "Chinstrap" "Chinstrap" "Gentoo"   
 [7] "Gentoo"    "Gentoo"    "Gentoo"    "Gentoo"

Факторы:

typeof(as_factor(ten_penguins$species))
[1] "integer"
class(as_factor(ten_penguins$species))
[1] "factor"
as.integer(as_factor(ten_penguins$species))
 [1] 1 2 1 3 2 3 1 1 3 1
sort(as_factor(ten_penguins$species)[c(6:10,1:5)])
 [1] Gentoo    Gentoo    Gentoo    Gentoo    Gentoo    Adelie    Adelie   
 [8] Chinstrap Chinstrap Chinstrap
Levels: Gentoo Adelie Chinstrap

Что такое уровни фактора? 🏠

levels(ten_penguins$species)
NULL
levels(as_factor(ten_penguins$species))
[1] "Gentoo"    "Adelie"    "Chinstrap"
as_factor(ten_penguins$species)[1:3]
[1] Gentoo Adelie Gentoo
Levels: Gentoo Adelie Chinstrap

Создание уровней фактора 🏠

some_penguins <- c(ten_penguins$species, "Adelii")

penguin_species_list <- c("Chinstrap", "Adelie", "Gentoo", "Emperor")

some_penguins <- factor(some_penguins, levels = penguin_species_list)
some_penguins
 [1] Gentoo    Adelie    Gentoo    Chinstrap Adelie    Chinstrap Gentoo   
 [8] Gentoo    Chinstrap Gentoo    <NA>     
Levels: Chinstrap Adelie Gentoo Emperor

Убрать лишнее 🏠

fct_drop()

Удалить неиспользуемые уровни фактора: Emperor.

fct_drop(some_penguins)
 [1] Gentoo    Adelie    Gentoo    Chinstrap Adelie    Chinstrap Gentoo   
 [8] Gentoo    Chinstrap Gentoo    <NA>     
Levels: Chinstrap Adelie Gentoo

fct_na_value_to_level()

Приписать отсутствующим уровням (NA) явное название.

fct_na_value_to_level(some_penguins) # level задает название нового уровня
 [1] Gentoo    Adelie    Gentoo    Chinstrap Adelie    Chinstrap Gentoo   
 [8] Gentoo    Chinstrap Gentoo    <NA>     
Levels: Chinstrap Adelie Gentoo Emperor <NA>

Упорядочить уровни фактора 🏠

  • fct_inorder()- упорядочить уровни фактора в порядке встречаемости в векторе
  • fct_infreq() - упорядочить уровни фактора по частоте встречаемости.
  • fct_rev() - изменить порядок уровней на обратный.
  • fct_shuffle() - перемешать уровни

Отсортировать уровни фактора в зависимости от других значений 🏠

  • Функции fct_reorder() и fct_reorder2()

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

penguins %>% 
  mutate(species = fct_reorder(
    species, flipper_length_mm, .fun = median, na.rm = TRUE, .desc = TRUE)) %>%
  ggplot(aes(x = species, y = flipper_length_mm, fill = species)) +
  geom_boxplot() +
  theme(legend.position = "none")

fct_reorder() 🏠

penguins %>% 
  group_by(species) %>% 
  summarise(mean_fl_length = mean(flipper_length_mm, na.rm = TRUE)) %>% 
  arrange(desc(mean_fl_length))
# A tibble: 3 × 2
  species   mean_fl_length
  <chr>              <dbl>
1 Gentoo              217.
2 Chinstrap           196.
3 Adelie              190.
penguins$species %>% 
  fct_reorder(penguins$flipper_length_mm, .fun = mean, na.rm = TRUE, .desc = TRUE) %>%
  levels()
[1] "Gentoo"    "Chinstrap" "Adelie"

Объединить редкие категории в одну 🏠

  • Функции семейства fct_lump():
    • fct_lump_n() - оставляет n самых многочисленных уровней, остальные объединяет
    • fct_lump_min() - оставляет уровни, встречающиеся не реже указанного числа раз
    • fct_lump_prop() - оставляет уровни, встречающиеся не чаще указанной частоты
    • fct_lump_lowfreq() - максимально наполняет группу Other так, чтобы она все равно оставалась самой малопредставленной

forcats - что почитать 🏠

fct_infreq()

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

penguins %>% 
  mutate(species = fct_infreq(species)) %>% 
  ggplot(aes(x = species, fill = species)) +
  geom_bar() +
  theme(legend.position = "none")

fct_reorder()

Хотим изобразить среднюю длину крыла и расположить столбцы по убыванию их высоты.

penguins %>% 
  mutate(species = fct_reorder(
    species, flipper_length_mm, .fun = mean, na.rm = TRUE, .desc = TRUE)) %>%
  ggplot(aes(x = species, y = flipper_length_mm, fill = species)) +
  geom_bar(stat = "summary", fun = "mean") +
  theme(legend.position = "none")

fct_reorder() 🏠

penguins %>% 
  group_by(species) %>% 
  summarise(mean_fl_length = mean(flipper_length_mm, na.rm = TRUE)) %>% 
  arrange(desc(mean_fl_length))
# A tibble: 3 × 2
  species   mean_fl_length
  <chr>              <dbl>
1 Gentoo              217.
2 Chinstrap           196.
3 Adelie              190.
penguins$species %>% 
  fct_reorder(penguins$flipper_length_mm, .fun = mean, na.rm = TRUE, .desc = TRUE) %>%
  levels()
[1] "Gentoo"    "Chinstrap" "Adelie"

Круговая диаграмма

ggplot(penguin_species, 
       aes(
         x = "",
         y = n,
         fill = species)) +
  geom_col() +
  coord_polar("y", start=0) +
  theme_void()

Трансформация датафреймов

Что подставить в aes?

Хотим на одном графике показать распределение длины и ширины клювов пингвинов.

aes переменные - из разных столбцов

Широкий формат

  • Несколько переменных, похожих по своей природе, в разных столбцах
# A tibble: 344 × 3
   species bill_length_mm bill_depth_mm
   <chr>            <dbl>         <dbl>
 1 Adelie            39.1          18.7
 2 Adelie            39.5          17.4
 3 Adelie            40.3          18  
 4 Adelie            NA            NA  
 5 Adelie            36.7          19.3
 6 Adelie            39.3          20.6
 7 Adelie            38.9          17.8
 8 Adelie            39.2          19.6
 9 Adelie            34.1          18.1
10 Adelie            42            20.2
# ℹ 334 more rows

Длинный формат

  • Столбец с названиями переменных
  • Столбец со значениями переменных
# A tibble: 688 × 3
   species measurement       mm
   <chr>   <chr>          <dbl>
 1 Adelie  bill_length_mm  39.1
 2 Adelie  bill_depth_mm   18.7
 3 Adelie  bill_length_mm  39.5
 4 Adelie  bill_depth_mm   17.4
 5 Adelie  bill_length_mm  40.3
 6 Adelie  bill_depth_mm   18  
 7 Adelie  bill_length_mm  NA  
 8 Adelie  bill_depth_mm   NA  
 9 Adelie  bill_length_mm  36.7
10 Adelie  bill_depth_mm   19.3
# ℹ 678 more rows

Длинный и широкий формат

  • Из широкого в длинный - pivot_longer()
  • Из длинного в широкий - pivot_wider()

pivot_longer()

  • Из широкого в длинный формат
penguins %>% 
  # Превращаем в длинный формат
  pivot_longer(
    cols = c(bill_length_mm, bill_depth_mm), 
    names_to = "measurement", values_to = "mm") %>% 
  ggplot() +
  geom_boxplot(aes(x = species, y = mm, fill = measurement))

pivot_wider()

  • Из длинного в широкий формат
penguins %>% 
  select(species, bill_length_mm, bill_depth_mm) %>% 
  mutate(id = 1:n()) %>% 
  # Превращаем в длинный формат
  pivot_longer(
    cols = c(bill_length_mm, bill_depth_mm), 
    names_to = "measurement", values_to = "mm") %>%
  # Обратно в широкий
  pivot_wider(names_from = measurement, values_from = mm)
# A tibble: 344 × 4
   species    id bill_length_mm bill_depth_mm
   <chr>   <int>          <dbl>         <dbl>
 1 Adelie      1           39.1          18.7
 2 Adelie      2           39.5          17.4
 3 Adelie      3           40.3          18  
 4 Adelie      4           NA            NA  
 5 Adelie      5           36.7          19.3
 6 Adelie      6           39.3          20.6
 7 Adelie      7           38.9          17.8
 8 Adelie      8           39.2          19.6
 9 Adelie      9           34.1          18.1
10 Adelie     10           42            20.2
# ℹ 334 more rows

Порисуем!

Конец!