Динамика современных экосистем в голоцене: методы, конференции, литература. 2019 г.

Использование R в палеоэкологических исследованиях

R это программа для статистической обработки и графического представления данных. Но не только! Это ещё и ясный логичный язык, с помощью которого можно общаться с компьютером и фактически безгранично расширять свои возможности. Здесь, мне кажется, уместна такая аналогия: потратив пару часов на изучение нотной грамоты, вы сможете наиграть одним пальцем простенькую мелодию на синтезаторе. Занимаясь по нескольку часов каждый день в течение нескольких лет, станете профессиональным музыкантом. И только от вашего желания и мотивации зависит, в какой точке этого широкого диапазона вы в итоге окажетесь. Также и с R: регулярные занятия приносят обильные плоды.

Для начала необходимо установить программу с сайта http://cran.r-project.org/. Там же можно найти и информацию о дополнительных пакетах программ (packages), которые значительно расширяют возможности пользователя. После установки и запуска программы вы окажетесь один на один с командной строкой, которая позволяет вызывать функции, создавать объекты, загружать данные, словом работать. Для того чтобы облегчить начальный период освоения этой сложной программы, в этом разделе сайта мы рассмотрим несколько примеров применения R в палеоэкологических исследованиях.

В качестве основного русскоязычного руководства по R можно рекомендовать вот эту превосходную книгу:

Шипунов А. Б., Балдин Е. М., Волкова П. А. и др. 2012. Наглядная статистика. Используем R!. М.: ДМК Пресс. 298 с.

Очень много полезного можно найти здесь http://herba.msu.ru/shipunov/software/r/r-ru.htm.

Обработка и визуализация данных изучения отложений

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

1. Подготовка данных

Будем исходить из того, что абсолютное большинство из нас для ввода и редактирования своих данных использует программу Excel. Для дальнейшей работы в R надо создать в Excel вот такую таблицу, в которой в первой колонке записана глубина образца, а в остальных — результаты анализа. Ограничим этот текст примером из спорово-пыльцевого анализа. Тогда в остальных колонках следует записать количество пыльцевых зёрен обнаруженных таксонов.

 

Depth  |  Var.1  |   Var.2 |  и т.д.        

 

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

Важное замечание: Порядок таксонов в таблице задаёт порядок таксонов в спорово-пыльцевой диаграмме.

Для последующей работы в R таблицу Excel следует сохранить как текстовый файл с разделителями табуляции (tab delimited). Удобно поместить этот файл в особую директорию вместе с несколькими другими файлами, необходимыми для работы (см. ниже).

2. Пакет rioja

Для анализа палеоэкологических данных в R, в том числе и построения спорово-пыльцевых диаграмм, был создан пакет rioja. Знакомство с R начнём с установки и инициализации этого пакета. В окне R есть панель с пунктом «Пакеты». Щёлкнув по ней мышкой, выберите пункт «Установить пакет». Программа сначала предложит выбрать зеркало (если вы находитесь в России, выбирайте Russia), а потом выдаст список пакетов в алфавитном порядке (выбирайте rioja). На этом этапе можно убедиться в том, что пакетов для решения разнообразных задач написано уже очень много! Установку пакета надо провести всего один раз, а вот инициализировать его надо будет при каждом новом запуске R. Для этого надо вызвать специальную функцию, т. е. напечатать её имя в командной строке и нажать Enter:

library(rioja)

Важное замечание: Управление R осуществляется функциями. Функции имеют имя, в данном случае library, за которым идут скобки. В скобках указываются аргументы функции, в данном случае имя требуемого пакета. Ниже мы рассмотрим некоторые другие функции.

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

 

3. Загрузка данных

Данные, набранные в Excel, уже сохранены в текстовом файле, хранящемся в директории. Теперь их надо «скормить» R. Для этого удобно сначала указать путь к этой директории:

 

setwd("D://Data")

 

В таком написании функция setwd устанавливает в качестве рабочей директории директорию Data в корневом каталоге диска D. Обратите внимание на особенности синтаксиса! Теперь файлы будут загружаться из этой директории, а вновь созданные сохраняться туда же. Для лучшего понимания текста я предлагаю примеры электронных таблиц, а также файл с полезными функциями в виде архива для скачивания example.zip (см. внизу страницы). Его лучше разархивировать в отдельную директорию. Путь к этой директории и следует указать функцией setwd. Для того чтобы загрузить файл надо напечатать:

data.main<-read.table ("data.txt",dec=",",h=T)

Функция read.table создаёт новый объект из файла data.txt. Стрелка показывает, что результат надо записать под названием data.main (пользователь может придумать любое название, только без пробелов, скобок, слешей и т. д. Понятно, что называть объекты надо коротко, но ясно!). В качестве аргументов этой функции мы указали полное название файла с расширением, а также dec="," и h=T. Первый нужен вот почему. В России десятичные знаки принято отделять запятыми, тогда как R по умолчанию использует точку. Второй указывает, что первую строку таблицы надо воспринимать как имена переменных. Аргументы функции (любой!) разделяются запятыми. В архиве для скачивания под именем data.txt содержится пример электронной таблицы с данными подсчёта пыльцы.

Таким образом, мы загрузили в R нашу электронную таблицу. В R такой тип объекта называется data frame. Для дальнейшей работы надо научиться ориентироваться в ней. Прежде всего, укажем на такую функцию:

 

names(data.main)

 

Она перечисляет имена всех переменных этой таблицы. К каждой переменной в отдельности можно обратиться, например вот так:

 

data.main$Depth

 

После знака $ указано имя переменной, перед — имя таблицы. Кстати, давайте создадим новую переменную, которая нам в дальнейшем понадобится не один раз:

 

Depth<-data.main$Depth

 

Теперь под именем Depth у нас существует отдельный объект. Для работы с таблицей также очень полезны квадратные скобки:

 

data.main[3,5]

 

Эта команда обращается к третьей строке и пятому столбцу таблицы. Можно обратиться и ко всему столбцу:

 

data.main[,1]

 

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

 

data.main[1:10,1:5]

 

Эта команда выдаёт первые десять строк и пять столбцов таблицы.

 

4. Построение диаграммы

Напомню, что в исходном текстовом файле содержатся результаты подсчёта пыльцевых зёрен в слоях отложения. Для начала их надо преобразовать в процентное содержание таксонов, т. е. для каждого слоя вычислить сумму пыльцевых зёрен всех таксонов и разделить количество зёрен данного таксона на величину суммы. Для этого я написал функцию percent, в качестве аргумента которой надо указать имя таблицы. Но сначала надо загрузить эту функцию. Она и несколько других полезных функций находятся в файле pollen.r. Этот файл содержится в архиве для скачивания вместе с примерами электронных таблиц. Чтобы сделать эти функции доступными, печатаем:

 

source("pollen.r")

 

Очевидно, что этот файл должен находиться в рабочей директории. Теперь можно считать проценты:

 

data.p<-percent(data.main[,-1])

 

Обратите внимание, что мы указали не просто имя нашей таблицы, а приделали к нему некое указание в квадратных скобках. Знак - перед числом 1 указывает, что первый столбец надо исключить. В самом деле, в этом столбце содержится информация о глубине слоёв, которую совсем не нужно переводить в проценты! Под именем data.p мы создали новую таблицу (data frame) уже с процентными данными. Возможно, вам захочется её сохранить, только сначала приделать к ней обратно столбец с глубиной:

 

data.p<-data.frame(Depth,data.p)

 

В этой строчке кода мы соединяем в одну таблицу вектор Depth, созданный ранее, и таблицу с процентными данными.

Теперь её можно записать в рабочую директорию:

 

write.table (data.p,"pcent.txt",dec=",",row.names=FALSE)

 

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

Обратимся к построению пыльцевой диаграммы. Для этого в пакете rioja есть функция strat.plot(). Эта функция вызывает новое окно, в котором и располагается рисунок. Получающаяся в результате диаграмма представляет собой довольно сложный рисунок, поэтому при вызове этой функции в скобках надо указать значения довольно большого количества аргументов. Вот примерный вид:

 

strat.plot(data.p[,-1],yvar=data.p$Depth,y.rev=TRUE,plot.line=FALSE,lwd.bar=5,cex.xlabel=0.6,srt.xlabel=60,
y.tks=c(10,20,30,40,50,60,70,80,90),scale.percent=TRUE,yTop=0.7,xRight=0.9)

 

Давайте разберёмся. На первом месте имя таблицы, как и в случае с вычислением процентов без первой колонки, так как там глубина. На втором месте как раз глубина, такое обращение к переменной вам уже знакомо. Аргумент y.rev указывает, какие значения глубины должны располагаться сверху: если выбрано TRUE, то меньшее значение будет вверху, если FALSE, то наоборот. plot.line указывает, надо ли рисовать ломаную кривую (TRUE) или представить данные в виде прямоугольников (FALSE), lwd.bar задаёт ширину прямоугольников. Следующие два аргумента задают параметры написания названий таксонов: cex.xlabel задаёт размер, а srt.xlabel угол наклона надписи. Для y.tks надо указать вектор числовых значений глубины отложения, в которых следует поставить tick marks. Аргумент scale.percent указывает, надо ли нарисовать все таксоны в одном масштабе (TRUE) или нет (FALSE). Последние два аргумента задают область страницы, которую должен занимать рисунок. Для более тонкой настройки к ним можно прибавить yBottom и xLeft. Поэкспериментируйте со значениями аргументов, чтобы увидеть, как меняется вид рисунка. При этом удобно изменять значения в отдельном текстовом файле, из которого копировать всю команду целиком. Расширив окно с диаграммой до полного размера, увидим меню, которое позволит сохранить рисунок как растровый или векторный файл.

Важное замечание: Сложные функции содержат большое число аргументов. При вызове функции не обязательно указывать их все. Значения аргументов, не указанных при вызове функции, будут установлены по умолчанию. Чтобы посмотреть на список всех аргументов, их значения по умолчанию и другие детали использования функции, надо набрать команду такого вида (знак вопроса + имя функции):

 

?strat.plot

 

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

 

5. Обработка данных

В рассмотренном примере мы построили пыльцевую диаграмму, используя процентное содержание, вычисленное от суммы всех таксонов. Однако во многих случаях стоит попробовать посчитать отдельно пыльцу деревьев и травянистых растений, или наземных и водных таксонов и т. д. Для упрощения этой задачи я предлагаю функцию extract из файла pollen.r. Для работы с ней нам понадобится ещё один файл, что-то вроде словаря, который по названию таксона определит его принадлежность к той или иной группе. В архиве для скачивания пример такого файла назван dict.txt. Этот файл должен быть организован следующим образом: в первой колонке названия таксонов, а в последующих — текстовые идентификаторы группы. Например, во второй колонке файла dict.txt всем древесным таксонам присвоен идентификатор AP (arboreal pollen), всем травянистым — NAP, всем споровым — Spore. Таких классификаций в этом словаре может быть несколько. Отмечу, что по разнообразию таксонов словарь может быть гораздо больше данного конкретного отложения, т. е. один и тот же словарь можно использовать для разных случаев. Важно только, чтобы все таксоны из данного разреза были представлены в словаре. Излишне говорить, что написание названий таксонов в основной таблице с данными и в словаре должны строго соответствовать друг другу. Вообще, R не прощает грамматических ошибок и отличает прописные и строчные буквы. Загрузим словарь:

 

dict<-read.table ("dict.txt",h=T)

 

Теперь можно вызвать функцию extract:

 

tree<-extract(data.main,dict,2,"AP")

 

Разберёмся с аргументами. На первом месте имя основной таблицы с результатами подсчётов пыльцевых зёрен. На втором — имя объекта, который был создан функцией read.table из файла со словарём. За ним идёт номер колонки, в которой в словаре содержится информация для классификации, и, наконец, значение идентификатора в кавычках, согласно которому надо провести отбор таксонов. Эта строчка, например, помещает в новую таблицу, которую мы назвали tree, только те таксоны, которые помечены во второй колонке словаря идентификатором AP. Создадим ещё две новые таблицы:

 

herb<-extract(data.main,dict,2,"NAP")

spore<-extract(data.main,dict,2,"Spore")

 

Теперь у нас есть отдельные таблицы с деревьями, травянистыми и споровыми растениями. В каждой из них результаты подсчёта пыльцы. Нам надо их перевести в процентные данные:

 

tree.p<-percent(tree)

herb.p<-percent(herb)
spore.p<-percent(spore)

 

Теперь у нас процентные соотношения таксонов деревьев подсчитаны относительно суммы всех древесных таксонов. То же верно и для травянистых и споровых растений. Объединим всё в одну таблицу и не забудем приделать к ней глубину слоёв:

 

newdata.p<-data.frame(Depth,tree.p,herb.p,spore.p)

 

Вы уже знаете, как сохранить эту таблицу и построить для неё пыльцевую диаграмму! Не забудьте исключить при построении диаграммы первый столбец таблицы. Возможно, что для анализа данных будет интересна ещё и такая информация — каковы процентные соотношения между деревьями, травами и споровыми в целом? Для этого нам надо для каждого слоя подсчитать число пыльцевых зёрен деревьев, трав и споровых, которые содержатся в трёх отдельных таблицах:

 

AP<-rowSums(tree)

NAP<-rowSums(herb)
SPR<-rowSums(spore)

 

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

 

6. Концентрация

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

Чтобы перевести исходные данные подсчётов пыльцевых зёрен (таблица data.main) в данные о концентрации пыльцы, понадобится ещё один дополнительный файл. Если определение концентрации проводилось первым способом, то этот файл должен быть организован следующим образом: первый столбец глубина; второй — конечный объём ёмкости, из которой отбирается материал для изготовления препарата; третий — количество капель, взятых для изготовления препарата; четвёртый — объём одной капли; пятый — начальный объём (или масса) образца. В архиве для скачивания пример такого файла назван drop.txt. Для расчётов концентрации используем функцию conc_drop из файла pollen.r:

 

drop<-read.table ("drop.txt",dec=",",h=T)

data.c<-conc_drop(data.main[,-1],drop)

 

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

Если концентрация определялась с использованием экзотической пыльцы, то дополнительный файл должен включать в себя следующие столбцы: глубина, количество экзотической пыльцы, количество добавленных таблеток, число пыльцевых зёрен в одной таблетке, исходный объём (или масса) образца. Пример такого файла в архиве для скачивания назван spike.txt, а концентрация рассчитывается с помощью функции conc_spike из файла pollen.r:

 

spike<-read.table ("spike.txt",dec=",",h=T)

data.c<-conc_spike(data.main[,-1],spike)

 

В обоих случаях порядок столбцов во вспомогательных файлах должен строго соблюдаться. Данные о концентрации пыльцы теперь содержатся в таблице data.c. Вы можете и сохранить её, прикрепив к ней вектор с глубинами образцов, и построить для неё новую диаграмму. В этом случае при вызове функции strat.plot стоит использовать аргумент scale.percent=FALSE.

 

7. Выделение зон

Теперь у нас есть три таблицы, в которых исходные данные обработаны разными способами. Это, во-первых, data.p, в которой процентные соотношения подсчитаны от суммы всех таксонов. Во-вторых, это newdata.p, в которой проценты подсчитаны для деревьев, травянистых и споровых растений отдельно. И, наконец, это data.c — данные о концентрации пыльцы. Какую из трёх таблиц использовать для выделения спорово-пыльцевых зон решать вам. Я лишь покажу, как это делается в R.

Прежде всего, отметим, что помимо визуального анализа есть и численный метод. Этот метод называется CONISS или Stratigraphically Constrained Cluster Analysis. Для того чтобы провести этот анализ, нам понадобятся ещё две функции:

 

new<-dist(data.p)

coniss.data<-chclust(new,method="coniss")

 

В первой строке мы вызвали встроенную в R функцию dist и назвали созданный ею объект new. По умолчанию эта функция рассчитывает эвклидовы дистанции между образцами, однако можно воспользоваться и другими методами (см. ?dist). Во второй строке объект new был обработан функцией chcluct из пакета rioja. Новый объект coniss.data понадобится нам при построении диаграммы:

 

plot.data<-strat.plot(data.p[,-1],yvar=data.p$Depth,y.rev=TRUE,plot.line=FALSE,lwd.bar=5, cex.xlabel=0.6,srt.xlabel=60,
y.tks=c(10,20,30,40,50,60,70,80,90),scale.percent=TRUE,yTop=0.7,xRight=0.9,clust=coniss.data)

 

Это уже знакомая нам функция strat.plot. Обратите внимание вот на какие моменты. Во-первых, мы вставили новый аргумент clust=coniss.data, т. е. аргумент clust берёт информацию из объекта coniss.data, созданного чуть раньше. Во-вторых, мы не просто вызвали функцию strat.plot, а записали её результат в новый объект plot.data. Это нужно для того, чтобы нарисовать границы пыльцевых зон:

 

addZone(plot.data,20)

 

Эта функция пакета rioja проводит верхнюю границу зоны на глубине 20 см. Обратите внимание — для того, чтобы эта функция работала, окно с рисунком диаграммы должно быть открыто!

 

8. Заключение

Мы рассмотрели примеры работы с данными спорово-пыльцевого анализа отложений. Осталось отметить ещё одно важное обстоятельство. Теперь вы можете не только построить спорово-пыльцевую диаграмму. Теперь ваши данные находятся в R, и вы можете применить к ним все те практически безграничные возможности статистической обработки, которыми обладает эта программа. Только не забывайте, что R – это плод многолетней работы большого коллектива, и не забывайте ссылаться на соответствующие публикации. Посмотреть, как это сделать, можно с помощью команды citation().

Файл для скачивания:

example.zip  – содержит текстовые файлы с примерами электронных таблиц и файл pollen.r

pollen.r – содержит функции, полезные при анализе палеоэкологических данных. Должен быть помещён в рабочую директорию и загружен функцией source("pollen.r")

 

 

 

 

Булат Хасанов

Лаборатория биогеоценологии и исторической экологии

Институт проблем экологии и эволюции РАН

bulatfk@gmail.com