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

Анализ дифференциальной экспрессии генов и функциональная аннотация результатов

Задания

В этом домашнем задании вы будете работать с данными РНК-секвенирования из статьи Li, L., Freitag, J., Asbrand, C. et al. Multi-omics profiling of collagen-induced arthritis mouse model reveals early metabolic dysregulation via SIRT1 axis. Sci Rep 12, 11830 (2022). https://doi.org/10.1038/s41598-022-16005-9.

В этой статье изучали прогрессирование ревматоидного артирита, используя модель коллаген-индуцированного артрита (CIA) у мышей. Секвенирование РНК проводили на нескольких временных точках в контрольной группе (Ctrl) и группе мышей с артритом (CIA).

Задание 0

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

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

Информация о скольки генах и скольки образцах содержится в таблице?

library(tidyverse)

counts <- read_csv("https://ftp.ncbi.nlm.nih.gov/geo/series/GSE193nnn/GSE193335/suppl/GSE193335%5FGEO%5FPorcessed%5FData%5FRaw%5FCount%5FCIA.csv.gz")

counts$geneid <- str_remove(counts$geneid, ".\\d+$")

dim(counts)

Задание 1

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

• временная точка (например, 3 weeks),
• группа мышей с артритом (CIA) или контрольная группа (Ctrl),
• идентификатор мыши.

Создайте датафрейм, в котором эти три элемента будут распределены по трем столбцам. Имена строк в этом датафрейме должны соответствовать исходным названиям образцов.

Задание 2

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

Из матрицы с каунтами удалите те гены, для которых средняя экспрессия по всем отобранным образцам меньше 10.

Задание 3

Создайте DESeq-объект, используя датафрейм и матрицу с каунтами, полученными на предыдущем шаге. В формуле для модели укажите переменную, описывающую наличие/отсутствие заболевания у мыши. Визуализируйте консистентность реплик с помощью графика PCA.

При использовании функции plotPCA() вы можете указать параметр returnData = TRUE, тогда вы получите датафрейм со всеми значениями PC1/PC2, по которым строится график. Этот датафрейм вы можете использовать для создания графика PCA с помощью ggplot2 с нуля. Также вы можете просто модифицировать получившийся после использования функции plotPCA() график с помощью ggplot2-функций.

Задание 4

Проведите анализ дифференциальной экспрессии генов, сравнив группу мышей с артритом против контрольной группы. Сколько дифференциально-экспрессируемых генов нашлось, если установить следующие пороги: поправленное значение p-value (padj) меньше 0.05 и значение log2FoldChange больше 1 или меньше -1?

Задание 5

Для визуализации результатов анализа дифференциальной экспрессии генов постройте график volcano plot. Цветом подсветите дифференциально-экспрессируемые гены. Подпишите символьные названия топ-5 ДЭ генов, повысивших и понизивших экспрессию.

Задание 6

Проведите анализ обогащения GO-категориями (GSEA). Визуализируйте результат с помощью дотплота.

Для визуализации вы можете воспользоваться функцией dotplot() из пакета clusterProfiler, которая на вход принимает объект, получаемый в результате работы функции gseGO() (и аналогичных), либо преобразовать объект с результатами GSEA в датафрейм и построить график самостоятельно с помощью ggplot2. Второй вариант будет удобнее для отрисовки двух графиков рядом, которые бы показывали активированные и подавленные клеточные процессы.

В случае GSEA к “активированным” генным категориям будут относится те категории, у которых значение NES (normalized enrichment score) положительное. Гены, ассоциированные с “активированными” категориями оказались перепредставлены в начале ранжированного списка генов - среди генов с высокими положительными log2FoldChange. В свою очередь, отрицательные значения NES показывают, что категория “подавлена”, то есть гены, входящие в нее, оказались в конце ранжированного списка - среди генов с отрицательными log2FoldChange.