Математик среди биологов: заметки с тегом кластерный анализ

Добавляем рисунки к листьям дендрограммы в R: используем ggtree

Антон Лях — Sun, 26 Nov 2023 11:31:42 +0300

Дендрограмма — это результат объединения данных при помощи иерархического кластерного анализа. Обычно на листьях дендрограммы написаны имена переменных — названия объектов. Рассказываю, как вместо названий показать изображения.

В конце заметки — финальный код на R ↓

Эта красочная дендрограмма содержит 38 листьев и разделена на 19 классов. На листьях дендрограммы написаны названия видов плоских паразитических червей из класса моногеней. Собираюсь показать вместо них изображения прикрепительных органов этих паразитов

Для работы понадобятся

Знание языка Эр.

Установленный пакет ggtree.

Создаю изображения объектов

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

Тестовые изображения
Тестовые изображения содержат номера, идущие по порядку. Они нужны для того, чтобы было понятно, куда именно разработанный код помещает картинки.

Вот код для генерации 99-ти пнг-файлов размером 20×20 пикселей. Имя файла с изображением совпадает с нарисованым номером.

imgs_path = "img_indexed/"

for (i in 1:99) {
  # Создаю новый пнг-файл img_{i}.png с нулевыми отступами
  png(file = paste0(imgs_path, "img_", i, ".png"), width = 20, height = 20)
  par(mar = c(0, 0, 0, 0))
  
  # Созданию рисунок числа
  plot(0, 0, type = "n", xlim = c(0, 1), ylim = c(0, 1), xlab = "", ylab = "", axes = FALSE)
  text(0.5, 0.5, i, cex = 1)
  
  # Сохраняю пнг-файл
  dev.off()
}

Вот сами изображения.

Вод код для их загрузки.

imgs_indexed = c()
imgs_path = "imgs_indexed/"

for(i in 1:75) {
  imgs_indexed = c(imgs_indexed, paste0(imgs_path, "img_", i, ".png"))
}

Изображения объектов

Изображение объекта — это уменьшенный рисунок прикрепительного крючка паразитического червя (моногенеи). Вот они все:

Общая картинка создана с помощью Имедж-меджик: montage -tile 20x2 -geometry +0+0 *.png all-anchors.png

Название файла с крючком совпадает с наванием вида. Эти названия храню списком в текстовом файле labels.txt. Для загрузки изображений объектов использую следующий код.

imgs_path = "img_objects/"
# Загружаем имена меток
img_objects <- readLines("labels.txt")
# и добавляем расширение '.png'
img_objects <- paste0(imgs_path, img_objects, ".png")

Изображения готовы, перейдем к дереву.

Строю дерево с картинками, используя `ggtree`

В книге «Data integration, manipulation and visualization of phylogenetic trees» описано, как прикрутить картинки к дендрограмме с помощью ggtree.

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

data — это матрица данных: в ее строках — объекты, в столбцах — признаки.

# Выполняю иерархический кластерный анализ
dist <- dist(data)
hc   <- hclust(dist, method = "complete") 

library(dendextend)

# Преобразую результат в дендрограмму
dend <- as.dendrogram(hc)

Из дендрограммы делаю филодерево.

library(ggtree)

phylo <- as.phylo(dend)

Рисую дерево.

ggtree(phylo)

И получаю пустой каркас.

R: ggtree(∙)

Строю дендрограмму с текстовыми листьями

Добавлю к листьям каркаса текстовые подписи при помощи функции geom_tiplab(∙). Чтобы подписи поместились на канве, оставлю пустое место справа при помощи xlim(∙).

ggtree(phylo) + xlim(0, 0.7) + geom_tiplab(geom="label", size=2.2)

R: ggtree(∙) + xlim(∙) + geom_tiplab(label)

Получилось так себе: текст мелкий и бесит рамка вокруг него, но как от нее избавиться — не понял.

Метки мелкие и иногда перекрывают друг на друга. От перекрытия меток избавит пакет ggrepel. Но в результате получается каша.

Рисую картинки на дендрограмме

У функции geom_tiplab(∙) есть параметр geom, задающий тип метки:

one of ’text’, ’label’, ’shadowtext’, ’image’ and ’phylopic’.

При помощи опции ’image’, можно на месте листьев нарисовать картинки.

На нашей дендрограмме 38 листьев. Значит передадим ей массив из 38 тестовых изображений: imgs_indexed[1:38].

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

ggtree(phylo) + geom_tiplab(geom="image", aes(image=imgs_indexed[1:38]), size=.03)

В результате выполнения кода получил ошибку:

Error in `label_geom()`:
! Problem while computing aesthetics.
i Error occurred in the 3rd layer.
Caused by error in `check_aesthetics()`:
! Aesthetics must be either length 1 or the same as the data (75)
x Fix the following mappings: `image`

После длительного анализа определил, что массив изображений должен состоять из 75 элементов, что написано в сообщении об ошибке. (75 = 38 листьев * 2 − 1.) Но код все равно использует только 38 первых изображений. (Это хорошо.)

ggtree(phylo) + geom_tiplab(geom="image", aes(image=imgs_indexed[1:75]), size=.03)

Вуаля.

R: ggtree(∙) + xlim(∙) + geom_tiplab(image)

Подставлю сюда изображения объектов.

ggtree(phylo) + geom_tiplab(geom="image", aes(image=img_objects[1:75]), size=.03)

Цель достигнута.

Финальный код

library(ggplot2)
library(ggtree)

# Формируем массив с названиями файлов с рисунками объектов.
imgs_path = "img_objects/"
img_objects <- readLines("labels.txt")
img_objects <- paste0(imgs_path, img_objects, ".png")

# Выполняю иерархический кластерный анализ
dist <- dist(data)
hc   <- hclust(dist, method = "complete") 

# Преобразую результат в филодерево
phylo <- as.phylo(dend)

# Рисую дерево с картинками
ggtree(phylo) + geom_tiplab(geom="image", aes(image=img_objects[1:75]), size=.03)

Итог

Результат не нравится.

Дендрограмма вверху заметки — красивая, а эта — средненькая. Поэтому в следующий раз покажу, как сделать красивую дендрограмму с картинками при помощи пакета dendextend.

Если после прочтения этой прекрасной заметки вам вдруг непреодолимо захотелось меня поблагодарить, переведите мне 200 руб. на круасан с чаем. А я пока напишу что-то новенькое.

Вырезаем из большой дендрограммы маленькую веточку с помощью R

Антон Лях — Fri, 03 Nov 2023 12:05:56 +0300

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

Нужно оставить только отмеченные ветви

Примечание. Размер шрифта на денрограмме уменьшил при помощи R-кода ↓.

# Размер подписей на оси
par(cex.axis=0.8)
# Размер подписей листьев
dend <- set(dend, "labels_cex", 0.9)
# Рисует дендрограмму
plot(dend)

В R (вроде бы пока еще) нет специальной процедуры для вырезания ветвей. Зато в пакете dendextend есть функция prune.

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

dend_pruned <- prune(dend, c("L_grandis__da", "L_kedahensis__TG__da", "L_chabaudi__da", "L_fenestrum__KM221913__da"))
plot(dend_pruned)

Отмеченные ветви удалены

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

«Пусть потеет машина», гласит дизайнерская мудрость. Поручим компьютеру сделать необходимые перестроения.

Для этого предлагаю функцию invert_prune. Она получает на вход список листьев, которые нужно оставить в дендрограмме, и возвращает список листьев, которые должна обрезать функция prune.

invert_prune <- function(dend, remained_leaves) {
  setdiff( labels(dend), remained_leaves )
}

Пользуемся ею так:

dend_pruned <- prune(dend, 
invert_prune(dend, c("L_grandis__da", "L_kedahensis__TG__da", "L_chabaudi__da", "L_fenestrum__KM221913__da"))
)

plot(dend_pruned)

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

Учите R.

—-

Если вам вдруг непреодолимо захотелось меня поблагодарить, переведите мне 200 руб. на чашеку кофе.

Стандартизация числовых данных

Антон Лях — Thu, 01 Jun 2023 11:27:27 +0300

Это продолжение серии заметок об анализе биологических данных. В прошлый раз мы разобрали понятие объекта и его признаков.

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

Допустим нам нужно статистически сравнить или сгруппировать несколько наблюдений по схожести признаков (например, выполнив кластерный анализ). Каждое наблюдение — это ряд числовых значений, к примеру, биомасс разных видов организмов. Ряды наблюдений упорядочены: первой всегда идет биомасса вида А, затем вида Б и так далее.

У значений есть особенность: в пределах ряда они могут отличаются в несколько раз, потому что мы изучаем как крупные, так и мелкие организмы. В результате биомасса одного или нескольких видов часто в несколько раз превышает остальные. Как на графике.

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

Такие наблюдения сравнивать некорректно. Виды, доминирующие по биомассе, будут сильнее всего влиять на результаты сравнения. Поэтому получится, что мы будем сравнивать не биомассы всех видов, а только нескольких доминирующих (в частном случае — одного).

Теперь выполним стандартизацию по видам (по столбцам).

Данные стандартизованные по признакам (столбцам)

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

Формула стандартизации

Для стандартизации ряда значений (строки или столбца) необходимо вычислить среднее \( \mu \) и стандартное отклонение \( \sigma \) значений ряда, затем отнять среднее от каждого элемента ряда и поделить разницу на стандартное отклонение:

\( y_i = (x_i — \mu) / \sigma \).

После таких манипуляций среднее ряда станет равным нулю, а стандартное отклонение — единице.

Полученные величины \( y_i \) будут выражать значения в пропорции к стандартному отклонению: то есть 1 будет означать одно стандартное отклонение, 2 — два и т. д.

Код на R

Если вы работаете в R, для стандартизации используйте функцию scale(). Она стандартизирует столбцы переданной матрицы.

standardized_columns <- scale(matrix)

Если необходимо стандартизировать строки, транспонируйте t матрицу, а затем транспонируйте результат.

standardized_rows <- t( scale( t(matrix) ) )

Когда применять стандартизацию

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