- 講義資料PDF(2016.07.12版; 約7MB)
- hoge.zip(2016.05.23版; 約3MB)
- 数値行列 (スライド5)
「サンプルデータ」の例題41は、以下と同じ。
Blekhman et al., Genome Res., 2010の
Supplementary Table1で提供されているエクセルファイル(http://genome.cshlp.org/content/suppl/2009/12/16/gr.099226.109.DC1/suppTable1.xls; 約4.3MB)
からカウントデータのみ抽出し、きれいに整形しなおしたものがここでの出力ファイルになります。
20,689 genes×36 samplesのカウントデータ(sample_blekhman_36.txt)です。
#in_f <- "http://genome.cshlp.org/content/suppl/2009/12/16/gr.099226.109.DC1/suppTable1.xls"
in_f <- "suppTable1.xls"
out_f <- "sample_blekhman_36.txt"
hoge <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
dim(hoge)
data <- cbind(
hoge$R1L4.HSF1, hoge$R4L2.HSF1, hoge$R2L7.HSF2, hoge$R3L2.HSF2, hoge$R8L1.HSF3, hoge$R8L2.HSF3,
hoge$R1L1.HSM1, hoge$R5L2.HSM1, hoge$R2L3.HSM2, hoge$R4L8.HSM2, hoge$R3L6.HSM3, hoge$R4L1.HSM3,
hoge$R1L2.PTF1, hoge$R4L4.PTF1, hoge$R2L4.PTF2, hoge$R6L6.PTF2, hoge$R3L7.PTF3, hoge$R5L3.PTF3,
hoge$R1L6.PTM1, hoge$R3L3.PTM1, hoge$R2L8.PTM2, hoge$R4L6.PTM2, hoge$R6L2.PTM3, hoge$R6L4.PTM3,
hoge$R1L7.RMF1, hoge$R5L1.RMF1, hoge$R2L2.RMF2, hoge$R5L8.RMF2, hoge$R3L4.RMF3, hoge$R4L7.RMF3,
hoge$R1L3.RMM1, hoge$R3L8.RMM1, hoge$R2L6.RMM2, hoge$R5L4.RMM2, hoge$R3L1.RMM3, hoge$R4L3.RMM3)
colnames(data) <- c(
"R1L4.HSF1", "R4L2.HSF1", "R2L7.HSF2", "R3L2.HSF2", "R8L1.HSF3", "R8L2.HSF3",
"R1L1.HSM1", "R5L2.HSM1", "R2L3.HSM2", "R4L8.HSM2", "R3L6.HSM3", "R4L1.HSM3",
"R1L2.PTF1", "R4L4.PTF1", "R2L4.PTF2", "R6L6.PTF2", "R3L7.PTF3", "R5L3.PTF3",
"R1L6.PTM1", "R3L3.PTM1", "R2L8.PTM2", "R4L6.PTM2", "R6L2.PTM3", "R6L4.PTM3",
"R1L7.RMF1", "R5L1.RMF1", "R2L2.RMF2", "R5L8.RMF2", "R3L4.RMF3", "R4L7.RMF3",
"R1L3.RMM1", "R3L8.RMM1", "R2L6.RMM2", "R5L4.RMM2", "R3L1.RMM3", "R4L3.RMM3")
rownames(data)<- rownames(hoge)
dim(data)
tmp <- cbind(rownames(data), data)
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)
- データの正規化 (スライド8)
「書籍 | トランスクリプトーム解析 | 3.3.2 データの正規化(基礎編)」と説明内容は同じです。
2262/3
(2262/6)/3
(2262/3000)/(3/500)
2262/3000
3/500
2262*(1000/3000)
3*(1000/500)
colSums(data)
- RPM正規化 (スライド14)
「正規化 | 基礎 | RPM or CPM (総リード数補正)」をベースに作成。
入力は、20,689 genes×36 samplesのカウントデータ(sample_blekhman_36.txt)です。
教科書p134。
in_f <- "sample_blekhman_36.txt"
out_f <- "hoge1.txt"
param1 <- 1000000
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
colSums(data)
nf <- param1/colSums(data)
data <- sweep(data, 2, nf, "*")
colSums(data)
tmp <- cbind(rownames(data), data)
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)
- RPKM正規化 (スライド18)
行名が"ENSG00000000971"で、1列目のカウント値(2262)に対するRPKM値の計算例を示す。
data["ENSG00000000971", 1]
nf[1]
2262 * nf[1]
2262 * (1000000/1665987)
1000/3000
data["ENSG00000000971", 1] * (1000/3000)
2262 * (1000000/1665987) * (1000/3000)
- クラスタリング(サンプル間) (スライド21)
「解析 | クラスタリング | サンプル間 | TCC(Sun_2013)」の例題7は以下と同じ。
入力は、20,689 genes×36 samplesのカウントデータ(sample_blekhman_36.txt)です。
in_f <- "sample_blekhman_36.txt"
out_f <- "hoge7.png"
param_fig <- c(700, 400)
library(TCC)
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
dim(data)
out <- clusterSample(data, dist.method="spearman",
hclust.method="average", unique.pattern=TRUE)
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()
- サブセット抽出と整形 (スライド25)
「サンプルデータ」の例題42は以下と同じ。
入力は、Blekhman et al., Genome Res., 2010の
Supplementary Table1で提供されているエクセルファイル(http://genome.cshlp.org/content/suppl/2009/12/16/gr.099226.109.DC1/suppTable1.xls; 約4.3MB)
ですが、hogeフォルダ中にあります。例題41とは違って、technical replicatesの2列分のデータは足して1列分のデータとしています。
出力は、20,689 genes×18 samplesのカウントデータ(sample_blekhman_18.txt)です。
#in_f <- "http://genome.cshlp.org/content/suppl/2009/12/16/gr.099226.109.DC1/suppTable1.xls"
in_f <- "suppTable1.xls"
out_f <- "sample_blekhman_18.txt"
hoge <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
dim(hoge)
data <- cbind(
hoge$R1L4.HSF1 + hoge$R4L2.HSF1, hoge$R2L7.HSF2 + hoge$R3L2.HSF2, hoge$R8L1.HSF3 + hoge$R8L2.HSF3,
hoge$R1L1.HSM1 + hoge$R5L2.HSM1, hoge$R2L3.HSM2 + hoge$R4L8.HSM2, hoge$R3L6.HSM3 + hoge$R4L1.HSM3,
hoge$R1L2.PTF1 + hoge$R4L4.PTF1, hoge$R2L4.PTF2 + hoge$R6L6.PTF2, hoge$R3L7.PTF3 + hoge$R5L3.PTF3,
hoge$R1L6.PTM1 + hoge$R3L3.PTM1, hoge$R2L8.PTM2 + hoge$R4L6.PTM2, hoge$R6L2.PTM3 + hoge$R6L4.PTM3,
hoge$R1L7.RMF1 + hoge$R5L1.RMF1, hoge$R2L2.RMF2 + hoge$R5L8.RMF2, hoge$R3L4.RMF3 + hoge$R4L7.RMF3,
hoge$R1L3.RMM1 + hoge$R3L8.RMM1, hoge$R2L6.RMM2 + hoge$R5L4.RMM2, hoge$R3L1.RMM3 + hoge$R4L3.RMM3)
colnames(data) <- c(
"HSF1", "HSF2", "HSF3", "HSM1", "HSM2", "HSM3",
"PTF1", "PTF2", "PTF3", "PTM1", "PTM2", "PTM3",
"RMF1", "RMF2", "RMF3", "RMM1", "RMM2", "RMM3")
rownames(data)<- rownames(hoge)
dim(data)
tmp <- cbind(rownames(data), data)
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)
- クラスタリング(サンプル間) (スライド27)
「解析 | クラスタリング | サンプル間 | TCC(Sun_2013)」の例題8は以下と同じ。
入力は、20,689 genes×18 samplesのカウントデータ(sample_blekhman_18.txt)です。
in_f <- "sample_blekhman_18.txt"
out_f <- "hoge8.png"
param_fig <- c(700, 400)
library(TCC)
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
dim(data)
out <- clusterSample(data, dist.method="spearman",
hclust.method="average", unique.pattern=TRUE)
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()
- クラスタリングとDEG数の関係 (スライド32)
原著論文はTang et al., BMC Bioinformatics, 2015です。
- Tips: cex.lab (スライド33)
cex.labは、この場合縦軸ラベル情報("Height")の文字の大きさをデフォルトの何倍にして表示するかを指定するオプション。
例えばデフォルトの2.0倍にしたいときは、cex.lab=2.0とすればよい。
in_f <- "sample_blekhman_18.txt"
out_f <- "hoge.png"
param_fig <- c(600, 400)
library(TCC)
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
dim(data)
out <- clusterSample(data, dist.method="spearman",
hclust.method="average", unique.pattern=TRUE)
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=2.0,
cex=1.3, main="", ylab="Height")
dev.off()
- Tips: mar (スライド34)
marの左部分の余白を4から5行分に変更することで、"Height"の文字が切れなくなる。
marはmargin (マージン)の意味です。
in_f <- "sample_blekhman_18.txt"
out_f <- "hoge.png"
param_fig <- c(600, 400)
library(TCC)
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
dim(data)
out <- clusterSample(data, dist.method="spearman",
hclust.method="average", unique.pattern=TRUE)
png(out_f, pointsize=13, width=param_fig[1], height=param_fig[2])
par(mar=c(0, 5, 1, 0))
plot(out, sub="", xlab="", cex.lab=2.0,
cex=1.3, main="", ylab="Height")
dev.off()
- ReCount:Frazee et al., BMC Bioinformatics, 2011 (スライド36)
マッピング済みの行列形式のカウントデータセットを多数提供しています。
getwd()
list.files(pattern="bodymap")
- bodymapデータ解析 (スライド39)
「書籍 | トランスクリプトーム解析 | 3.3.3 クラスタリング」の「p143-144の網掛け部分」は、以下と同じ。
in_f1 <- "bodymap_count_table.txt"
in_f2 <- "bodymap_phenodata.txt"
data <- read.table(in_f1, header=TRUE, row.names=1, sep="\t", quote="")
phenotype <- read.table(in_f2, header=TRUE, row.names=1, sep=" ", quote="")
colnames(data) <- phenotype$tissue.type
obj <- rowSums(data) != 0
hoge <- unique(data[obj,])
data.dist <- as.dist(1 - cor(hoge, method = "spearman"))
out <- hclust(data.dist, method = "average")
plot(out)
- TCCで実行すると... (スライド41)
「書籍 | トランスクリプトーム解析 | 3.3.3 クラスタリング」の「p143-144の網掛け部分」を、
「解析 | クラスタリング | サンプル間 | TCC(Sun_2013)」
のように書くと以下のような感じになります。
in_f1 <- "bodymap_count_table.txt"
in_f2 <- "bodymap_phenodata.txt"
out_f <- "hoge.png"
param_fig <- c(600, 400)
library(TCC)
data <- read.table(in_f1, header=TRUE, row.names=1, sep="\t", quote="")
phenotype <- read.table(in_f2, header=TRUE, row.names=1, sep=" ", quote="")
phenotype
colnames(data) <- phenotype$tissue.type
out <- clusterSample(data, dist.method="spearman",
hclust.method="average", unique.pattern=TRUE)
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()
- ReCountのヒトデータは (スライド42)
bodymapを含む、ReCountで提供されているヒトデータは同じ遺伝子数
dim(data)
head(rownames(data))
- giladデータ解析 (スライド45)
getwd()
list.files(pattern="gilad")
- giladデータ解析 (スライド46)
giladデータのサンプル間クラスタリングを実行。phenotype情報は入力ファイルによってフォーマットが異なる点に注意。
in_f1 <- "gilad_count_table.txt"
in_f2 <- "gilad_phenodata.txt"
out_f <- "hoge_gilad.png"
param_fig <- c(600, 400)
library(TCC)
data <- read.table(in_f1, header=TRUE, row.names=1, sep="\t", quote="")
phenotype <- read.table(in_f2, header=TRUE, row.names=1, sep=" ", quote="")
phenotype
colnames(data) <- paste(phenotype$gender, rownames(phenotype), sep="_")
colnames(data)
out <- clusterSample(data, dist.method="spearman",
hclust.method="average", unique.pattern=TRUE)
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()
- ReCountのヒトデータは (スライド53)
giladを含む、ReCountで提供されているヒトデータは同じ遺伝子数
dim(data)
head(rownames(data))
- bodymap + gilad (スライド55)
2つのデータセットを読み込んで、マージ(列方向で結合)して、サンプル間クラスタリングを実行。
out_f <- "hoge_merge.png"
param_fig <- c(700, 400)
library(TCC)
in_f1 <- "bodymap_count_table.txt"
in_f2 <- "bodymap_phenodata.txt"
data <- read.table(in_f1, header=TRUE, row.names=1, sep="\t", quote="")
phenotype <- read.table(in_f2, header=TRUE, row.names=1, sep=" ", quote="")
phenotype
colnames(data) <- phenotype$tissue.type
colnames(data)
dim(data)
data_bodymap <- data
in_f1 <- "gilad_count_table.txt"
in_f2 <- "gilad_phenodata.txt"
data <- read.table(in_f1, header=TRUE, row.names=1, sep="\t", quote="")
phenotype <- read.table(in_f2, header=TRUE, row.names=1, sep=" ", quote="")
phenotype
colnames(data) <- paste(phenotype$gender, rownames(phenotype), sep="_")
colnames(data)
dim(data)
data_gilad <- data
data <- cbind(data_bodymap, data_gilad)
colnames(data)
dim(data)
out <- clusterSample(data, dist.method="spearman",
hclust.method="average", unique.pattern=TRUE)
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()
- HS vs. RM (スライド64)
「解析 | 発現変動 | 2群間 | 対応なし | 複製あり | Blekhmanデータ | TCC(Sun_2013)」
の例題1は以下と同じ。ヒト2サンプル(G1群:HSF1とHSM1) vs. アカゲザル2サンプル(G2群:RMF1とRMM1)です。
1, 4, 13, 16
列目のデータのみ抽出しています。
in_f <- "sample_blekhman_18.txt"
out_f1 <- "hoge1.txt"
out_f2 <- "hoge1.png"
param_subset <- c(1, 4, 13, 16)
param_G1 <- 2
param_G2 <- 2
param_FDR <- 0.05
param_fig <- c(430, 350)
param_mar <- c(4, 4, 0, 0)
library(TCC)
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
data <- data[,param_subset]
data.cl <- c(rep(1, param_G1), rep(2, param_G2))
tcc <- new("TCC", data, data.cl)
dim(data)
head(data)
tcc <- calcNormFactors(tcc, norm.method="tmm", test.method="edger",
iteration=3, FDR=0.1, floorPDEG=0.05)
normalized <- getNormalizedData(tcc)
tcc <- estimateDE(tcc, test.method="edger", FDR=param_FDR)
result <- getResult(tcc, sort=FALSE)
sum(tcc$stat$q.value < param_FDR)
tmp <- cbind(rownames(tcc$count), normalized, result)
tmp <- tmp[order(tmp$rank),]
write.table(tmp, out_f1, sep="\t", append=F, quote=F, row.names=F)
png(out_f2, pointsize=13, width=param_fig[1], height=param_fig[2])
par(mar=param_mar)
plot(tcc, FDR=param_FDR, xlim=c(-2, 17), ylim=c(-10, 10),
cex=0.8, cex.lab=1.2,
cex.axis=1.2, main="",
xlab="A = (log2(G2) + log2(G1))/2",
ylab="M = log2(G2) - log2(G1)")
legend("topright", c(paste("DEG(FDR<", param_FDR, ")", sep=""), "non-DEG"),
col=c("magenta", "black"), pch=20, cex=1.2)
dev.off()
sum(tcc$stat$q.value < 0.05)
sum(tcc$stat$q.value < 0.10)
sum(tcc$stat$q.value < 0.20)
sum(tcc$stat$q.value < 0.30)
- DEG数の見積もり (スライド71)
2489*(1 - 0.05)
3121*(1 - 0.10)
4049*(1 - 0.20)
4785*(1 - 0.30)
- DEG検出結果1位 (スライド77)
「解析 | 発現変動 | 2群間 | 対応なし | 複製あり | Blekhmanデータ | TCC(Sun_2013)」
の例題2は以下と同じ。ランキング1位の"ENSG00000208570"をハイライトさせるテクニックです。
1位のy軸の値(m.value)が11.29なので、y軸の範囲を[-10.0, 11.5]の範囲に変更する操作も行っています。
in_f <- "sample_blekhman_18.txt"
out_f1 <- "hoge2.txt"
out_f2 <- "hoge2.png"
param_subset <- c(1, 4, 13, 16)
param_G1 <- 2
param_G2 <- 2
param_FDR <- 0.05
param_fig <- c(400, 310)
param_mar <- c(4, 4, 0, 0)
param_geneid <- "ENSG00000208570"
param_col <- "magenta"
param_cex <- 2
param_pch <- 20
library(TCC)
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
data <- data[,param_subset]
data.cl <- c(rep(1, param_G1), rep(2, param_G2))
tcc <- new("TCC", data, data.cl)
dim(data)
head(data)
tcc <- calcNormFactors(tcc, norm.method="tmm", test.method="edger",
iteration=3, FDR=0.1, floorPDEG=0.05)
normalized <- getNormalizedData(tcc)
tcc <- estimateDE(tcc, test.method="edger", FDR=param_FDR)
result <- getResult(tcc, sort=FALSE)
sum(tcc$stat$q.value < param_FDR)
tmp <- cbind(rownames(tcc$count), normalized, result)
tmp <- tmp[order(tmp$rank),]
write.table(tmp, out_f1, sep="\t", append=F, quote=F, row.names=F)
png(out_f2, pointsize=13, width=param_fig[1], height=param_fig[2])
par(mar=param_mar)
plot(tcc, FDR=param_FDR, xlim=c(-2, 17), ylim=c(-10.0, 11.5),
cex=0.8, cex.lab=1.2,
cex.axis=1.2, main="",
xlab="A = (log2(G2) + log2(G1))/2",
ylab="M = log2(G2) - log2(G1)")
legend("topright", c(paste("DEG(FDR<", param_FDR, ")", sep=""), "non-DEG"),
col=c("magenta", "black"), pch=20, cex=1.2)
obj <- is.element(tcc$gene_id, param_geneid)
points(result$a.value[obj], result$m.value[obj],
pch=param_pch, cex=param_cex, col=param_col)
dev.off()
sum(tcc$stat$q.value < 0.05)
sum(tcc$stat$q.value < 0.10)
sum(tcc$stat$q.value < 0.20)
sum(tcc$stat$q.value < 0.30)
- DEG検出結果2位 (スライド78)
「解析 | 発現変動 | 2群間 | 対応なし | 複製あり | Blekhmanデータ | TCC(Sun_2013)」
の例題3は以下と同じ。ランキング2位の"ENSG00000220191"をハイライトさせるテクニックです。
param_cexやparam_pchを変更しています。
in_f <- "sample_blekhman_18.txt"
out_f1 <- "hoge3.txt"
out_f2 <- "hoge3.png"
param_subset <- c(1, 4, 13, 16)
param_G1 <- 2
param_G2 <- 2
param_FDR <- 0.05
param_fig <- c(400, 310)
param_mar <- c(4, 4, 0, 0)
param_geneid <- "ENSG00000220191"
param_col <- "magenta"
param_cex <- 2.3
param_pch <- 15
library(TCC)
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
data <- data[,param_subset]
data.cl <- c(rep(1, param_G1), rep(2, param_G2))
tcc <- new("TCC", data, data.cl)
dim(data)
head(data)
tcc <- calcNormFactors(tcc, norm.method="tmm", test.method="edger",
iteration=3, FDR=0.1, floorPDEG=0.05)
normalized <- getNormalizedData(tcc)
tcc <- estimateDE(tcc, test.method="edger", FDR=param_FDR)
result <- getResult(tcc, sort=FALSE)
sum(tcc$stat$q.value < param_FDR)
tmp <- cbind(rownames(tcc$count), normalized, result)
tmp <- tmp[order(tmp$rank),]
write.table(tmp, out_f1, sep="\t", append=F, quote=F, row.names=F)
png(out_f2, pointsize=13, width=param_fig[1], height=param_fig[2])
par(mar=param_mar)
plot(tcc, FDR=param_FDR, xlim=c(-2, 17), ylim=c(-10.0, 11.5),
cex=0.8, cex.lab=1.2,
cex.axis=1.2, main="",
xlab="A = (log2(G2) + log2(G1))/2",
ylab="M = log2(G2) - log2(G1)")
legend("topright", c(paste("DEG(FDR<", param_FDR, ")", sep=""), "non-DEG"),
col=c("magenta", "black"), pch=20, cex=1.2)
obj <- is.element(tcc$gene_id, param_geneid)
points(result$a.value[obj], result$m.value[obj],
pch=param_pch, cex=param_cex, col=param_col)
dev.off()
sum(tcc$stat$q.value < 0.05)
sum(tcc$stat$q.value < 0.10)
sum(tcc$stat$q.value < 0.20)
sum(tcc$stat$q.value < 0.30)
- DEG検出結果3位 (スライド79)
「解析 | 発現変動 | 2群間 | 対応なし | 複製あり | Blekhmanデータ | TCC(Sun_2013)」
の例題4は以下と同じ。ランキング3位の"ENSG00000106366"をハイライトさせるテクニックです。
param_cexやparam_pchを変更しています。
in_f <- "sample_blekhman_18.txt"
out_f1 <- "hoge4.txt"
out_f2 <- "hoge4.png"
param_subset <- c(1, 4, 13, 16)
param_G1 <- 2
param_G2 <- 2
param_FDR <- 0.05
param_fig <- c(400, 310)
param_mar <- c(4, 4, 0, 0)
param_geneid <- "ENSG00000106366"
param_col <- "magenta"
param_cex <- 3.1
param_pch <- 18
library(TCC)
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
data <- data[,param_subset]
data.cl <- c(rep(1, param_G1), rep(2, param_G2))
tcc <- new("TCC", data, data.cl)
dim(data)
head(data)
tcc <- calcNormFactors(tcc, norm.method="tmm", test.method="edger",
iteration=3, FDR=0.1, floorPDEG=0.05)
normalized <- getNormalizedData(tcc)
tcc <- estimateDE(tcc, test.method="edger", FDR=param_FDR)
result <- getResult(tcc, sort=FALSE)
sum(tcc$stat$q.value < param_FDR)
tmp <- cbind(rownames(tcc$count), normalized, result)
tmp <- tmp[order(tmp$rank),]
write.table(tmp, out_f1, sep="\t", append=F, quote=F, row.names=F)
png(out_f2, pointsize=13, width=param_fig[1], height=param_fig[2])
par(mar=param_mar)
plot(tcc, FDR=param_FDR, xlim=c(-2, 17), ylim=c(-10.0, 11.5),
cex=0.8, cex.lab=1.2,
cex.axis=1.2, main="",
xlab="A = (log2(G2) + log2(G1))/2",
ylab="M = log2(G2) - log2(G1)")
legend("topright", c(paste("DEG(FDR<", param_FDR, ")", sep=""), "non-DEG"),
col=c("magenta", "black"), pch=20, cex=1.2)
obj <- is.element(tcc$gene_id, param_geneid)
points(result$a.value[obj], result$m.value[obj],
pch=param_pch, cex=param_cex, col=param_col)
dev.off()
sum(tcc$stat$q.value < 0.05)
sum(tcc$stat$q.value < 0.10)
sum(tcc$stat$q.value < 0.20)
sum(tcc$stat$q.value < 0.30)
- DEG検出結果2,489位 (スライド80)
「解析 | 発現変動 | 2群間 | 対応なし | 複製あり | Blekhmanデータ | TCC(Sun_2013)」
の例題5は以下と同じ。ランキング2,489位の"ENSG00000139445"をハイライトさせるテクニックです。
param_col, param_cex, param_pchなどを変更しています。
in_f <- "sample_blekhman_18.txt"
out_f1 <- "hoge5.txt"
out_f2 <- "hoge5.png"
param_subset <- c(1, 4, 13, 16)
param_G1 <- 2
param_G2 <- 2
param_FDR <- 0.05
param_fig <- c(400, 310)
param_mar <- c(4, 4, 0, 0)
param_geneid <- "ENSG00000139445"
param_col <- "skyblue"
param_cex <- 2.5
param_pch <- 17
library(TCC)
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
data <- data[,param_subset]
data.cl <- c(rep(1, param_G1), rep(2, param_G2))
tcc <- new("TCC", data, data.cl)
dim(data)
head(data)
tcc <- calcNormFactors(tcc, norm.method="tmm", test.method="edger",
iteration=3, FDR=0.1, floorPDEG=0.05)
normalized <- getNormalizedData(tcc)
tcc <- estimateDE(tcc, test.method="edger", FDR=param_FDR)
result <- getResult(tcc, sort=FALSE)
sum(tcc$stat$q.value < param_FDR)
tmp <- cbind(rownames(tcc$count), normalized, result)
tmp <- tmp[order(tmp$rank),]
write.table(tmp, out_f1, sep="\t", append=F, quote=F, row.names=F)
png(out_f2, pointsize=13, width=param_fig[1], height=param_fig[2])
par(mar=param_mar)
plot(tcc, FDR=param_FDR, xlim=c(-2, 17), ylim=c(-10.0, 11.5),
cex=0.8, cex.lab=1.2,
cex.axis=1.2, main="",
xlab="A = (log2(G2) + log2(G1))/2",
ylab="M = log2(G2) - log2(G1)")
legend("topright", c(paste("DEG(FDR<", param_FDR, ")", sep=""), "non-DEG"),
col=c("magenta", "black"), pch=20, cex=1.2)
obj <- is.element(tcc$gene_id, param_geneid)
points(result$a.value[obj], result$m.value[obj],
pch=param_pch, cex=param_cex, col=param_col)
dev.off()
sum(tcc$stat$q.value < 0.05)
sum(tcc$stat$q.value < 0.10)
sum(tcc$stat$q.value < 0.20)
sum(tcc$stat$q.value < 0.30)
- 分布やモデル (スライド82)
「解析 | 発現変動 | 2群間 | 対応なし | 複製あり | Blekhmanデータ | TCC(Sun_2013)」
の例題6は以下と同じ。3種類のFDR閾値(0.0001, 0.05, 0.4)でのM-A plotを出力させています。
in_f <- "sample_blekhman_18.txt"
out_f1 <- "hoge6_1.png"
out_f2 <- "hoge6_2.png"
out_f3 <- "hoge6_3.png"
param_subset <- c(1, 4, 13, 16)
param_G1 <- 2
param_G2 <- 2
param_FDR <- c(0.0001, 0.05, 0.4)
param_fig <- c(375, 350)
param_mar <- c(4, 4, 0, 0)
library(TCC)
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
data <- data[,param_subset]
data.cl <- c(rep(1, param_G1), rep(2, param_G2))
tcc <- new("TCC", data, data.cl)
dim(data)
head(data)
tcc <- calcNormFactors(tcc, norm.method="tmm", test.method="edger",
iteration=3, FDR=0.1, floorPDEG=0.05)
normalized <- getNormalizedData(tcc)
tcc <- estimateDE(tcc, test.method="edger")
result <- getResult(tcc, sort=FALSE)
sum(tcc$stat$q.value < param_FDR[1])
sum(tcc$stat$q.value < param_FDR[2])
sum(tcc$stat$q.value < param_FDR[3])
png(out_f1, pointsize=13, width=param_fig[1], height=param_fig[2])
par(mar=param_mar)
plot(tcc, FDR=param_FDR[1], xlim=c(-2, 17), ylim=c(-10.0, 11.5),
cex=0.8, cex.lab=1.2,
cex.axis=1.2, main="",
xlab="A = (log2(G2) + log2(G1))/2",
ylab="M = log2(G2) - log2(G1)")
legend("topright", c(paste("DEG(FDR<", param_FDR[1], ")", sep=""), "non-DEG"),
col=c("magenta", "black"), pch=20, cex=1.2)
dev.off()
png(out_f2, pointsize=13, width=param_fig[1], height=param_fig[2])
par(mar=param_mar)
plot(tcc, FDR=param_FDR[2], xlim=c(-2, 17), ylim=c(-10.0, 11.5),
cex=0.8, cex.lab=1.2,
cex.axis=1.2, main="",
xlab="A = (log2(G2) + log2(G1))/2",
ylab="M = log2(G2) - log2(G1)")
legend("topright", c(paste("DEG(FDR<", param_FDR[2], ")", sep=""), "non-DEG"),
col=c("magenta", "black"), pch=20, cex=1.2)
dev.off()
png(out_f3, pointsize=13, width=param_fig[1], height=param_fig[2])
par(mar=param_mar)
plot(tcc, FDR=param_FDR[3], xlim=c(-2, 17), ylim=c(-10.0, 11.5),
cex=0.8, cex.lab=1.2,
cex.axis=1.2, main="",
xlab="A = (log2(G2) + log2(G1))/2",
ylab="M = log2(G2) - log2(G1)")
legend("topright", c(paste("DEG(FDR<", param_FDR[3], ")", sep=""), "non-DEG"),
col=c("magenta", "black"), pch=20, cex=1.2)
dev.off()
- HS vs. PT (スライド85)
「解析 | 発現変動 | 2群間 | 対応なし | 複製あり | Blekhmanデータ | TCC(Sun_2013)」
の例題7は以下と同じ。1, 4, 7, 10
列目のデータのみ抽出して、ヒト2サンプル(G1群:HSF1とHSM1) vs. チンパンジー2サンプル(G2群:PTF1とPTM1)の2群間比較を行っています。
in_f <- "sample_blekhman_18.txt"
out_f1 <- "hoge7.txt"
out_f2 <- "hoge7.png"
param_subset <- c(1, 4, 7, 10)
param_G1 <- 2
param_G2 <- 2
param_FDR <- 0.05
param_fig <- c(430, 350)
param_mar <- c(4, 4, 0, 0)
library(TCC)
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
data <- data[,param_subset]
data.cl <- c(rep(1, param_G1), rep(2, param_G2))
tcc <- new("TCC", data, data.cl)
colnames(data)
tcc <- calcNormFactors(tcc, norm.method="tmm", test.method="edger",
iteration=3, FDR=0.1, floorPDEG=0.05)
normalized <- getNormalizedData(tcc)
tcc <- estimateDE(tcc, test.method="edger", FDR=param_FDR)
result <- getResult(tcc, sort=FALSE)
sum(tcc$stat$q.value < param_FDR)
tmp <- cbind(rownames(tcc$count), normalized, result)
tmp <- tmp[order(tmp$rank),]
write.table(tmp, out_f1, sep="\t", append=F, quote=F, row.names=F)
png(out_f2, pointsize=13, width=param_fig[1], height=param_fig[2])
par(mar=param_mar)
plot(tcc, FDR=param_FDR, xlim=c(-2, 17), ylim=c(-10.0, 11.5),
cex=0.8, cex.lab=1.2,
cex.axis=1.2, main="",
xlab="A = (log2(G2) + log2(G1))/2",
ylab="M = log2(G2) - log2(G1)")
legend("topright", c(paste("DEG(FDR<", param_FDR, ")", sep=""), "non-DEG"),
col=c("magenta", "black"), pch=20, cex=1.2)
dev.off()
sum(tcc$stat$q.value < 0.05)
sum(tcc$stat$q.value < 0.10)
sum(tcc$stat$q.value < 0.20)
sum(tcc$stat$q.value < 0.30)
- サンプル間類似度RM (スライド92)
unique関数周辺については、教科書p143辺りで解説しています。
in_f <- "sample_blekhman_18.txt"
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
dim(data)
data <- unique(data)
dim(data)
cor(data$HSM1, data$HSF1, method="spearman")
cor(data$RMM1, data$RMF1, method="spearman")
cor(data$HSM1, data$RMM1, method="spearman")
- サンプル間類似度PT (スライド93)
unique関数周辺については、教科書p143辺りで解説しています。
in_f <- "sample_blekhman_18.txt"
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
dim(data)
data <- unique(data)
dim(data)
cor(data$HSM1, data$HSF1, method="spearman")
cor(data$PTM1, data$PTF1, method="spearman")
cor(data$HSM1, data$PTM1, method="spearman")
- HS1 vs. HS2 (スライド97)
「解析 | 発現変動 | 2群間 | 対応なし | 複製あり | Blekhmanデータ | TCC(Sun_2013)」
の例題8は以下と同じ。1, 4, 2, 5
列目のデータのみ抽出して、ヒト2サンプル(G1群:HSF1とHSM1) vs. ヒト2サンプル(G2群:HSF2とHSM2)の2群間比較を行います。
in_f <- "sample_blekhman_18.txt"
out_f1 <- "hoge8.txt"
out_f2 <- "hoge8.png"
param_subset <- c(1, 4, 2, 5)
param_G1 <- 2
param_G2 <- 2
param_FDR <- 0.05
param_fig <- c(430, 350)
param_mar <- c(4, 4, 0, 0)
library(TCC)
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
data <- data[,param_subset]
data.cl <- c(rep(1, param_G1), rep(2, param_G2))
tcc <- new("TCC", data, data.cl)
colnames(data)
tcc <- calcNormFactors(tcc, norm.method="tmm", test.method="edger",
iteration=3, FDR=0.1, floorPDEG=0.05)
normalized <- getNormalizedData(tcc)
tcc <- estimateDE(tcc, test.method="edger", FDR=param_FDR)
result <- getResult(tcc, sort=FALSE)
sum(tcc$stat$q.value < param_FDR)
tmp <- cbind(rownames(tcc$count), normalized, result)
tmp <- tmp[order(tmp$rank),]
write.table(tmp, out_f1, sep="\t", append=F, quote=F, row.names=F)
png(out_f2, pointsize=13, width=param_fig[1], height=param_fig[2])
par(mar=param_mar)
plot(tcc, FDR=param_FDR, xlim=c(-2, 17), ylim=c(-10.0, 11.5),
cex=0.8, cex.lab=1.2,
cex.axis=1.2, main="",
xlab="A = (log2(G2) + log2(G1))/2",
ylab="M = log2(G2) - log2(G1)")
obj <- as.logical(tcc$stat$q.value < param_FDR)
points(result$a.value[obj], result$m.value[obj],
pch=20, cex=1.5, col="magenta")
legend("topright", c(paste("DEG(FDR<", param_FDR, ")", sep=""), "non-DEG"),
col=c("magenta", "black"), pch=20, cex=1.2)
dev.off()
sum(tcc$stat$q.value < 0.05)
sum(tcc$stat$q.value < 0.10)
sum(tcc$stat$q.value < 0.20)
sum(tcc$stat$q.value < 0.30)
- PT1 vs. PT2 (スライド98)
「解析 | 発現変動 | 2群間 | 対応なし | 複製あり | Blekhmanデータ | TCC(Sun_2013)」
の例題9は以下と同じ。7, 10, 8, 11列目のデータのみ抽出して、
チンパンジー2サンプル(G1群:PTF1とPTM1) vs. チンパンジー2サンプル(G2群:PTF2とPTM2)の2群間比較を行います。
in_f <- "sample_blekhman_18.txt"
out_f1 <- "hoge9.txt"
out_f2 <- "hoge9.png"
param_subset <- c(7, 10, 8, 11)
param_G1 <- 2
param_G2 <- 2
param_FDR <- 0.05
param_fig <- c(430, 350)
param_mar <- c(4, 4, 0, 0)
library(TCC)
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
data <- data[,param_subset]
data.cl <- c(rep(1, param_G1), rep(2, param_G2))
tcc <- new("TCC", data, data.cl)
colnames(data)
tcc <- calcNormFactors(tcc, norm.method="tmm", test.method="edger",
iteration=3, FDR=0.1, floorPDEG=0.05)
normalized <- getNormalizedData(tcc)
tcc <- estimateDE(tcc, test.method="edger", FDR=param_FDR)
result <- getResult(tcc, sort=FALSE)
sum(tcc$stat$q.value < param_FDR)
tmp <- cbind(rownames(tcc$count), normalized, result)
tmp <- tmp[order(tmp$rank),]
write.table(tmp, out_f1, sep="\t", append=F, quote=F, row.names=F)
png(out_f2, pointsize=13, width=param_fig[1], height=param_fig[2])
par(mar=param_mar)
plot(tcc, FDR=param_FDR, xlim=c(-2, 17), ylim=c(-10.0, 11.5),
cex=0.8, cex.lab=1.2,
cex.axis=1.2, main="",
xlab="A = (log2(G2) + log2(G1))/2",
ylab="M = log2(G2) - log2(G1)")
obj <- as.logical(tcc$stat$q.value < param_FDR)
points(result$a.value[obj], result$m.value[obj],
pch=20, cex=1.5, col="magenta")
legend("topright", c(paste("DEG(FDR<", param_FDR, ")", sep=""), "non-DEG"),
col=c("magenta", "black"), pch=20, cex=1.2)
dev.off()
sum(tcc$stat$q.value < 0.05)
sum(tcc$stat$q.value < 0.10)
sum(tcc$stat$q.value < 0.20)
sum(tcc$stat$q.value < 0.30)
- RM1 vs. RM2 (スライド99)
「解析 | 発現変動 | 2群間 | 対応なし | 複製あり | Blekhmanデータ | TCC(Sun_2013)」
の例題10は以下と同じ。13, 16, 14, 17列目のデータのみ抽出して、
アカゲザル2サンプル(G1群:RMF1とRMM1) vs. アカゲザル2サンプル(G2群:RMF2とRMM2)の2群間比較を行います。
in_f <- "sample_blekhman_18.txt"
out_f1 <- "hoge10.txt"
out_f2 <- "hoge10.png"
param_subset <- c(13, 16, 14, 17)
param_G1 <- 2
param_G2 <- 2
param_FDR <- 0.05
param_fig <- c(430, 350)
param_mar <- c(4, 4, 0, 0)
library(TCC)
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
data <- data[,param_subset]
data.cl <- c(rep(1, param_G1), rep(2, param_G2))
tcc <- new("TCC", data, data.cl)
colnames(data)
tcc <- calcNormFactors(tcc, norm.method="tmm", test.method="edger",
iteration=3, FDR=0.1, floorPDEG=0.05)
normalized <- getNormalizedData(tcc)
tcc <- estimateDE(tcc, test.method="edger", FDR=param_FDR)
result <- getResult(tcc, sort=FALSE)
sum(tcc$stat$q.value < param_FDR)
tmp <- cbind(rownames(tcc$count), normalized, result)
tmp <- tmp[order(tmp$rank),]
write.table(tmp, out_f1, sep="\t", append=F, quote=F, row.names=F)
png(out_f2, pointsize=13, width=param_fig[1], height=param_fig[2])
par(mar=param_mar)
plot(tcc, FDR=param_FDR, xlim=c(-2, 17), ylim=c(-10.0, 11.5),
cex=0.8, cex.lab=1.2,
cex.axis=1.2, main="",
xlab="A = (log2(G2) + log2(G1))/2",
ylab="M = log2(G2) - log2(G1)")
obj <- as.logical(tcc$stat$q.value < param_FDR)
points(result$a.value[obj], result$m.value[obj],
pch=20, cex=1.5, col="magenta")
legend("topright", c(paste("DEG(FDR<", param_FDR, ")", sep=""), "non-DEG"),
col=c("magenta", "black"), pch=20, cex=1.2)
dev.off()
sum(tcc$stat$q.value < 0.05)
sum(tcc$stat$q.value < 0.10)
sum(tcc$stat$q.value < 0.20)
sum(tcc$stat$q.value < 0.30)
- RM1 vs. RM3 (スライド100)
「解析 | 発現変動 | 2群間 | 対応なし | 複製あり | Blekhmanデータ | TCC(Sun_2013)」
の例題11は以下と同じ。13, 16, 15, 18列目のデータのみ抽出して、
アカゲザル2サンプル(G1群:RMF1とRMM1) vs. アカゲザル2サンプル(G2群:RMF3とRMM3)の2群間比較を行います。
in_f <- "sample_blekhman_18.txt"
out_f1 <- "hoge11.txt"
out_f2 <- "hoge11.png"
param_subset <- c(13, 16, 15, 18)
param_G1 <- 2
param_G2 <- 2
param_FDR <- 0.05
param_fig <- c(430, 350)
param_mar <- c(4, 4, 0, 0)
library(TCC)
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
data <- data[,param_subset]
data.cl <- c(rep(1, param_G1), rep(2, param_G2))
tcc <- new("TCC", data, data.cl)
colnames(data)
tcc <- calcNormFactors(tcc, norm.method="tmm", test.method="edger",
iteration=3, FDR=0.1, floorPDEG=0.05)
normalized <- getNormalizedData(tcc)
tcc <- estimateDE(tcc, test.method="edger", FDR=param_FDR)
result <- getResult(tcc, sort=FALSE)
sum(tcc$stat$q.value < param_FDR)
tmp <- cbind(rownames(tcc$count), normalized, result)
tmp <- tmp[order(tmp$rank),]
write.table(tmp, out_f1, sep="\t", append=F, quote=F, row.names=F)
png(out_f2, pointsize=13, width=param_fig[1], height=param_fig[2])
par(mar=param_mar)
plot(tcc, FDR=param_FDR, xlim=c(-2, 17), ylim=c(-10.0, 11.5),
cex=0.8, cex.lab=1.2,
cex.axis=1.2, main="",
xlab="A = (log2(G2) + log2(G1))/2",
ylab="M = log2(G2) - log2(G1)")
obj <- as.logical(tcc$stat$q.value < param_FDR)
points(result$a.value[obj], result$m.value[obj],
pch=20, cex=1.5, col="magenta")
legend("topright", c(paste("DEG(FDR<", param_FDR, ")", sep=""), "non-DEG"),
col=c("magenta", "black"), pch=20, cex=1.2)
dev.off()
sum(tcc$stat$q.value < 0.05)
sum(tcc$stat$q.value < 0.10)
sum(tcc$stat$q.value < 0.20)
sum(tcc$stat$q.value < 0.30)
- 統計的手法とは (スライド103)
「解析 | 発現変動 | 2群間 | 対応なし | 複製あり | Blekhmanデータ | TCC(Sun_2013)」
の例題12は以下と同じ。1, 4, 13, 16列目のデータのみ抽出して、
ヒト2サンプル(G1群:HSF1とHSM1) vs. アカゲザル2サンプル(G2群:RMF1とRMM1)の2群間比較を行っています。
q.value = 1の"ENSG00000125844"と
"ENSG00000115325"をハイライトさせ、明らかにDEGでないものがどのあたりに位置するのかを示しています。
in_f <- "sample_blekhman_18.txt"
out_f1 <- "hoge12.txt"
out_f2 <- "hoge12.png"
param_subset <- c(1, 4, 13, 16)
param_G1 <- 2
param_G2 <- 2
param_FDR <- 0.05
param_fig <- c(400, 310)
param_mar <- c(4, 4, 0, 0)
param_geneid <- c("ENSG00000125844", "ENSG00000115325")
param_col <- "skyblue"
param_cex <- 2
param_pch <- 20
library(TCC)
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
data <- data[,param_subset]
data.cl <- c(rep(1, param_G1), rep(2, param_G2))
tcc <- new("TCC", data, data.cl)
dim(data)
head(data)
tcc <- calcNormFactors(tcc, norm.method="tmm", test.method="edger",
iteration=3, FDR=0.1, floorPDEG=0.05)
normalized <- getNormalizedData(tcc)
tcc <- estimateDE(tcc, test.method="edger", FDR=param_FDR)
result <- getResult(tcc, sort=FALSE)
sum(tcc$stat$q.value < param_FDR)
tmp <- cbind(rownames(tcc$count), normalized, result)
tmp <- tmp[order(tmp$rank),]
write.table(tmp, out_f1, sep="\t", append=F, quote=F, row.names=F)
png(out_f2, pointsize=13, width=param_fig[1], height=param_fig[2])
par(mar=param_mar)
plot(tcc, FDR=param_FDR, xlim=c(-2, 17), ylim=c(-10.0, 11.5),
cex=0.8, cex.lab=1.2,
cex.axis=1.2, main="",
xlab="A = (log2(G2) + log2(G1))/2",
ylab="M = log2(G2) - log2(G1)")
legend("topright", c(paste("DEG(FDR<", param_FDR, ")", sep=""), "non-DEG"),
col=c("magenta", "black"), pch=20, cex=1.2)
obj <- is.element(tcc$gene_id, param_geneid)
points(result$a.value[obj], result$m.value[obj],
pch=param_pch, cex=param_cex, col=param_col)
dev.off()
sum(tcc$stat$q.value < 0.05)
sum(tcc$stat$q.value < 0.10)
sum(tcc$stat$q.value < 0.20)
sum(tcc$stat$q.value < 0.30)
- 結果の比較(倍率変化とFDR) (スライド105)
「解析 | 発現変動 | 2群間 | 対応なし | 複製あり | Blekhmanデータ | TCC(Sun_2013)」
の例題13は以下と同じ。デフォルトのparam_subsetは、1, 4, 13, 16列目のデータのみ抽出して、
ヒト2サンプル(G1群:HSF1とHSM1) vs. アカゲザル2サンプル(G2群:RMF1とRMM1)の2群間比較を行っていますが、適宜変更してください。
「FDR閾値を満たすもの」と「fold-change閾値を満たすもの」それぞれのM-A plotを作成しています。
in_f <- "sample_blekhman_18.txt"
out_f1 <- "hoge13.txt"
out_f2 <- "hoge13_FDR.png"
out_f3 <- "hoge13_FC.png"
param_subset <- c(1, 4, 13, 16)
#param_subset <- c(1, 4, 7, 10)
#param_subset <- c(1, 4, 2, 5)
#param_subset <- c(7, 10, 8, 11)
#param_subset <- c(13, 16, 14, 17)
#param_subset <- c(13, 16, 15, 18)
param_G1 <- 2
param_G2 <- 2
param_FDR <- 0.05
param_FC <- 2
param_fig <- c(400, 310)
param_mar <- c(4, 4, 0, 0)
library(TCC)
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
data <- data[,param_subset]
data.cl <- c(rep(1, param_G1), rep(2, param_G2))
tcc <- new("TCC", data, data.cl)
dim(data)
head(data)
tcc <- calcNormFactors(tcc, norm.method="tmm", test.method="edger",
iteration=3, FDR=0.1, floorPDEG=0.05)
normalized <- getNormalizedData(tcc)
tcc <- estimateDE(tcc, test.method="edger", FDR=param_FDR)
result <- getResult(tcc, sort=FALSE)
sum(tcc$stat$q.value < param_FDR)
tmp <- cbind(rownames(tcc$count), normalized, result)
tmp <- tmp[order(tmp$rank),]
write.table(tmp, out_f1, sep="\t", append=F, quote=F, row.names=F)
png(out_f2, pointsize=13, width=param_fig[1], height=param_fig[2])
par(mar=param_mar)
plot(tcc, FDR=param_FDR, xlim=c(-2, 17), ylim=c(-10.0, 11.5),
cex=0.8, cex.lab=1.2,
cex.axis=1.2, main="",
xlab="A = (log2(G2) + log2(G1))/2",
ylab="M = log2(G2) - log2(G1)")
legend("topright", c(paste("DEG(FDR<", param_FDR, ")", sep=""), "non-DEG"),
col=c("magenta", "black"), pch=20, cex=1.2)
dev.off()
sum(tcc$stat$q.value < 0.05)
sum(tcc$stat$q.value < 0.10)
sum(tcc$stat$q.value < 0.20)
sum(tcc$stat$q.value < 0.30)
M <- getResult(tcc)$m.value
hoge <- rep(1, length(M))
hoge[abs(M) > log2(param_FC)] <- 2
cols <- c("black", "magenta")
png(out_f3, pointsize=13, width=param_fig[1], height=param_fig[2])
par(mar=param_mar)
plot(tcc, col=cols, col.tag=hoge, xlim=c(-2, 17), ylim=c(-10.0, 11.5),
cex=0.8, cex.lab=1.2,
cex.axis=1.2, main="",
xlab="A = (log2(G2) + log2(G1))/2",
ylab="M = log2(G2) - log2(G1)")
legend("topright", c(paste("DEG(", param_FC, "-fold)", sep=""), "non-DEG"),
col=c("magenta", "black"), pch=20)
dev.off()
sum(abs(M) > log2(16))
sum(abs(M) > log2(8))
sum(abs(M) > log2(4))
sum(abs(M) > log2(2))
- giladデータでDEG同定 (スライド109)
ReCountのgiladデータ(gilad_count_table.txt)で反復あり2群間比較を行う。
「解析 | 発現変動 | 2群間 | 対応なし | 複製あり | TCC(Sun_2013)」例題1をベースに作成。
M-A plot周辺はオプションなどを多少変更しています。
in_f <- "gilad_count_table.txt"
out_f1 <- "hoge1.txt"
out_f2 <- "hoge1.png"
param_G1 <- 3
param_G2 <- 3
param_FDR <- 0.05
param_fig <- c(430, 350)
library(TCC)
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
data.cl <- c(rep(1, param_G1), rep(2, param_G2))
tcc <- new("TCC", data, data.cl)
tcc <- calcNormFactors(tcc, norm.method="tmm", test.method="edger",
iteration=3, FDR=0.1, floorPDEG=0.05)
tcc <- estimateDE(tcc, test.method="edger", FDR=param_FDR)
result <- getResult(tcc, sort=FALSE)
sum(tcc$stat$q.value < param_FDR)
tmp <- cbind(rownames(tcc$count), tcc$count, result)
write.table(tmp, out_f1, sep="\t", append=F, quote=F, row.names=F)
png(out_f2, pointsize=13, width=param_fig[1], height=param_fig[2])
par(mar=c(4, 4, 0, 0))
plot(tcc, FDR=param_FDR, main="",
xlab="A = (log2(G2) + log2(G1))/2",
ylab="M = log2(G2) - log2(G1)")
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)
sum(tcc$stat$q.value < 0.10)
sum(tcc$stat$q.value < 0.20)
sum(tcc$stat$q.value < 0.30)
- クラスタリング(maqc) (スライド116)
ReCountのmaqcデータ(
maqc_count_table.txtと
maqc_phenodata.txt)
を入力としてサンプル間クラスタリング。
in_f1 <- "maqc_count_table.txt"
in_f2 <- "maqc_phenodata.txt"
out_f <- "hoge.png"
param_fig <- c(600, 400)
library(TCC)
data <- read.table(in_f1, header=TRUE, row.names=1, sep="\t", quote="")
phenotype <- read.table(in_f2, header=TRUE, row.names=1, sep=" ", quote="")
phenotype
colnames(data) <- phenotype$tissue
out <- clusterSample(data, dist.method="spearman",
hclust.method="average", unique.pattern=TRUE)
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()
- クラスタリング(blekhman) (スライド118)
「解析 | クラスタリング | サンプル間 | TCC(Sun_2013)」の例題7をベースに作成。
入力は、20,689 genes×36 samplesのカウントデータ(sample_blekhman_36.txt)です。
縦軸の範囲を[0, 0.28]にするやり方を誰か教えて下さいm(_ _)m
(「○○氏提供情報」とさせていただきます。)
in_f <- "sample_blekhman_36.txt"
out_f <- "hoge7.png"
param_fig <- c(700, 400)
param_yrange <- c(0, 0.28)
library(TCC)
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
dim(data)
out <- clusterSample(data, dist.method="spearman",
hclust.method="average", unique.pattern=TRUE)
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",
ylim=param_yrange)
dev.off()
- maqcデータでDEG同定 (スライド120)
ReCountのmaqcデータ(
maqc_count_table.txt)で反復あり2群間比較を行う。
「解析 | 発現変動 | 2群間 | 対応なし | 複製あり | TCC(Sun_2013)」例題1をベースに作成。
M-A plot周辺はオプションなどを多少変更しています。
in_f <- "maqc_count_table.txt"
out_f1 <- "hoge1.txt"
out_f2 <- "hoge1.png"
param_G1 <- 7
param_G2 <- 7
param_FDR <- 0.05
param_fig <- c(430, 350)
library(TCC)
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
data.cl <- c(rep(1, param_G1), rep(2, param_G2))
tcc <- new("TCC", data, data.cl)
tcc <- calcNormFactors(tcc, norm.method="tmm", test.method="edger",
iteration=3, FDR=0.1, floorPDEG=0.05)
tcc <- estimateDE(tcc, test.method="edger", FDR=param_FDR)
result <- getResult(tcc, sort=FALSE)
sum(tcc$stat$q.value < param_FDR)
tmp <- cbind(rownames(tcc$count), tcc$count, result)
write.table(tmp, out_f1, sep="\t", append=F, quote=F, row.names=F)
png(out_f2, pointsize=13, width=param_fig[1], height=param_fig[2])
par(mar=c(4, 4, 0, 0))
plot(tcc, FDR=param_FDR, main="",
xlab="A = (log2(G2) + log2(G1))/2",
ylab="M = log2(G2) - log2(G1)")
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)
sum(tcc$stat$q.value < 0.10)
sum(tcc$stat$q.value < 0.20)
sum(tcc$stat$q.value < 0.30)
- DEG数の見積もり(maqcデータ) (スライド124)
sum(tcc$stat$q.value < 0.05)
sum(tcc$stat$q.value < 0.10)
sum(tcc$stat$q.value < 0.20)
sum(tcc$stat$q.value < 0.30)
7618*(1 - 0.05)
7843*(1 - 0.10)
8114*(1 - 0.20)
8338*(1 - 0.30)
- mergeの基本形 (スライド127)
technical replicatesのデータを足すだけ。
in_f <- "maqc_count_table.txt"
param_G1 <- 7
param_G2 <- 7
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
data.cl <- c(rep(1, param_G1), rep(2, param_G2))
head(data[, data.cl == 1], n=3)
head(data[, data.cl == 2], n=3)
data_G1 <- rowSums(data[, data.cl == 1])
data_G2 <- rowSums(data[, data.cl == 2])
data <- cbind(data_G1, data_G2)
dim(data)
head(data)
dim(unique(data))
sum(rowSums(data) > 0)
- maqc (pooled) (スライド133)
getwd()
list.files(pattern="maqc_pool")
- maqc (pooled)でDEG同定 (スライド135)
「解析 | 発現変動 | 2群間 | 対応なし | 複製なし | TCC(Sun_2013)」
の例題5は以下と同じ。ReCountのmaqc (pooled)データ(
maqc_pooledreps_count_table.txt)です。
52,580 genes×2 samplesのカウントデータ(G1群1サンプル vs. G2群1サンプル)です。
in_f <- "maqc_pooledreps_count_table.txt"
out_f1 <- "hoge5.txt"
out_f2 <- "hoge5.png"
param_G1 <- 1
param_G2 <- 1
param_FDR <- 0.05
param_fig <- c(400, 380)
library(TCC)
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
data.cl <- c(rep(1, param_G1), rep(2, param_G2))
tcc <- new("TCC", data, data.cl)
tcc <- calcNormFactors(tcc, norm.method="deseq2", test.method="deseq2",
iteration=3, FDR=0.1, floorPDEG=0.05)
normalized <- getNormalizedData(tcc)
tcc <- estimateDE(tcc, test.method="deseq2", FDR=param_FDR)
result <- getResult(tcc, sort=FALSE)
tmp <- cbind(rownames(tcc$count), normalized, result)
tmp <- tmp[order(tmp$rank),]
write.table(tmp, out_f1, sep="\t", append=F, quote=F, row.names=F)
png(out_f2, pointsize=13, width=param_fig[1], height=param_fig[2])
par(mar=c(4, 4, 0, 0))
plot(tcc, FDR=param_FDR, main="",
xlab="A = (log2(G2) + log2(G1))/2",
ylab="M = log2(G2) - log2(G1)")
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)
sum(tcc$stat$q.value < 0.10)
sum(tcc$stat$q.value < 0.20)
sum(tcc$stat$q.value < 0.30)
- 反復あり or なし(TCC) (スライド143)
TCCパッケージ(Sun et al., BMC Bioinformatics, 2013)のやり方です。
in_f <- "maqc_count_table.txt"
param_G1 <- 7
param_G2 <- 7
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
data.cl <- c(rep(1, param_G1), rep(2, param_G2))
library(TCC)
tcc <- new("TCC", data, data.cl)
tcc <- calcNormFactors(tcc, norm.method="tmm", test.method="edger",
iteration=3, FDR=0.1, floorPDEG=0.05)
tcc <- estimateDE(tcc, test.method="edger")
p.value <- tcc$stat$p.value
p.value[is.na(p.value)] <- 1
ranking <- rank(p.value)
ranking_TCC_ari <- ranking
q.value <- tcc$stat$q.value
sum(q.value < 0.05)
sum(q.value < 0.10)
sum(q.value < 0.20)
sum(q.value < 0.30)
sum(q.value < 0.50)
sum(q.value < 0.70)
sum(q.value < 0.80)
in_f <- "maqc_pooledreps_count_table.txt"
param_G1 <- 1
param_G2 <- 1
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
data.cl <- c(rep(1, param_G1), rep(2, param_G2))
library(TCC)
tcc <- new("TCC", data, data.cl)
tcc <- calcNormFactors(tcc, norm.method="deseq2", test.method="deseq2",
iteration=3, FDR=0.1, floorPDEG=0.05)
tcc <- estimateDE(tcc, test.method="deseq2")
p.value <- tcc$stat$p.value
p.value[is.na(p.value)] <- 1
ranking <- rank(p.value)
ranking_TCC_nasi <- ranking
q.value <- tcc$stat$q.value
sum(q.value < 0.05)
sum(q.value < 0.10)
sum(q.value < 0.20)
sum(q.value < 0.30)
sum(q.value < 0.50)
sum(q.value < 0.70)
sum(q.value < 0.80)
- 反復あり or なし(edgeR) (スライド145)
edgeRパッケージ(Robinson et al., Bioinformatics, 2010)のやり方です。
in_f <- "maqc_count_table.txt"
param_G1 <- 7
param_G2 <- 7
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
data <- as.matrix(data)
data.cl <- c(rep(1, param_G1), rep(2, param_G2))
library(edgeR)
d <- DGEList(counts=data, group=data.cl)
d <- calcNormFactors(d)
d <- estimateCommonDisp(d)
d <- estimateTagwiseDisp(d)
out <- exactTest(d)
p.value <- out$table$PValue
ranking <- rank(p.value)
ranking_edgeR_ari <- ranking
q.value <- p.adjust(p.value, method="BH")
sum(q.value < 0.05)
sum(q.value < 0.10)
sum(q.value < 0.20)
sum(q.value < 0.30)
sum(q.value < 0.50)
sum(q.value < 0.70)
sum(q.value < 0.80)
in_f <- "maqc_pooledreps_count_table.txt"
param_G1 <- 1
param_G2 <- 1
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
data <- as.matrix(data)
data.cl <- c(rep(1, param_G1), rep(2, param_G2))
library(edgeR)
d <- DGEList(counts=data, group=data.cl)
d <- calcNormFactors(d)
d <- estimateGLMCommonDisp(d, method="deviance", robust=TRUE, subset=NULL)
out <- exactTest(d)
p.value <- out$table$PValue
ranking <- rank(p.value)
ranking_edgeR_nasi <- ranking
q.value <- p.adjust(p.value, method="BH")
sum(q.value < 0.05)
sum(q.value < 0.10)
sum(q.value < 0.20)
sum(q.value < 0.30)
sum(q.value < 0.50)
sum(q.value < 0.70)
sum(q.value < 0.80)
- 反復あり or なし(DESeq2) (スライド146)
DESeq2パッケージ(Love et al., Genome Biol., 2014)のやり方です。
in_f <- "maqc_count_table.txt"
param_G1 <- 7
param_G2 <- 7
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
data <- as.matrix(data)
data.cl <- c(rep(1, param_G1), rep(2, param_G2))
library(DESeq2)
colData <- data.frame(condition=as.factor(data.cl))
d <- DESeqDataSetFromMatrix(countData=data, colData=colData, design=~condition)
d <- DESeq(d)
tmp <- results(d)
p.value <- tmp$pvalue
p.value[is.na(p.value)] <- 1
ranking <- rank(p.value)
ranking_DESeq2_ari <- ranking
q.value <- tmp$padj
q.value[is.na(q.value)] <- 1
#q.value <- p.adjust(p.value, method="BH")
sum(q.value < 0.05)
sum(q.value < 0.10)
sum(q.value < 0.20)
sum(q.value < 0.30)
sum(q.value < 0.50)
sum(q.value < 0.70)
sum(q.value < 0.80)
in_f <- "maqc_pooledreps_count_table.txt"
param_G1 <- 1
param_G2 <- 1
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
data <- as.matrix(data)
data.cl <- c(rep(1, param_G1), rep(2, param_G2))
library(DESeq2)
colData <- data.frame(condition=as.factor(data.cl))
d <- DESeqDataSetFromMatrix(countData=data, colData=colData, design=~condition)
d <- DESeq(d)
tmp <- results(d)
p.value <- tmp$pvalue
p.value[is.na(p.value)] <- 1
ranking <- rank(p.value)
ranking_DESeq2_nasi <- ranking
q.value <- tmp$padj
q.value[is.na(q.value)] <- 1
#q.value <- p.adjust(p.value, method="BH")
sum(q.value < 0.05)
sum(q.value < 0.10)
sum(q.value < 0.20)
sum(q.value < 0.30)
sum(q.value < 0.50)
sum(q.value < 0.70)
sum(q.value < 0.80)
- 反復あり or なし(DESeq2) (スライド148)
DESeq2パッケージ
(Love et al., Genome Biol., 2014)のやり方です。
DESeq2オリジナルのq-valueではなく、p-valueからp.adjust関数を用いてq.valueを計算しています。
in_f <- "maqc_count_table.txt"
param_G1 <- 7
param_G2 <- 7
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
data <- as.matrix(data)
data.cl <- c(rep(1, param_G1), rep(2, param_G2))
library(DESeq2)
colData <- data.frame(condition=as.factor(data.cl))
d <- DESeqDataSetFromMatrix(countData=data, colData=colData, design=~condition)
d <- DESeq(d)
tmp <- results(d)
p.value <- tmp$pvalue
p.value[is.na(p.value)] <- 1
ranking <- rank(p.value)
ranking_DESeq2_ari <- ranking
#q.value <- tmp$padj
#q.value[is.na(q.value)] <- 1
q.value <- p.adjust(p.value, method="BH")
sum(q.value < 0.05)
sum(q.value < 0.10)
sum(q.value < 0.20)
sum(q.value < 0.30)
sum(q.value < 0.50)
sum(q.value < 0.70)
sum(q.value < 0.80)
in_f <- "maqc_pooledreps_count_table.txt"
param_G1 <- 1
param_G2 <- 1
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
data <- as.matrix(data)
data.cl <- c(rep(1, param_G1), rep(2, param_G2))
library(DESeq2)
colData <- data.frame(condition=as.factor(data.cl))
d <- DESeqDataSetFromMatrix(countData=data, colData=colData, design=~condition)
d <- DESeq(d)
tmp <- results(d)
p.value <- tmp$pvalue
p.value[is.na(p.value)] <- 1
ranking <- rank(p.value)
ranking_DESeq2_nasi <- ranking
#q.value <- tmp$padj
#q.value[is.na(q.value)] <- 1
q.value <- p.adjust(p.value, method="BH")
sum(q.value < 0.05)
sum(q.value < 0.10)
sum(q.value < 0.20)
sum(q.value < 0.30)
sum(q.value < 0.50)
sum(q.value < 0.70)
sum(q.value < 0.80)