RNA-seqその9、DESeq2を使った正規化
以前edgeRでリードカウントの正規化をした
同じ結果をDESeq2を使って正規化してみる
正規化の方法が違うらしい
DESeq2のマニュアル
これは読むべき
Analyzing RNA-seq data with DESeq2
マニュアルはDESeq2をインストール後に
browseVignettes("DESeq2")
でも表示できる
DESeq2はインストール済み
mecobalamin.hatenablog.com
バージョンが上がっているのでちょっと修正して再インストール
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("DESeq2", version = "3.9")
ライブラリの読み込みは
library(DESeq2)
DESeq2以外にマニュアルに出てきたライブラリのインストール
まずパッケージマネージャーを指定?する
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
次にそれぞれのライブラリをインストールする
# vsn BiocManager::install("vsn") # hexbin BiocManager::install("hexbin") # BiocParallel BiocManager::install("BiocParallel")
hexbinはbioconductorに載っている以下の方法は
Rのバージョンによっては使えないっぽい
source("http://bioconductor.org/biocLite.R") biocLite("hexbin")
この方法ではRのバージョン3.6ではエラーが出たので
BiocManagerを使ってインストールした
BiocParallelで呼び出すMulticoreParam()はwindowsでは使用できなかった
SnowParamは使えるがDESeq()のparallelに対応していないっぽい
インストールはしておく
DESeq2を読み込ませるとエラーが出たので追加でパッケージをインストール
install.packages("cluster", dependencies = TRUE)
マニュアル以外には以下のサイトを参考にさせてもらった
RNA-seqによる発現量解析
RNA-seqによる発現量解析 - Palmsonntagmorgen
サルマップ2018 (3) DESeq2による標準化と可視化まで
サルマップ2018 (3) DESeq2による標準化と可視化まで - ノンコーディングRNAネオタクソノミ
二群間比較(DESeq2)
二群間比較(DESeq2) | 一般化線形モデルによる発現変動遺伝子の検出
準備ここまで
リードカウントのデータは以前使ったのを再利用する
mecobalamin.hatenablog.com
コードも必要なところをedgeRからDESeq2のコマンドに書き変えて再利用
とりあえず計算結果
MA-plotはなんかずれてる?
分散もマニュアルとは雰囲気が違うな。。。
ばらつきが平坦にならない
variance stabilizing transformationとregularized log transformationも試したが
こっちも平坦にならない
variance stabilizing transformation
regularized log transformation
biological replicateではないのでしょうがないのか
リードカウントの正規化はそれっぽいけど。。。
acc.egfr iso.egfr H1 H2 H3 J1 J2 J3 1 NM_005228 a 8.375 11.061 11.400 0 0 0 2 NM_201282 b 6.379 9.113 10.042 0 0 0 3 NM_201283 c 5.287 8.695 9.071 0 0 0 4 NM_201284 d 6.402 9.119 10.056 0 0 0 5 NM_001346897 e 7.240 10.072 10.977 0 0 0 6 NM_001346898 f 7.266 10.072 10.981 0 0 0 7 NM_001346899 g 8.363 11.061 11.398 0 0 0 8 NM_001346900 h 8.348 11.034 11.389 0 0 0 9 NM_001346941 i 8.293 11.007 11.234 0 0 0
一連の操作ができていそうなので使ったコードをまとめておく
rm(list = ls(all = T)) library(DESeq2) library(vsn) # Paralellのオプションを使う場合に必要、windowsは不可 #library(BiocParallel) #register(MulticoreParam(4)) FindGeneInterest <- function(d){ #find genes of interest acc.egfr <- c("NM_005228", "NM_201282", "NM_201283", "NM_201284", "NM_001346897", "NM_001346898", "NM_001346899", "NM_001346900", "NM_001346941") iso.egfr <- c("a", "b", "c", "d", "e", "f", "g", "h", "i") var.egfr <- t(rbind(acc.egfr, iso.egfr)) egfr <- c() for(i in var.egfr[,1]){ egfr <- rbind(egfr, d[i,]) } print(data.frame(var.egfr, egfr)) } WriteGraphs <- function(res, vsd, rld, ntd){ #MA-plot # plot(log(res$baseMean),res$log2FoldChange) png(sprintf("DESeq2_MA-Plot.png"), width = 600, height = 450) plotMA(res) dev.off() #volcano plot png(sprintf("DESeq2_Volcano-Plot.png"), width = 600, height = 450) plot(res$log2FoldChange,-log(res$padj)) dev.off() # meanSdPlot png(sprintf("DESeq2_MeanSD-Plot_ntd.png"), width = 600, height = 450) meanSdPlot(assay(ntd)) dev.off() png(sprintf("DESeq2_MeanSD-Plot_vsd.png"), width = 600, height = 450) meanSdPlot(assay(vsd)) dev.off() png(sprintf("DESeq2_MeanSD-Plot_rld.png"), width = 600, height = 450) meanSdPlot(assay(rld)) dev.off() } DESeqNormalization <- function(Count.Data){ group <- data.frame(con = factor(c("H", "H", "H", "J", "J", "J"))) # チルダは必要、Rのformula dds <- DESeqDataSetFromMatrix(countData = Count.Data, colData = group, design = ~ con) dds <- DESeq(dds, parallel = FALSE) res <- results(dds, parallel = FALSE) res$gene_id <- row.names(res) # This gives log2(n + 1) ntd <- normTransform(dds) print(head(assay(ntd), 3)) # variance stabilizing transformation vsd <- vst(dds, blind = FALSE) print(head(assay(vsd), 3)) # regularized log transformation rld <- rlog(dds, blind = FALSE) print(head(assay(rld), 3)) WriteGraphs(res, vsd, rld, ntd) return(assay(ntd)) } options(digits = 4, width = 120) # ファイルの読み込み setwd("E:\\pathto\\countfile") Results.FC <- read.table("counts.txt", header = TRUE, skip = 1, sep="\t") # カラムに名前をつける Names.ID <- c("Geneid", "Chr", "Start", "End", "Strand", "Length", "H1", "H2", "H3", "J1", "J2", "J3") colnames(Results.FC) <- Names.ID rownames(Results.FC) <- sub("\\.[0-9]{1,}$", "", as.character(Results.FC[, 1])) # リードカウントの列を抜き出して行列にする Count.Data <- as.matrix(Results.FC[, 7:12]) # Filtering, リードカウントが0以上を残す Count.Keep <- rowSums(Count.Data) > 0 Count.Data <- Count.Data[Count.Keep, ] Data.Norm <- DESeqNormalization(Count.Data) FindGeneInterest(Data.Norm)
あとおまけ
計算時間の測定について
Rコードの実行時間を測るのに便利です「tictoc」パッケージ - Qiita
28 April 2020 追記
setwd("E:\\pathto\\countfile")
この行はcounts.txtのあるディレクトリを記入する
例えばcドライブのtmpディレクトリにcounts.txtがある場合は
setwd("c:\\tmp")
となる
counts.txtはリードカウントの結果の入ったファイル
追記ここまで