RのedgeRパッケージを使って発現変動遺伝子を抽出する方法

RのedgeRパッケージを使って、RNA-seqデータでグループ間で発現が変動している遺伝子の抽出方法です。

データセット

GEOデータベースのGSE143213のデータを使用しました。
元論文は
Kared H, Tan SW, Lau MC, Chevrier M et al. Immunological history governs human stem cell memory CD4 heterogeneity via the Wnt signaling pathway. Nat Commun 2020 Feb 10;11(1):821. PMID: 32041953
加齢によるCD4＋T細胞の遺伝子発現レベルの変化を、バルクのRNA-seqとscRNA-seqを使って明らかにしようとしています。ここではバルクのRNA-seqデータであるGSE143213_read_count_bulk.txt.gzをダウンロードして、解凍したものを使います。
データの内訳は、Naive T cell, recent thymic emigrants, memory T cell with naive phenotype, stem-cell memory T lymphocyteの４種類のT細胞で、それぞれ若年５サンプル、老年４サンプル。

edgeR package

発現変動遺伝子の抽出にはRのedgeRパッケージを使いました。参考文献は
Lun, ATL, Chen, Y, and Smyth, GK. It’s DE-licious: a recipe for differential expression analyses of RNA-seq experiments using quasi-likelihood methods in edgeR.　Methods in　Molecular Biology 2016 1418, 391–416.
R version 3.6.3 windows 64-bit版で実行しています

# edgeRのインストール
BiocManager::install("edgeR")
library(edgeR)
# データの読み込み
data_set <- read.table("GSE143213_read_count_bulk.txt", sep = "\t", header = T)
# 列名をensembl gene IDに変更
rownames(data_set) <- data_set$X
# カウント値のみを抽出
counts <- x[,3:38]
# グループの定義
# ラベルはYoungがY, OldがO, NaiveがN, Recent thymic emigrantsがR, T with naive phenotypeがT, Stem-cell memory T lymphocyteがSの組み合わせで、"Y_N", "Y_R", "Y_T", "Y_S", "O_N", "O_R", "O_T", "O_S"とした 
group <- factor(c("Y_N", "Y_R", "Y_T", "Y_S", rep(c("O_N", "O_R", "O_T", "O_S"), 4), rep("Y_N", "Y_R", "Y_T", "Y_S"), 4)))
# design行列の定義
design <- model.matrix(~group)
# データリストを設定
x <- DGEList(counts = counts, group = group)
# dispersion推定
x <- estimateDisp(x, design)
x <- estimateGLMCommonDisp(x, design)
x <- estimateGLMTrendedDisp(x, design)
x <- estimateGLMTagwiseDisp(x, design)
# モデルのフィッティング
fit <- glmQLFit(x, design)
# Young TS とOld TS間でF-Test
qlf.3vs7 <- glmQLFTest(fit, contrast = c(0,0,1, 0, 0, 0, -1, 0))
# 結果の抽出
table <- topTags(qlf.3vs7, nrow(counts))
# FDR < 0.05でDifferentially Expressed Genes (DEG) を抽出
DEG <- subset(table[[1]], table[[1]]$FDR < 0.05)
# 結果のまとめ
DEG_genes <- as.matrix(rownames(DEG))
colnames(DEG_genes) <- "X"
result <- merge(DEG_genes, data_set, by = "X")
table[[1]] <- cbind(table[[1]], X = rownames(table[[1]]))
result <- merge(result, table[[1]], by = "X")
# DEGをプロット
par(cex = 0.7)
plot(result$logFC, result$PValue, pch = 20, xlab = "logFC", ylab = "p-value")
text(result$logFC, result$PValue + 0.0000001, result$gene_name)

結果

データの元論文のSUPPLEMENTARY FIGURE3-Cと比較すると、発現変動遺伝子として抽出された遺伝子が少ないことが解ります。判定が厳しいのかもしれません。

＃参考
データセット
Kared H, Tan SW, Lau MC, Chevrier M et al. Immunological history governs human stem cell memory CD4 heterogeneity via the Wnt signaling pathway. Nat Commun 2020 Feb 10;11(1):821. PMID: 32041953

edgeR
Lun, ATL, Chen, Y, and Smyth, GK. It’s DE-licious: a recipe for differential expression analyses of RNA-seq experiments using quasi-likelihood methods in edgeR. Methods in Molecular Biology 2016 1418, 391–416.