ex803 - nibb-gitc/gitc2021mar-rnaseq GitHub Wiki
ここでは、演習80１で作成したアノテーションを用いてGOエンリッチメント解析を行う。
S. eubayanus ゲノムにはマルトースを取り込んで利用する遺伝子は存在するものの、実際にはその能力がないことが知られているが、この研究では、S. cerevisiaeのある遺伝子を添加することでその能力を回復できることを示している。このデータは、その遺伝子を添加した状態でのグルコース、およびマルトースの存在下での遺伝子発現を比較したものであり、マルトース存在下でマルトースの代謝系の発現が上昇しているかどうかがポイントとなっている。このことを念頭に、以下の解析を行ってみよう。
Rを用いるので、ローカル環境で作業する。ディレクトリ~/gitc/data/IU/yeast以下のファイルを用いる。ここで用いるのは以下の4ファイルである。
| file | contents | remarks |
|---|---|---|
| topTags.transcript.txt | EdgeRの結果（トランスクリプトレベル） | |
| topTags.gene.txt | EdgeRの結果（遺伝子レベル） | |
| blast_bbh.tab または blast_top.tab | BLASTベストヒット | 演習801のStep 4で作成 |
| seub_genes.emapper.go | EggNOG-mapperによるGOアサインメント | 演習801のStep 7で作成 |
このうち、後者２つは演習801のStep4, Step7でそれぞれ作成しているので、それらをbias5からローカルに転送するか、演習801におけるコマンドをローカルで実行して作成しておくこと。ただし、blast_bbh.tabの代わりにblast_top.tabを用いても良い。以下は、bias5上で作業ディレクトリex1を作って作業した場合の転送例。
% scp [email protected]:gitc/data/IU/ex1/blast_bbh.tab .
% scp [email protected]:gitc/data/IU/ex1/seub_genes.emapper.go .
差し当たり、トランスクリプトレベルの方が簡単なので、こちらを使おう。R を立ち上げてedgeRの結果を読み込むみ、headで確認しよう。
> etab.transcript <- read.table("topTags.transcript.txt")
> head(etab.transcript)
logFC logCPM F PValue FDR
XM_018364781.1 6.872644 7.169865 3160.089 2.479345e-07 0.0006804321
XM_018366251.1 5.288273 7.587097 2970.877 2.831039e-07 0.0006804321
XM_018368189.1 6.916983 7.706345 2361.930 4.633663e-07 0.0006804321
logFC>0がマルトース存在下で発現上昇した遺伝子である。演習802のStep1, Step4を参照して、解析対象とする発現上昇遺伝子リスト(upregGenes)、全遺伝子リスト(allGenes)、発現スコアリスト(expScore)を作成しておこう。
blast_bbh.tab は、S. eubayanusの各遺伝子について、双方向ベストヒット基準を用いて S. cerevisiaeにおけるオーソログを同定した結果である。(このファイルの代わりにblast_top.tab（クエリから見た一方向のベストヒット）を用いてもよいこととするが、その場合はパラログがいくらか混入するおそれがある代わりに、アノテーション可能な遺伝子数は増える。)
DI49_1142 NP_009669.1 93.421 228 15 0 8 235 4 231 9.83e-162 444 ADP-ribose diphosphatase [Saccharomyces cerevisiae S288C]
DI49_2764 NP_012252.1 88.604 1404 156 3 1 1401 1 1403 0.0 2560 ATP-binding cassette multidrug transporter PDR11 [Saccharomyces cerevisiae S288C]
MSTRG.1000.1 NP_013901.1 68.702 1425 423 7 1 1416 1 1411 0.0 2005 Ecm5p [Saccharomyces cerevisiae S288C]
MSTRG.1039.1 NP_013956.1 86.944 1731 219 2 1 1726 1 1729 0.0 3026 Rrp5p [Saccharomyces cerevisiae S288C]
MSTRG.105.1 NP_012826.1 79.747 79 8 1 1 79 1 71 2.61e-26 91.3 Cwp2p [Saccharomyces cerevisiae S288C]
この結果に基づいて、S. eubayanus の遺伝子とS. cerevisiaeを結びつけることによって、公開データベースのアノテーションが使えるようになる。まずはこの方針で解析してみよう。
R に読み込む。ここで読み込みにread.tableを用いることもできるが、その場合は sep="\t"（タブで分割）に加えてquote="" （引用記号を解釈しない）を指定する必要がある（検索結果のタイトル中にシングルクオーテーション(')が入っているため）。read.delimであればこれらは必要がない。
> blasthit <- read.delim("blast_bbh.tab", header=FALSE)
> head(blasthit)
V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12
1 DI49_1142 NP_009669.1 93.421 228 15 0 8 235 4 231 9.83e-162 444.0
2 DI49_2764 NP_012252.1 88.604 1404 156 3 1 1401 1 1403 0.00e+00 2560.0
3 MSTRG.1000.1 NP_013901.1 68.702 1425 423 7 1 1416 1 1411 0.00e+00 2005.0
4 MSTRG.1039.1 NP_013956.1 86.944 1731 219 2 1 1726 1 1729 0.00e+00 3026.0
２カラム目のIDの末尾についている .1 などの番号は、異なるバージョンを区別するバージョン番号で、REFSEQのアクセション番号には含まれない。アクセション番号と対応を取るにはこれらを除かなければならない。これは、正規表現置換を行う sub 関数を用いて行える。以下、まずblasthitの最初の2カラムをortho_pairという変数にコピーし、sub関数を用いて２カラム目のバージョン番号を除く。
> ortho_pair <- blasthit[,1:2]
> ortho_pair[,2] <- sub("\\.[0-9]+", "", ortho_pair[,2])
> head(ortho_pair)
V1 V2
1 DI49_1142 NP_009669
2 DI49_2764 NP_012252
3 MSTRG.1000.1 NP_013901
S. cerevisiaeのアノテーションデータベースorg.Sc.sgd.dbを使ってGOをアサインしよう。ここで使われているNP_009669などのIDはREFSEQのアクセション番号なので、keytype="REFSEQ"を指定して、以下のようにデータを取得できる。
> library(org.Sc.sgd.db)
> select(org.Sc.sgd.db, keys="NP_009669", keytype="REFSEQ", columns="GO")
'select()' returned 1:many mapping between keys and columns
REFSEQ GO EVIDENCE ONTOLOGY
1 NP_009669 GO:0005634 IDA CC
2 NP_009669 GO:0005737 IDA CC
3 NP_009669 GO:0005739 IDA CC
そこで、まずオーソログテーブルortho_pairに含まれるすべてのREFSEQ IDについて、org.Sc.sgd.dbに格納されているGOの対応を抜き出す。
> go_annot_sce <- select(org.Sc.sgd.db, keys=ortho_pair[,2], keytype="REFSEQ", columns="GO")
このテーブルとさきほどのオーソログテーブルをつなげよう。これはmerge関数で行える。この関数は２つのテーブルを引数にとり、共通に存在するカラム(by.x, by.yで指定)を使って同じ値を持つ行を結び付けて一つのテーブルにまとめる。リレーショナルデータベースにおけるjoinの機能である。
> go_annot_merged <- merge(ortho_pair, go_annot_sce, by.x=2, by.y=1)
> head(go_annot_merged)
V2 V1 GO EVIDENCE ONTOLOGY
1 NP_001018030 MSTRG.2869.1 GO:0003941 IDA MF
2 NP_001018030 MSTRG.2869.1 GO:0003941 IEA MF
3 NP_001018030 MSTRG.2869.1 GO:0004794 IBA MF
これで遺伝子とGOの対応関係はついたので、あとはclusterProfilerで扱えるようにデータの形を整える。
> go2gene <- unique(go_annot_merged[,c(3,2)])
> names(go2gene) <- c('term', 'gene')
> head(go2gene)
term gene
1 GO:0003941 MSTRG.2869.1
3 GO:0004794 MSTRG.2869.1
6 GO:0005739 MSTRG.2869.1
これでGOとの対応表はできた。clusterProfilerをロードし、Case 2-1の「解析の下準備」以下を行って解析を進めよう。
演習80１ではEggNOG mapperを使ってGOアサインメントを行った。その結果 seub_genes.emapper.go を用いてエンリッチメント解析を行おう。中身は以下のようになっている。
% head seub_genes.emapper.go
DI49_1142 GO:0003674,GO:0003824,GO:0005575,GO:0005622,GO:0005623,GO:0005634,GO:0005737,GO:0005739,GO:0006139,...
DI49_2764 GO:0003674,GO:0003824,GO:0005215,GO:0005319,GO:0005575,GO:0005622,GO:0005623,GO:0005737,GO:0005739,...
MSTRG.1000.1 GO:0000118,GO:0000228,GO:0000785,GO:0000790,GO:0000812,GO:0003674,GO:0003824,GO:0005575,GO:0005622,...
MSTRG.1000.3 GO:0000118,GO:0000228,GO:0000785,GO:0000790,GO:0000812,GO:0003674,GO:0003824,GO:0005575,GO:0005622,...
clusterProfilerには、read.gmt関数が用意されており、以下のGMT形式でGOと遺伝子の対応表を読み込むことができる。
GO-1 GOname1 gone1 gene2 gene3 ...
GO-2 GOname2 gone4 gene5 gene6 ...
残念ながら、GOとgeneの関係が逆になっているが、形式的には同じ形をしているので、とりあえず強引に読み込んで、あとでGOとgeneを入れ替えることにする。まず、2カラム目にGO nameを入れ、タブ区切りに変換する。ここではGO nameは使わないのでNA（欠損値）を挿入している。GO IDとGO nameの対応表は後ほど別に作成する。なお、"\t"はタブを意味するが、MacOSのsedは"\t"をタブと認識しないので、<Control-v><TAB>で入力すること。
% sed 's/\t/\tNA\t/' seub_genes.emapper.go | sed 's/,/\t/g' > seub_genes.emapper.gmt
% head seub_genes.emapper.gmt
DI49_1142 NA GO:0003674 GO:0003824 GO:0005575 GO:0005622 GO:0005623 GO:0005634
DI49_2764 NA GO:0003674 GO:0003824 GO:0005215 GO:0005319 GO:0005575 GO:0005622
MSTRG.1000.1 NA GO:0000118 GO:0000228 GO:0000785 GO:0000790 GO:0000812 GO:0003674
MSTRG.1000.3 NA GO:0000118 GO:0000228 GO:0000785 GO:0000790 GO:0000812 GO:0003674
Rで、read.gmtを用いて読み込む。
% R
> library(clusterProfiler)
> gene2go <- read.gmt("seub_genes.emapper.gmt")
> head(gene2go)
term gene
1 DI49_1142 GO:0003674
2 DI49_1142 GO:0003824
3 DI49_1142 GO:0005575
名前を付け直した上で、1列目と2列目の値を交換する。
> colnames(gene2go) <- colnames(gene2go)[c(2,1)]
> go2gene <- gene2go[, c(2,1)]
> head(go2gene)
term gene
1 GO:0003674 DI49_1142
2 GO:0003824 DI49_1142
3 GO:0005575 DI49_1142
clusterProfilerでは、ユーザアノテーションを用いる場合は、GOに特化したインターフェイスであるenrichGOではなく、汎用の低レベルインターフェイスenricherを使う必要がある。このため、いくつかの処理を自前で行う必要がある。まず、GO階層を遡って、子ノードにアサインされた遺伝子を親ノードにもアサインする処理を、自前で行わないといけない。これは、buildGOmap関数で行う。
> go2gene.expand <- buildGOmap(go2gene)
GO ID と名前の対応表も自前で用意する必要がある。これは、OrgDbなどと同じAnnotationDbiインターフェイスで実装されたGOデータベースのパッケージ、GO.dbを用いて作成できる。GO term を３つの種類、BP/MF/CC に分けることを想定して、その情報も付加する。
> library(GO.db)
> columns(GO.db)
[1] "DEFINITION" "GOID" "ONTOLOGY" "TERM"
> gonames <- select(GO.db, keys=keys(GO.db), columns=c("TERM", "ONTOLOGY"))
先ほど作成したgo2geneは、すべてのGO termが含まれているため、これをBP/MF/CCに分けて解析できるように分割する。まず、先ほど作成したgonamesをオントロジーの種類ごとに分けて、その情報を使って先ほど拡張したgo2gene.expandの情報を分割する。以下は BP を抽出する場合。
> gonames.bp <- subset(gonames, ONTOLOGY=="BP")
> go2gene.expand.bp=subset(go2gene.expand, GO %in% gonames.bp$GOID)
これで準備は整った。以下、低レベルインターフェイスを用いてエンリッチメント解析を行う。
低レベルインターフェイスを用いてFisher's exact test を実行する。
> cprof.up.bp.fisher2 <- enricher(gene=upregGenes, TERM2GENE=go2gene.expand.bp, TERM2NAME=gonames)
低レベルインターフェイスを用いてGene Set Enrichment Analysis を実行する。
> cprof.bp.gsea2 <- GSEA(geneList=expScore, TERM2GENE=go2gene.expand.bp, TERM2NAME=gonames)
topGOでは、ユーザ定義アノテーションを取り込んだ解析がよりシンプルに行える。readMapping 関数が定義されており、これを用いて遺伝子とGOの対応表を読み込める。想定されているのは以下のような形式で、区切り文字はオプションで変更できる。
Gene1 GO1, GO2, GO3
Gene2 GO4, GO5
seub_genes.emapper.goは、まさにこの形式なので、そのまま読み込むことができる。
> gene2go.topGO <- readMappings("seub_genes.emapper.go")
これを用いて解析するには、topGOdata オブジェクト作成時にアノテーション用関数としてannFUN.gene2GOを用い、引数gene2GOにこの変換テーブルを指定する。
> topgoData <- new("topGOdata", ontology="BP", allGenes=expScore, geneSel=selfun, annot=annFUN.gene2GO, gene2GO=gene2go.topGO)
あとは、runTestで解析を実行すれば良い。
実際の解析では、トランスクリプト単位ではなく、遺伝子単位でアイソフォームを集計したカウントマトリクスを使うことも多い。アノテーションはトランスクリプト配列に対してつけているので、その場合は遺伝子とトランスクリプトの対応表を作る必要がある。一般に、一つの遺伝子に複数のトランスクリプトが対応づくので、その中で代表を一つ選んで対応づけるのが簡単である。代表の選び方としては発現量の大きいものをとるという考え方もあるが、アノテーションの目的であれば最も長いものを取るというのもリーズナブルである。この場合、gffread を使って以下のように行える。
% gffread --table @geneid,@id,@covlen stringtie_merged.gtf > gene2transcript_all.txt
% sort -k 1,1 -k 3,3nr gene2transcript_all.txt | sort -k 1,1 -s -u > gene2transcript_maxlen.txt
% head gene2transcript_maxlen.txt
DI49_0001 XM_018363258.1 1251
DI49_0054 DI49_0054 74
DI49_0066 XM_018363322.1 582
DI49_0073 XM_018363328.1 267
これが遺伝子ID（１カラム目）とトランスクリプトID（２カラム目）の対応表になる。この結果と、遺伝子レベルのedgeRの結果(topTags.gene)をRで読み込んでジョインする。遺伝子IDを共通フィールドとしてmergeを実行する。
> etab.gene <- read.table("topTags.gene.txt")
> head(etab.gene)
logFC logCPM F PValue FDR
MSTRG.3686 6.791698 10.403534 2043.018 6.215580e-09 8.791816e-06
MSTRG.1820 6.915789 7.706148 2038.774 6.254992e-09 8.791816e-06
MSTRG.2 5.288289 7.586859 1961.359 7.035849e-09 8.791816e-06
> gene2transcript <- read.delim("gene2transcript_maxlen.txt", header=FALSE)
> head(gene2transcript)
V1 V2 V3
1 DI49_0001 XM_018363258.1 1251
2 DI49_0054 DI49_0054 74
3 DI49_0066 XM_018363322.1 582
> etab2.gene <- merge(etab.transcript, gene2transcript, by.x=0, by.y=1)
> head(etab2.gene)
Row.names logFC logCPM F PValue FDR
1 DI49_0001 0.16578946 0.2879501 0.12955874 0.731036019 0.81856500
2 DI49_0079 1.54214299 -1.3284993 11.75625001 0.013655489 0.04530080
3 DI49_0855 -0.59553762 0.3485453 1.88125843 0.218547353 0.34622968
V2 V3
1 XM_018363258.1 1251
2 XM_018363333.1 1842
3 XM_018364071.1 1761
7カラム目にトランスクリプトIDがついたテーブルができた。これを、これまでのedgeR出力と同じ形式で、遺伝子IDをトランスクリプトIDに置き換えたものに変換したい。これは以下のようにしてに行える。
> row.names(etab2.gene) <- etab2.gene[,7]
> etab2.gene <- etab2.gene[,2:6]
> head(etab2.gene)
logFC logCPM F PValue FDR
XM_018363258.1 0.16578946 0.2879501 0.12955874 0.731036019 0.81856500
XM_018363333.1 1.54214299 -1.3284993 11.75625001 0.013655489 0.04530080
XM_018364071.1 -0.59553762 0.3485453 1.88125843 0.218547353 0.34622968