Три метода дифференциального экспрессионного анализа для секвенирования РНК: limma, EdgeR, DESeq2

Резюме

Приведен подробный протокол методов дифференциального экспрессионного анализа для секвенирования РНК: limma, EdgeR, DESeq2.

Аннотация

Секвенирование РНК (RNA-seq) является одной из наиболее широко используемых технологий в транскриптомике, поскольку оно может выявить связь между генетическим изменением и сложными биологическими процессами и имеет большое значение в диагностике, прогностике и терапии опухолей. Дифференциальный анализ данных RNA-seq имеет решающее значение для выявления аберрантных транскрипций, а limma, EdgeR и DESeq2 являются эффективными инструментами для дифференциального анализа. Однако дифференциальный анализ RNA-seq требует определенных навыков владения языком R и умения выбирать соответствующий метод, чего не хватает в учебной программе медицинского образования.

Здесь мы предоставляем подробный протокол для идентификации дифференциально экспрессированных генов (DEG) между холангиокарциномой (CHOL) и нормальными тканями через limma, DESeq2 и EdgeR, соответственно, и результаты показаны на графиках вулканов и диаграммах Венна. Три протокола limma, DESeq2 и EdgeR похожи, но имеют разные этапы среди процессов анализа. Например, линейная модель используется для статистики в лимме, в то время как отрицательное биномиальное распределение используется в edgeR и DESeq2. Кроме того, нормализованные данные о количестве РНК-seq необходимы для EdgeR и limma, но не нужны для DESeq2.

Здесь мы предоставляем подробный протокол для трех методов дифференциального анализа: limma, EdgeR и DESeq2. Результаты применения трех методов частично перекрываются. Все три метода имеют свои преимущества, и выбор метода зависит только от данных.

Введение

РНК-секвенирование (RNA-seq) является одной из наиболее широко используемых технологий в транскриптомике со многими преимуществами (например, высокой воспроизводимостью данных) и значительно расширило наше понимание функций и динамики сложных биологических процессов^1,2. Идентификация аберратных транскриптов в различных биологических контекстах, которые также известны как дифференциально экспрессированные гены (ДЭГ), является ключевым шагом в анализе РНК-seq. RNA-seq позволяет получить глубокое понимание молекулярных механизмов и биологических функций, связанных с патогенезом. Поэтому дифференциальный анализ был расценен как ценный для диагностики, прогностики и терапии опухолей^3,4,5. В настоящее время для дифференциального экспресс-анализа РНК-seq разработано больше пакетов R/Bioconductor с открытым исходным кодом, в частности limma, DESeq2 и EdgeR^1,6,7. Однако дифференциальный анализ требует определенных навыков владения языком R и умения выбирать подходящий метод, чего не хватает в учебной программе медицинского образования.

В этом протоколе, основанном на данных о количестве РНК-seq холангиокарциномы (CHOL), извлеченных из Атласа генома рака (TCGA), три наиболее известных метода (limma^8,EdgeR⁹ и DESeq2¹⁰⁾были проведены, соответственно, программой R¹¹ для идентификации DEG между CHOL и нормальными тканями. Три протокола limma, EdgeR и DESeq2 похожи, но имеют разные этапы среди процессов анализа. Например, нормализованные данные о количестве РНК-seq необходимы для EdgeR и limma^8,9, тогда как DESeq2 использует свои собственные библиотечные расхождения для исправления данных вместо нормализации^10. Кроме того, edgeR специально подходит для данных RNA-seq, в то время как limma используется для микрочипов и RNA-seq. Линейная модель принята Limma для оценки DEG^12,в то время как статистика в edgeR основана на отрицательных биномиальных распределениях, включая эмпирическую оценку Байеса, точные тесты, обобщенные линейные модели и квазивероятностные тесты^9.

Таким образом, мы предоставляем подробные протоколы дифференциального экспрессивного анализа RNA-seq с использованием limma, DESeq2 и EdgeR соответственно. Ссылаясь на эту статью, пользователи могут легко выполнить дифференциальный анализ RNA-seq и выбрать подходящие методы дифференциального анализа для своих данных.

протокол

ПРИМЕЧАНИЕ: Откройте программу R-studio и загрузите R файл "DEGs.R", файл можно получить из Дополнительных файлов/Скриптов.

1. Загрузка и предварительная обработка данных

1. Загрузите данные о количестве высокопроизводительного секвенирования (HTSeq) холангиокарциномы (CHOL) из Атласа генома рака (TCGA). Этот шаг может быть легко достигнут с помощью следующего кода R.
   1. Нажмите кнопку Выполнить, чтобы установить пакеты R.
   2. Нажмите кнопку Выполнить, чтобы загрузить пакеты R.
      if(!requireNamespace("BiocManager", quietly=TRUE))
      + install.packages("BiocManager")
      BiocManager::install(c("TCGAbiolinks", "SummarizedExperiment"))
   3. Задайте рабочий каталог.
      библиотека (TCGAbiolinks)
      библиотека(Обобщенныйопыт)
      setwd("C:/Пользователи/LIUSHIYI/Рабочий стол")
   4. Выберите тип рака.
      < рака - "TCGA-CHOL"
   5. Запустите код R из файла "GDCquery.R", чтобы загрузить данные. Файл "GDCquery.R" можно получить из Дополнительных файлов/Скриптов:
      source("Дополнительные файлы/Скрипты/GDCquery.R")
      голова (cnt)
      ##TCGA-3X-AAVA-01A-11R-A41I-07
      no #ENSG00000000003 4262
      No #ENSG00000000005 1
      no #ENSG00000000419 1254
      No #ENSG00000000457 699
      No #ENSG00000000460 239
      No #ENSG00000000938 334
      ПРИМЕЧАНИЕ: После выполнения данные подсчета CHOLHTSeq будут загружены и названы "cnt", где строки представляют идентификаторы генов ансамбля, а столбцы представляют идентификаторы образцов. Пожалуйста, обратите внимание на цифры на позициях 14-15 в идентификаторах образцов; числа в диапазоне от 01 до 09 указывают на опухоли и в диапазоне от 10 до 19 указывают на нормальные ткани.
2. Преобразуйте идентификаторы генов ансамбля в символы генов.
   1. Импортируйте файл аннотации в R в соответствии с его путем к хранилищу. Файл аннотации (gencode.v22.annotation.gtf) можно получить из дополнительных файлов.
      gtf_v22 <- rtracklayer::import('Дополнительные файлы/gencode.v22.annotation.gtf')
   2. Запустите код R из gtf_v22. R", который можно получить из Дополнительных файлов/Скриптов:
      source("Дополнительные файлы/Скрипты/gtf_v22. R")
   3. Примените функцию "ann" для преобразования идентификаторов генов ансамбля в символы генов.
      cnt=ann(cnt;gtf_v22)
3. Фильтрация низко экспрессированных генов
   1. Нажмите кнопку Выполнить, чтобы установить пакет R "edgeR".
      BiocManager::install("edgeR")
   2. Нажмите кнопку Выполнить, чтобы загрузить пакет R "edgeR".
      библиотека(edgeR)
   3. Выполните следующий код R, чтобы сохранить гены со значениями количества на миллион (CPM) больше, чем один из по крайней мере двух выборок.
      сохранить <- rowSums(cpm(cnt)>1)>=2
      cnt <- as.matrix(cnt[keep,])
      ПРИМЕЧАНИЕ: Значение счетчиков на миллион (CPM) используется вместо счетчика считывания для устранения отклонения, вызванного различными глубинами последовательности.

2. Дифференциальный экспрессионный анализ через "limma"

1. Нажмите кнопку Выполнить, чтобы установить пакет R "limma".
   BiocManager::install("limma")
2. Нажмите кнопку Выполнить, чтобы загрузить пакеты R "limma", "edgeR".
   библиотека(лимма)
   библиотека(edgeR)
3. Выполните следующий код R, чтобы создать матрицу проектирования.
   группа <- substring(colnames(cnt),14,15) # Extract group information
   группа [группа %in% "01"] <- "Cancer" # set '01' as tumor tissue
   группа [группа %in% "11"] <- "Normal" # set '11' as normal tissue
   группа <- factor (group, levels = c("Normal","Cancer"))
   1. Создайте матрицу проектирования.
      проектирование <- model.matrix (~group)
      rownames(design) <- colnames(cnt)
   2. Создайте объект DGEList.
      dge <- DGEList(counts = cnt, group = group)
   3. Нормализуйте данные.
      dge <- calcNormFactors(dge, method = "TMM")
   4. Выполните следующий код R для выполнения анализа дифференциальных выражений на основе метода limma-trend.
      дгэ
      ##An объект класса "DGEList"
      ##$counts
      ##TCGA-3X-AAVA-01A-11R-A41I-07
      no #TSPAN6 4262
      #DPM1 1254
      No #SCYL3 699
      No #C1orf112 239
      No #FGR 334
   5. Рассчитайте значение CPM.
      logdge <- cpm(dge, log=TRUE, prior.count=3)
   6. Нажмите кнопку Выполнить, чтобы соответствовать линейной модели для прогнозирования данных или вывода взаимосвязи между переменными.
      подходит <- lmFit (бревно, дизайн)
   7. Рассчитайте значение T, значение F и логарифмы коэффициентов на основе байесовского значения.
      подходит <- eBayes(fit, trend=TRUE)
   8. Извлеките таблицу результатов.
      res_limma<- as.data.frame(topTable(fit;n=Inf))
      
      голова(res_limma)
      ## logFC AveExpr t P.Value adj. П.Вал Б
      ##RP11-252E2.2 -4.899493 -2.488589 -20.88052 2.386656e-25 4.931786e-21 47.28823
      ##BX842568.1 -4.347930 -2.595205 -20.14532 1.082759e-24 1.118706e-20 45.83656
      ##CTC-537E7.3 -5.154894 -2.143292 -19.59571 3.452354e-24 2.216114e-20 44.72001
      ##RP11-468N14.3 -6.532259 -2.029714 -19.49409 4.289807e-24 2.216114e-20 44.51056
      ##AP006216.5 -4.507051 -2.670915 -19.25649 7.153356e-24 2.956339e-20 44.01704
      ##RP11-669E14.4 -4.107204 -2.828311 -18.93246 1.448209e-23 4.987633e-20 43.33543
      #The результат дифференциального экспрессионного анализа сохраняется в «res_limma», который включает в себя идентификатор гена, значение log2-кратного изменения (logFC), средний уровень экспрессии log2 гена в эксперименте (AveExpr), модифицированную статистику t (t), значение relavent p (P.Value), ложное значение обнаружения (FDR), скорректированное значение p (adj. P.Val) и логарифмических шансов дифференциально экспрессированных генов (B)
      ПРИМЕЧАНИЕ: Функция "calcNormFactors()" "edgeR" использовалась для нормализации данных с целью устранения влияния, вызванного подготовкой образцов или построением библиотеки и секвенированием. При построении проектной матрицы необходимо сопоставить экспериментальный проект (например, тип ткани: нормальные или опухолевые ткани) с образцом идентификаторов матрицы. limma-trend подходит для данных, глубина секвенирования которых одинакова, в то время как limma-voom подходит: (i) когда размер библиотеки выборки отличается; ii) данные, не нормализованные ТММ; iii) в данных много "шума". Положительный logFC означает, что ген регулируется в эксперименте, в то время как отрицательное число означает, что ген регулируется понижением.
   9. Определите DEG.
      res_limma$sig <- as.factor(
      ifelse(res_limma$adj. P.Val < 0.05 & abs(res_limma$logFC) > 2,
      ifelse(res_limma$logFC > 2 ,'up','down'),'not')) # Значение adj.p < 0.05 и |log2FC| >= 2 — пороговые значения для идентификации ДЭГ
      резюме(res_limma$sig)
      ##down не вверх
      ##1880 ​17341 1443
   10. Выведем таблицу результатов в файл.
       запись.csv(res_limma, файл = 'result_limma.csv')
   11. Нажмите кнопку Выполнить, чтобы установить пакет R "ggplot2".
       install.packages("ggplot2")
   12. Нажмите кнопку Выполнить, чтобы загрузить пакет R "ggplot2".
       библиотека(ggplot2)
   13. Запустите код R из «вулкана. R" для создания сюжета вулкана. Файл «Вулкан. R" можно получить из дополнительных файлов.
       source("Дополнительные файлы/Скрипты/вулкан. R")
       вулкан(res_limma,"logFC","adj. П.Вал",2,0.05)
       ПРИМЕЧАНИЕ: Гены могут быть отображены в различные позиции в соответствии с их значениями log2FC и adj-p, регулируемые вверх DEG окрашены в красный цвет, а deG с пониженной регулией окрашены в зеленый цвет.
   14. Нажмите «Экспорт», чтобы сохранить график вулкана.
       ПРИМЕЧАНИЕ: Графики вулканов могут быть сгенерированы и загружены в различных форматах (например, pdf, TIFF, PNG, JPEG формат). Гены могут быть отображены в различные позиции в соответствии с их значениями log2FC и adj p, регулируемые deG (log2FC > 2, adj p < 0,05) окрашены в красный цвет, а deG с пониженной регулизией (log2FC < -2, adj p < 0,05) окрашены в зеленый, не-DEG окрашены в серый.

3. Дифференциальный анализ выражений через "edgeR"

1. Нажмите кнопку Выполнить, чтобы загрузить пакет R "edgeR".
   библиотека(edgeR)
2. Выполните следующий код R, чтобы создать матрицу проектирования.
   группа <-подстрока(colnames(cnt),14,15)
   группа [группа %in% "01"] <- "Рак"
   группа [группа %in% "11"] <- "Нормальный"
   group=factor(группа, уровни = c("Нормальный","Рак"))
   проектирование <-model.matrix(~group)
   rownames(design) = colnames(cnt)
3. Нажмите кнопку Выполнить, чтобы создать объект DGEList.
   dge <- DGEList(counts=cnt)
4. Нормализуйте данные.
   dge <- calcNormFactors(dge, method = "TMM")
5. Нажмите кнопку Выполнить, чтобы оценить дисперсию значений экспрессии генов.
   dge <- estimateDisp(dge, design, robust = T)
6. Нажмите кнопку Выполнить, чтобы подогнать модель для подсчета данных.
   подходит <- glmQLFit(dge, дизайн)
7. Проведите статистический тест.
   подходит <- glmQLFTest(fit)
8. Извлеките таблицу результатов. Результат сохраняется в "res_edgeR", который включает в себя значение изменения сгиба журнала, CPM журнала, значение F, p и скорректированное значение fDR p.
   res_edgeR=as.data.frame(topTags(fit, n=Inf))
   голова(res_edgeR)
   ## logFC logCPM F PValue FDR
   ##GCDH -3.299633 5.802700 458.5991 1.441773e-25 2.979280e-21
   ##MSMO1 -3.761400 7.521111 407.0416 1.730539e-24 1.787993e-20R
   ##CL1 -3.829504 5.319641 376.5043 8.652474e-24 5.516791e-20
   ##ADI1 -3.533664 8.211281 372.6671 1.067904e-23 5.516791e-20
   ##KCNN2 -5.583794 3.504017 358.6525 2.342106e-23 9.679455e-20
   ##GLUD1 -3.287447 8.738080 350.0344 3.848408e-23 1.194406e-19
   #The результат сохраняется в "res_edgeR", который включает в себя значение изменения сгиба журнала (logFC), CPM журнала, F, значение p и скорректированное значение FDR p
9. Определите DEG.
   res_edgeR$sig = as.factor(
   ifelse(res_edgeR$FDR < 0.05 & abs(res_edgeR$logFC) > 2,
   ifelse(res_edgeR$logFC > 2 , 'вверх', 'вниз'),'не'))
   резюме(res_edgeR$sig)
   ##down не вверх
   ##1578 15965 3121
10. Выведем таблицу результатов в файл.
    запись.csv(res_edgeR, файл = 'res_edgeR.csv')
11. Создайте сюжет вулкана.
    вулкан(res_edgeR,"logFC","FDR",2,0.05)
12. Нажмите «Экспорт», чтобы сохранить график вулкана.

4. Дифференциальный анализ выражений с помощью "DESeq2"

1. Нажмите кнопку Выполнить, чтобы установить пакеты R "DESeq2".
   BiocManager::install("DESeq2")
2. Нажмите кнопку Выполнить, чтобы загрузить пакеты R "DESeq2".
   библиотека(DESeq2)
3. Выполните следующий код R, чтобы определить фактор группировки.
   группа <-подстрока(colnames(cnt),14,15)
   группа [группа %in% "01"] <- "Рак"
   группа [группа %in% "11"] <- "Нормальный"
   group=factor(группа, уровни = c("Нормальный","Рак"))
4. Создайте объект DESeqDataSet.
   dds <-DESeqDataSetFromMatrix (cnt, DataFrame(group), design = ~group)
   дд
   ##class: DESeqDataSet
   ##dim: 20664 45
   ##metadata(1): версия
   ##assays(1): количество
   ##rownames(20664): TSPAN6 DPM1 ... РП11-274Б21.13 ЛИНК01144
   ##rowData имена(0):
   ##colnames(45): TCGA-3X-AAVA-01A-11R-A41I-07 ...
   ##colData имена(1): группа
5. Выполните анализ.
   dds <- DESeq(dds)
6. Создайте таблицу результатов.
   res_DESeq2 <- data.frame(results(dds))
   
   голова(res_DESeq2)
   ## baseМеанский log2FoldChange lfcSE stat pvalue padj
   ##TSPAN6 4704.9243 -0.8204515 0.3371667 -2.433370 1.495899e-02 2.760180e-02
   ##DPM1 1205.9087 -0.3692497 0.1202418 -3.070894 2.134191e-03 4.838281e-03
   ##SCYL3 954.9772 0.2652530 0.2476441 1.071106 2.841218e-01 3.629059e-01
   ##C1orf112 277.7756 0.7536911 0.2518929 2.992109 2.770575e-03 6.101584e-03
   ##FGR 345.8789 -0.6423198 0.3712729 -1.730047 8.362180e-02 1.266833e-01
   ##CFH 27982.3546 -3.8761382 0.5473363 -7.081823 1.422708e-12 1.673241e-11
   ПРИМЕЧАНИЕ: Результат сохраняется в "res_DESeq2", который включает в себя среднее значение нормализованного счетчика чтения (baseMean), значение изменения сгиба журнала (log2FoldChange), стандартную ошибку изменения сгиба журнала (lfcSE), статистику Wald (stat), исходное значение p (pvalue) и исправленное значение p (padj)
7. Идентификация DEG.
   res_DESeq2$sig = as.factor(
   ifelse(res_DESeq2$padj < 0.05 & abs(res_DESeq2$log2FoldChange) > 2,
   ifelse(res_DESeq2$log2FoldChange > 2 ,'вверх', 'вниз'),'не'))
   резюме(res_DESeq2$sig)
   ##down не вверх
   ##1616 16110 2938
8. Выведем таблицу результатов в файл.
   запись.csv(res_DESeq2, файл = 'res_DESeq2.csv')
9. Создайте сюжет вулкана.
   вулкан(res_DESeq2,"log2FoldChange","padj",2,0.05)
10. Нажмите «Экспорт», чтобы сохранить график вулкана.

5. Диаграмма Венна

1. Нажмите кнопку Выполнить, чтобы установить пакет R "VennDiagram".
   install.packages("VennDiagram")
2. Нажмите кнопку Выполнить, чтобы загрузить пакет R "VennDiagram".
   библиотека (ВеннДиаграмма)
3. Составьте диаграмму Венна из регулируемых DEG.
   grid.newpage()
   grid.draw(venn.diagram(list(Limma=rownames(res_
   limma[res_limma$sig=="up",]),
   edgeR=имена строк(res_edgeR[res_edgeR$sig=="up",]),
   DESeq2=rownames(res_DESeq2[res_DESeq2$sig==
   "вверх",])),
   NULL,высота = 3,ширина = 3,единицы = "in",
   col="черный",lwd=0.3,fill=c("#FF6666","#FFFF00",
   "#993366"),
   alpha=c(0.5, 0.5, 0.5),main = "Up-регулируемые DEG"))
4. Нажмите кнопку Экспорт, чтобы сохранить диаграмму Венна.
5. Составьте диаграмму Венна с пониженными регулируемыми DEG.
   grid.newpage()
   grid.draw(venn.diagram(list(Limma=rownames(res_
   limma[res_limma$sig=="down",]),
   edgeR=имена строк(res_edgeR[res_edgeR$sig==
   "вниз",]),
   DESeq2=rownames(res_DESeq2[res_DESeq2$sig=="down",])),
   NULL,высота = 3,ширина = 3,единицы = "in",
   col="черный",lwd=0.3,fill=c("#FF6666","#FFFF00",
   "#993366"),
   alpha=c(0.5, 0.5, 0.5),main = "Пониженные регулируемые DEG"))
6. Нажмите кнопку Экспорт, чтобы сохранить диаграмму Венна.

Результаты

Существуют различные подходы к визуализации результата дифференциального экспрессионного анализа, среди которых особенно используются график вулкана и диаграмма Венна. Лимма идентифицировала 3323 ДЭГ между CHOL и нормальными тканями с |logFC|≥2 и adj. P.Val <0,05 в качестве пороговых значений, ср?...

Обсуждение

Обильные аберратные транскрипты при раке могут быть легко идентифицированы с помощью дифференциального анализа RNA-seq^5. Однако применение дифференциального экспрессивного анализа RNA-seq часто ограничено, поскольку оно требует определенных навыков владения языком R и способ?...

Материалы

Name	Company	Catalog Number	Comments
R		version 3.6.2	free software
Rstudio			free software