This repository has been archived by the owner on Apr 26, 2024. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathDiscordant_concordant_pairs.R
116 lines (103 loc) · 3.62 KB
/
Discordant_concordant_pairs.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
setwd("~/Documents/IFN_enhancer/R/")
library(GenomicRanges)
library(GenomicFeatures)
library(gplots)
library(rtracklayer)
# Discordant/concordant enhancer list
# Load RPKM matrix
load("Data/eRNA_expression_GM_unfiltered.RData")
gene.rpkm <- read.table("Data/Expression_matrix.symbol.csv")
rpkm <- rbind(eRNA.rpkm, gene.rpkm)
#################################################################################################################
# Heatmaps
plot_heatmap <- function(l_dc, l_cc) {
gene_list <- c(l_dc, l_cc)
heat.data <- rpkm[gene_list, ]
# heat.data <- log2((heat.data+0.25) /(heat.data[,1]+0.25))
# heat.data <- (heat.data - rowMeans(heat.data)) / apply(heat.data, 1, sd)
heat.data <- (heat.data - heat.data[, 1]) / apply(heat.data, 1, sd)
heat.data <- as.matrix(heat.data)
for (i in c(1:length(gene_list))) {
rownames(heat.data)[i] <- substr(gene_list[i], 14, 1000L)
}
for (i in c(1:ncol(heat.data))) {
colnames(heat.data)[i] <- substr(colnames(heat.data)[i], 3, 1000L)
}
my_palette <- colorRampPalette(c("blue", "white", "yellow"), bias = 1)(n = 100)
breaks <- c(seq(-1, 3, length.out = 101))
heatmap.2(heat.data,
col = my_palette, breaks = breaks,
Colv = F, Rowv = F, dendrogram = "none", labRow = "",
RowSideColors = c(rep("darkgrey", length(l_dc)), rep("lightgrey", length(l_cc))),
keysize = 1.2, density.info = "density", key.ylab = "", key.xlab = "", key.title = "",
trace = "none", scale = "none"
)
}
plot_heatmap(dc[, 2], cc[, 2])
amplitude_index <- function(x, abs = T) {
fc <- log2((x[, -1] + 0.1) / (x[, 1] + 0.1))
if (abs) {
return(apply(abs(fc), 1, max))
}
else {
lab <- apply(abs(fc), 1, function(x) {
xx <- which(x == max(x))
return(xx[1])
})
tp <- NULL
for (i in 1:length(lab)) {
tp <- c(tp, fc[i, lab[i]])
}
names(tp) <- rownames(x)
return(tp)
}
}
get_idx <- function(pairs, method = 1) {
# Method 1:
if (method == 1) {
res <- c()
for (i in 1:nrow(pairs)) {
res <- c(res, cor(t(rpkm[pairs[i, 1], ]), t(rpkm[pairs[i, 2], ]), method = "spearman"))
}
return(res)
}
# Method 2:
if (method == 2) {
res <- c()
for (i in 1:nrow(pairs)) {
res <- c(res, max(as.vector(rpkm[pairs[i, 2], 5:9])) / max(as.vector(rpkm[pairs[i, 2], 10:12])))
}
return(res)
}
# Method 3:
if (method == 3) {
res <- c()
for (i in 1:nrow(pairs)) {
res <- c(res, apply(rpkm[pairs[i, 2], 5:9], 1, mean) / apply(rpkm[pairs[i, 2], 10:12], 1, mean))
}
return(res)
}
}
idx_dc <- get_idx(dc, 1)
idx_cc <- get_idx(cc, 1)
hist(idx_dc, col = rgb(0.7, 0, 0, 0.5), xlim = range(c(idx_dc, idx_cc)), breaks = 20)
hist(idx_cc, col = rgb(0, 0, 0.7, 0.5), add = T, breaks = 60)
pairs <- read.csv(file = "Data/Induced_EP.csv", colClasses = "character")[, c("symbol", "id_b")]
pairs <- pairs[!(is.na(pairs[, 1])), ]
ai_e <- amplitude_index(rpkm[pairs[, 2], 1:5], F)
ai_g <- amplitude_index(rpkm[pairs[, 1], 1:5], F)
pairs <- pairs[(ai_e > 1) & (ai_g > 1), ]
idx <- get_idx(pairs, 1)
dc <- pairs[idx < quantile(idx, 0.3), ]
cc <- pairs[idx > quantile(idx, 0.7), ]
plot_heatmap(dc[, 1], cc[, 1])
plot_heatmap(dc[, 2], cc[, 2])
idx[pairs[, 1] == "IFNB1"]
# dc <- read.table("~/Documents/eRNA/link/IFN project/discordant_pairs_gene-control-enhancer.pairs.txt", header = F, colClasses = 'character')
# cc <- read.table("~/Documents/eRNA/link/IFN project/concordant_pairs_gene-control-enhancer.pairs.txt", header = F, colClasses = 'character')
pdf("Figure/Discordant_concordant_pairs_gene.pdf")
plot_heatmap(dc[, 1], cc[, 1])
dev.off()
pdf("Figure/Discordant_concordant_pairs_eRNA.pdf")
plot_heatmap(dc[, 2], cc[, 2])
dev.off()