Bash on Ubuntu on WindowsでSingle-cell RNA sequence解析 ⑥Rによる発現解析

※記事を読んでくださるみなさまへ

・対象はプログラミングもバイオインフォ解析も全く未経験のバイオ研究者です
・筆者は基礎研究を初めて約1年3ヶ月の大学院生です
・いわゆる”DRY”な研究を始めたのは約5ヶ月前からです
・そのため記事の内容にはとんでもない間違いやずっと良い別の方法が存在するかもしれません
・また、筆者と同じレベルの人が0から始めるときの参考になるよう、説明は結構まだるっこいです
・ご指摘の度に記事をブラッシュアップしていきたいです、よろしくお願いします

（新規注意）
・今回、統計学的内容やRの説明がやや薄いかもしれません。筆者がもともと臨床研究でRをある程度使い慣れているのが原因です。下記のサイトやGoogle検索を駆使して、コードの記載法や各種関数の機能、統計について適宜調べてください。

Google 無敵のGoogleさん。"R　〇〇"という形で検索。
Rseek 本家本元R専用検索サイト。
seekR こちらのほうがやや日本語ページがヒットしやすい印象。
RjpWiki 日本語のR情報まとめ。必要な情報がやや見つけにくい。一番下に検索バーがある。
統計ソフトRの「困った」を解決する12（+α）の方法 素晴らしいRの情報源まとめ記事。

過去記事一覧
①環境構築まで
②NGSデータのダウンロードとクオリティチェック
③NGSデータのトリミング
④Pipeline (hisat2, samtools, stringtie)
⑤番外編：コードのブラッシュアップ
⑥Rによる発現解析

RStudioの導入、必要なパッケージの導入

Rとはなにか、RStudioとはなにかについては先人の素晴らしいまとめがありますのでそちらを参考にしてRの環境を導入しましょう。準備ができたらbashの画面より

mv ballgown/ERR* /mnt/c/Users/(ユーザー名)/Document/analyze/ballgown

などのようにして必要なデータをWindows側の領域に持ってきておきましょう。事前に解析用のanalyzeのようなフォルダを作成しておいてください。なぜかはよくわかりませんが、Linux側のballgownフォルダをまるごとWindows側に持ってこようとしてもうまく行きませんでした。

同様にしてLinux側のChrX_dataディレクトリのgeuvadis_phenodata.csvファイルも持ってきてanalyzeフォルダにおきましょう。

次にRStudioを起動し、左上のメニューから

File→New project→Existing directory

と選択し、Windows側のanalyzeフォルダを選択しましょう。このプロジェクトの作成により、コードや解析結果が自動的に保存されるので便利です。

ここでR本体に加え、解析に必要なパッケージをインストール、可用化します。

install.packages("devtools")
devtools::install_github("vqv/ggbiplot")

source("https://bioconductor.org/biocLite.R")
biocLite("genefilter")

library(ggbiplot)
library(devtools)
library(ballgown)
library(RSkittleBrewer)
library(genefilter)
library(dplyr)

ballgownファイルの作成

次に必要なデータの読み込みと、ballgownファイルを作成します。

# フェノタイプデータの読み込み、オブジェクト化（pheno_data）
pheno_data = read.csv("ChrX_data/geuvadis_phenodata.csv")

# ballgownデータ（bg_chrX）の作成
# dataDirは各サンプルのgtf用フォルダが入っている、最初に作ったフォルダ
# samplePatternは各サンプルのフォルダの共通する頭文字
# pDataはさっき読み込んだフェノタイプデータ（""をつけない）
bg_chrX = ballgown(dataDir = "ballgown", samplePattern = "ERR", pData=pheno_data)

フェノタイプデータは人種、性別などの情報が入っています。自前データの解析のときにはエクセルでサンプルのデータ（臓器なのか、組織なのか、培養細胞なのかなど）を作成してcsvで保存しておくと良いでしょう。Rではいちいち元データを参照するのではなく、”オブジェクト”という属性情報なども含めた形でデータを扱っていきます。

個別の解析

次に低発現な遺伝子を除外したデータセットのオブジェクトを作成します。

bg_chrX_filt = subset(bg_chrX,"rowVars(texpr(bg_chrX)) >1",genomesubset=TRUE)

# rowVars()が無いと怒られるときはlibrary(matrixStats)を加える

subset()関数は本来データフレームから条件に合う部分を抜き出して新たなデータフレームを作る関数。ballgownパッケージに全く同名称の関数として用意されており、genomesubset=TRUEにより条件にあった遺伝子・転写産物のみを取り出す。texpr()はballgownファイルから転写産物の発現量（デフォルトでFPKM）を取り出す。rowVarsは行の分散を取り出す。なので統計的にはサンプル間のバラツキが少ない遺伝子を除去していると思うのだが、元論文には低発現の遺伝子を除去、と書いている…。また、本来オブジェクトを作るときは<-が推奨だと思うのだが…。とりあえず元論文に従ってやりましょう。Nature様に逆らってはいけません。

次にフェノタイプデータに書いてある人種を交絡因子とした、性別により発現に差がある転写産物および遺伝子を抽出します。

# 先の条件を満たす転写産物および遺伝子の抽出,statttest()関数
# getFCはFold changeも計算する、というオプション
# measは発現量の指標は〇〇にしなさい、というオプション
results_transcripts = stattest(bg_chrX_filt, feature="transcript",covariate="sex",adjustvars = c("population"), getFC=TRUE, meas="FPKM")
results_genes = stattest(bg_chrX_filt, feature="gene", covariate="sex", adjustvars = c("population"), getFC=TRUE, meas="FPKM")

# 遺伝子名、遺伝子IDなどと結合した表を作る
results_transcripts = data.frame(geneNames=ballgown::geneNames(bg_chrX_filt), geneIDs=ballgown::geneIDs(bg_chrX_filt), results_transcripts)

# p値ごとに並び替える
results_transcripts = arrange(results_transcripts,pval)
results_genes = arrange(results_genes,pval)

p値のヒストグラムを眺めることで、一様分布であれば全くのランダム、0付近に一部固まるようであればサンプルに差がある、ということがわかります。これは後に出てくるq値（Tukey法）の理論的基礎でもあります。

結果のファイルをcsv出力してみましょう。

write.csv(results_transcripts, "chrX_transcript_results.csv", row.names=FALSE)
write.csv(results_genes, "chrX_gene_results.csv", row.names=FALSE)

chrX_transcript_results.csvなどは出力するファイル名です。row.namesオプションは気にせず書かれているとおりにつけてください。

では、発現変動が統計学的に有意な遺伝子をq値という指標を使って取り出します。バイオインフォでは、何千何万という数の遺伝子や転写産物について統計検定が行われます。統計に詳しい方は知っているかもしれませんが、このような多重検定において通常のp<0.05という判断基準は意味を持ちません。そのため、普段の検定と同じように使える指標としてq値というものが開発されています。詳細はここにわかりやすく書かれていますので参考にしてください。ちなみにこの部分は元論文からちょっと改変してあります。

significant_results_transcripts<-subset(results_transcripts,results_transcripts$qval<0.05)
significant_results_genes<-subset(results_genes,results_genes$qval<0.05)

print(significant_results_transcripts)
print(significant_results_genes)

また、元論文ではサンプルごとのfpkmの分布を見ています。

# fpkmを抜き出し、扱いやすいログに変換
fpkm = texpr(bg_chrX, meas="FPKM")
fpkm = log2(fpkm+1)

# 箱ひげ図の作図
boxplot(fpkm,col=as.numeric(pheno_data$sex),las=2,ylab='log2(FPKM+1)')

細かい作図のTipsについてはRグラフィックスクックブックなどを見てください。
さらにNM_012227という転写産物につき、発現量の男女差をみています。

# リストの12番めを確認
ballgown::transcriptNames(bg_chrX)[12]
ballgown::geneNames(bg_chrX)[12]

# 性別ごとの箱ひげ図を作成
plot(fpkm[12,] ~ pheno_data$sex, border=c(1,2),
main=paste(ballgown::geneNames(bg_chrX)[12],' : ',
ballgown::transcriptNames(bg_chrX)[12]),pch=19, xlab="Sex", ylab='log2(FPKM+1)')
points(fpkm[12,] ~ jitter(as.numeric(pheno_data$sex)),
col=as.numeric(pheno_data$sex))

最後にERR188234のサンプルでのXIST遺伝子転写産物の構造と発現量一覧を可視化しています。

plotTranscripts(ballgown::geneIDs(bg_chrX)[1729], bg_chrX, main=c('Gene XIST in
sample ERR188234'), sample=c('ERR188234'))
plotMeans('MSTRG.56', bg_chrX_filt,groupvar="sex",legend=FALSE)

以上がサンプルデータを使った可視化の練習でした。ですが、実際に自分の解析に活かすには今回の記事だけでは不十分かと思います。次回以降は筆者のSingle-cell RNA seqデータをどのように解析したかの概略を示そうかと思います。