###################################################### ### Ichihashi et al., Plant Cell Physiol., 2018 ### https://pubmed.ncbi.nlm.nih.gov/29281058/ ### から取得可能な"pcp-2017-e-00401-File010.csv" ### が作業ディレクトリ上にある前提で記載しています。 ### setwd("~/2020/データサイエンス講習会R2/講習会R2/day1") ###################################################### in_f <- "pcp-2017-e-00401-File010.csv" #入力ファイル名を指定 param_G1 <- 6 #G1群のサンプル数を指定(Haustria) param_G2 <- 6 #G2群のサンプル数を指定(Roots) #必要なパッケージをロード library(TCC) #パッケージの読み込み library(edgeR) #パッケージの読み込み #パッケージのバージョン情報 packageVersion("TCC") #1.28.0 (私の確認時) packageVersion("edgeR") #3.30.3 (私の確認時) #入力ファイルの読み込みとラベル情報の作成 data <- read.csv(in_f, row.names="Transcript_ID")#ファイルの読み込み data <- as.matrix(data) #データの型をmatrixにしている data.cl <- c(rep(1, param_G1), rep(2, param_G2))#G1群を1、G2群を2としたベクトルdata.clを作成 group <- factor(data.cl) #factor形式に変更 ############################################ ### Case 1: TCC (ver. 1.28.0)の手順 ############################################ tcc <- new("TCC", data, data.cl) #TCCクラスオブジェクトtccを作成 tcc <- calcNormFactors(tcc) #正規化を実行した結果をtccに格納 tcc <- estimateDE(tcc) #DEG検出を実行した結果をtccに格納 p.value <- tcc$stat$p.value #p値をp.valueに格納 q.value <- p.adjust(p.value, method="BH")#q値をq.valueに格納 sum(q.value < 0.05) #q-value < 0.05を満たす遺伝子数は2,225個 ############################################ ### Case 2: edgeR (ver. 3.30.3)の手順 ### User's Guide中の1.4 Quick startに準拠していますが、 ### Case 1と条件を揃えるために低発現遺伝子のフィルタリング ### のところをコメントアウトしています。 ############################################ d <- DGEList(counts=data,group=group) #DGEListオブジェクトを作成してdに格納 #keep <- filterByExpr(d) #低発現遺伝子のフィルタリング #d <- d[keep,,keep.lib.sizes=FALSE] #低発現遺伝子のフィルタリング d <- calcNormFactors(d) #TMM正規化を実行 design <- model.matrix(~group) #デザイン行列を作成 d <- estimateDisp(d, design) #モデル構築 fit <- glmQLFit(d, design) #モデル構築 out <- glmQLFTest(fit, coef=2) #検定 p.value <- out$table$PValue #p値をp.valueに格納 q.value <- p.adjust(p.value, method="BH")#q値をq.valueに格納 sum(q.value < 0.05) #q-value < 0.05を満たす遺伝子数は1,764個 ############################################ ### Case 3: ### 「(Rで)塩基配列解析」で提供している ### 2015年頃のedgeRのUser's Guideに準拠したコードです。 ### Case 1と類似した結果(2,213個)であるが、 ### Case 2とは違っていることが分かります。 ############################################ d <- DGEList(counts=data,group=group) #DGEListオブジェクトを作成してdに格納 d <- calcNormFactors(d) #TMM正規化を実行 d <- estimateCommonDisp(d) #モデル構築 d <- estimateTagwiseDisp(d) #モデル構築 out <- exactTest(d) #検定 p.value <- out$table$PValue #p値をp.valueに格納 q.value <- p.adjust(p.value, method="BH")#q値をq.valueに格納 sum(q.value < 0.05) #q-value < 0.05を満たす遺伝子数は2,213個