はじめに
前回はRNA-seqデータのダウンロードからQC、そしてSTARによる参照ゲノムへのマッピングとRSEMによる発現定量までを行いました！
はじめに RNA-seq（RNAシーケンシング）は、サンプル中の遺伝子の発現レベルを分析するために使用される方法です！ 近年ではRNA-seqが安価で行えるようになり、非常に盛んに行われる手法になりました！ それに伴ってSeq[…]
今回はedgeRというRのパッケージを使用して発現変動解析（DEG解析）を行なっていきます！
次回は抽出した発現変動遺伝子（DEG）を使用して、パスウェイ解析によってその遺伝子がどのような機能を持っているのかを調べます！
まずは必要なライブラリーをインポートします！
library(edgeR)
library(tidyverse)
library(tximport)
library(plotly)
library(gplots)
library(viridis)
今回発現変動解析はDexamethasone (Dex)を投与したサンプル vs Dexを投与していないサンプルで比較を行います！
使用したデータセットについてはPart1で解説してありますので前回の記事をご参照ください！
はじめに RNA-seq（RNAシーケンシング）は、サンプル中の遺伝子の発現レベルを分析するために使用される方法です！ 近年ではRNA-seqが安価で行えるようになり、非常に盛んに行われる手法になりました！ それに伴ってSeq[…]
使用したサンプルのより詳細を知りたい方はこちらのページをご参考ください！(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE52778 )
まずはこれまでに扱ってきたデータがどのようなグループに分かれているのかと、そのファイルパスがまとまったデータフレーム(sampleinfo)を作成します！
sampleinfo <- tibble(run=c("SRR1039509", "SRR1039513", "SRR1039517", "SRR1039521", "SRR1039508", "SRR1039512", "SRR1039516", "SRR1039520"),
Dex=c("trt", "trt", "trt", "trt", "untrt", "untrt", "untrt", "untrt")) %>%
mutate(group=recode(Dex, trt="case", untrt="control"),
genes_results_path=str_c(run, ".genes.results"),
isoforms_results_path=str_c(run, ".isoforms.results"))
sampleinfoAirwayデータセットはDexだけではなく、albuterolの投与・非投与や細胞種などの情報も含まれていますが、今回はシンプルな2群比較を行うためにDexの投与・非投与のみに注目して解析を行なっています！
tximportによるカウントデータの読み込み
まずはRSEMで発現定量を行った結果である、カウントデータを読み込みます！
これには「tximport」というRのライブラリーを使用します！
「tximport」はRSEMやSalmon、kallistoの出力ファイルを読み込んで、edgeRやDEseq2のインプットとして使用できるように処理してくれるめちゃくちゃ便利なライブラリーです！
つまりは、発現定量するソフトウェア（RSEMやSalmon、kallisto）と発現変動解析をするソフトウェア（edgeRやDEseq2）の架け橋となってくれるということです！
詳細な使用方法はtximportのvignettes(https://bioconductor.org/packages/release/bioc/vignettes/tximport/inst/doc/tximport.html) をご参照ください！
（Salmonなどのファイルの読み込み方,edgeR/DEseq2への渡し方がまとめられていて必読です！）
さて早速RSEMの出力ファイルを読み込んでいきます！
rsem_output_dir = file.path("result", "RSEM")
files <- file.path(rsem_output_dir, sampleinfo$genes_results_path)
names(files) <- paste0("sample", 1:8)
txi <- tximport(files, type = "rsem", txIn = FALSE, txOut = FALSE,
countsFromAbundance = "scaledTPM")
countdata <- txi$counts
head(countdata)
低発現遺伝子のフィルタリング
まずは、あまりにも発現量の低い遺伝子を取り除きます！
なぜ低発現遺伝子を取り除く必要があるかというと、のちに行う発現変動解析の邪魔になってしまうからです。全てのサンプルを通して発現量が非常に少ないものはFold change(FC）が大きくなりやすく、本来誤差であるはずの差でさえも発現変動であると勘違いをしてしまうリスクがあります。また、遺伝子はタンパク質に翻訳されたりして生物学的に意味を持つために、最低限のレベルで発現している必要があります。
上記の理由から低発現遺伝子のフィルタリングは生物学的に意味のある発現変動を取り出すという視点から非常に重要な工程です！
edgeRパッケージで低発現遺伝子のフィルタリングを自動で行ってくれるのがfilterByExpr です！
y <- DGEList(countdata, group=sampleinfo$group)
keep <- filterByExpr(y)
table(keep)
min.count
y <- y[keep, , keep.lib.sizes=FALSE]ほんの数行でフィルタリングができました！
filterByExprについて詳しく説明します！
filterByExpr は具体的にどんなことをしてるの？この論文(https://f1000research.com/articles/5-1438)で記載されている低発現遺伝子のフィルタリングの処理を行っています！
具体的には「n個のサンプルでk以上の発現量を持つ遺伝子を保持する」という戦略で遺伝子のフィルタリングを行っています！
ここでの「n個」は最小のグループサンプルサイズ（そのグループに含まれるサンプルの数、今回は4個）であり、
「k以上の発現量」とはfilterByExpr のmin.count 引数で規定されるカウント値です！（デフォルトは10です！） CPMではなくカウント値でフィルタリングする点にご注意ください。
正規化（Normalization）
次に正規化を行います！
端的にいうと「生物学的なエフェクトとは無関係な外的要因による発現量の変化」を無視できるようにするために「各サンプルの遺伝子発現量の分布を揃える」ことです！
RNA-seqのデータは、サンプル調整やシーケンシングの過程で意図せず遺伝子発現量に影響を与えられることがあります。例えば、サンプルAとサンプルBは実験が行われた場所が違ったために、サンプルの置かれた温度に差があったなどです。
このあと行う発現変動解析は「全てのサンプルで発現量の分布がほぼ同じ」であることを前提にしているため、正規化を行う必要があります！
edgeRではcalcNormFactors という関数が正規化を行ってくれます！実際にやってみましょう！
y <- calcNormFactors(y)
y$samples$norm.factors
# [1] 1.0216003 0.9467672 0.9722918 0.9562999 1.0587144 0.9899555 1.0312374 1.0288023
たったこれだけです！これだけで正規化は終了です！
calcNormFactors は何をやってるの？Trimmed mean of M values (TMM) 正規化という手法で正規化を行っています！
TMM正規化とは、各サンプルのライブラリーサイズによるバイアスを取り除き、サンプル間での遺伝子発現量の比較を正確に行うための正規化手法です！
TMM正規化はハウスキーピング遺伝子などをはじめとして、「ほとんどの遺伝子が発現変動しない」ことを前提とした正規化方法でほとんどのRNA-seqのデータに対して推奨される方法です！
少しだけ詳細を解説すると、発現変動していない遺伝子の発現変動量の対数比（M値）は0であると想定されますが、上記で述べた「外的要因による発現量の変化」があったとすると、このM値が0から遠く離れてしうことがあります！TMM正規化は、ライブラリーサイズに応じた係数（スケーリングファクター）を計算し、それを各サンプルの遺伝子発現量にかけることで、遺伝子間のM値の期待値が0に近づくように調整する方法です！
このライブラリーサイズに応じた係数（スケーリングファクター）がcalcNormFactors の計算してくれるnorm.factors です！
TMM正規化については以下のサイトが非常に参考になります！
https://bi.biopapyrus.jp/rnaseq/analysis/normalizaiton/tmm
発現変動遺伝子（DEG）の抽出
さて！ここからが本題の発現変動遺伝子の抽出です！
デザインマトリックスの定義
まずはデザインマトリックスという実験デザインを表現するための行列を定義します！
design <- model.matrix(~group)
design
デザインマトリックスは、異なる条件下でのサンプルのグループ分けを指定するためのものです！デザインマトリックスの各行は、RNA-Seqデータセット内の各サンプルを表し、各列はそのサンプルがどのグループまたは条件に属するかを示します！具体的には、各列はサンプルの説明変数（因子）に対応し、各因子の異なるレベルに対応する値が含まれます！
上記は二群比較のためのデザインマトリックスなのでDexを投与したグループはInterceptのみが1となり、Dexを投与しなかったたグループはInterceptに加えてsampleinfo$groupconrtolの列が1になっています！
参考のために三郡比較のデザインマトリックスを以下に示します！
example_group <- factor(c("A", "A", "A", "B", "B", "B", "C", "C"))
model.matrix(~example_group)
edgeRによる発現変動遺伝子の検定
さてデザインマトリックスが定義できたところで、早速、発現変動の比較を行っていきましょう！
発現変動遺伝子を抽出する手法はいくつも報告されていますが、edgeRではシンプルな二群比較のみに対応した「exact tests」と、一般化線型モデルを使用して複雑な実験デザインにも対応した「quasi-likelihood F-tests」「likelihood ratio tests」の三つの手法が提供されています！
今回の解析はシンプルな二群比較なので三つのどれを使用しても問題ないのですが、edgeR user guideのQuick startでも使用されている「quasi-likelihood F-tests」を使用していきます！
edgeRは遺伝子の発現量が負の二項分布（Negative Binomial）に従うと仮定して、モデル化を行います！そのため、まずは個々の遺伝子の分散パラメータを推定します！
y <- estimateDisp(y, design)
分布の推定ができたら、一般化線形モデルを適合していって、異なるグループ間での発現量の比較を行います！
fit <- glmQLFit(y,design)
qlf <- glmQLFTest(fit,coef=2)
top.tag <- topTags(qlf, n=Inf)
top.tag
そして最後にtopTags 関数を使って発現変動している遺伝子群を抽出してきて終了です！
topTagsはデフォルトでn=10が指定されていて、10個の発現変動遺伝子しか返してくれないのでn=Inf として全ての遺伝子の結果を返してもらっています！
estimateDisp 、glmQLFit、glmQLFTest が何を行っているかの詳細についてはedgeR user guideをご参照ください！
発現変動解析の結果をcsvファイルとして保存します！
top.tag$table %>%
rownames_to_column("ensemble_id") %>%
write_csv("DE_result.csv")
可視化
さて！発現変動遺伝子（DEG）の抽出ができました！
ここで、取れてきたDEGの全体的な傾向を掴むために幾つかの可視化をしてみましょう！
Volcano plot
発現変化の大きさと統計的有意性を同時に可視化することができるのがボルケーノプロットです！
ボルケーノプロットは、X軸に発現変化の対数値（例えば対照群と実験群の比）を、Y軸に対数変換された調整済みp値（例えばBenjamini-Hochberg法による多重検定補正後の値）をプロットします。
たとえば、大きな発現変化と高い統計的有意性を持つ遺伝子は、プロット上部の右側に配置されます。逆に、小さな発現変化と低い統計的有意性を持つ遺伝子は、プロットの下部に配置されます。
RNA-seq実験では、何千もの遺伝子が同時に解析されるため、膨大な量のデータが得られます。ボルケーノプロットを用いることで、大量の遺伝子の中から統計的に有意な遺伝子を特定することが容易になります！
今回はFDRの基準値を0.01, logFCの基準値を1として、upregurated geneを赤、down-regulated geneを青で表示してみました！
また、 ggplot2で描画した上で、plotlyというパッケージのggplotly関数を使用することでインターラクティブに操作可能なボルケーノプロットを作成しました！遺伝子のプロットにマウスオーバーすることでどんな遺伝子がどのくらいのFCとFDRがあるのかをすぐに確認できます！
FDR.threshold <- 0.01
logFC.threshold <- 1
upregulated.genes.color <- "red"
downregulated.genes.color <- "blue"
ns.genes.color <- "grey"
p <- top.tag$table %>%
rownames_to_column(var = "geneid") %>%
mutate(color = ifelse(FDR >= FDR.threshold , ns.genes.color,
ifelse(logFC >= logFC.threshold, upregulated.genes.color,
ifelse(logFC <= -logFC.threshold, downregulated.genes.color, ns.genes.color)))) %>%
ggplot() +
geom_point(aes(x = logFC, y = -log10(FDR), color = color, text = geneid), size = 1) +
labs(x = "log2 (Fold Change)", y = "-log10 (FDR)") +
scale_color_identity() +
geom_hline(yintercept = -log10(FDR.threshold), size = 0.2, color = "dark green") +
geom_vline(xintercept = -logFC.threshold, size = 0.2, color = "dark green") +
geom_vline(xintercept = logFC.threshold, size = 0.2, color = "dark green") +
theme(legend.position = "none", panel.grid = element_blank())
ggplotly(p)
M-A plot
次は2つのサンプル間の遺伝子発現量の差と平均値の比較を可視化するM-A plotです！
M-A plotの横軸は平均発現量（Average expression）を表し、縦軸は発現変動（logFC）を表します！
edgeRではplotSmear で可視化できます！
deg.ids <- top.tag$table %>%
filter(FDR < 0.05) %>%
rownames()
plotSmear(et, de.tags=deg.ids, cex=0.3)
Heatmap
次は似た発現パターンを持った遺伝子同士やサンプル同士のクラスタリングも同時に可視化することのできるHeatmapです！
RNA-seq解析でヒートマップを書く場合、一般的には正規化された発現量を用いることが推奨されているので、カウント値をlogCPMに変換して描画します！
Heatmapを描画するライブラリはいくつかあるので、その中でも特におすすめの二つ「gplots」のheatmap.2 関数と「Heatplus」のregHeatmap を使用します！
FDR.threshold <- 0.01
logFC.threshold <- 1
deg.ids <- top.tag$table %>%
filter((FDR < FDR.threshold) & (abs(logFC) > logFC.threshold)) %>%
rownames()
lcpm <- cpm(y$counts, log=TRUE)
data.for.heatmap <- lcpm %>%
as.data.frame() %>%
rownames_to_column("ensemble_id") %>%
filter(ensemble_id %in% deg.ids) %>%
column_to_rownames("ensemble_id") %>%
as.matrix()
heatmap.2(data.for.heatmap, col=viridis, trace="none", scale="row")
# BiocManager::install("Heatplus")
library(Heatplus)
regHM <- regHeatmap(data.for.heatmap)
plot(regHM)
Box plot（遺伝子ごとに）
最後に遺伝子別に発現量をプロットするBoxplotです！
ggplotでBoxplotを作成するコードが以下です！
target.gene.ids <- top.tag$table %>%
head(n=3) %>%
rownames()
sample.to.group <- tibble(
sample=paste0("sample", 1:8),
group=c("case", "case", "case", "case", "control", "control", "control", "control")
)
lcpm <- cpm(y$counts, log=TRUE)
data.for.boxplot <- lcpm %>%
as.data.frame() %>%
rownames_to_column("ensembl_id") %>%
filter(ensembl_id %in% target.gene.ids) %>%
pivot_longer(cols = !ensembl_id, values_to = "cpm", names_to = "sample") %>%
left_join(sample.to.group, by = "sample")
ggplot(data.for.boxplot, aes(x = group, y = cpm, fill = group)) +
geom_boxplot(outlier.colour = NA) +
geom_jitter(aes(col = group)) +
facet_wrap(~ensembl_id)
次回予告
次回は抽出した発現変動遺伝子（DEG）を使用して、パスウェイ解析によってその遺伝子がどのような機能を持っているのかを調べます！
参考
https://bioconductor.org/packages/release/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf
https://bioconductor.riken.jp/packages/3.9/workflows/vignettes/RNAseq123/inst/doc/limmaWorkflow.html
https://bi.biopapyrus.jp/rnaseq/analysis/normalizaiton/tmm
https://bi.biopapyrus.jp/rnaseq/analysis/de-analysis/2g-edger.html