Введение в R: часть 2

Pavel Polishchuk, 2014


Содержание

  1. Работа с текстовыми и бинарными файлами (чтение/запись).
  2. Создание пользовательских функций, их сохранение и повторное использование.
  3. Манипуляции с данными (векторы, data.frames).
  4. Простейшие функции для работы с графикой.
  5. Семейство функций apply.

Чтение/запись файлов в R

Чтение текстовых файлов

В базовом пакете base есть функция read.table, которая отвечает за чтение текстовых файлов. Она имеет множество параметров, что позволяет гибко управлять процессом чтения. Эта функция возвращает data.frame.

y1 <- read.table("data/sol_y1.txt", 
                header=TRUE, sep="\t", as.is=TRUE, 
                check.names=FALSE, comment.char="", 
                row.names=1)

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

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

y1 <- read.table("D:/Teaching/R/introduction/md/data/sol_y1.txt", 
                header=TRUE, sep="\t", as.is=TRUE, 
                check.names=FALSE, comment.char="", 
                row.names=1)

Для представления пути используются либо символы слэша (/) либо удвоенные символы обратного слэша (\). Это обусловлено тем, что обратный слэш в R это спецсимвол и добавление второго слэша его экранирует и указываает R, что следующий за ним символ надо вопринимать как обычный.
Т.е. обе эти записи пути эквивалетны: "D:/Teaching/R/introduction/md/data/sol_y1.txt" или "D:\\Teaching\\R\\introduction\\md\\data\\sol_y1.txt"


Посмотрим на данные, которые были загружены. Сделать это можно несколькими способами, вот некоторые из них

head(y1)
            Sol
sol_10001 -3.18
sol_10002 -2.64
sol_10003 -3.84
sol_10004 -3.74
sol_10005 -3.55
sol_10006 -3.10
head(y1, 10)
            Sol
sol_10001 -3.18
sol_10002 -2.64
sol_10003 -3.84
sol_10004 -3.74
sol_10005 -3.55
sol_10006 -3.10
sol_10007 -3.30
sol_10008 -4.53
sol_10009 -3.85
sol_10010 -5.24
tail(y1)
            Sol
sol_10795 -3.33
sol_10796 -3.31
sol_10797 -2.73
sol_10798 -0.12
sol_10799 -0.47
sol_10800  0.06

Функция read.table является не очень эффективной с точки зрения производительности. Файлы, содержащие большое число столбцов считываются очень медленно. Если попробовать прочитать файл sol_x1.txt, который содержит 4058 дескрипторов для 800 соединений, то это займет заметное время.
Можете попробовать выполнить следующий код

x1 <- read.table("data/sol_x1.txt", 
                 header=TRUE, sep="\t", as.is=TRUE, 
                 check.names=FALSE, comment.char="", 
                 row.names=1)

Для преодоления этого недостатка есть несколько вариантов, вот два из них:
1. Использовать экспериментальную функцию fread из пакета data.table.
2. Конвертировать текстовые файлы в бинарный формат и в дальнейшем работать с этими файлами.


Чтение текстовых файлов с использованием функции fread (пакет data.table)

x1 <- data.table::fread("data/sol_x1.txt", sep="\t", header=TRUE)

Для просмотра загруженной таблицы в RStudio можно использовать команду

View(x1)

или кликнуть по имени переменной x1 в списке alt text

Результат в обоих случаях будет одинаковым alt text

Проверим объект какого класса получился при загрузке таблицы с использованием функции fread

class(x1)
[1] "data.table" "data.frame"

Приведем объект x1 к типу data.frame

x1 <- as.data.frame(x1)
class(x1)
[1] "data.frame"

Отметим, что первая колонка содержит названия соединений. Переместим значения этой колонки в поле rownames

# создание имен строк
rownames(x1) <- x1[,1]

# удаление первой колонки
x1 <- x1[,-1]

Проверим, что получилось. alt text


Сохранение данных в бинарном формате

Любые объекты созданные в R можно сохранить в бинарном формате в виде файлов .RData. Это позволяет осуществлять быструю загрузку и доступ к сохраненным данным.

Таким образом можно сохранить данные из однажды прочитанного текстового файла в бинарном формате и в дальнейшем загружать эти данные из него.

Чтобы сохранить содержимое переменной x1 используем следующий вызов (указание ключевого слова file является обязательным)

save(x1, file="data/sol_x1.RData")

Очистка рабочей области от загруженных и используемых переменных

Перед загрузкой данных, только что сохраненных в файл sol_x1.RData, сперва очистим рабочую область (удалим все загруженные переменные). Сделать это можно либо функцией rm. Приведенный вызов удалит все загруженные объекты.

rm(list = ls())

Либо кликнув по кнопке Clear. alt text

Объекты также можно удалять выборочно, указывая их названия.

rm(x1)

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

gc()
          used (Mb) gc trigger (Mb) max used (Mb)
Ncells  342570  9.2     597831   16   497955 13.3
Vcells 2397623 18.3    6411868   49  6411635 49.0

Чтение данных в бинарном формате

Загрузим ранее сохраненный файл sol_x1.RData

load("data/sol_x1.RData")

Как вы могли заметить в списке загруженных переменных появилась новая с именем x1 как у ранее сохраненного data.frame.
Однако если в разных файлах сохранены объекты с одинаковым названием, то для их одновременной загрузки эти объекты необходио переименовать. Это можно во время загрузки данных, использовав следующий набор команд:

x1.1 <- local(get(load("data/sol_x1.RData")))

Проверим загруженный объект x1.1 на идентичность с x1

all.equal(x1, x1.1) 
[1] TRUE

Сохранение данных в текстовом формате

Функция write.table сохраняет данные в текстовом формате. Она имеет много параметров для гибкой настройки вида сохраняемого файла.

write.table(y1, "data/new_y.txt", quote=FALSE, sep="\t", row.names=TRUE, col.names=TRUE)

Другие функции для работы с текстовыми файлами

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

read.csv(file, header = TRUE, sep = ",", quote = "\"",
         dec = ".", fill = TRUE, comment.char = "", ...)
read.csv2(file, header = TRUE, sep = ";", quote = "\"",
          dec = ",", fill = TRUE, comment.char = "", ...)
read.delim(file, header = TRUE, sep = "\t", quote = "\"",
           dec = ".", fill = TRUE, comment.char = "", ...)
read.delim2(file, header = TRUE, sep = "\t", quote = "\"",
            dec = ",", fill = TRUE, comment.char = "", ...)
write.csv(...)
write.csv2(...)

Другие способы ввода-вывода

Помимо текстовых и бинарных данных R может читать данные из различных баз данных, а также файлов MS Excel, для чего используется подключение дополнительных пакетов. В частности пакет xlsx позволяет работать с данными, сохраненными в одноименном формате.


Создание пользовательских функций

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

Создадим простейшую функцию рассчитывающую разность двух векторов

f1 <- function(x, y) {
  z <- x - y
  return(z)
}

Применим ее

f1(1:5, 10:6)
[1] -9 -7 -5 -3 -1
a <- 1:5
b <- 2:6
f1(a, b)
[1] -1 -1 -1 -1 -1
f1(a, a)
[1] 0 0 0 0 0

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

f1(a, b)
[1] -1 -1 -1 -1 -1
f1(b, a)
[1] 1 1 1 1 1
f1(y=b, x=a)
[1] -1 -1 -1 -1 -1

Продемонстрируем области видимости переменных на примере собственных функций.

f2 <- function() {
  x <- x + 1
  return(x)
}

Создадим вектор x и вызовем созданную нами функцию. Результат функции это измененный вектор x, однако сам вектор x не изменился.

x <- 1:5
f2()
[1] 2 3 4 5 6
x
[1] 1 2 3 4 5

Таким образом все что передается в функцию попадает в локальную области видимости только этой функции как и все изменения.


Реализуем собственную функцию загрузки бинарных файлов

local.load <- function(fname) {
  local(get(load(fname)))
}

Ключевое слово return в завершении функции писать не обязательно. По умолчанию функция возвращает результат последней операции.

Для того чтобы сохранить созданную функцию текст функции посещается в (новый) файл R script, который можно создать вызвав меню File - New file - R script или нажав комбинацию Ctrl + Shift + N. После чего пишется текст функции и файл сохраняется с расширением “.R”.

alt text

Для загрузки функции из файла используется функция source. Но предварительно очистим содержимое рабочей области.

rm(list = ls())
source("scripts/load.R")

После чего для загрузки бинарного файла достаточно вызвать

x1 <- local.load("data/sol_x1.RData")

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

x1 <- read.table("data/sol_x1.txt", 
                 header=TRUE, sep="\t", as.is=TRUE, 
                 check.names=FALSE, comment.char="", 
                 row.names=1)

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

read_descriptors_1 <- function(fname) {
  df <- read.table(fname, header=TRUE, sep="\t", as.is=TRUE, 
                   check.names=FALSE, comment.char="", row.names=1)
  return(df)
}

Второй вариант

read_descriptors_2 <- function(fname, header=TRUE, sep="\t", 
                             as.is=TRUE, check.names=FALSE, 
                             comment.char="", row.names=1, ...) {
  df <- read.table(fname, header=header, sep=sep, as.is=as.is, 
                   check.names=check.names, comment.char=comment.char, 
                   row.names=row.names, ...)
  return(df)
}

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

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

x1 <- read_descriptors_1("data/sol_x1.txt")

или

x1 <- read_descriptors_2("data/sol_x1.txt")

Какие достоинства и недостатки у каждого из вариантов?


Задания

  1. Создать функцию, которая будет считывать текстовый файл, содержащий дескрипторы, с использованием функции fread, и возвращать data.frame.
  2. Прочитать с использованием функции из п.1 и сохранить в бинарном формате файл sol_x2.txt.
  3. Прочитать файл sol_y1.txt, конвертировать загруженные данные в именованный вектор и сохранить его в бинарном формате. Создать для этого соответствующую функцию.
  4. Повторить операции п.3 для файла sol_y2.txt

Манипуляции с данными

Векторы

Загрузим файл sol_y1.RData и определим среднее значение растворимости в выборке, медианное значение, максимальное и минимальное значения.

y1 <- local.load("data/sol_y1.RData")
mean(y1)
[1] -2.591
median(y1)
[1] -2.295
min(y1)
[1] -11.62
max(y1)
[1] 1.58
range(y1)
[1] -11.62   1.58

Все приведенные функции векторизированы, поэтому расчет происходит очень быстро и эффективно.

Функция summary возвращает всю вышеприведенную статистику в виде одного вектора значений.

summary(y1)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -11.60   -3.77   -2.29   -2.59   -1.06    1.58 

Помимо этих есть еще множество функций рассчета различных статистических характеристик.


Data.frames

Продемонстрируем работу с data.frames на примере.
Есть два набора соединений с рассчитанными дескрипторами и их необходимо объединить в один data.frame.

Загрузим оба имеющихся файла данных:

x1 <- local.load("data/sol_x1.RData")
x2 <- local.load("data/sol_x2.RData")

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

dim(x1)
[1]  800 4058
dim(x2)
[1]  233 3518

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

Функция setdiff возвращает элементы первого вектора, которые отсутствуют во втором

var.names1 <- setdiff(colnames(x1), colnames(x2))
var.names2 <- setdiff(colnames(x2), colnames(x1))

var.names1 содержит имена дескрипторов, которые присутствуют в x1 и отсутствуют в x2. Аналогично для var.names2.

head(var.names1, 10)
 [1] "S_A(chg)/A_A_A_C/1_2s,3_4a/3"     
 [2] "S_A(chg)/A_A_A_D/1_4d,2_4a,3_4a/5"
 [3] "S_A(chg)/A_A_A_D/1_4s,2_3a,2_4s/6"
 [4] "S_A(chg)/A_A_A_D/2_3a,2_4s/4"     
 [5] "S_A(chg)/A_A_A_D/2_4d,3_4a/4"     
 [6] "S_A(chg)/A_A_A_D/2_4s,3_4t/4"     
 [7] "S_A(chg)/A_A_B_B/1_2d,3_4a/3"     
 [8] "S_A(chg)/A_A_B_B/1_2t,3_4s/3"     
 [9] "S_A(chg)/A_A_B_B/1_3d,2_4d,3_4s/6"
[10] "S_A(chg)/A_A_B_B/1_3s,2_4s/3"     
head(var.names2, 10)
 [1] "S_A(chg)/A_A_A_B/1_2s,3_4s/3"     
 [2] "S_A(chg)/A_A_A_C/1_2a,3_4a/3"     
 [3] "S_A(chg)/A_A_A_C/2_4s,3_4d/4"     
 [4] "S_A(chg)/A_A_A_D/1_2s,3_4t/3"     
 [5] "S_A(chg)/A_A_A_D/1_4a,2_4a,3_4a/5"
 [6] "S_A(chg)/A_A_B_B/2_3s,2_4s/4"     
 [7] "S_A(chg)/A_A_B_C/1_2a,1_4s,3_4a/6"
 [8] "S_A(chg)/A_A_B_C/1_2s,2_4s,3_4d/6"
 [9] "S_A(chg)/A_A_B_C/1_3d,2_4a/3"     
[10] "S_A(chg)/A_A_B_C/1_3s,2_4a/3"     

Посмотрим сколько таких дескриптров в каждом случае

length(var.names1)
[1] 923
length(var.names2)
[1] 383

Т.е. первый набор данных содержит 923 дескриптора, которые отсутствуют вов втором наборе, а второй - 383, которые отсутствуют в первом.

Поскольку мы используем фрагментные дескрипторы, то их отсутствие фактически означает, что число дескриторов данного вида равно нулю. Следовательно мы должны добавить отсутствующие дескрипторы в каждый набор данных и приравнять все их значения нулю.
Чтобы добавить нулевые значения этим дескрипторам мы вызываем их для соответствующего data.frame и присваиваем значение 0.

x1[,var.names2] <- 0
x2[,var.names1] <- 0

Проверям размерность полученных данных

dim(x1)
[1]  800 4441
dim(x2)
[1]  233 4441

Число колонок (дескрипторов) теперь идентично. Проверим порядок следования дескриптров в каждом наборе, используя следующую конструкцию

all(colnames(x1) == colnames(x2))
[1] FALSE

Выражение в скобках возвращает вектор логических значений, который содержит TRUE в случае если на одинаковых позициях находятся одинаковые значения и FALSE в противном случае. Функция all возвращает TRUE, если все значения логического вектора равны TRUE.

Проверим действительно ли имена дескрипторов в обоих data.frames совпадают и мы ничего не упустили. Используем для этого функцию %in%, которая принимает два вектора, и если значение первого вектора присутствует во втором, то для этого элемента возвращается значение TRUE. Таким образом выражение записанное ниже можно прочитать как “все ли названия колонок первого набора данных присутствуют во втором наборе”

all(colnames(x1) %in% colnames(x2))
[1] TRUE

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

all(colnames(x2) %in% colnames(x1))
[1] TRUE

Две эти проверки можно заменить одной, если предварительно отсортировать имена колонок в обоих наборах почле чего сравнить их

all(sort(colnames(x1)) == sort(colnames(x2)))
[1] TRUE

Теперь перед объединением строк обоих data.frames необходимо расположить колонки в одинаковом порядке. Воспользуемся для этого свойством индексов и расположим колонки в x2 в той же последовательности что и в x1

x2 <- x2[,colnames(x1)]

Теперь колонки расположены в одной последовательности. Проверим

all(colnames(x1) == colnames(x2))
[1] TRUE

Полученные data.frames теперь можно объединить в один с использованием функции rbind - объединение данных по строкам (row bind).

x <- rbind(x1, x2)
dim(x)
[1] 1033 4441

Аналогично когда надо объединить данные по колонкам используется функция cbind.


Загрузим и объединим значения свойства для двух наборов данных

y1 <- local.load("data/sol_y1.RData")
y2 <- local.load("data/sol_y2.RData")
y <- c(y1, y2)
head(y)
sol_10001 sol_10002 sol_10003 sol_10004 sol_10005 sol_10006 
    -3.18     -2.64     -3.84     -3.74     -3.55     -3.10 

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

summary(y)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -11.60   -3.89   -2.41   -2.70   -1.32    1.58 

Простейшие функции для работы с графикой

Базовый пакет R имеет в своем составе много функций для визуализации данных. Остановимся только на некоторых, которые могут быть полезны при первом знакомстве с данными.

Функция hist - строит гистограмму распределения какой-либо величины. Может помочь составить первое впечатление о нормальности распределения данных. Применим эту функцию к свойству y.

hist(y)

plot of chunk unnamed-chunk-54

Если стандартное отображение не устраивает, то можно его изменить, например, сделав интервалы меньше.

hist(y, seq(-12, 2, 0.5))

plot of chunk unnamed-chunk-55

hist(y, 50)

plot of chunk unnamed-chunk-55

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

qqnorm(y)
qqline(y, col=2)

plot of chunk unnamed-chunk-56


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

1, R Base Graphics: An Idiot's Guide http://rpubs.com/SusanEJohnston/7953
2. Различные примеры работы с графикой (и не только) http://flowingdata.com/category/tutorials/


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

apply

Семейство функций apply является более удобным аналогом цикла for.

Рассмотрим пример - необходимо найти среднее значение каждой колонки в data.frame.

df <- data.frame(var1=c(11,21,31), var2=c(12,22,32), var3=c(13,23,33), var4=c(14,24,34), row.names=c("case1", "case2", "case3"))

df
      var1 var2 var3 var4
case1   11   12   13   14
case2   21   22   23   24
case3   31   32   33   34

Решение с использованием цикла for

result <- vector()
for (i in 1:ncol(df)) {
  result[i] <- mean(df[,i])
}
result
[1] 21 22 23 24

Решение с использованием функции apply

result <- apply(df, 2, mean)
result
var1 var2 var3 var4 
  21   22   23   24 

Правда короче и проще? Кроме того вектор с результатами содержит названия переменных.

А вот выражение, которое вычисляет среднее значение по строкам

result <- apply(df, 1, mean)
result
case1 case2 case3 
 12.5  22.5  32.5 

На самом деле для операции нахождения среднего и суммы по строкам/столбцам есть отдельные векторизированные функции rowMeans, rowSums и т.д.

colMeans(df)
var1 var2 var3 var4 
  21   22   23   24 

Общий вид вызова функции apply

apply(X, MARGIN, FUN, ...)

X - матрица, массив или data.frame
MARGIN - порядковый номер размерноси, к которой будет применяться функция FUN (1 для строк, 2 для колонок)
FUN - применяемая функция

Результатом будет являться вектор (если используемая функция возвращает одно значение), массив (если функция возвращает вектор значений) или список (если функция возвращает результат в виде более сложной структуры данных, например data.frame или матрица).

ВАЖНО! При применении функции apply к data.frame, data.frame неявно приводится к матрице. Матрицы как мы помним содержат данные только одного типа. Это означеат что если в data.frame 9 числовых колонок и 1 текстовая, то будет произведена конвертация в текстовую матрицу, и следовательно все действия будут производиться над текстовой матрицей.
Пример

df2 <- data.frame(v1=1:3, v2=3:5, v3=c("a", "b", "c"))
df2
  v1 v2 v3
1  1  3  a
2  2  4  b
3  3  5  c
apply(df2, 2, is.character)
  v1   v2   v3 
TRUE TRUE TRUE 

Чуть более сложный пример - надо посчитать для каждой колонки среднее значение и извлечь из него квадратный корень. выполним это с использованием функции apply

result <- apply(df, 2, function(v) {sqrt(mean(v))})
result
 var1  var2  var3  var4 
4.583 4.690 4.796 4.899 

Предложите альтернативный вариант расчета.


sapply & lapply

Помимо функции apply есть еще функции sapply и lapply, отличие которых состоит в том, что на вход они могут принимать вектор или список и возвращают вектор/матрицу (sapply) или список (lapply).

Вспомним, что data.frame является списком векторов-столбцов. Тогда предыдущий пример можно переписать как

sapply(df, function(v) sqrt(mean(v)))
 var1  var2  var3  var4 
4.583 4.690 4.796 4.899 

Использование lapply вернет уже список

lapply(df, function(v) sqrt(mean(v)))
$var1
[1] 4.583

$var2
[1] 4.69

$var3
[1] 4.796

$var4
[1] 4.899

Возведем значения каждой колонки в степень, соответствующую номеру этой колонки - результатом будет новая матрица.

sapply(1:4, function(i) {df[,i]^i})
     [,1] [,2]  [,3]    [,4]
[1,]   11  144  2197   38416
[2,]   21  484 12167  331776
[3,]   31 1024 35937 1336336

Обратите внимание, что в этом случае в качестве первого параметра передается не data.frame, а вектор индексов колонок, по которым происходит итерация.

sapply(1:ncol(df), function(i) {df[,i]^i})
     [,1] [,2]  [,3]    [,4]
[1,]   11  144  2197   38416
[2,]   21  484 12167  331776
[3,]   31 1024 35937 1336336
sapply(seq_along(df), function(i) {df[,i]^i})
     [,1] [,2]  [,3]    [,4]
[1,]   11  144  2197   38416
[2,]   21  484 12167  331776
[3,]   31 1024 35937 1336336

Функция seq_along возвращает вектор индексов с первого до последнего элемента вектора/списка.

seq_along(df)
[1] 1 2 3 4

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

Например, надо определить класс каждой колонки в data.frame

sapply(df2, class)
       v1        v2        v3 
"integer" "integer"  "factor" 

Как видим в отличие от apply фугкция sapply не производит неявной конвертации data.frame в матрицу и типы данных в колонках остаются неизменными.

Задание

Проверить совпадает ли порядок следования имен соединений в x и y.

Домашнее задание

Создать функцию, которая читала бы формат файлов дескрипторов dat/cds. Подсказка: можно использовать функцию readBin.

Создать функцию, которая будет объединять два набора фрагментных дескрипторов.