超初心者向け！！RNA-seq解析シリーズ⑥ Ballgownで発現差解析

Q: Ballgownがインストールできない

install.packages("ballgown") ではインストールできません。BallgownはCRANではなくBioconductorのパッケージなので、BiocManager::install("ballgown") を使います。それでも失敗する場合は、RのバージョンとBioconductorのバージョンが対応しているかを確認してください。古いRのままだと、対応する古いパッケージが入るか、そもそも解決に失敗します。

Q: stattest の結果はどう読めばいい？

qval（FDR補正後のq値）でまず絞り込み、次に fc（fold change）で効果量を見ます。pval で絞ってはいけません。数万の転写産物を同時に検定しているため、p < 0.05 の中には偶然有意になったものが大量に含まれます。qval < 0.05 かつ fold change が十分大きいものを候補にするのが標準的な進め方です。

Q: 低発現の転写産物はなぜ削除するの？

発現量がほとんどゼロの転写産物は、わずかなカウントのゆらぎで極端なfold changeを示すため、偽陽性の温床になります。また検定数が減ることで多重比較補正が緩くなり、本当に見たいシグナルが埋もれにくくなります。このシリーズでは分散でフィルタしていますが、最低発現量でフィルタする方法もあります。

Q: BallgownとDESeq2・edgeRはどう使い分ける？

Ballgownは、StringTieが推定したFPKMなどの連続値を入力として、転写産物レベルの発現差を検定します。新規転写産物やアイソフォーム単位の解析に向いています。一方DESeq2やedgeRは、遺伝子ごとの生カウント値を入力とし、負の二項分布を仮定して検定します。遺伝子レベルの発現差解析では現在こちらが標準的で、サンプル数が少ない場合の頑健性でも優れています。目的がアイソフォームでなく遺伝子レベルの発現変動であれば、featureCountsなどでカウントを取ってDESeq2に渡す構成が一般的です。

Q: サンプル情報csvが読み込めない

ballgown() の dataDir と samplePattern を確認してください。samplePattern はサンプルフォルダ名の頭にある共通文字列で、この例では "SRR" です。また data_info.csv の ids 列の並び順とフォルダ名が対応している必要があります。作業ディレクトリがballgownフォルダの親になっているかも、つまずきやすいところです。 Google Colaboratoryを使えば、Macやターミナルの環境構築なしで同じ解析を進められます。そちらの手順は別シリーズにまとめています。

RNA-seq解析シリーズです。
前回までに、
Hisat2によるマッピング、
StringTieによるカウントが終了しましたので、
今回はBallgownによって発現差解析データを作成していきます。

これまでは、MacOSのターミナルやLinuxのコマンドラインを用いて行いましたが、
BallgownからはRを用いていきます。

環境

RはRStudioで使っています。
Rのバージョンは以下のとおりです。
Rのインストールがまだの方はインストールしてください。
すでにインストール済の方もできるだけ最新版にしましょう。

version

platform       x86_64-apple-darwin17.0     
arch           x86_64                      
os             darwin17.0                  
system         x86_64, darwin17.0          
status                                     
major          4                           
minor          0.2                         
year           2020                        
month          06                          
day            22                          
svn rev        78730                       
language       R                           
version.string R version 4.0.2 (2020-06-22)
nickname       Taking Off Again

Ballgownのインストール

BiocManagerを使ってBallgownをインストールします。

#R 3.5以降の方はBiocManagerを使う
install.packages("BiocManager")
BiocManager::install("ballgown")

#途中でUpdate all/some/none?と聞かれたら、「a」を入力してenter

install.packages('metaMA')
library(metaMA)

install.packages('tidyr')
library(tidyr)

install.packages('dplyr')
library(dplyr)

以下の入力してエラーが出なければ、無事インストール完了です。

library(ballgown)

サンプル情報ファイルの作成

サンプル情報の確認

前回までの記事の内容で、StringTieでBallgown用のファイルを作成しました。
Ballgownを実行し、発現解析データを取得するためには、サンプル情報を与えて上げる必要があります。
サンプルの情報は、RNA-seqデータのダウンロード元のDDBJのサイトを見れば大体書いてあります。

csvファイルの作成と保存

今回のデータは、SRR1571967SRR1571969がControl、
SRR1571970SRR1571972がHIF1 mutantでしたので、
↓のような表をcsvで作ります。

ids	type
SRR1571967	Control
SRR1571968	Control
SRR1571969	Control
SRR1571970	Hif1 mutant
SRR1571971	Hif1 mutant
SRR1571972	Hif1 mutant

サンプル情報ファイルの保存場所ですが、
StringTieで作ったballgownフォルダの直下に作ります。
↓のような構造となっているかと思います。

ballgown
├ SRR1571967
│     ├ e_data.ctab
│     ├ e2t.ctab
│     ├ i_data.ctab
│     ├ i2t.ctab
│     ├ SRR1571967.gtf
│     └ t_data.ctab
├ SRR1571968
├ SRR1571969
├ SRR1571970
├ SRR1571971
└ SRR1571972
#SRR1571967だけ全てのファイルを示していますが、 
#全てのフォルダで同様のファイルがあります。

ここにdata_info.csvという名前のサンプル情報ファイルを作成します。

ballgown
├ data_info.csv #このファイルを作る
├ SRR1571967
├ SRR1571968
├ SRR1571969
├ SRR1571970
├ SRR1571971
└ SRR1571972

↑こんな感じになる予定です。
今回はRを使って作ってみます。

作業ディレクトリの変更

まずは作業ディレクトリを変更します。
ballgownフォルダの直下にサンプル情報ファイルを保存するため、
ballgownフォルダに移動します。
作業ディレクトリの移動はsetwd(パス)を用います。

getwd()   # 現在の作業ディレクトリの確認
# [1] "/Users/***********/"   # 現在の作業ディレクトリが表示される
setwd("***********/ballgown")      # ballgownフォルダを作業ディレクトリにする
getwd()     # 作業ディレクトリが変更できたことを確認
# [1] "/Users/***********/ballgown"

ファイル作成&保存

#ids列
ids <- c('SRR1571967', 'SRR1571968', 'SRR1571969', 'SRR1571970', 'SRR1571971', 'SRR1571972')

#type列
type <- c('Control', 'Control', 'Control', 'Hif1 mutant', 'Hif1 mutant', 'Hif1 mutant')

#DataFrameの作成
data <- data.frame('ids'=ids, 'type'=type)

#csvファイルとして保存
write.csv(data, 'data_info.csv', row.names = FALSE)

Ballgownの実行

実際にballgownを実行し、発現差データを取得します。

作業ディレクトリの移動とサンプル情報の読み込み

ballgownフォルダの親フォルダに移動します。
ballgownフォルダの一つ上のフォルダであることを気をつけてください。

setwd("ballgownフォルダの親フォルダ")
data_info <- read.csv("ballgown/data_info.csv")

ballgown()の実行と読み込み結果

ballgwonを実行していきます。
ballgown()で実行できます。
引数の説明ですが
dataDir ballgownフォルダの名前
pData サンプル情報オブジェクトを指定
samplePatter サンプル名の頭3文字を指定

bg <- ballgown(dataDir = "ballgown", pData=data_info, samplePattern = "SRR")
bg
ballgown instance with 117257 transcripts and 6samples

6サンプル、117257の転写産物が読み込まれました。

bg_cut <- subset(bg,"rowVars(texpr(bg)) >1",genomesubset=TRUE)
bg_cut
ballgown instance with 12789 transcripts and 6 samples

低発現の転写産物を削除すると、117257 → 12789 と、およそ1/10まで絞り込まれました。

rowVars(texpr(bg)) > 1 は「サンプル間で発現量の分散が1より大きい転写産物だけを残す」という条件です。ほとんど発現していない転写産物や、全サンプルでほぼ同じ値の転写産物を先に除いておくと、この後の検定で多重比較の負担が減り、意味のあるシグナルが埋もれにくくなります。

transcript_data <- stattest(bg_cut, 
                                feature = 'transcript', 
                                covariate = 'type',
                                getFC = 'TRUE',
                                meas = 'FPKM')

transcript_table <- data.frame(gene = geneNames(bg_cut),
                                  gene_id = geneIDs(bg_cut),
                                  transcript_data)

gene_data <- stattest(bg_cut, 
                        feature = 'gene', 
                        covariate = 'type',
                        getFC = 'TRUE',
                        meas = 'FPKM')

stattest() の引数は次の意味です。

feature 検定の単位。'transcript' なら転写産物ごと、'gene' なら遺伝子ごとに検定します
covariate 群分けに使う列名。ここでは data_info.csv の type 列（Control / Hif1 mutant）
getFC = 'TRUE' fold change を結果に含める
meas = 'FPKM' 検定に使う発現量の指標

結果の見方

stattest() の出力には次の列が入ります。

列	意味
`feature`	検定の単位（transcript / gene）
`id`	転写産物IDまたは遺伝子ID
`fc`	fold change（`getFC = 'TRUE'` のときのみ）
`pval`	検定のp値
`qval`	多重比較補正後のq値（FDR）

絞り込みに使うのは pval ではなく qval です。数万の転写産物を同時に検定しているため、p < 0.05 だけでは偶然の当たりが大量に混ざります。慣例的には qval < 0.05 を目安にし、さらに fc で効果量の大きいものに絞ります。

write.csv(transcript_table, "transcript_results.csv", row.names=FALSE)
write.csv(gene_data, "gene_results.csv", row.names=FALSE)

これで、長きに渡ったRNA-seq解析記事は終了となります。

よくある質問

Q. Ballgownがインストールできない install.packages("ballgown") ではインストールできません。BallgownはCRANではなくBioconductorのパッケージなので、BiocManager::install("ballgown") を使います。それでも失敗する場合は、RのバージョンとBioconductorのバージョンが対応しているかを確認してください。古いRのままだと、対応する古いパッケージが入るか、そもそも解決に失敗します。

Q. stattest の結果はどう読めばいい？ qval（FDR補正後のq値）でまず絞り込み、次に fc（fold change）で効果量を見ます。pval で絞ってはいけません。数万の転写産物を同時に検定しているため、p < 0.05 の中には偶然有意になったものが大量に含まれます。qval < 0.05 かつ fold change が十分大きいものを候補にするのが標準的な進め方です。

Q. 低発現の転写産物はなぜ削除するの？ 発現量がほとんどゼロの転写産物は、わずかなカウントのゆらぎで極端なfold changeを示すため、偽陽性の温床になります。また検定数が減ることで多重比較補正が緩くなり、本当に見たいシグナルが埋もれにくくなります。このシリーズでは分散でフィルタしていますが、最低発現量でフィルタする方法もあります。

Q. BallgownとDESeq2・edgeRはどう使い分ける？ Ballgownは、StringTieが推定したFPKMなどの連続値を入力として、転写産物レベルの発現差を検定します。新規転写産物やアイソフォーム単位の解析に向いています。一方DESeq2やedgeRは、遺伝子ごとの生カウント値を入力とし、負の二項分布を仮定して検定します。遺伝子レベルの発現差解析では現在こちらが標準的で、サンプル数が少ない場合の頑健性でも優れています。目的がアイソフォームでなく遺伝子レベルの発現変動であれば、featureCountsなどでカウントを取ってDESeq2に渡す構成が一般的です。

Q. サンプル情報csvが読み込めない ballgown() の dataDir と samplePattern を確認してください。samplePattern はサンプルフォルダ名の頭にある共通文字列で、この例では "SRR" です。また data_info.csv の ids 列の並び順とフォルダ名が対応している必要があります。作業ディレクトリがballgownフォルダの親になっているかも、つまずきやすいところです。

Google Colaboratoryを使えば、Macやターミナルの環境構築なしで同じ解析を進められます。そちらの手順は別シリーズにまとめています。