Язык R

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

Лекция 4

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

27 сентября 2024

Функции

Зачем нужны функции?

  • Если вы заметили, что несколько раз используете один и тот же код, то запишите его в функцию.
    • Если вы задали функции говорящее имя (и с помощью комментариев пояснили, как функция работает), ваш код станет более понятным.
    • Если нужно будет что-то изменить в вашей функции, это нужно будет сделать единожды в коде функции.
    • Будет меньше вероятность ошибиться при дублировании кода.
    • Одну и ту же функцию можно использовать в разных проектах.

Пример функции

  • Создали функцию: перевод температуры из градусов Фаренгейта в Цельсия
fahrenheit2celsius <- function(x) {
  (x - 32) * 5 / 9
}

Пример функции

  • Применили функцию
nhtemp
Time Series:
Start = 1912 
End = 1971 
Frequency = 1 
 [1] 49.9 52.3 49.4 51.1 49.4 47.9 49.8 50.9 49.3 51.9 50.8 49.6 49.3 50.6 48.4
[16] 50.7 50.9 50.6 51.5 52.8 51.8 51.1 49.8 50.2 50.4 51.6 51.8 50.9 48.8 51.7
[31] 51.0 50.6 51.7 51.5 52.1 51.3 51.0 54.0 51.4 52.7 53.1 54.6 52.0 52.0 50.9
[46] 52.6 50.2 52.6 51.6 51.9 50.5 50.9 51.7 51.4 51.7 50.8 51.9 51.8 51.9 53.0
fahrenheit2celsius(nhtemp)
Time Series:
Start = 1912 
End = 1971 
Frequency = 1 
 [1]  9.944444 11.277778  9.666667 10.611111  9.666667  8.833333  9.888889
 [8] 10.500000  9.611111 11.055556 10.444444  9.777778  9.611111 10.333333
[15]  9.111111 10.388889 10.500000 10.333333 10.833333 11.555556 11.000000
[22] 10.611111  9.888889 10.111111 10.222222 10.888889 11.000000 10.500000
[29]  9.333333 10.944444 10.555556 10.333333 10.944444 10.833333 11.166667
[36] 10.722222 10.555556 12.222222 10.777778 11.500000 11.722222 12.555556
[43] 11.111111 11.111111 10.500000 11.444444 10.111111 11.444444 10.888889
[50] 11.055556 10.277778 10.500000 10.944444 10.777778 10.944444 10.444444
[57] 11.055556 11.000000 11.055556 11.666667

Функции

Нужно придумать:

  • имя функции. Оно не должно совпадать с именами функций из базового R или пакетов, которые вы используете. В идеале оно отражает смысл вашей функции.
# не делайте так!
mean <- function(x) { sum(x) }

  • список параметров, которые функция принимает на вход. Например, function(x, y, z).

  • сам код, выполняющий работу, который вы записываете в тело функции внутри {...}.

smart_name <- function(input1, input2, param3) {
  ...
  body
  ...
}

Что на вход и выход?

  • У функции может быть 0 или несколько параметров.
  • Функция может возвращать максимум 1 объект.
# ничего не требует
say_hello <- function() {
  print("Привет!")
}


# ничего не возвращает
save_res <- function(df) {
  df = df[df$pval < 0.05, c(1,3:5)]
  write.csv(df, "path-to-file.csv")
}
say_hello()
[1] "Привет!"
save_res(multiple_tests)

return()

  • Функция возвращает результат последнего выражения либо то, что указано как return(...).
fahrenheit2celsius <- function(x) {
  y = (x - 32) * 5 / 9
  return(y)
}

fahrenheit2celsius(54)
[1] 12.22222
fahrenheit2celsius <- function(x) {
  return(x)
  y = (x - 32) * 5 / 9
  y
}

fahrenheit2celsius(54)
[1] 54

Выполнение кода по условию 🏠

if (condition) {
  # что делать, когда condition = TRUE
} else {
  # что делать, когда condition = FALSE
}

Логическое выражение condition должно возвращать либо TRUE, либо FALSE.

Несколько условий 🏠

if (this) {
  # делай это, когда первое условие выполняется
} else if (that) {
  # делай что-то другое, когда первое условие не выполняется, а второе выполняется
} else {
  # делай что-то третье, когда ни одно из условий выше не выполняется
}

Не путайте else if () {...} с ifelse().

Параметры 🏠

Через параметры на вход функции передаются данные или какие-то детали. Обычно данные передаются первому параметру. В таком случае эту функцию будет легко использовать с конвейером (%>% или |>).

Для параметров можно задать значение по умолчанию:

fahrenheit2celsius <- function(x = 32) {
  (x - 32) * 5 / 9
}

fahrenheit2celsius()
[1] 0
fahrenheit2celsius(86)
[1] 30

Названия параметров 🏠

Идеи для названий параметров:

  • x, y, z - вектора,
  • w - вектор весов,
  • df - датафрейм,
  • i, j - индексы (строки и столбцы),
  • n - длина или число строк,
  • p - число столбцов.

Провека формата входных данных 🏠

В каком случае нужно остановиться.

if() + stop()

fahrenheit2celsius <- function(x) {
  if (!is.numeric(x)) {        
    stop("`x` must be numeric") 
    }
  (x - 32) * 5 / 9
}

fahrenheit2celsius("23")
Error in fahrenheit2celsius("23"): `x` must be numeric

stopifnot()

fahrenheit2celsius <- function(x) {
  stopifnot(is.numeric(x))
  (x - 32) * 5 / 9
}

fahrenheit2celsius("23")
Error in fahrenheit2celsius("23"): is.numeric(x) не TRUE

Multiple returns 🏠

Чтобы функция возвращала несколько объектов, нужно эти объекты возвращать в виде списка.

return_two_and_four <- function(){
  list(2, 4)
}

return_two_and_four()
[[1]]
[1] 2

[[2]]
[1] 4

Локальные переменные 🏠

x внутри функции (в ее среде) и вне функции (в глобальной среде) существуют независимо.

x <- 1000

add_ten <- function(x){
  x + 10
}

add_ten(32)
[1] 42
x
[1] 1000

Глобальные переменные 🏠

Изнутри функции можно переписать глобальную переменную с помощью оператора <<-.

x <- 1000

add_ten <- function(x){
  x <<- 32
  x + 10
}

add_ten(32)
[1] 42
x
[1] 32

Глобальные переменные 🏠

Если R не нашел переменную в среде функции, то он будет искать ее в глобальной среде.

y <- 500

add_ten <- function(){
  y + 10
}

add_ten()
[1] 510

Функции, работающие с датафреймами 🏠

  • Хотим название столбца использовать как один из параметров
library(tidyverse)
summary6 <- function(df, var) {
  summarize(
    df,
    min = min(var, na.rm = TRUE),
    mean = mean(var, na.rm = TRUE),
    sd = sd(var, na.rm = TRUE),
    median = median(var, na.rm = TRUE),
    iqr = IQR(var, na.rm = TRUE),
    max = max(var, na.rm = TRUE),
    .groups = "drop")
}

summary6(mtcars, var = mpg)
Error in `summarize()`:
ℹ In argument: `min = min(var, na.rm = TRUE)`.
Caused by error in `FUN()`:
! only defined on a data frame with all numeric-alike variables

Функции, работающие с датафреймами 🏠

  • Нужно использовать {{ ... }} - embracing
  • то же с функциями, работающими с ggplot2
summary6 <- function(df, var) {
  summarize(
    df,
    min = min({{ var }}, na.rm = TRUE),
    mean = mean({{ var }}, na.rm = TRUE),
    sd = sd({{ var }}, na.rm = TRUE),
    median = median({{ var }}, na.rm = TRUE),
    iqr = IQR({{ var }}, na.rm = TRUE),
    max = max({{ var }}, na.rm = TRUE),
    .groups = "drop")
}

summary6(mtcars, var = mpg)
   min     mean       sd median   iqr  max
1 10.4 20.09062 6.026948   19.2 7.375 33.9

Импорт функции из файла 🏠

  • В R скрипте записаны одна или несколько функций
# temperature_conversion.R

fahrenheit2celsius <- function(x) {
  (x - 32) * 5 / 9
}

celsius2kelvin <- function(x) {
  x + 273.15
}
  • Импорт функций из R скрипта в текущую R сессию
source("temperature_conversion.R")

Выполнение R скрипта в командной строке 🏠

# say_hello.R

say_hello <- function(name = "Anya") {
  print(paste0("Hello, ", name, "!"))
}

say_hello()


Rscript say_hello.R
[1] "Hello, Nastya!"

Передача параметров скрипту из командной строки 🏠

# say_hello_to_anyone.R

args <- commandArgs(trailingOnly=TRUE)

say_hello <- function(name = "Nastya") {
  print(paste0("Hello, ", name, "!"))
}

if (length(args)==0) {
  say_hello()
} else {
  name = args[1]
  say_hello(name)
}


Rscript say_hello_to_anyone.R "Anya"
[1] "Hello, Nastya!"

Выполнение R скрипта в фоне 🏠

  • Запустить выполнение скрипта в текущей R сессии:
source("my_script.R")
  • Чтобы запустить выполнение скрипта в другой (фоновой) R сессии:

Применяем функции

Векторизация

Многие функции в R приспособлены к работе с векторами и итерации по элементам вектора. Благодаря этому нам не нужно лишний раз писать циклы.

1:5 * 2
[1]  2  4  6  8 10
1:5 > 3
[1] FALSE FALSE FALSE  TRUE  TRUE
out_vec <- c()

for (i in 1:ncol(airquality)) {
  out_vec[i] = typeof(airquality[, i])
}

out_vec
[1] "integer" "integer" "double"  "integer" "integer" "integer"

Как избегать циклов?

Можно использовать функции из семейства apply:

  • apply()
  • lapply()
  • sapply()
  • tapply()

Как избегать циклов?

Если нужно применить функцию к элементам некоторого вектора или списка, то используйте map из пакета purrr (входит в tidyverse).

library(tidyverse)
# library(purrr)

# вместо
for (i in 1:3){
  f(i)
}

# или
list(f(1), f(2), f(3))

# нужно всего лишь...
map(1:3, f)

Семейство функций map

cube <- function(x) x ** 3

map(1:3, cube)
[[1]]
[1] 1

[[2]]
[1] 8

[[3]]
[1] 27

Разные map_

  • Простой map() всегда возвращает list.
  • Если вы уверены, что ваш результат подходит под определение вектора (данные одного типа), используйте map_:
    • map_chr()
    • map_lgl()
    • map_int()
    • map_dbl()
  • Для других типов данных, которые можно записать в вектор (даты, факторы, …), можно использовать map_vec().

Разные map_

map_chr(mtcars, typeof)
     mpg      cyl     disp       hp     drat       wt     qsec       vs 
"double" "double" "double" "double" "double" "double" "double" "double" 
      am     gear     carb 
"double" "double" "double" 
map_lgl(mtcars, is.double)
 mpg  cyl disp   hp drat   wt qsec   vs   am gear carb 
TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE 
map_int(mtcars, n_distinct)
 mpg  cyl disp   hp drat   wt qsec   vs   am gear carb 
  25    3   27   22   22   29   30    2    2    3    6 
map_vec(c("a", "b", "c"), factor)
[1] a b c
Levels: a b c

Анонимные функции

lambda functions

  • Посчитаем число уникальных значений в каждом столбце mtcars (забыли про n_distinct() из dplyr)
# Полная запись:
map_dbl(mtcars, function(x) length(unique(x)))
 mpg  cyl disp   hp drat   wt qsec   vs   am gear carb 
  25    3   27   22   22   29   30    2    2    3    6 
# Лень писать так много. Так проще:
map_dbl(mtcars, \(x) length(unique(x)))
 mpg  cyl disp   hp drat   wt qsec   vs   am gear carb 
  25    3   27   22   22   29   30    2    2    3    6 
# Устаревший вариант записи:
map_dbl(mtcars, ~ length(unique(.x)))
 mpg  cyl disp   hp drat   wt qsec   vs   am gear carb 
  25    3   27   22   22   29   30    2    2    3    6

Задание 1

Возьмите встроенный набор данных - msleep.

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

map_df

Удобно использовать для чтения нескольких таблиц одного формата в единый датафрейм.

  • map_dfr - это map() + bind_rows()
  • map_dfc - это map() + bind_cols()
  • bind_rows() и bind_cols() - аналоги rbind() и cbind() из базового R

map_df

  • 3 файла, содержащими таблицы с одинаковыми названиями и порядком столбцов
input_files <- c("file1.csv", "file2.csv", "file3.csv")
  • Муторный вариант:
file1 <- read_csv("file1.csv")
file2 <- read_csv("file2.csv")
file3 <- read_csv("file3.csv")

file <- bind_rows(file1, file2, file3)
  • Оптимальный вариант:
file <- map_dfr(input_files, read_csv)

Вместо map_df

  • С версии purrr 1.0.0 для чтения файлов в один датафрейм рекомендуется использовать map() + list_cbind() или map() + list_rbind()
file <- input_files %>% 
  map(read_csv) %>%     # получился список из датафремов
  list_rbind()          # объединили датафреймы из списка

map ❤️ списки

  • Плюс новые функции из purrr 1.0.0: list_flatten(), list_simplify(), list_c()
x <- list(
  list(-1, x = 1, y = c(2), z = "a"),
  list(-2, x = 4, y = c(5, 6), z = "b"),
  list(-3, x = 8, y = c(9, 10, 11)))

map_dbl(x, "x") # по имени элемента
[1] 1 4 8
map_dbl(x, list("y", 1)) # по позиции
[1] 2 5 9
map_chr(x, "z", .default = NA) 
[1] "a" "b" NA

Аргументы функции

x <- list(1:5, c(1:10, NA))
  • Вариант не очень:
map_dbl(x, \(x) mean(x, na.rm = TRUE))
[1] 3.0 5.5
  • Вариант получше:
map_dbl(x, mean, na.rm = TRUE)
[1] 3.0 5.5

map2

  • Если нужно итерировать и по элементам списка, и по вектору аргумента функции.
xs <- map(1:8, \(x) runif(10))
ws <- map(1:8, \(x) rpois(10, 5) + 1)

map2_dbl(xs, ws, weighted.mean, na.rm = TRUE)
[1] 0.6281366 0.5560616 0.3529188 0.3790625 0.4429240 0.4149723 0.5539007
[8] 0.5637695

pmap

  • Когда не хватает map2, а нужен map3 или даже map4…
  • Нужно подать список всех аргументов функции.
  • map2(x, y, f) - то же, что и pmap(list(x, y), f).
pmap_dbl(list(xs, ws), weighted.mean)
[1] 0.6281366 0.5560616 0.3529188 0.3790625 0.4429240 0.4149723 0.5539007
[8] 0.5637695

imap

  • imap(x, f) - то же, что и map2(x, seq_along(x), f) или map2(x, names(x), f)
# x - элемент списка
# idx - это название (индекс) элемента списка
x <- map(1:6, \(x) sample(1000, 10))
imap_chr(x,  \(x, idx) paste0("The highest value of ", idx, " is ", max(x)))
[1] "The highest value of 1 is 820" "The highest value of 2 is 896"
[3] "The highest value of 3 is 938" "The highest value of 4 is 987"
[5] "The highest value of 5 is 945" "The highest value of 6 is 942"

Задание 2

Возьмите встроенный набор данных - msleep.

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

walk

  • Для функций, которые не возращают никакие результаты в R сессию
ggplots <- list(gg1, gg2, gg3)
output_files <- c("plot1.png", "plot2.png", "plot3.png")

walk2(output_files, ggplots, ggsave)

Объединение датафреймов по ключу

  • Как объединить два датафрейма по совпадающим значениям (“ключам”) в одном из столбцов (или нескольких столбцах)?
x <- tribble(
  ~key, ~val_x,
     1, "x1",
     2, "x2",
     3, "x3")
y <- tribble(
  ~key, ~val_y,
     1, "y1",
     2, "y2",
     4, "y3")


inner_join(x, y, by = "key")
# A tibble: 2 × 3
    key val_x val_y
  <dbl> <chr> <chr>
1     1 x1    y1   
2     2 x2    y2

Объединение датафреймов

Объедиение происходит по одному или нескольким столбцам с “ключами”.

  • inner_join() - сохранить только строчки с общими “ключами”
  • full_join() - сохранить все строчки из обеих таблиц
  • left_join() - сохранить все строчки из “первой” таблицы
  • right_join() - сохранить все строчки из “второй” таблицы
full_join(x, y, by = "key")
# A tibble: 4 × 3
    key val_x val_y
  <dbl> <chr> <chr>
1     1 x1    y1   
2     2 x2    y2   
3     3 x3    <NA> 
4     4 <NA>  y3

Семейство функций reduce

  • reduce() берет на вход вектор длины n и возвращает вектор длины 1, применяя функцию к элементам вектора попарно.
  • Удобно, если надо объединить несколько датафреймов.
# вместо
f(f(f(1, 2), 3), 4)

# нужно всего лишь...
reduce(1:4, f)

reduce2()

Если нужно, например, объединить несколько датафреймов, но использовать разные переменные, по которым делать объединение.

Отслеживать и ловить ошибки 🏠

  • safely() возвращает список из двух элементов:

    • result нужный результат или NULL если была ошибка,
    • error error object или NULL, если ошибки не было.

  • possibly() позволяет использовать default value, если возникает ошибка.

  • quietly() выдает result, output, messages и warnings.

safely() 🏠

  • Можно воспринимать safely() как наречие, а функцию, к которой оно применяется, как глагол
safe_log <- safely(log)
str(safe_log(10))
List of 2
 $ result: num 2.3
 $ error : NULL
str(safe_log("a"))
List of 2
 $ result: NULL
 $ error :List of 2
  ..$ message: chr "нечисловой аргумент для математической функции"
  ..$ call   : language .Primitive("log")(x, base)
  ..- attr(*, "class")= chr [1:3] "simpleError" "error" "condition"

safely() 🏠

  • safely() удобно использовать вместе с map
x <- list(1, 10, "a")
y <- x %>% map(safely(log))
str(y)
List of 3
 $ :List of 2
  ..$ result: num 0
  ..$ error : NULL
 $ :List of 2
  ..$ result: num 2.3
  ..$ error : NULL
 $ :List of 2
  ..$ result: NULL
  ..$ error :List of 2
  .. ..$ message: chr "нечисловой аргумент для математической функции"
  .. ..$ call   : language .Primitive("log")(x, base)
  .. ..- attr(*, "class")= chr [1:3] "simpleError" "error" "condition"

safely() 🏠

Можно отделить результаты от ошибок:

  • Транспонируем список:
y <- transpose(y)
str(y)
List of 2
 $ result:List of 3
  ..$ : num 0
  ..$ : num 2.3
  ..$ : NULL
 $ error :List of 3
  ..$ : NULL
  ..$ : NULL
  ..$ :List of 2
  .. ..$ message: chr "нечисловой аргумент для математической функции"
  .. ..$ call   : language .Primitive("log")(x, base)
  .. ..- attr(*, "class")= chr [1:3] "simpleError" "error" "condition"
  • Затем находим элементы, которые вызвали ошибку:
is_ok <- y$error %>% map_lgl(is_null)
x[!is_ok]
[[1]]
[1] "a"

Подробный текст ошибок у map 🏠

  • С версии purrr 1.0.0 текст ошибок map стал более подробным: указывается индекс элемента на котором все сломалось.
  • Но выполнение кода из-за ошибки останавливается.
map(x, log)
Error in `map()`:
ℹ In index: 3.
Caused by error:
! нечисловой аргумент для математической функции

possibly() 🏠

Позволяет использовать default value, если возникает ошибка.

x <- list(1, 10, "a")
x %>% map_dbl(possibly(log, otherwise = NA_real_))
[1] 0.000000 2.302585       NA

От purrr к furrr

  • тот же map, но с возможностью параллелизации процессов
# install.packages("furrr")
library(furrr)

# План того, как код должен быть запущен
plan(multisession, workers = 2)

# Код распараллелизован
future_map_dbl(1:3, cube)
[1]  1  8 27

Что почитать про функции и purrr

Конец!