FY2024_生命情報学@静岡大学農学部

RNA-seq：群間比較_GOエンリッチメント解析

Licenses

CC BY 4.0

はじめに

RNA-seq解析では、実験群間を統計的に比較して、抽出した発現変動遺伝子群 (Differentially Expressed Genes: DEGs) を一つのセットとして解析することで、生物学的なプロセスを推定することが多い．edgeRは、遺伝子発現が負の二項分布に従うと仮定した検定方法であり、正規化はTMM (Trimmed mean of M values) で行う．edgeRでは、一般化線形モデルを利用した検定法が利用できるため、二群間比較のみならず、二群間二因子比較や多群間比較など複雑な実験デザインにも対応可能である．edgeRの他には、DESeq2やvoomが発現変動遺伝子群の検出によく利用される．本セクションでは、edgeRを用いた発現変動遺伝子群 (DEGs)の検出法について説明する．

遺伝子オントロジー（gene ontology; GO）は、遺伝子の生物的プロセス、細胞の構成要素および分子機能に着目して、遺伝子に付けられるアノテーションである．ある遺伝子に付けられた GO を調べることによって、その遺伝子の機能や細胞内局在を推定できる場合がある．RNA-seq による発現変動解析では、データセットによっては、数百もの発現変動遺伝子を同定できる場合がある．このような大量な遺伝子群に対して、その機能を解明する方法の一つとして遺伝子オントロジーのエンリッチメント解析がある．GOエンリッチメント解析では、詳細な分子経路を理解することはできないが、おおまかな生物学的機能を推定することができる．

edgeRによる発現変動遺伝子群(DEGs)の検出

Countデータの読み込み&サンプル情報付与

まずは解析に必要なedgeRパッケージをインストールする．

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

BiocManager::install("edgeR")

edgeRパッケージを読み込む． 解析で出力されるファイルを格納するディレクトリ(./edgeR)もdir.create関数により作成しておく．

#load library
library(edgeR)
library(tidyverse)

#make dir
dir.create("edgeR")

遺伝子発現カウントデータを読み込む (RSEMによる出力ファイル)． edgeRによる解析への入力は、正規化されたTPMやFPKMデータではなく、countデータを用いる必要があるので注意． (edgeRでは、TMM正規化を実施する．)

本演習では、チャ冬芽を対象に、12月と3月にサンプリングした検体(n=5)のトランスクリプトームデータについて解析を行い、季節によって(統計的有意に)変動する遺伝子群を検出する．

#input count data
count <- read.csv("CountData/Count_Demo.csv", header=T, row.names = 1)

低count遺伝子については、精度が悪いので、除去しておく． ここでは、rowMeans関数を用いて、全サンプル平均count < 10の遺伝子を除去する．

#Genes with average counts <10 were excluded from the analysis
count <- count %>% 
  filter(rowMeans(count)>10)

このあと実施するMDSプロットのために、各サンプルの名前(列名)をわかりやすい情報に変えておく． 12月に採取した検体(n=5)を「Dec」、3月に採取した検体(n=5)を「Mar」とラベルした．

#colnames
colnames(count) <- c("Dec_01", "Dec_02", "Dec_03", "Dec_04", "Dec_05", "Mar_01", "Mar_02", "Mar_03", "Mar_04", "Mar_05")

edgeRによるDEGs検出解析を始める．まずは、統計処理する際の実験群のグルーピングを行う． ここでは、12月 (n=5)、3月 (n=5)のデータなので、以下のように設定する． countデータのサンプル順のメタ情報に対応させること．

#Label for experiment
group <- factor(c(rep("Dec",5), rep("Mar",5)))

わかりやすくするとこんな感じに、各countデータのメタ情報を付与するイメージ

#Label for experiment
group <- factor(c("Dec","Dec","Dec","Dec","Dec","Mar","Mar","Mar","Mar","Mar"))

正確確率検定による比較(二群間比較のみ)

ここでは、正確確率検定による二群間比較を行う． DGEList関数を用いて、edgeRの解析フォーマットにcountデータ、groupデータを格納する． また、TMM正規化等を行う．

#Objectification of table
DGEList <- DGEList(counts = count, group = group)

#TMM normalization
DGEList <- calcNormFactors(DGEList, method="TMM") %>%
  estimateCommonDisp() %>% 
  estimateTagwiseDisp()

plotMDS関数で、MDS (多次元尺度構成法) プロットを描いて、遺伝子発現プロファイルのサンプル間におけるバラツキを確認してみる．

#plotMDS
plotMDS(DGEList)

上図のような結果が得られる．12月 (Dec) と3月 (Mar) でプロットが分離していることがわかる． ここで、実験群毎に遺伝子発現プロファイルがある程度分離していないと、統計検定の結果が妥当なものとして判定しにくいので、解析データのプロファイルを可視化して確認しておくことは重要．

exactTest関数で正確確率検定を実施する．

#extract test
stat.result <- exactTest(DGEList)

ここでは、stat.result変数に統計解析結果が格納されている． 解析結果をわかりやすくまとめるために、topTags関数を用いて、logFC、P-value、FDRなどが整理されたtableにする．

#conversion to dataframe with topTags and Output
res.stat.table <- as.data.frame(topTags(stat.result, n = nrow(count),adjust.method="BH", sort.by="PValue"))
#bind GeneID
res.stat.table <- bind_cols(rownames(res.stat.table), res.stat.table)
colnames(res.stat.table)[1] <- "GeneID"

head関数でTableを確認してみる．

head(res.stat.table, 30)

FC (Fold Change) とは、2群間でどの程度変化があったかを表す指標であり、発現量の比を意味する． RNA-seqデータは、ダイナミックレンジが広いので、対数変換することが多く、ここでもlogFCを算出している． P-valueは帰無仮説が成立する確率なので、値が低ければ低いほど、統計的には確からしいことを意味する． つまり、任意に設定した有意水準(5%など)を満たす遺伝子群が、発現変動遺伝子群(DEGs)と考えられるが、一つ落とし穴がある． 多数の遺伝子発現データを扱うRNA-seqデータでは、多重比較の問題を考慮する必要がある． 例えば、有意水準5%(p <0.05)で100回の検定(多重比較)を行うと、5回の偽陽性(False Positive)が起こる． つまり、発現変動遺伝子検出(DEGs)について、10000個の遺伝子に対して、有意水準5%(p <0.05)を満たすものを調べるとすると、500個の遺伝子は、本当はDEGsでないにも関わらずDEGsと判断されてしまう． このような検定結果は信用できないので、補正を加える必要がある．その一つに、FDR (False Discovery Rate) を使って補正する方法があり、RNA-seqデータ解析では、FDRを適用する場合がほとんどである．

このTableを解析結果として、出力しておく．

#output
write.table(res.stat.table, "edgeR/DEGs_edgeR.txt", col.names = T, row.names = F, sep = "\t")

DGEsとして検出された遺伝子について、アノテーション情報(Gene onthology, KEGGなど)を統合する． ここでは、チャの遺伝子アノテーション情報がまとまったTeaCSS.emapper.annotations_partial.csvファイルを用いる． 詳細は割愛するが、このファイルは機能アノテーションを行うWebツールであるeggNOG-mapperで得られた出力ファイルを整理したものである．

##combine gene annotation
#input annotation data
GeneAno <- read.csv("GeneAno/TeaCSS.emapper.annotations_partial.csv", header=T)

#check
head(GeneAno, 30)

アノテーション情報のGeneAnoデータフレームとDEGs情報のres.stat.tableデータフレームを、left_join関数で統合する． このとき、2つのデータフレームで重複するGeneID列の情報で統合している．

#combine
res.stat.table.ano <- left_join(res.stat.table, GeneAno, by="GeneID", na_matches="never") 

#head
head(res.stat.table.ano)

今回のデータでは、例えば「CSS0022212」、「CSS0043596」、「CSS0006318」などの遺伝子が上位にきている． これら遺伝子のアノテーション情報を確認すると、「Belongs to the cytochrome P450 family」、「Belongs to the expansin family」、「Beta-amylase」であり、他、GOやKEGGのIDも紐づいていることがわかる． これら遺伝子リストのアノテーション情報から、どういった遺伝子が生命現象に関与するかを推定することもできる．例えば、今回のデータでは、expansin fammilyのCSS0043596が12月(冬)に比べて、3月(春)で発現増加していることがわかる．expansinは細胞伸長に関わる遺伝子であり(細胞壁を緩ませる作用)、季節に応じた植物の成長特性を反映していることが推察できる．

次に、FDRで閾値を有意水準を決めて、統計的有意に変動した発現変動遺伝子群(DEGs)を抽出する． ここでは、FDR < 0.05としたが、解析結果や目的に応じて、FDR < 0.01やFDR < 0.001など閾値を変えてみて結果を比較してみるとよい． 統計の閾値(有意水準)の数字が重要というわけではなく、データの中身を向き合うようことが大切である． (FDR =0.051の遺伝子が意味がないかと言われればそうではないということ．)

dplyr::filter関数で、FDR列を<0.05でフィルタリングしている． また、logFC列が0より大きいか小さいかで、発現増加(up-regulated)と発現減少(down-regulated)の遺伝子群に分けている． フィルタリング後のファイルについては、write.table関数で出力．

##filtered significant
#FDR only
#up
res.stat.table.upDEG <- res.stat.table.ano %>% dplyr::filter(FDR<0.05, logFC>0)
#up
res.stat.table.downDEG <- res.stat.table.ano %>% dplyr::filter(FDR<0.05, logFC<0)

#output
write.table(res.stat.table.upDEG, "edgeR/DEGs_edgeR_up_FDR005.txt", col.names = T, row.names = F, sep = "\t")
write.table(res.stat.table.downDEG, "edgeR/DEGs_edgeR_down_FDR005.txt", col.names = T, row.names = F, sep = "\t")

統計的有意に発現変動した遺伝子がどの程度存在しているのかをMA-plotで可視化してみる． ここでは、FDR < 0.05を満たす遺伝子群をDEGsとして抽出している．

##MA-plot
#Data arrangement for MA-plot
is.DEG <- as.logical(res.stat.table$FDR < 0.05) 
DEG.names <- rownames(res.stat.table)[is.DEG] 

#plot
plotSmear(stat.result, de.tags = DEG.names, cex=0.3)

上図のような結果が得られる．赤プロットが発現変動した遺伝子群を意味する． 縦軸がlogFC(発現変動比)、横軸がlogCPM(count per million；正規化された発現量を表す一つの指標)を意味する． 本データセットでも、ある程度発現変動している遺伝子群が存在することがわかる．

topGOによるGO enrichment解析

RパッケージのtopGOを用いてGO enrichment解析を実施する．

まずは解析に必要なedgeRパッケージをインストールする．

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

BiocManager::install("topGO")

topGOパッケージを読み込む． 解析で出力されるファイルを格納するディレクトリ(./topGO)もdir.create関数により作成しておく．

#library
library(topGO)

#make dir
dir.create("topGO")

上述したeggNOG mapperにより出力した全遺伝子のGO term情報を読み込む． ここでは、チャ遺伝子のGO term情報がまとまったorf2go.mapファイルを用いる．

##read output of eggNOG mapper: Go information of all genes
geneID2GO <- readMappings("GeneAno/orf2go.map", sep = "\t", IDsep = ",")
#gene name
geneNames <- names(geneID2GO)
#check
str(head(geneID2GO))

## List of 6
##  $ #query      : chr "GOs"
##  $ CSS0037220.1: chr [1:57] "GO:0005575" "GO:0005622" "GO:0005623" "GO:0005634" ...
##  $ CSS0024626.1: chr "-"
##  $ CSS0019345.1: chr "-"
##  $ CSS0024537.1: chr [1:157] "GO:0000003" "GO:0000988" "GO:0000989" "GO:0001076" ...
##  $ CSS0046910.1: chr "-"

解析対象としたい遺伝子群を読み込む． ここでは、上述で同定した発現増加したDEGs (res.stat.table.upDEGs) を対象とする．

#read query genes
TargetGene <- res.stat.table.upDEG$GeneID %>% as.data.frame()
colnames(TargetGene) <- "GeneID"

対象遺伝子群を、topGOで解析するための形式にフォーマットする．

#resetting gene names corresponding to GO table
QueryGenes <- str_c(as.character(TargetGene[,1]), ".1")
#read as list
list <- unlist(as.list(QueryGenes))
#count query genes number in all genes
geneList <- factor(as.integer(geneNames %in% list)) 
names(geneList) <- geneNames

全遺伝子群と対象遺伝子群の情報を含めて、topGOで解析するための形式にフォーマットする． またここで、解析を実施するGOカテゴリーを選択する (BP or CC or MF)．

BP: Biological Process (生物学的なプロセス) CC: Cellular Component (細胞の構成要素) MF: Molecular Function (分子機能)

ここでは、BPを対象に解析を行う．

#Construct topGOdata
GOdata <- new("topGOdata", 
              ontology = "BP", #BP or CC or MF
              allGenes = geneList,
              annot = annFUN.gene2GO, 
              gene2GO = geneID2GO, 
              nodeSize = 20) 

###Descrition of parameter
##annFUN: 
#gene2GO this function is used when the annotations are provided as a gene-to-GOs mapping.
##nodeSize: 
#an integer larger or equal to 1. 
#This parameter is used to prune the GO hierarchy from the termswhich have less than nodeSize annotated genes 
#(after the true path rule is applied)

topGO::runTest関数により、enrichment解析を実行する． アルゴリズム (algorithm) と統計手法 (statistic) をオプションで選択できる． これらの手法の詳細については、公式マニュアルの記事を確認するとよい． ここでは、algorithm = "elim"、statistic = "fisher"に設定している(基本的にはこれで大丈夫)．

#Enrichment test
resultFis <- topGO::runTest(GOdata, 
                            algorithm = "elim", #classic or elim or weight01(default)
                            statistic = "fisher") 

enrichment解析結果について、上記の出力のtopGOフォーマットからテーブル形式に出力する． GenTableに解析結果が出力される．

#GenTable
num_nodes=length(GOdata@graph@nodes) 

GenTable <- GenTable(GOdata, 
                     Pvalue = resultFis,
                     topNodes = num_nodes)

#check
head(GenTable, 10)

続いて、多重検定を行う (多重検定については、上述した通り)． ここでは、Benjamini-Hochberg 法によるFDRで統計処理． 多重検定結果をそのまま、GenTableデータフレームに結合している．

Adj.Pvalueの列カラムが追加されており、これが多重検定による補正後のp値(FDR)である．

#Multiple test
GenTable <- bind_cols(GenTable, p.adjust(GenTable$Pvalue, method = "BH"))
colnames(GenTable)[7] <- "Adj.Pvalue"

#check
head(GenTable, 10)

Adj.Pvalue(FDR) < 0.05を有意水準として、統計的有意に濃縮(enrich)されたGO termを抽出する．

#fileter adj.pvalue < 0.05
GenTable.Sig <- GenTable %>% filter(Adj.Pvalue < 0.05)

#check
head(GenTable.Sig, 30)

上表の通り、本データセットでは、統計的有意に濃縮(enrich)されたGO termは15個検出された． 例えば、“response to cytokinin(GO:0009735)”、“jasmonic acid biosynthetic process GO:0009695)”、“auxin effluxGO:0010315)”のGO termが濃縮されている．これは、サイトカイニン、ジャスモン酸、オーキシン等の植物ホルモン関連のGO termであり、冬芽の季節変化において植物ホルモン応答が関与することを示唆する結果となっている． このように、GO enrichment解析では、対象とする遺伝子群がどのような機能を有しているかを推定することで、対象とする生命現象に関連する分子応答を大まかに理解することできる． ただし、GO enrichment解析のみでは、分子メカニズムを語ることは難しいので、あくまでデータ駆動研究における仮説立てに役立てると良い．

最後に、統計処理前後のGOenrichment解析結果について、出力しておく．

write.csv(GenTable, "topGO/DEGs_up_GOenrich.csv", row.names = FALSE)
write.csv(GenTable.Sig, "topGO/DEGs_up_GOenrich_FDR005.csv", row.names = FALSE)