EBSeqの使い方（三群以上の多群比較） - アメリエフの技術ブログ

こんにちは。
本日はバイオインフォマティクスの中でもメジャーなRNA-seqデータの発現比較解析において、ちょっとユニークな解析手法をご紹介します。

EBSeqはRNA-seqデータの遺伝子またはアイソフォームの発現比較解析のためのRパッケージです。
経験ベイズモデルを利用して発現変動遺伝子（DEG）を検出します。３群以上の多群比較におけるDEGの検出にも使用できます。

ベイズモデルを利用しているため、

とくに多群比較でどの群間に差があるのかを1回の計算で判定できるのが長所
サンプル数が大きいとDESeqなどよりDEGの精度が良いという報告もある
一方、p-valueといったお馴染みの値が出てこないのが短所でしょうか。

この記事では結果のイメージや使い方を説明します。

多群比較の結果イメージ
インストール
実行
参考

多群比較の結果イメージ

多群比較の場合、各遺伝子の平均発現レベルが群間で差があるか・そして差がある群のパターンを分類します。
群間のFold Change(FC) や、群ごとの平均値も求めることができます。

「差がある群のパターン」については棒グラフでイメージをご覧ください。

例1. 平均が同じだろう（three conditions have the same latent mean expression level）

例2. A群とB&C群に差があるだろう。BとCは平均が同じだろう。
※ この場合、A、B、Cを入れ替えたパターンで計3パターンあります。

例3. 全ての群間で差があるだろう。

インストール

Bioconductor からインストールします。

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("EBSeq")

使用したバージョンは次のとおりです。

R 4.1.2
Bioconductor version 3.14
EBSeq 1.7.0

実行

EBSeq 公式ドキュメントに従い、３群以上の遺伝子発現変動解析を行います。

テストデータMultiGeneMat は 500遺伝子 × 6サンプルの発現マトリクスです。
サンプルは"C1", "C2", "C3" の3群、各 n=2です。

> data(MultiGeneMat)
  
> class(MultiGeneMat)
[1] "matrix" "array" 
> dim(MultiGeneMat)
[1] 500   6
> head(MultiGeneMat)
        [,1] [,2] [,3] [,4] [,5] [,6]
Gene_1   411  588  395  425 1515 1585
Gene_3   268  240  240  292  956 1049
Gene_5   768  860  698  820 2538 2077
Gene_7  1853 1841 2077 1912  478  500
Gene_9   878  634  853  679  253  204
Gene_11  844  791 1122  724  271  235

サンプルの群名を設定し、パターンを計算します。

Conditions=c("C1","C1","C2","C2","C3","C3")  
PosParti=GetPatterns(Conditions)

パターンの中身は、異なる符号間に差があることを示しています。前述の「差がある群のパターン」をイメージしてください。
3群だと全部で5パターンになります。 ※ 発現の上下はここでは考慮しないので、あとでFold Changeから求めましょう。

> PosParti
         C1 C2 C3
Pattern1  1  1  1
Pattern2  1  1  2
Pattern3  1  2  1
Pattern4  1  2  2
Pattern5  1  2  3

サンプル間の正規化（のためのスケーリング値の作成）行います。
正規化にはmedian-by-ratio normalization法 *1 を使います。

MultiSize = MedianNorm(MultiGeneMat)

formal test を行い、次に各遺伝子がどの「差がある群のパターン」に属するかの事後確率を求めます。
平たく言うと、どの「差がある群のパターン」の可能性が最も高いかを調べることです。

MultiOut = EBMultiTest(MultiGeneMat, NgVector=NULL,Conditions= Conditions, AllParti=PosParti, sizeFactors=MultiSize, maxround=5)
MultiPP = GetMultiPP(MultiOut)

⇒ 結果は list になっていて、 PP 列には事後確率（posterior probability）、 MAPには最もありそうなパターンが入っています。
Patterns は、上で作った PosParti と同じです。
遺伝子を１つ１つ見てみると、確かにPPが最も大きいPattern に判定されていますね。

> names(MultiPP)
[1] "PP"       "MAP"      "Patterns"

> MultiPP$PP[1:5,]
            Pattern1  Pattern2      Pattern3      Pattern4   Pattern5
Gene_1  8.868983e-94 0.4178697  1.998138e-55  5.377997e-72 0.58213030
Gene_3 9.641153e-164 0.9717510 3.065620e-116 6.723412e-109 0.02824904
Gene_5  6.226920e-26 0.9382824  1.183473e-18  6.413616e-20 0.06171764
Gene_7  0.000000e+00 0.5762613 1.335965e-294  0.000000e+00 0.42373874
Gene_9  4.999072e-16 0.9477071  6.553113e-16  2.030393e-15 0.05229290

> MultiPP$MAP[1:5]
    Gene_1     Gene_3     Gene_5     Gene_7     Gene_9 
"Pattern5" "Pattern2" "Pattern2" "Pattern2" "Pattern2" 

> MultiPP$Patterns
         C1 C2 C3
Pattern1  1  1  1
Pattern2  1  1  2
Pattern3  1  2  1
Pattern4  1  2  2
Pattern5  1  2  3

集計してみると、Pattern1 すなわち"群間で平均が同じだろう"（いわゆる、有意差がない）遺伝子の個数が最も多く、一般的にイメージされる結果に近いと思われます。

> table(MultiPP$MAP)

Pattern1 Pattern2 Pattern3 Pattern4 Pattern5 
     323       56        3       53       65

例えば、Pattern5 すなわち "全ての群間で差がある" 遺伝子のリストを得るには次のようにします。

> Genes = rownames(MultiGeneMat)
> Genes[MultiPP$MAP == "Pattern5"]
 [1] "Gene_1"   "Gene_29"  "Gene_31"  "Gene_69"  "Gene_91"  "Gene_105"
 [7] "Gene_123" "Gene_125" "Gene_165" "Gene_167" "Gene_201" "Gene_203"
...
[61] "Gene_515" "Gene_605" "Gene_789" "Gene_845" "Gene_911"

群間のFold Change(FC) や、平均値を求めるには次の関数を使います。

MultiFC = GetMultiFC(MultiOut)

結果はリストになっています。その中身は以下の通りです。
Fold Changeは、組み合わせの数だけ（3C2 = 3通り）列があります。

FCMat　群間のFold Change。生データから算出したもの。
Log2FCMat　群間のlog2(Fold Change)
PostFCMat　群間のFold Change。正規化を考慮した事後FC (posteriot FC)。
Log2PostFCMat　群間のlog2(Fold Change)
CondMeans　正規化を考慮した、各群または各遺伝子の平均値
ConditionOrder　群名の順番

> str(MultiFC)
List of 6
 $ FCMat         : num [1:500, 1:3] 1.217 0.951 1.069 0.923 0.983 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:500] "Gene_1" "Gene_3" "Gene_5" "Gene_7" ...
  .. ..$ : chr [1:3] "C1OverC2" "C1OverC3" "C2OverC3"
 $ Log2FCMat     : num [1:500, 1:3] 0.2828 -0.0724 0.0969 -0.1151 -0.0251 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:500] "Gene_1" "Gene_3" "Gene_5" "Gene_7" ...
  .. ..$ : chr [1:3] "C1OverC2" "C1OverC3" "C2OverC3"
 $ PostFCMat     : num [1:500, 1:3] 1.216 0.951 1.069 0.923 0.983 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:500] "Gene_1" "Gene_3" "Gene_5" "Gene_7" ...
  .. ..$ : chr [1:3] "C1OverC2" "C1OverC3" "C2OverC3"
 $ Log2PostFCMat : num [1:500, 1:3] 0.2819 -0.072 0.0967 -0.115 -0.0251 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:500] "Gene_1" "Gene_3" "Gene_5" "Gene_7" ...
  .. ..$ : chr [1:3] "C1OverC2" "C1OverC3" "C2OverC3"
 $ CondMeans     : num [1:500, 1:3] 499 253 813 1843 753 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:500] "Gene_1" "Gene_3" "Gene_5" "Gene_7" ...
  .. ..$ : chr [1:3] "C1" "C2" "C3"
 $ ConditionOrder: Named chr [1:3] "C1" "C2" "C3"
  ..- attr(*, "names")= chr [1:3] "Condition1" "Condition2" "Condition3"


# Fold Change  
> head(MultiFC$PostFCMat)
         C1OverC2  C1OverC3  C2OverC3
Gene_1  1.2157601 0.3252506 0.2675286
Gene_3  0.9513238 0.2557711 0.2688581
Gene_5  1.0692988 0.3550779 0.3320661
Gene_7  0.9233929 3.7881693 4.1024459
Gene_9  0.9827845 3.3027632 3.3606179
Gene_11 0.8836688 3.2320519 3.6575378

それでは、よいデータ解析をエンジョイしてくださいませ〜

参考

Leng, N., J.A. Dawson, J.A. Thomson, V. Ruotti, A.I. Rissman, B.M.G. Smits, J.D. Haag, M.N. Gould, R.M. Stewart, and C. Kendziorski. EBSeq: An empirical Bayes hierarchical model for inference in RNA-seq experiments, Bioinformatics, 2013.

Bioconductor - EBSeq

*1:Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11(10):R106. doi: 10.1186/gb-2010-11-10-r106. Epub 2010 Oct 27. PMID: 20979621; PMCID: PMC3218662.