############################################# ### サンプル間クラスタリング ### http://www.iu.a.u-tokyo.ac.jp/~kadota/r_seq.html#analysis_clustering_samples_TCC ### の例題7をテンプレートとしている ############################################# in_f <- "SRP001540_23_13.txt" #入力ファイル名を指定してin_fに格納 out_f <- "hoge7.png" #出力ファイル名を指定してout_fに格納 param_fig <- c(700, 400) #ファイル出力時の横幅と縦幅を指定(単位はピクセル) #必要なパッケージをロード library(TCC) #パッケージの読み込み #入力ファイルの読み込み data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み dim(data) #オブジェクトdataの行数と列数を表示 #本番 out <- clusterSample(data, dist.method="spearman",#クラスタリング実行結果をoutに格納 hclust.method="average", unique.pattern=TRUE)#クラスタリング実行結果をoutに格納 #ファイルに保存 png(out_f, pointsize=13, width=param_fig[1], height=param_fig[2])#出力ファイルの各種パラメータを指定 par(mar=c(0, 4, 1, 0)) #下、左、上、右の順で余白（行）を指定 plot(out, sub="", xlab="", cex.lab=1.2,#樹形図(デンドログラム)の表示 cex=1.3, main="", ylab="Height") #樹形図(デンドログラム)の表示 dev.off() #おまじない ############################################# ### Silhouette scores(シルエットスコア)を算出 ### http://www.iu.a.u-tokyo.ac.jp/~kadota/r_seq.html#analysis_general_silhouette ### の例題1をテンプレートとしている ### スコアは0.01873837と限りなくゼロに近い値 ############################################# in_f <- "SRP001540_23_13.txt" #入力ファイル名を指定してin_fに格納 param_G1 <- 23 #G1群のサンプル数を指定 param_G2 <- 13 #G2群のサンプル数を指定 #必要なパッケージをロード library(cluster) #パッケージの読み込み #入力ファイルの読み込みとラベル情報の作成 data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み data.cl <- c(rep(1, param_G1), rep(2, param_G2))#G1群を1、G2群を2としたベクトルdata.clを作成 dim(data) #オブジェクトdataの行数と列数を表示 #前処理(フィルタリング) obj <- as.logical(rowSums(data) > 0) #条件を満たすかどうかを判定した結果をobjに格納 data <- unique(data[obj,]) #objがTRUEとなる行のみ抽出し、ユニークパターンのみにした結果をdataに格納 dim(data) #オブジェクトdataの行数と列数を表示 #本番(AS値の計算) d <- as.dist(1 - cor(data, method="spearman"))#サンプル間の距離を計算し、結果をdに格納 AS <- mean(silhouette(data.cl, d)[, "sil_width"])#サンプルごとのSilhouette scoreを計算し、その平均値をASに格納 AS #AS値を表示 ############################################# ### 解析 | 発現変動 | 2群間 | 対応なし | 複製あり | TCC(Sun_2013) ### http://www.iu.a.u-tokyo.ac.jp/~kadota/r_seq.html#analysis_deg_2_unpaired_ari_TCC ### の例題1をテンプレートとしている ### 10%FDRを満たす遺伝子数は85。ないとは言えない程度か...。 ############################################# in_f <- "SRP001540_23_13.txt" #入力ファイル名を指定してin_fに格納 out_f1 <- "hoge1.txt" #出力ファイル名を指定してout_f1に格納 out_f2 <- "hoge1.png" #出力ファイル名を指定してout_f2に格納 param_G1 <- 23 #G1群のサンプル数を指定 param_G2 <- 13 #G2群のサンプル数を指定 param_FDR <- 0.05 #false discovery rate (FDR)閾値を指定 param_fig <- c(400, 380) #ファイル出力時の横幅と縦幅を指定(単位はピクセル) #必要なパッケージをロード library(TCC) #パッケージの読み込み #入力ファイルの読み込み data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み #前処理(TCCクラスオブジェクトの作成) data.cl <- c(rep(1, param_G1), rep(2, param_G2))#G1群を1、G2群を2としたベクトルdata.clを作成 tcc <- new("TCC", data, data.cl) #TCCクラスオブジェクトtccを作成 #本番(正規化) tcc <- calcNormFactors(tcc, norm.method="tmm", test.method="edger",#正規化を実行した結果をtccに格納 iteration=3, FDR=0.1, floorPDEG=0.05)#正規化を実行した結果をtccに格納 #本番(DEG検出) tcc <- estimateDE(tcc, test.method="edger", FDR=param_FDR)#DEG検出を実行した結果をtccに格納 result <- getResult(tcc, sort=FALSE) #p値などの計算結果をresultに格納 sum(tcc$stat$q.value < param_FDR) #条件を満たす遺伝子数を表示 #ファイルに保存(テキストファイル) tmp <- cbind(rownames(tcc$count), tcc$count, result)#入力データの右側にDEG検出結果を結合したものをtmpに格納 write.table(tmp, out_f1, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存 #ファイルに保存(M-A plot) png(out_f2, pointsize=13, width=param_fig[1], height=param_fig[2])#出力ファイルの各種パラメータを指定 plot(tcc, FDR=param_FDR) #param_FDRで指定した閾値を満たすDEGをマゼンタ色にして描画 legend("topright", c(paste("DEG(FDR<", param_FDR, ")", sep=""), "non-DEG"),#凡例を作成 col=c("magenta", "black"), pch=20)#凡例を作成 dev.off() #おまじない sum(tcc$stat$q.value < 0.05) #FDR < 0.05を満たす遺伝子数を表示 sum(tcc$stat$q.value < 0.10) #FDR < 0.10を満たす遺伝子数を表示 sum(tcc$stat$q.value < 0.20) #FDR < 0.20を満たす遺伝子数を表示 sum(tcc$stat$q.value < 0.30) #FDR < 0.30を満たす遺伝子数を表示