clusterProfiler | RNASeq に対応した遺伝子オントロジーエンリッチメント解析

clusterProfiler は R/Bioconductor のパッケージの一つで、発現変動遺伝子にアノテーションされた GO と KEGG パスウェイなどを統計的に有意のあるものを検出する機能を提供してる。ここでは GO 解析の部分を紹介する。

サンプルに用いるデータは ReCount データベースの katz.mouse を利用する。このデータは 2 群間比較のデータであり、ここで、まず二つの群を比較し発現変動遺伝子を同定する。つづいて、発現変動遺伝子に有意な GO term を検出する作業を行う。

まず、データベースからデータをダウンロードする。

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 を利用して発現変動遺伝子を検出する。

library(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)))

clusterProfiler を利用して GO 解析を行う。この際に発現変動遺伝子を FDR < 0.01 の遺伝子とする。ヒトやマウスなどの場合は Ensembl ID ではなく、 Entrez ID を利用するので、GO 解析の前に遺伝子名の変更を行う。

library(clusterProfiler)
library(org.Mm.eg.db)   # マウスのアノテーションデータ

all.genes <- rownames(table)
is.degs <- all.genes[table$FDR < 0.01]

all.genes.entrez <- bitr(all.genes, fromType="ENSEMBL", toType="ENTREZID", annoDb="org.Mm.eg.db")
is.degs.entrez <- bitr(is.degs, fromType="ENSEMBL", toType="ENTREZID", annoDb="org.Mm.eg.db")


ego <- enrichGO(gene = is.degs.entrez[, 2],          # 発現変動遺伝子
                universe = all.genes.entrez[, 2],    # 全遺伝子
                organism = "mouse",                  # 生物種
                ont = "BP",                          # オントロジー
                pAdjustMethod = "BH", pvalueCutoff = 0.1, qvalueCutoff = 0.1,
                readable = TRUE)
res <- summary(ego)

解析結果は res オブジェクトに保存される。これをファイルに保存すればよい。head 関数で確認すると以下のようになる。（※ 8 列目が非常に長いため、ここでは 8 列目を除いて表示させて見やすいようにしている。）

head(res[, -8])
##                    ID                                    Description GeneRatio
## GO:0044707 GO:0044707          single-multicellular organism process   210/443
## GO:0032501 GO:0032501               multicellular organismal process   211/443
## GO:0009653 GO:0009653             anatomical structure morphogenesis   119/443
## GO:0007275 GO:0007275           multicellular organismal development   172/443
## GO:0048731 GO:0048731                             system development   157/443
## GO:0051239 GO:0051239 regulation of multicellular organismal process   115/443
##              BgRatio       pvalue     p.adjust       qvalue Count
## GO:0044707 1926/6672 9.828877e-18 3.738905e-14 3.086267e-14   210
## GO:0032501 1972/6672 8.060899e-17 1.533183e-13 1.265561e-13   211
## GO:0009653  938/6672 1.496734e-13 1.897858e-10 1.566581e-10   119
## GO:0007275 1589/6672 2.168149e-13 2.011456e-10 1.660350e-10   172
## GO:0048731 1402/6672 2.643870e-13 2.011456e-10 1.660350e-10   157
## GO:0051239  910/6672 5.970621e-13 3.785374e-10 3.124625e-10   115