-
Notifications
You must be signed in to change notification settings - Fork 0
/
analys_DESeq2genes.r
117 lines (91 loc) · 4.45 KB
/
analys_DESeq2genes.r
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
# source("https://bioconductor.org/biocLite.R")
# biocLite(c("edgeR"))
library(DESeq2)
#library(yaml)
tb <- "TLX3vsRAGvsTAP_featCounts_genes.txt"
counts <- read.table(tb, header=TRUE, row.names=1)
counts <- counts[ ,6:ncol(counts)]
counts <- as.matrix(counts)
condition <- factor(c(rep("rag", 3), rep("tlx", 3), rep("tap", 3)))
coldata <- data.frame(row.names=colnames(counts), condition)
# Run the DESeq pipeline
dds <- DESeqDataSetFromMatrix(countData=counts, colData=coldata, design=~condition)
dds <- DESeq(dds)
#plotDispEsts(dds, main="Dispersion plot")
## === Get differential expression results ===
# res <- results(dds)
# === TLX3 vs RAG ===
res_TLX_RAG <- results(dds, contrast=c("condition","tlx","rag"))
#table(res$padj<0.05)
## Order by adjusted p-value
res_TLX_RAG <- res_TLX_RAG[order(res_TLX_RAG$padj), ]
## Merge with normalized count data
res_TLX_RAGout <- merge(as.data.frame(res_TLX_RAG), as.data.frame(counts(dds, normalized=TRUE)[,1:6]), by="row.names", sort=FALSE)
names(res_TLX_RAGout)[1] <- "Gene"
## Write results
write.table(res_TLX_RAGout, file="TLX3vsRAG-results_genes.txt", row.names = FALSE, quote = FALSE, sep = "\t")
# === TAP vs RAG ===
res_TAP_RAG <- results(dds, contrast=c("condition","tap","rag"))
#table(res$padj<0.05)
## Order by adjusted p-value
res_TAP_RAG <- res_TAP_RAG[order(res_TAP_RAG$padj), ]
## Merge with normalized count data
res_TAP_RAGout <- merge(as.data.frame(res_TAP_RAG), as.data.frame(counts(dds, normalized=TRUE)[,c(1,2,3,7,8,9)]), by="row.names", sort=FALSE)
names(res_TAP_RAGout)[1] <- "Gene"
## Write results
write.table(res_TAP_RAGout, file="TAPvsRAG-results_genes.txt", row.names = FALSE, quote = FALSE, sep = "\t")
# === TAP vs TLX3
res_TAP_TLX3 <- results(dds, contrast=c("condition","tap","tlx"))
## Order by adjusted p-value
res_TAP_TLX3 <- res_TAP_TLX3[order(res_TAP_TLX3$padj), ]
## Merge with normalized count data
res_TAP_TLX3out <- merge(as.data.frame(res_TAP_TLX3), as.data.frame(counts(dds, normalized=TRUE)[,4:9]), by="row.names", sort=FALSE)
names(res_TAP_TLX3out)[1] <- "Gene"
## Write results
write.table(res_TAP_TLX3out, file="TAPvsTLX3-results_genes.txt", row.names = FALSE, quote = FALSE, sep = "\t")
# cooks = assays(dds)[["cooks"]]
###################### ===================== ##################################################################
#resdata["ENSMUST00000103276", "RAGS.RAGZ"]
#res <- results(dds)
#hist(res$pvalue, breaks=50, col="grey")
# library(gplots)
# library(RColorBrewer)
# library("pheatmap")
#
# # Regularized log transformation for clustering/heatmaps, etc
# rld <- rlogTransformation(dds)
# mycols <- brewer.pal(8, "Dark2")[1:length(unique(condition))]
# # Sample distance heatmap
# sampleDists <- as.matrix(dist(t(assay(rld))))
# heatmap.2(as.matrix(sampleDists), key=F, trace="none",
# col=colorpanel(100, "black", "white"),
# ColSideColors=mycols[condition], RowSideColors=mycols[condition],
# margin=c(10, 10), main="Sample Distance Matrix")
#
# plotPCA(rld, intgroup=c("condition"))
#
#
# DESeq2::plotMA(dds, ylim=c(-15,15), cex=1)
# volcanoplot <- function (res, lfcthresh=2, sigthresh=0.05, main="Volcano Plot", legendpos="bottomright", labelsig=TRUE, textcx=1, ...) {
# with(res, plot(log2FoldChange, -log10(pvalue), pch=20, main=main, ...))
# with(subset(res, padj<sigthresh ), points(log2FoldChange, -log10(pvalue), pch=20, col="red", ...))
# with(subset(res, abs(log2FoldChange)>lfcthresh), points(log2FoldChange, -log10(pvalue), pch=20, col="orange", ...))
# with(subset(res, padj<sigthresh & abs(log2FoldChange)>lfcthresh), points(log2FoldChange, -log10(pvalue), pch=20, col="green", ...))
# if (labelsig) {
# require(calibrate)
# #with(subset(res, padj<sigthresh & abs(log2FoldChange)>lfcthresh), textxy(log2FoldChange, -log10(pvalue), labs=Gene, cex=textcx, ...))
# }
# legend(legendpos, xjust=1, yjust=1, legend=c(paste("FDR<",sigthresh,sep=""), paste("|LogFC|>",lfcthresh,sep=""), "both"), pch=20, col=c("red","orange","green"))
# }
# #png("diffexpr-volcanoplot.png", 1200, 1000, pointsize=20)
# volcanoplot(resdata, lfcthresh=1, sigthresh=0.05, textcx=.8, xlim=c(-15, 15))
#dev.off()
# select <- order(rowMeans(counts(dds,normalized=TRUE)), decreasing=TRUE)[1:55]
# nt <- normTransform(dds) # defaults to log2(x+1)
# df <- as.data.frame(colData(dds))
#
# pheatmap(assay(nt)[select,], cluster_rows=FALSE, show_rownames=TRUE, cluster_cols=FALSE, annotation_col=df)
# ===== EdgeR =====
# library(edgeR)
#
# y <- DGEList(counts=counts, group=condition)