RNA-seq解析シリーズです。
前回までに、
Hisat2によるマッピング、
StringTieによるカウントが終了しましたので、
今回はBallgownによって発現差解析データを作成していきます。
これまでは、MacOSのターミナルやLinuxのコマンドラインを用いて行いましたが、
BallgownからはRを用いていきます。
環境
RはRStudioで使っています。
Rのバージョンは以下のとおりです。
Rのインストールがまだの方はインストールしてください。
すでにインストール済の方もできるだけ最新版にしましょう。
version
<div></div>
'''
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")
<div></div>
#途中でUpdate all/some/none?と聞かれたら、「a」を入力してenter
<div></div>
install.packages('metaMA')
library(metaMA)
<div></div>
install.packages('tidyr')
library(tidyr)
<div></div>
install.packages('dplyr')
library(dplyr)以下の入力してエラーが出なければ、無事インストール完了です。
library(ballgown)サンプル情報ファイルの作成
サンプル情報の確認
前回までの記事の内容で、StringTieでBallgown用のファイルを作成しました。
Ballgownを実行し、発現解析データを取得するためには、サンプル情報を与えて上げる必要があります。
サンプルの情報は、RNA-seqデータのダウンロード元のDDBJのサイトを見れば大体書いてあります。
csvファイルの作成と保存
今回のデータは、SRR1571967SRR1571969がControl、SRR1571972がHIF1 mutantでしたので、
SRR1571970
↓のような表を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_ballgown.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 6samples6サンプル、117257の遺伝子が読み込まれました。
bg_cut <- subset(bg,"rowVars(texpr(bg)) >1",genomesubset=TRUE)
bg_cut
ballgown instance with 12789 transcripts and 6 samples低発現の転写物を削除すると、117257→127891と転写物の数が1/10になりました。
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')write.csv(transcript_table, "transcript_results.csv", row.names=FALSE)
write.csv(gene_data, "gene_results.csv", row.names=FALSE)これで、長きに渡ったRNA-seq解析記事は終了となります。
需要があるかわかりませんが、気が向いたらGoogle colaboratoryでマッピングを行う方法でも記事にしたいと思います。