・お題：遺伝子発現の多寡を解析するパッケージとして、EdgeR、DESeq2、limmaなどいろいろあるらしい。今回もDESeq2をやってみたい。
・私は素人なので、tutorialをなぞる感じで学んでいきたい。今回のtutorialは以下。正しいことは元サイトを見て頂きたい。
・データセットは以下の文献で得られたものを使用する。酵母の野生型株48 サンプルと、Snf2ノックアウト株の48 サンプルのRNA-Seq で構成されているらしい。
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4754627/pdf/btv425.pdf
・データを読み込む。countTable にカウントデータを、phenoTable にサンプル情報を読み込む。
>library(tidyverse)#とりあえずtidyverseを使えるようにしておく。
> countTable <- read.table("http://jvanheld.github.io/stats_avec_RStudio_EBA/practicals/yeast_2x48_replicates/data/counts.txt", header=TRUE, row.names=1)
> countTable %>% head()
> countTable %>% str()
> phenoTable <- read.table("http://jvanheld.github.io/stats_avec_RStudio_EBA/practicals/yeast_2x48_replicates/data/expDesign.txt", header=TRUE, row.names=1)
・各遺伝子型に色を着ける。
>strainColor <- c("WT"="green","Snf2"="orange")
>phenoTable$color[phenoTable$strain == "Snf"] <- strainColor["Snf2"]
>phenoTable$color[phenoTable$strain == "WT"] <- strainColor["WT"]
> phenoTable %>% head()
> phenoTable %>% str()
・分布をboxplotで見る。
> boxplot(log2(countTable + 1), col=phenoTable$color, pch=".",
+ horizontal=TRUE, cex.axis=0.5,
+ las=1, ylab="Samples", xlab="log2(counts +1)")
・density plotでも見てみる。affyというパッケージがデータを探索的に視覚化するのに便利らしいので、使ってみる（以下を参照。インストールは割愛）。
・描いてみる。
> library(affy)
> plotDensity(log2(countTable + 1), lty=1, col=phenoTable$color, lwd=2)
> grid()
> legend("topright", legend=names(strainColor), col=strainColor, lwd=2)
・すべてのサンプルで共通に発現が低い遺伝子が結構ある。今回はすべてのサンプルで発現ゼロの遺伝子を除く。
> countTable <- countTable[rowSums(countTable)!=0,]
> countTable %>% str()
・6887個の遺伝子が残ったらしい。次に、正規化する。DESeq2の機能をそのまま使わせてもらう。
> dds.norm <- estimateSizeFactors(dds0)
・せっかくなので、正規化前後で分布を比較してみる。左半分が正規化前、右半分が正規化後、上半分で箱ひげ図を、下半分で密度分布を示してみた。
> boxplot(log2(counts(dds.norm, normalized=TRUE)+1), col=col.pheno.selected, cex.axis=0.7,
+ las=1, xlab="log2(normalized counts)", horizontal=TRUE, main="Normalized counts")
> par(mfrow=c(2,2),cex.lab=0.7)
> boxplot(log2(counts(dds.norm)+1), col=col.pheno.selected, cex.axis=0.7,
+ las=1, xlab="log2(counts+1)", horizontal=TRUE, main="Raw counts")
> boxplot(log2(counts(dds.norm, normalized=TRUE)+1), col=col.pheno.selected, cex.axis=0.7,
+ las=1, xlab="log2(normalized counts)", horizontal=TRUE, main="Normalized counts")
> plotDensity(log2(counts(dds.norm)+1), col=col.pheno.selected,
+ xlab="log2(counts+1)", cex.lab=0.7, panel.first=grid())
> plotDensity(log2(counts(dds.norm, normalized=TRUE)+1), col=col.pheno.selected,
+ xlab="log2(normalized counts)", cex.lab=0.7, panel.first=grid())
・確かになんだか揃っている。
つづく。
