RNA-seq解析シリーズです。
前回までに、
Hisat2によるマッピング、
StringTieによるカウントが終了しましたので、
今回はBallgownによって発現差解析データを作成していきます。
これまでは、MacOSのターミナルやLinuxのコマンドラインを用いて行いましたが、
BallgownからはRを用いていきます。
↓↓これまでの記事はこちら↓↓
lifesciencehack-ai.hatenablog.com
環境
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のサイトを見れば大体書いてあります。
過去記事:超初心者向け！！RNA-seq解析シリーズ③公共データベースからRNA-seqデータをダウンロード
をご参照ください。
csvファイルの作成と保存
今回のデータは、SRR1571967~SRR1571969がControl、
SRR1571970~SRR1571972が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_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 6samples
6サンプル、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でマッピングを行う方法でも記事にしたいと思います。
