【Ballgownの使い方】RNA-Seq解析における発現変動遺伝子の検出
はじめに
次世代シーケンサーを用いてRNA-Seq解析を行うと、FASTQファイルと呼ばれる生データが得られます。生データを参照配列にマッピングした後に、それぞれの遺伝子にマッピングされたリードをカウントすることによって遺伝子発現量を定量していきます。複数のサンプルの遺伝子発現量の定量結果をもとに、グループ間比較を行うことで発現変動遺伝子の検出が行われます。
本ページでは、発現変動遺伝子を検出するソフトウェアであるBallgownの使い方を説明します。RNA-Seqデータ解析の全体の流れを知りたい方はこちらをご覧ください。
インストール
まずはじめに、R がない場合には、R をインストールします。
$ brew install r
Rを起動して、以下を実行します。
> if (!requireNamespace("BiocManager", quietly=TRUE))
> install.packages("BiocManager")
> BiocManager::install("ballgown")
以下を実行して、エラーが表示されなければインストール成功です。
> library(ballgown)
ちなみに、私の環境では上記では以下のようなエラーが出てしまい、インストールできませんでした。
Error: package or namespace load failed for ‘RCurl’ in dyn.load(file, DLLpath = DLLpath, ...)
condaの影響で失敗しているようです。
$ conda deactivate
としてから、再度上記コマンドを実行すると無事インストールできました。
データの準備
StringTieを用いて、Ballgownに読み込むための遺伝子発現量の定量データを作成します。StringTieの詳しい説明はこちらでご覧ください。
上記リンクの手順で解析を行うと、結果は以下のようなフォルダ構成となっています。gtfファイルも含まれていますが、Ballgownの解析で使用されるファイルは*.ctabのみです。
.
└── result
├── sample1
│ ├── e2t.ctab
│ ├── e_data.ctab
│ ├── i2t.ctab
│ ├── i_data.ctab
│ ├── sample1.gtf
│ └── t_data.ctab
├── sample2
│ ├── e2t.ctab
│ ├── e_data.ctab
│ ├── i2t.ctab
│ ├── i_data.ctab
│ ├── sample2.gtf
│ └── t_data.ctab
├── sample3
│ ├── e2t.ctab
│ ├── e_data.ctab
│ ├── i2t.ctab
│ ├── i_data.ctab
│ ├── sample3.gtf
│ └── t_data.ctab
└── sample4
├── e2t.ctab
├── e_data.ctab
├── i2t.ctab
├── i_data.ctab
├── sample4.gtf
└── t_data.ctab
データの読み込み
以下を実行してデータを読み込みます。
> bg = ballgown(dataDir='./result', samplePattern='sample', meas='all')
dataDirで結果を保存してあるディレクトリを指定して、samplePatternでサンプルごとのディレクトリ名を正規表現で指定します。
transcriptごとのFPKMを表示してみます。
> texpr(bg)
FPKM.sample1 FPKM.sample2 FPKM.sample3 FPKM.sample4
1 0.000000 0.000000 0.000000 0.000000
2 0.000000 0.000000 0.000000 0.000000
3 0.032410 0.157659 0.075953 0.000000
4 2.594643 2.855274 0.694965 1.707786
5 0.005810 0.163933 0.846295 0.000000
6 0.000000 0.000000 1.071161 0.065238
...
遺伝子ごとのFPKMを表示してみます。
> gexpr(bg)
FPKM.sample1 FPKM.sample2 FPKM.sample3 FPKM.sample4
gene:ENSG00000000005 0.00000000 0.0000000 0.0000000 0.0000000
gene:ENSG00000000938 0.00000000 0.0000000 0.0000000 0.0000000
gene:ENSG00000001561 0.44382000 0.5909544 0.8668932 0.8514636
gene:ENSG00000002587 0.09133599 0.3305672 1.1535695 0.4534936
gene:ENSG00000002726 0.00000000 0.0000000 0.0000000 0.0000000
gene:ENSG00000002745 0.41149328 0.4678598 0.0000000 0.3937639
...
発現変動遺伝子の抽出
今回は、２群間での比較を行っていきます。以下のようにグループの情報を付加します。
> pData(bg) <- data.frame(id=sampleNames(bg), group=c(0,0,1,1))
以下で検定を実施します。
> stat_results = stattest(bg, feature="transcript", meas="FPKM", covariate="group")
p値でソートして表示してみます。
> stat_results[order(stat_results$pval),]
以下のように結果が得られました。
feature id pval qval
71628 transcript 71628 9.481568e-06 0.4774591
230911 transcript 230911 3.129703e-05 0.4774591
50751 transcript 50751 3.244529e-05 0.4774591
226104 transcript 226104 4.594516e-05 0.4774591
59948 transcript 59948 4.630268e-05 0.4774591
60911 transcript 60911 5.690807e-05 0.4774591
215748 transcript 215748 6.540641e-05 0.4774591
246032 transcript 246032 7.733410e-05 0.4774591
22717 transcript 22717 8.017697e-05 0.4774591
97454 transcript 97454 8.172791e-05 0.4774591
...
論文に必要な解析が簡単にできるRNA-Seqデータ解析ツール
RNA-Seqデータ解析ツールを利用すれば、外部委託や共同研究者への依頼は必要ありません。高スペックなコンピュータの準備やLinuxコマンドの操作も不要ですので、いますぐにご自身で解析できるようになります。
遺伝子発現量の定量、発現変動遺伝子抽出(DEG解析)、Volcano plot描画、MAプロット描画、ヒートマップ描画、GO解析、パスウェイ解析等 を簡単に実施できます。
