mRNA-Seq入門【第四回】

mRNA-Seq 解析の流れをざっくりと説明してみた mRNA-Seq 解析 発現量解析編 2021/12/03 ⽔産⽣物環境学（九州⼤学） ⾼井優⽣

前回がちゃがちゃ⾔いましたが 発現量解析も同じです

マニュアルを しっかりと 読んでください

解析⼿法によって結果は⼤きく変わってきます 特に なので発現量解析は です

でもマニュアルって英語で書かれとるしな、、、

ちょっと調べるだけで ⽇本語で解析⽅法を説明しているブログが たくさん出てきます （実際僕も参考にさせていただいてます）

が

きちんと実⾏内容を理解して解析に臨みましょう

訳もわからぬまま ただひたすらに

ブログからの コピペ＆コピペ

そして発覚する 盛⼤な解析ミス

阿⿐叫喚 荒れるゼミ 地獄絵図

こんな悲しい物語を 避けたければ

マニュアルを しっかりと 読んでください

というわけで、今回の説明も あくまで参考であることを 脳みそに刻み込み

実際にご⾃⾝で解析をする際は

マニュアルを

しっかりと

読んで ください

mRNA-Seq 解析の流れをざっくりと説明してみた mRNA-Seq 解析 発現量解析編 2021/12/03 ⽔産⽣物環境学（九州⼤学） ⾼井優⽣

R Project の作成 R（4.0.3） とりあえず Rstudio を起動して、解析⽤の R project を作りましょう これは必須ではないけども、R Project 単位で管理した⽅がいろいろと楽です

R Project の作成 R（4.0.3） R Project ファイルを保存するための ディレクトリを指定します ← 新しくディレクトリを作る場合 ← 既存のディレクトリを使う場合 あとはぽちぽちしていってください

データセット作成 edgeR

基本的な⼆群間⽐較
edgeR（3.32.1）

edgeR を読み込みます
> library(edgeR)

カウントデータを読み込んで、head で読み込んだデータを確認します
> count <- read.table("count_data.tsv", header = T, row.names = 1)
> head(count)
C1 C2 C3 C4 A1 A2 A3 A4
ENSOJAG00000000068 0 0 0 0 0 0 0 0
ENSOJAG00000000073 0 0 0 0 0 0 1 0
ENSOJAG00000000082 0 0 0 0 0 0 0 0
ENSOJAG00000000089 0 0 0 0 0 0 0 0
ENSOJAG00000000092 4 4 0 0 0 2 5 0
ENSOJAG00000000110 0 0 0 0 0 0 1 0

解析⽤の識別ラベルを作ります（ここめっちゃ重要です、特に levels に気をつけて）
> group <- factor(rep(c('C', 'A'), c(4, 4)), levels = c('C', 'A'))

基本的な⼆群間⽐較 edgeR（3.32.1） DGEList 関数で edgeR が扱うためのデータセットを作ります（詳細はマニュアル参照） で、データセットの内の遺伝⼦数を確認します > y <- DGEList(counts = count, group = group) > nrow(y$counts) [1] 24179 filterByExp 関数で発現量の少ない遺伝⼦を除去します（詳細はマニュアル参照） で、解析に使⽤する遺伝⼦数を確認します > keep <- filterByExpr(y, group = group) > y <- y[keep, , keep.lib.sizes = FALSE] > nrow(y$counts) [1] 11128 今回の場合は 24,179 遺伝⼦中 11,128 遺伝⼦が解析に使⽤されるようです

基本的な⼆群間⽐較 edgeR（3.32.1） 正規化してデータの分散を推定します（詳細はマニュアル参照） > y <- calcNormFactors(y) > design <- model.matrix(~ group) > y <- estimateDisp(y, design) ここらで⼀度、サンプル間のばらつき具合を⾒るために、MDS プロットを描いてみましょう > plotMDS(y)

基本的な⼆群間⽐較 edgeR（3.32.1） 正規化してデータの分散を推定します（詳細はマニュアル参照） > y <- calcNormFactors(y) > design <- model.matrix(~ group) > y <- estimateDisp(y, design) ここらで⼀度、サンプル間のばらつき具合を⾒るために、MDS プロットを描いてみましょう > plotMDS(y) C グループと A グループが まあそれなりに分かれているので よしとしましょう

発現量⽐較 edgeR

発現量⽐較の⽅法 edgeR（3.32.1） 発現量⽐較は以下 3 つの⽅法を選べます ・正確確率検定（⼆群間⽐較） ・擬似尤度 F 検定（⼆群間⽐較、多群⽐較） ・尤度⽐検定（⼆群間⽐較、多群⽐較） マニュアルを読む感じでは edgeR は 擬似尤度 F 検定を推奨しているようです まあ、でも⾃分の主張したいことが⾔える検定があればそれを、、、

正確確率検定の場合 edgeR（3.32.1） 正確確率検定で⼆群を⽐較します（詳細はマニュアル参照） > et <- exactTest(y) > topTags(et) Comparison of groups: ENSOJAG00000018326 ENSOJAG00000005553 ENSOJAG00000011505 ENSOJAG00000004540 ENSOJAG00000009155 ENSOJAG00000011793 ENSOJAG00000001481 ENSOJAG00000018807 ENSOJAG00000012089 ENSOJAG00000020385 A-C logFC 4.719640 3.891036 4.018177 3.384223 3.346690 -7.423284 6.351456 3.066405 -2.578423 2.537101 logCPM 3.927181 4.765416 5.064951 6.664245 5.292157 4.604369 3.361985 2.419279 6.215812 3.757806 PValue 7.072340e-28 1.687091e-27 2.873377e-24 2.437048e-21 6.394592e-19 1.842783e-16 2.993140e-15 2.044566e-12 8.035441e-12 9.192831e-12 FDR 7.870100e-24 9.386976e-24 1.065831e-20 6.779869e-18 1.423180e-15 3.417748e-13 4.758238e-12 2.843992e-09 9.935377e-09 1.022978e-08

疑似尤度 F 検定の場合 edgeR（3.32.1） 疑似尤度 F 検定で⼆群を⽐較します（詳細はマニュアル参照） > fit <- glmQLFit(y, design) > qlf <- glmQLFTest(fit) > topTags(qlf) Coefficient: groupA logFC ENSOJAG00000005553 3.883698 ENSOJAG00000011505 4.016347 ENSOJAG00000004540 3.386678 ENSOJAG00000018326 4.726539 ENSOJAG00000009155 3.350300 ENSOJAG00000011793 -7.424337 ENSOJAG00000016468 2.141731 ENSOJAG00000012089 -2.578196 ENSOJAG00000001481 6.348479 ENSOJAG00000019017 2.268050 logCPM 4.765416 5.064951 6.664245 3.927181 5.292157 4.604369 3.375748 6.215812 3.361985 8.575967 F 201.39841 149.29239 143.65134 130.14953 112.05123 70.08246 69.26781 67.51540 66.54909 62.11718 PValue 2.553015e-08 1.175634e-07 1.428170e-07 2.346230e-07 4.949007e-07 4.818372e-06 5.092098e-06 5.746007e-06 6.149079e-06 8.488193e-06 FDR 0.0002840995 0.0005297560 0.0005297560 0.0006527213 0.0011014510 0.0076029944 0.0076029944 0.0076029944 0.0076029944 0.0091734060 多群⽐較にも応⽤可能です（詳細はマニュアル参照）

尤度⽐検定の場合 edgeR（3.32.1） 尤度⽐検定で⼆群を⽐較します（詳細はマニュアル参照） > fit <- glmFit(y, design) > lrt <- glmLRT(fit) > topTags(lrt) Coefficient: groupA logFC ENSOJAG00000018326 4.719640 ENSOJAG00000005553 3.891036 ENSOJAG00000011505 4.018177 ENSOJAG00000004540 3.384223 ENSOJAG00000009155 3.346690 ENSOJAG00000011793 -7.423284 ENSOJAG00000001481 6.351456 ENSOJAG00000018807 3.066405 ENSOJAG00000017959 7.258648 ENSOJAG00000012089 -2.578423 logCPM LR PValue FDR 3.9271815 118.74352 1.191857e-27 7.025014e-24 4.7654160 118.62917 1.262583e-27 7.025014e-24 5.0649509 103.75510 2.289350e-24 8.491961e-21 6.6642448 90.27790 2.069509e-21 5.757374e-18 5.2921566 79.15275 5.748870e-19 1.279468e-15 4.6043690 70.04172 5.806340e-17 1.076882e-13 3.3619854 63.34838 1.731977e-15 2.753348e-12 2.4192786 48.92766 2.655784e-12 3.694195e-09 0.3709419 48.15884 3.930537e-12 4.859890e-09 6.2158120 47.13966 6.610461e-12 7.356121e-09 多群⽐較にも応⽤可能です（詳細はマニュアル参照）

発現量の定量 edgeR

発現量の定量 edgeR（3.32.1） edgeR では⼀応発現量の定量もできます（CPM） > logcpm <- cpm(y, log = T) log2CPM で計算すればとりあえずのヒートマップもぽんっと描けます > logcpm_deg <- logcpm[row.names(table_deg),] > heatmap(logcpm_deg, labRow = F) 各遺伝⼦の塩基配列⻑データがあれば RPKM にも換算できます （マニュアル参照） でも最近は TPM を推奨してる⼈が多いから、論⽂とかにするときは ちょっと頑張って TPM にした⽅が良い気もします

データの整理 R、エクセル

データの書き出し ファイルに書き出します 解析されたデータをデータフレームにします > table_all <- as.data.frame(topTags(qlf, n = nrow(qlf))) ⼀般的に FDR < 0.05 が発現変動遺伝⼦の⽬安とされているので それに従って発現変動遺伝⼦を抽出してみます > table_deg <- table_all[table_all$FDR < 0.05,] 発現変動遺伝⼦の数を確認してみます > nrow(table_all) [1] 11128 > nrow(table_deg) [1] 794 ふーんと思いながら、解析結果を csv ファイルに書き出します > write.csv(table_all, “table_all.csv”)

データの書き出し ファイルに書き出します 解析されたデータをデータフレームにします > table_all <- as.data.frame(topTags(qlf, n = nrow(qlf))) ⼀般的に FDR < 0.05 が発現変動遺伝⼦の⽬安とされているので それに従って発現変動遺伝⼦を抽出してみます > table_deg <- table_all[table_all$FDR < 0.05,] 発現変動遺伝⼦の数を確認してみます > nrow(table_all) [1] 11128 > nrow(table_deg) [1] 794 ふーんと思いながら、解析結果を csv ファイルに書き出します > write.csv(table_all, “table_all.csv”) ここまでの流れは R Markdown（R Notebook）で 保存するのがおすすめです

データの書き出し ファイルに書き出します R Markdown（R Notebook）を使えば こんな感じで解析結果を保存できます

遺伝⼦情報（遺伝⼦名とか）の取得 BioMart BioMart から遺伝⼦の情報を抜き出してきます Ensembl のツールが使いやすい気がします（ensembl biomart でぐぐれば出てきます） で、解析結果（table_all.csv）の遺伝⼦ ID と BioMart から抜き出してきた遺伝⼦情報 （mart_export.csv）をエクセルの vlookup 関数で紐付けてあげます

vlookup 関数でデータの結合 エクセル BioMart（mart_export.csv） 解析結果（table_all.csv） BioMart の遺伝⼦ ID と解析結果の遺伝⼦ ID を vlookup 関数で 紐付けすれば、遺伝⼦の名前と発現量の変化を紐付けて⾒れます

vlookup 関数でデータの結合 エクセル 2 つのデータをエクセルシートにコピペして vlookup 関数で紐付けすると こんな感じになります わーっとスクロールしながら頑張って考察してください

パスウェイ解析、GO 解析 頑張ってください ⼀通りの作業が終わったら、パスウェイ解析や GO 解析をしていくと思います パスウェイ解析とか GO 解析は頑張れば R でも実⾏できますが、、、 ShinyGO というウェブツールがとても使いやすいです 今回は ここで 終わり de novo の解析については また次回やります