GOエンリッチメント解析の概説とR上での実践
概略
GOエンリッチメント解析は、遺伝子全体のうち特定のGOでアノテーションされた遺伝子の割合と、発現変動遺伝子のうち特定のGOでアノテーションされた遺伝子の割合を計算し、その割合を比べることによって、そのGOが発現変動遺伝子の中で有意に多く観測できるかを検定する方法です。
GOとは、遺伝子の生物的プロセス、細胞の構成要素および分子機能に着目して、遺伝子に付けられるアノテーションであり、言わば漢字辞典の部首索引みたいなものです。ある遺伝子がどのような機能を持つか、どのような代謝プロセスにかかわっているかなどの情報を、データベース上でその遺伝子に紐づけしておりその時に使われる用語がGene Ontology(GO)になります。
アノテーションはGOと遺伝子を紐づけることを指し、「遺伝子Aは～というGOでアノテーションされている」というような言い方をします。
GOエンリッチメント解析では変動発現遺伝子の中で特定のGOが有意に多く観測できるかを解析し、その変動発現遺伝子群の性質を見抜くものです。例えるなら、「全校生徒(すべての遺伝子)のうち自転車通学をしている(というGOがアノテーションされている)人を考えるとき、特定地域から通学している(変動発現遺伝子)人に絞ったときその割合は有意に増えるか？」という問題と似ています。もしその特定地域からの自転車通学者が特別多いなら「公共交通機関が発達してないのかな？」というような性質を見抜くヒントになりえます。
数理モデルとその解釈
まずはフィッシャーの正確確率検定に基づいて計算する。計算するものは「発現変動遺伝子(DEGs)」のうち、$${k}$$個の遺伝子が特定のGOアノテーションされる確率です。これはwith GOが$${m}$$個、without GOが$${N-m}$$個の合計$${N}$$個あり、$${n}$$個取り出した時にwith GOが$${k}$$個になる確率に言い換えられるので、超幾何分布の計算とみなすことが出来ます。
分母は$${N}$$個の中から$${n}$$個の遺伝子を選ぶので$${_NC_n}$$、分子はwith GOに対して$${m}$$個の中から$${k}$$個選ぶので$${_mC_k}$$、同様にwithout GOに対して$${N-m}$$個の中から$${n-k}$$個選ぶので$${_{N-m}C_{n-k}}$$。以上より求める確率変数は
$$
P(X=k)=\dfrac{_mC_k\cdot_{N-m}C_{n-k}}{_NC_n}
$$
サイトの方では二項展開の表記で書いてありますが同じことです。
この確率について$${k}$$以上になる事象のp値を求めれば良いことになります。すべてのGOについて行えば有意に増加した発現変動遺伝子を検出できる。（本当は累積分布関数でやればいいんだろうが、一般超幾何関数が登場するため初学者向けではない。$${k}$$が自然数なので一個一個足す方が感覚的に楽。計算もどうせプログラムがやってくれるので問題なし。)
注意点
サイトに記載されたものを引用
GO はツリー構造をしている。ツリーの下の層の GO でアノテーションされた遺伝子は、その上の層の GO でもアノテーションされていることになる。そのため、ツリーの上の層にある GO は検定で有意差が出にくい。また、上の層の GO にあるほど一般的な機能が記述されることが多く、下の層の GO ほど具体的な機能が記述される場合が多い。そのため、エンリッチメント解析を行うにあたって、上の層にある GO を除いてから行うこともある。さらに、GO によって、アノテーションされた遺伝子数が異なる。500 個の遺伝子もアノテーションされた GO と 5 個の遺伝子だけがアノテーションされた GO の重要度が異なる。そのため、あまりにも少ない遺伝子でしかアノテーションされていない GO を除いて解析する場合もある。例えば、遺伝子数 50-500 でアノテーションされた第 4 層以下の GO のみに対してエンリッチメント解析を行うのが一つの目安と言える。もちろん、この目安は生物種などによって調整する必要があり、試行錯誤が必要である。
実際にやってみる
使うデータはマウスでのコントロール群(2replicates)とCUG-BP1をノックダウンした筋芽細胞(2replicates)群の計4replicatesなっている。この遺伝子は筋ジストロフィーに関わっているとされるものである。まず二つの群を比較し発現変動遺伝子を同定する。つづいて、発現変動遺伝子に有意な GO term を検出する作業を行う。Rのverは4.3で、edgeRをBiocManagerを通してインストールしておく。
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("edgeR")
BiocManager::install("org.Mm.eg.db")
BiocManager::install("clusterProfiler")パッケージの読み込み
library(edgeR)
library(clusterProfiler)
library(org.Mm.eg.db) #マウスのアノテーションデータcounts <- read.table("http://bowtie-bio.sourceforge.net/recount/countTables/katzmouse_count_table.txt", header = TRUE, row.names = 1)
#データのダウンロード
#遺伝子のカウントデータのみなので準備もきっとしやすい
counts <- counts[rowSums(counts) > 10, ] # 低発現量のデータを取り除く
group <- c("CUGBP1", "control", "CUGBP1", "control")edgeR を利用して発現変動遺伝子を検出する。
d <- DGEList(counts = counts, group = group)
d <- calcNormFactors(d)
d <- estimateCommonDisp(d)
d <- estimateTagwiseDisp(d)
result <- exactTest(d)
table <- as.data.frame(topTags(result, n = nrow(counts)))関数の解説は別個ノートを書いたので別ウインドウを開いて参照してください。
ここから GO 解析を行う。この際に発現変動遺伝子を FDR < 0.01 の遺伝子とする。ヒトやマウスなどの場合は Ensembl ID ではなく、 Entrez ID を利用するので、GO 解析の前に遺伝子名の変更を行う。
all.genes <- rownames(table) #列名抜き出し
is.degs <- all.genes[table$FDR < 0.01] #FDR<0.01の遺伝子を抜き出す。
all.genes.entrez <- bitr(all.genes, fromType="ENSEMBL", toType="ENTREZID", OrgDb="org.Mm.eg.db")
#遺伝子名の変更
is.degs.entrez <- bitr(is.degs, fromType="ENSEMBL", toType="ENTREZID", OrgDb="org.Mm.eg.db")
ego <- enrichGO(gene = is.degs.entrez[, 2], # 発現変動遺伝子
universe = all.genes.entrez[, 2], # 全遺伝子
OrgDb = "org.Mm.eg.db", # 生物種
ont = "BP", # オントロジー
pAdjustMethod = "BH", pvalueCutoff = 0.1, qvalueCutoff = 0.1,
readable = TRUE)
res <- as.data.frame(ego)
write.csv(res,"GOenrichment_sample.csv")これで出力されるファイルにp値が小さいものから並びます。中身を見てみるとcell adhesion(細胞接着)やtissue morphogenesis(細胞形態形成)などのGOが変動していたため、筋ジストロフィーとの関連性が示唆されそうです。
実際に運用する場合はサンプルデータと同じようなデータを作って生物種を指定すれば良いです。