・お題：遺伝子発現の多寡を解析するパッケージとして、EdgeR、DESeq2、limmaなどいろいろあるらしい。今回もDESeq2をやってみたい。
・私は素人なので、tutorialをなぞる感じで学んでいきたい。今回のtutorialは以下。正しいことは元サイトを見て頂きたい。
・前回countデータの正規化までやったので、今回は比較をしていきたい。DESeq2では負の二項分布というやつを仮定するらしい。分散と平均の関係が厳密には直線的ではないらしく、ポアソン分布を仮定するのはよろしくないということで、負の二項分布を仮定するようになったみたい。
・ライブラリのアクチベートなどは前回に実施した通りなので割愛。正規化済みの解析用データをdds.normと名付けて処理を進める。
・まず、分散パラメータを推定する。次に、「実際の分散と発現量の関係」と「推定した分散と発現量の関係」を一つのグラフに図示してみる。
> dds.disp <- estimateDispersions(dds.norm)#分散パラメータを推定
> plotDispEsts(dds.disp)#図示
・fittingまで終わったので、この後、Wald検定で遺伝子発現変動の検定をしていくことになる。Wald検定は以下の記事を参照。統計の勉強が必要。。
・やってみる。
> alpha <- 0.0001
> waldTestResult <- nbinomWaldTest(dds.disp)
> resultDESeq2 <- results(waldTestResult, alpha=alpha, pAdjustMethod="BH")
> resultDESeq2
・DESeq2の最初の記事で一つの関数で算出されていた結果が出てきた。本来はこういうステップを経て出てくるものだったらしい。。
・せっかくなので、p値の分布をみてみる。
h1 <- hist(resultDESeq2$padj, breaks=20, col="grey", main="DESeq2 adj-p-value distribution", xlab="Adj-P-value", ylab="Number of genes")
・この後に、volcano plotを描く。volcano plotは、横軸にlog2(fold-change)を、縦軸に-log10(adjusted p-value)をとって、各遺伝子をプロットした散布図のこと。EnhancedVolcanoパッケージを使って描く。
>library(EnhancedVolcano)
> EnhancedVolcano(resultDESeq2,
+ lab = rownames(resultDESeq2),
+ x = "log2FoldChange",
+ y = "padj",
+ xlab = "Effect size: log2(fold-change)",
+ ylab = "-log10(adjusted p-value)",
+ pCutoff = 1*10^(-100),
+ FCcutoff = 1.0,
+ pointSize = 2.0)
Warning message:
One or more p-values is 0. Converting to 10^-1 * current lowest non-zero p-value...と出てきた。どうもp-valueが0のものがあるらしい。最小値の10分の1のp値を代入してくれるらしい。助かる。
・ラベルが被ってしまい、プロットが見づらい。ラベルサイズをゼロにしてみる。
> EnhancedVolcano(resultDESeq2,
+ lab = rownames(resultDESeq2),
+ x = "log2FoldChange",
+ y = "padj",
+ xlab = "Effect size: log2(fold-change)",
+ ylab = "-log10(adjusted p-value)",
+ pCutoff = 1*10^(-100),
+ FCcutoff = 1.0,
+ pointSize = 2.0,
+ labSize = 0)
・元サイトでは、この後個別の遺伝子の発現を棒グラフで見ていく感じだけれど、インタラクティブに遺伝子発現を確認できる方が便利だと思う。ということで、Glimmaパッケージでvolcano plotからインタラクティブに遺伝子発現を見れるようにしてみた。
> #CutOff
> pCutoff <- 1*10^(-100)
> FCcutoff <- 1.0
>
> #volcano plot color
> resultDESeq2$col<-0#0はグレー
> resultDESeq2$col[(resultDESeq2$log2FoldChange > FCcutoff)&(resultDESeq2$padj < pCutoff)] <- 1#1は赤
> resultDESeq2$col[(resultDESeq2$log2FoldChange < -FCcutoff)&(resultDESeq2$padj < pCutoff)] <- -1#-1は水色
>
> glXYPlot(x = resultDESeq2$log2FoldChange,
+ y = -log10(resultDESeq2$padj),
+ xlab = "log2FC",
+ ylab = "-log10(adjusted p-value)",
+ counts = counts(dds.norm, normalized=TRUE), groups = phenoTable$strain,
+ status = resultDESeq2$col,sample.cols = phenoTable$color#statusはvolvcano plotの色、samples.colsは遺伝子発現プロットの色。
+ )
・インタラクティブに遺伝子発現を確認できるようになったぽい。
・次に、MA plotを描いてみる。
> DESeq2::plotMA(resultDESeq2,
+ colSig = "red",
+ colNonSig = "blue",
+ main="MA plot")
> abline(h=c(-1:1), col="red")
・Heat mapも描いてみたい。Heat map用にデータを絞る。p-adj <=0.01かつp-adjがある遺伝子を釣ってくる。
> gene.kept <- rownames(resultDESeq2)[resultDESeq2$padj <= 0.01 & !is.na(resultDESeq2$padj)]
> gene.kept %>% head()
[1] "hra1" "icr1" "rna170" "rpr1" "ruf23" "ruf5-1"
・対応する遺伝子のlog2したカウントデータをとってくる。
> countTable.kept <- log2(countTable + 1)[gene.kept, ]
> dim(countTable.kept)
[1] 4374 96
・次に、Heat mapを描いてみる。
> heatmap.2(as.matrix(countTable.kept),
+ scale="row",
+ hclust=function(x) hclust(x,method="average"),
+ distfun=function(x) as.dist((1-cor(t(x)))/2),
+ trace="none",
+ density="none",
+ labRow="",
+ cexCol=0.7)
・ちょっと変えてみる。
> mycol <- colorpanel(1000, "blue", "white", "red") #color map
> heatmap.2(as.matrix(countTable.kept),
+ scale="row",
+ hclust=function(x) hclust(x,method="average"),
+ distfun=function(x) as.dist((1-cor(t(x)))/2),
+ col = mycol, #color map
+ ColSideColors=phenoTable$color, #strain color bar
+ trace="none",
+ density="none",
+ labRow="",
+ cexCol=0.7)
・気持ち見やすくなった…？
つづく。
