library(edgeR) #я использую другой протокол, он мне привычнее

urld <- "https://www.ncbi.nlm.nih.gov/geo/download/?format=file&type=rnaseq_counts" 
path <- paste(urld, "acc=GSE183464", "file=GSE183464_raw_counts_GRCh38.p13_NCBI.tsv.gz", sep="&");
GSE183464 <- as.matrix(data.table::fread(path, header=T, colClasses="integer"), rownames=1)
#эти три строки загрузили нам данные
groups <- c(rep("Aneu", 7), rep("Ctrl", 7))
y <- DGEList(counts=GSE183464,group=groups)
keep <- filterByExpr(y)
y <- y[keep,,keep.lib.sizes=FALSE]
y <- normLibSizes(y)
design <- model.matrix(~0+groups)
y <- estimateDisp(y, design, robust=TRUE)
y$common.dispersion
## [1] 0.2129857
# здесь можно посмотреть на распределение значений дисперсии после сжатия (всё нормально, продолжаем анализ)
plotBCV(y)

# здесь можно посмотреть результат снижения размерностей. Первая компонента хорошо разделяет образцы на 2 группы. 
plotMDS(y)

con <- makeContrasts(groupsAneu - groupsCtrl, levels=design)
fit <- glmQLFit(y, design)
qlf<-glmQLFTest(fit,contrast=con)
myy <- rownames(as.data.frame(topTags(qlf,n=nrow(y[["counts"]])))[as.data.frame(topTags(qlf,n=nrow(y[["counts"]])))$FDR < 0.005 & as.data.frame(topTags(qlf,n=nrow(y[["counts"]])))$logFC > 1,]) # Мой способ расчётов более строгий, поэтому пороги у меня мягче, чем в варианте GEO2R.

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

# вот как он начинается:
head(myy)
## [1] "1846"      "401463"    "28420"     "57105"     "101928399" "10018"
# вот сколько в нём генов
length(myy)
## [1] 1024
# Давайте сделаем его вид более "привычным"
library(org.Hs.eg.db)
myy <- na.omit(mapIds(org.Hs.eg.db, myy, "SYMBOL", "ENTREZID"))
# вот как он начинается:
head(myy)
##          1846        401463         28420         57105     101928399 
##       "DUSP4" "BHLHE22-AS1"    "IGHV3-53"     "CYSLTR2"   "LINC01679" 
##         10018 
##     "BCL2L11"
# вот сколько в нём генов
length(myy)
## [1] 1012

Давайте попробуем повторить enrichment

library(clusterProfiler)
Fig_1 <- enrichGO(myy, OrgDb = "org.Hs.eg.db",
keyType       = 'SYMBOL',
ont           = "BP",
pAdjustMethod = "BH",
pvalueCutoff  = 0.05,
qvalueCutoff  = 0.05)

dotplot(Fig_1)

Мы получили график, похожий на тот, что сделали с помощью web-инструмента, но, кажется, в нём много повторов. Попробуем улучшить наш график.

dotplot(simplify(Fig_1))

Какие ещё есть способы визуализации?

cnetplot(simplify(Fig_1))

library(ggplot2)
heatplot(simplify(Fig_1), showCategory = 10)+
theme(
axis.text.x = element_blank(),
axis.text.y = element_text(size = 14)
)

Теперь давайте добавим Background и посмотрим, что изменится. В этот раз нашим контролем будет не геном, а гены, мРНК которых есть в аорте в норме.

background <- na.omit(mapIds(org.Hs.eg.db, rownames(GSE183464[rowSums(cpm(GSE183464)[, c(8:14)]) > 10, ]), "SYMBOL", "ENTREZID")) #так мы получили список для фона
Fig_2 <- enrichGO(myy, OrgDb = "org.Hs.eg.db",
keyType       = 'SYMBOL',
ont           = "BP",
universe      = background, #наш фон
pAdjustMethod = "BH",
pvalueCutoff  = 0.05,
qvalueCutoff  = 0.05)
dotplot(simplify(Fig_2))

cnetplot(simplify(Fig_2))

heatplot(simplify(Fig_2), showCategory = 10)+
theme(
axis.text.x = element_blank(),
axis.text.y = element_text(size = 14)
)

Давайте сравним результаты с результатами GSEA

#Для начала нужно подготовить данные
GSEA_set <- as.data.frame(topTags(qlf,n=nrow(y[["counts"]])))$logFC
names(GSEA_set) <- rownames(as.data.frame(topTags(qlf,n=nrow(y[["counts"]]))))
GSEA_set <- sort(GSEA_set, decreasing = TRUE)
Fig_3 <- gseGO(geneList=GSEA_set,
ont ="BP",
keyType = "ENTREZID",
nPerm = 1000,
minGSSize = 3,
maxGSSize = 800,
pvalueCutoff = 0.05,
verbose = FALSE,
OrgDb = org.Hs.eg.db,
pAdjustMethod = "fdr")
dotplot(simplify(Fig_3)) # это привычная визуализация dotplot

ridgeplot(Fig_3, showCategory = 7) # визуализациия ridgeplot позволяет оценить распределение единицы онтологии в ранжированном ряду.

Существуют и другие онтологии: рассмотрим на примере KEGG.

Fig_4 <-  enrichKEGG(mapIds(org.Hs.eg.db, myy, "ENTREZID", "SYMBOL"), organism     = 'hsa', qvalueCutoff = 0.05, pvalueCutoff = 0.05)
dotplot(Fig_4)

cnetplot(Fig_4)

Что нового мы узнали?

Группа генов у которых есть что-то общее это кластер.

Чтобы узнать больше о кластере генов, можно использовать различные методы обогащения. Это можно делать благодаря разнообразным аннотациям. Кстати, что такое обогащение и аннотации?

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

То же самое касается аннотаций. Подумайте, почему так получилось?

Мы узнали новые виды графиков: dotplot, cnetplot, ridgeplot и enrichplot. Попробуйте поискать картинки с ними в интернете и разобраться с ними на новых, незнакомых примерах.