---
title: mRNA-Seq入門【第四回】
tags:  #バイオインフォマティクス #mrna-seq #edger  
author: [膝栗毛之男](https://convert.docswell.com/user/hizaaakurige)
site: [Docswell](https://www.docswell.com/)
thumbnail: https://bcdn.docswell.com/page/54JM956QEW.jpg?width=480
description: mRNA-Seq の発現量解析についてまとめたもの。edgeR を使って二群間比較をする。
published: December 10, 21
canonical: https://convert.docswell.com/s/hizaaakurige/K7E4GZ-mRNASeq-4
---
# Page. 1

![Page Image](https://bcdn.docswell.com/page/54JM956QEW.jpg)

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


# Page. 2

![Page Image](https://bcdn.docswell.com/page/ZPJRG3PKE9.jpg)

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


# Page. 3

![Page Image](https://bcdn.docswell.com/page/NPEXX435EX.jpg)

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


# Page. 4

![Page Image](https://bcdn.docswell.com/page/63EKWZYRED.jpg)

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


# Page. 5

![Page Image](https://bcdn.docswell.com/page/4L731D9GJ5.jpg)

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


# Page. 6

![Page Image](https://bcdn.docswell.com/page/V87DXPGNEG.jpg)

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


# Page. 7

![Page Image](https://bcdn.docswell.com/page/4VJPPN3NJ8.jpg)

が


# Page. 8

![Page Image](https://bcdn.docswell.com/page/L2EV2K4YEQ.jpg)

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


# Page. 9

![Page Image](https://bcdn.docswell.com/page/Z57GRQ1WEL.jpg)

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


# Page. 10

![Page Image](https://bcdn.docswell.com/page/V4EQVKDQJP.jpg)

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


# Page. 11

![Page Image](https://bcdn.docswell.com/page/DKJ4MPZYE1.jpg)

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


# Page. 12

![Page Image](https://bcdn.docswell.com/page/PLE189RN7G.jpg)

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


# Page. 13

![Page Image](https://bcdn.docswell.com/page/8GEWZQYME2.jpg)

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


# Page. 14

![Page Image](https://bcdn.docswell.com/page/W47Z1WXMJ3.jpg)

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


# Page. 15

![Page Image](https://bcdn.docswell.com/page/KYJ6LG45EV.jpg)

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


# Page. 16

![Page Image](https://bcdn.docswell.com/page/5GJ51VQG74.jpg)

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


# Page. 17

![Page Image](https://bcdn.docswell.com/page/V9E213PDJR.jpg)

マニュアルを


# Page. 18

![Page Image](https://bcdn.docswell.com/page/5D7YZ15MJM.jpg)

しっかりと


# Page. 19

![Page Image](https://bcdn.docswell.com/page/2VEN32N3E8.jpg)

読んで
ください


# Page. 20

![Page Image](https://bcdn.docswell.com/page/GY7991RP73.jpg)

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


# Page. 21

![Page Image](https://bcdn.docswell.com/page/5G7894WXJD.jpg)

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


# Page. 22

![Page Image](https://bcdn.docswell.com/page/ML7LW6NNER.jpg)

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


# Page. 23

![Page Image](https://bcdn.docswell.com/page/54EM95XQ7W.jpg)

データセット作成
edgeR


# Page. 24

![Page Image](https://bcdn.docswell.com/page/ZPERG3NK79.jpg)

基本的な⼆群間⽐較
edgeR（3.32.1）
edgeR を読み込みます
&gt; library(edgeR)
カウントデータを読み込んで、head で読み込んだデータを確認します
&gt; count &lt;- read.table(&quot;count_data.tsv&quot;, header = T, row.names = 1)
&gt; 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 に気をつけて）
&gt; group &lt;- factor(rep(c(&#039;C&#039;, &#039;A&#039;), c(4, 4)), levels = c(&#039;C&#039;, &#039;A&#039;))


# Page. 25

![Page Image](https://bcdn.docswell.com/page/NP7XX4N57X.jpg)

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


# Page. 26

![Page Image](https://bcdn.docswell.com/page/637KWZNR7D.jpg)

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


# Page. 27

![Page Image](https://bcdn.docswell.com/page/4LJ31DVGE5.jpg)

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


# Page. 28

![Page Image](https://bcdn.docswell.com/page/V8JDXP8N7G.jpg)

発現量⽐較
edgeR


# Page. 29

![Page Image](https://bcdn.docswell.com/page/4VEPPN8NE8.jpg)

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


# Page. 30

![Page Image](https://bcdn.docswell.com/page/L27V2KNY7Q.jpg)

正確確率検定の場合
edgeR（3.32.1）
正確確率検定で⼆群を⽐較します（詳細はマニュアル参照）
&gt; et &lt;- exactTest(y)
&gt; 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


# Page. 31

![Page Image](https://bcdn.docswell.com/page/Z5JGRQKW7L.jpg)

疑似尤度 F 検定の場合
edgeR（3.32.1）
疑似尤度 F 検定で⼆群を⽐較します（詳細はマニュアル参照）
&gt; fit &lt;- glmQLFit(y, design)
&gt; qlf &lt;- glmQLFTest(fit)
&gt; 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
多群⽐較にも応⽤可能です（詳細はマニュアル参照）


# Page. 32

![Page Image](https://bcdn.docswell.com/page/V47QVKNQEP.jpg)

尤度⽐検定の場合
edgeR（3.32.1）
尤度⽐検定で⼆群を⽐較します（詳細はマニュアル参照）
&gt; fit &lt;- glmFit(y, design)
&gt; lrt &lt;- glmLRT(fit)
&gt; 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
多群⽐較にも応⽤可能です（詳細はマニュアル参照）


# Page. 33

![Page Image](https://bcdn.docswell.com/page/DKE4MPGY71.jpg)

発現量の定量
edgeR


# Page. 34

![Page Image](https://bcdn.docswell.com/page/PL7189DNJG.jpg)

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


# Page. 35

![Page Image](https://bcdn.docswell.com/page/8G7WZQ8M72.jpg)

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


# Page. 36

![Page Image](https://bcdn.docswell.com/page/W4JZ1W8ME3.jpg)

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


# Page. 37

![Page Image](https://bcdn.docswell.com/page/KYE6LGP57V.jpg)

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


# Page. 38

![Page Image](https://bcdn.docswell.com/page/5GE51VKGJ4.jpg)

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


# Page. 39

![Page Image](https://bcdn.docswell.com/page/V97213WDER.jpg)

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


# Page. 40

![Page Image](https://bcdn.docswell.com/page/5DJYZ1LMEM.jpg)

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


# Page. 41

![Page Image](https://bcdn.docswell.com/page/2V7N324378.jpg)

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


# Page. 42

![Page Image](https://bcdn.docswell.com/page/GYJ991QPJ3.jpg)

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


