mecobalamin’s diary

人間万事塞翁が馬

https://help.hatenablog.com/entry/developer-option

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はリードカウントの結果の入ったファイル
追記ここまで