・お題：遺伝子発現の多寡を解析するパッケージとして、EdgeR、DESeq2、limmaなどいろいろあるらしい。今回はDESeq2をやってみたい。
・私は素人なので、tutorialをなぞる感じで学んでいきたい。今回のtutorialは以下。正しいことは元サイトを見て頂きたい。
・とりあえず、以下の公式サイトを参考に、DESeq2をインストールする（割愛）。
・次に、ライブラリを読み込む。
>library("tidyverse")
>library("DESeq2")
・今回は、以下のサイトで提供されているサンプルデータを用いる。airway_scaledcounts.csvairway_scaledcounts.csvというファイルをダウンロードし、カレントディレクトリに保存した。
・次に、データを読み込む。このファイルには、遺伝子のID、サンプル名およびリード数（？）が入っているらしい。
> countData <- read.csv('airway_scaledcounts.csv', header = TRUE, sep = ",")
> head(countData)
・サンプル情報も同様に、airway_metadata.csvをダウンロードし、読み込む。
> metaData <- read.csv('airway_metadata.csv', header = TRUE, sep = ",")
> metaData
・次に、これらの情報をDESeq2で取り扱うデータ形式に変換する。
> dds <- DESeqDataSetFromMatrix(countData=countData,
+ colData=metaData,
+ design=~dex, tidy = TRUE)
>dds
・次に、ddsをDESeqする。何をやっているのかというと、各サンプルの相対的なライブラリの深さを計算（正規化）、各遺伝子のカウントの分散を推定、サイズと分散の出力を用いて、Negative Binomial GLMにおける係数の有意性を計算、しているらしい。この一連の処理が一つの関数で済むとは。。
> dds <- DESeq(dds)
・次に、結果を見てみる。log2FCやpvalueなどを確認することができる。
> res <- results(dds)
> res_sorted <- res[order(res$padj),]
> res_sorted %>% head()
> res %>% summary()
・変動の激しかった遺伝子に関して、発現量をプロットしてみる。
・Volcano Plotを描きたいと思って調べていたら、EnhancedVolcanoなるパッケージを見つけた。インストールして使ってみた（インストールは割愛）。
>library(EnhancedVolcano)
> EnhancedVolcano(res,
+ lab = rownames(res),
+ x = 'log2FoldChange',
+ y = 'pvalue',
+ title = 'CTRL vs TREATMENT',
+ pCutoff = 10e-32,
+ FCcutoff = 2,
+ pointSize = 3.0,
+ labSize = 3.0)
・MDSプロットはGlimmaを使って可視化してみる。。
>library("Glimma")
> glMDSPlot(dds)
・もう少し他の例もやってみたいと思った。
とりあえず、おわり。
