【R】DEseq2を用いた RNA-seq解析データの発現変動遺伝子抽出【RNA-seq】

RNA-seqデータ解析では、健常者と疾患患者、コントロールと薬剤刺激など、様々な発現変動遺伝子（DEG）を取得するかと思います。しかし、多変量となる遺伝子から手作業でDEGを抽出するのは大変です。

ここでは、DEGの抽出方法としてDEseq2ライブラリを用いて行う方法を解説します。

これによりRNA−seqデータから簡単にDEGが抽出できるようになります。ぜひ挑戦してみましょう。

Tested Environment

macOS Monterey (12.4), Quad-Core Intel Core i7, Memory 32GB metafor 4.2-0

公共データを用いたSingle Cell RNA-seq解析に関する初心者向け技術書を販売中

プログラミング初心者でも始められるわかりやすい解説！
RとSeuratで始めるSingle Cell RNA-seq解析！

技術書ページへ

DEseq2とは

DESeq2は、RNAシーケンシング (RNA-seq) データの差異発現解析のためのRパッケージです。RNA-seqは、トランスクリプトーム全体の発現レベルを量定する技術であり、これを用いて異なる条件や時間点での遺伝子の発現の変動を調査することができます。

Gene-level differential expression analysis with DESeq2 より引用

DESeq2の主な特徴

正規化: サンプル間のリードの総数の違いを考慮して、遺伝子の発現量を正規化します。
差異発現解析: 2つ以上の条件間での遺伝子の発現の変動を統計的に評価します。
バリアンスの推定: データ内の遺伝子間でのバリアンスを共有することで、発現量の変動の推定を改善します。
複数検定補正: 多数の遺伝子を同時に評価するため、偽陽性のリスクを減少させるための補正が行われます。

DESeq2による差次的発現解析は、DESeq2はライブラリーの深さの違いを考慮するために正規化因子（サイズ因子）を使用して、生のカウントをモデル化します。次に、遺伝子ごとの分散を推定し、その推定値を縮小して、より正確な分散の推定値を生成し、カウントをモデル化します。最後に、DESeq2は負の二項モデルをフィットし、Wald検定または尤度比検定を用いて仮説検定を行います。

DEseq2の実行

今回のDEseq2で使うファイルは以下になります。それぞれをダウンロードしてデスクトップにおいてください。

counts_data.csv

sample_info.csv

必要なライブラリは以下になります。それぞれ読み込みます。

library(DESeq2)
library(tidyverse)

それぞれのデータを読み込んでください。

setwd("~/Desktop")

counts_data <- read.csv('counts_data.csv')
colData <- read.csv('sample_info.csv')

counts_data, sample_infoがどのようなものか見てみましょう。

counts_dataには行に遺伝子名、列にサンプル名が入っており、それぞれのリードカウント数が記載されたマトリクスファイルになっています。

> head(counts_data)
#                 SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517
# ENSG00000000003        679        448        873        408       1138       1047
# ENSG00000000005          0          0          0          0          0          0
# ENSG00000000419        467        515        621        365        587        799
# ENSG00000000457        260        211        263        164        245        331
# ENSG00000000460         60         55         40         35         78         63
# ENSG00000000938          0          0          2          0          1          0
#                 SRR1039520 SRR1039521
# ENSG00000000003        770        572
# ENSG00000000005          0          0
# ENSG00000000419        417        508
# ENSG00000000457        233        229
# ENSG00000000460         76         60
# ENSG00000000938          0          0

colData は各サンプルに関する実験的な変数や条件などのメタデータを提供します。例えば、あるサンプルが治療群か対照群か、サンプルの取得日時や場所、その他の実験的な条件や変数などの情報を含みます。

> head(colData)
#            cellLine dexamethasone
# SRR1039508   N61311     untreated
# SRR1039509   N61311       treated
# SRR1039512  N052611     untreated
# SRR1039513  N052611       treated
# SRR1039516  N080611     untreated
# SRR1039517  N080611       treated

この colData のサンプルの出力を見ると、サンプルのメタデータが2つの列で構成されていることがわかります：

cellLine: サンプルの細胞株を示す。この列は、異なる細胞株からのサンプルを区別するための識別子として機能します。例えば、N61311 や N052611 などの細胞株が示されています。
dexamethasone: この列は、サンプルがデキサメタゾンという薬物で処理されたかどうかを示す。具体的には、サンプルがデキサメタゾンで「treated」（処理された）か「untreated」（処理されていない）かを示します。

このデータを基に、DESeq2 は design = ~ dexamethasone の指定に従って、デキサメタゾンの処理が遺伝子の発現に与える影響を評価します。

それではDESeqDataSetFromMatrix 関数を使用して、DESeq2解析のためのデータセットオブジェクト（DESeqDataSet）を作成していきます。countDataとcolDataに先程読み込んだcounts_dataとcolDataを指定していきます。

dds <- DESeqDataSetFromMatrix(countData = counts_data,
                              colData = colData,
                              design = ~ dexamethasone)

designはこれは解析のデザイン式を指定します。この式は、データセットの変数（この場合はdexamethasone）を使用して、どの変数を基に統計的なモデルを構築するかを指定します。この例では、dexamethasone（デキサメタゾン処理の有無）に基づいて遺伝子の発現の変動を調べることを意味しています。

ちなみに、~ はなんだろう…?と思った方もいると思います。~はR の統計的モデリング関数で使用されるデザイン式の一部として標準的に使用される記号です。デザイン式は、統計モデルの構造や関係を定義するための簡潔な方法を提供します。

具体的には、~ の左側には応答変数（目的変数や従属変数とも呼ばれる）が、右側には説明変数（独立変数や予測変数とも呼ばれる）が来ることが一般的です。たとえば、y ~ x は、y（応答変数）がx（説明変数）によってどのように影響を受けるかをモデル化することを意味します。

しかし、DESeq2 や他の多くのバイオインフォマティクスのツールでは、応答変数は通常、遺伝子の発現量となるため、デザイン式では右側だけを指定します。この場合、~ は「以下の変数に基づいてモデル化する」という意味になります。

DESeq2を実行する前にプリフィルタリングステップをはさみます。プリフィルタリングステップは、多くの差異発現解析のワークフローで推奨されています。総リード数が10以上のもの（発現量が無いものをカット）に絞っていきます。低いカウントの遺伝子は、ノイズの影響を受けやすく、偽の結果を生む可能性が高まりますので、これらの遺伝子を除外することで、偽陽性のリスクを低減することができます。

keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]

それではDESeq2を実行していきます。終わったら結果を取得して中身を確認していきましょう。

dds <- DESeq(dds)

res <- results(dds)

DEseq2結果の確認

それでは実際にDEseq2結果の確認をしていきます。まずはresを見ていきましょう。

resにはDESeq2を使用して実行された差異発現解析の結果が入っています。各行は遺伝子を、各列は解析結果の異なる側面を示しています。

> res

# log2 fold change (MLE): dexamethasone untreated vs treated 
# Wald test p-value: dexamethasone untreated vs treated 
# DataFrame with 22369 rows and 6 columns
#                  baseMean log2FoldChange     lfcSE      stat    pvalue      padj
#                 <numeric>      <numeric> <numeric> <numeric> <numeric> <numeric>
# ENSG00000000003  708.6022      0.3788474  0.173165  2.187781 0.0286856  0.138465
# ENSG00000000419  520.2979     -0.2037622  0.100756 -2.022329 0.0431424  0.183141
# ENSG00000000457  237.1630     -0.0340364  0.126493 -0.269078 0.7878699  0.929861
# ENSG00000000460   57.9326      0.1171825  0.301588  0.388551 0.6976080  0.894183
# ENSG00000000971 5817.3529     -0.4409590  0.258783 -1.703975 0.0883858  0.297174
# ...                   ...            ...       ...       ...       ...       ...
# ENSG00000273483   2.68957      -0.600436  1.084493 -0.553656 0.5798145        NA
# ENSG00000273485   1.28645      -0.194034  1.346631 -0.144089 0.8854304        NA
# ENSG00000273486  15.45254       0.113349  0.426061  0.266040 0.7902086  0.930772
# ENSG00000273487   8.16323      -1.017777  0.575826 -1.767508 0.0771432  0.271676
# ENSG00000273488   8.58448      -0.218082  0.570756 -0.382094 0.7023919  0.896539

baseMean: 遺伝子の平均カウント値。これは、全てのサンプルを通じての遺伝子の平均発現量を示しています。
log2FoldChange: 対数2倍率の変化。この場合、dexamethasone untreated vs treatedの効果を示す値です。正の値は、treatedサンプルでの遺伝子の発現がuntreatedサンプルに比べて増加していることを示し、負の値は減少していることを示します。
lfcSE: 対数2倍率の変化の標準誤差。これは、推定された対数2倍率の変化の不確実性を示す指標です。
stat: Wald統計量。この値は、遺伝子の発現が治療間で有意に異なるかどうかを判断するための統計的検定に使用される値です。
pvalue: Wald検定のp値。この値は、観察されたデータが、指定されたモデル（この場合は、dexamethasone untreated vs treatedの効果がないという帰無仮説）の下で期待されるものと比べて、どれほど極端かを示します。
padj: 補正後のp値。これは、多重比較の誤りを考慮して補正されたp値です。この値を使用して、遺伝子が有意に差異発現しているかどうかを判断するのが一般的です。

遺伝子が差異発現しているかどうかを判断するためには、通常、補正後のp値（padj）と対数2倍率の変化の絶対値を使用します。

次にresのsummaryを見てみましょう。

> summary(res)

# out of 22369 with nonzero total read count
# adjusted p-value < 0.1
# LFC > 0 (up)       : 1502, 6.7%
# LFC < 0 (down)     : 1884, 8.4%
# outliers [1]       : 51, 0.23%
# low counts [2]     : 3903, 17%
# (mean count < 4)
# [1] see 'cooksCutoff' argument of ?results
# [2] see 'independentFiltering' argument of ?results

summary(res)の出力は、DESeq2によって生成された差異発現解析の結果の概要を示しています。この概要は、データセット内の遺伝子の発現動向や特性を理解するのに役立ちます。

out of 22369 with nonzero total read count:
- このデータセットには、ゼロでない総リードカウントを持つ22,369の遺伝子が含まれています。
adjusted p-value < 0.1:
- この後の行は、補正後のp値が0.1未満の遺伝子に関する情報を提供します。
LFC > 0 (up) : 1502, 6.7%:
- 対数2倍率（LFC）が0より大きい遺伝子は1,502あり、これは全体の6.7%に相当します。これは、これらの遺伝子がアップレギュレートされていることを示しています。
LFC < 0 (down) : 1884, 8.4%:
- 対数2倍率（LFC）が0より小さい遺伝子は1,884あり、これは全体の8.4%に相当します。これは、これらの遺伝子がダウンレギュレートされていることを示しています。
outliers [1] : 51, 0.23%:
- 51の遺伝子（全体の0.23%）が外れ値として識別されています。
low counts [2] : 3903, 17% (mean count < 4):
- 3,903の遺伝子（全体の17%）は、低いカウントを持っており、平均カウントが4未満です。