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 path4C_analysis_v2.R
180 lines (166 loc) · 6.08 KB
/
4C_analysis_v2.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
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
# Last update: Feb 6 2018
# source("https://bioconductor.org/biocLite.R")
# biocLite("BSgenome.Hsapiens.UCSC.hg19")
# biocLite("BSgenome.Hsapiens.UCSC.hg19.masked")
# biocLite("BSgenome.Mmusculus.UCSC.mm9")
# biocLite("BSgenome.Mmusculus.UCSC.mm9.masked")
# biocLite("r3Cseq")
setwd("~/Documents/IFN_enhancer/R/Data/4C_01292018/")
library(r3Cseq)
library(BSgenome.Hsapiens.UCSC.hg19)
library(RColorBrewer)
#######################################################
# Method_1
#######################################################
track_plot <- function(my_obj, tag, scope_size = 200, xaxt = FALSE, ylim = c(0, 4500)) {
viewpoint <- getViewpoint(my_obj)
anchor_point <- as.integer((end(viewpoint) + start(viewpoint)) / 2)
scope <- GRanges(viewpoint)
start(scope) <- anchor_point - scope_size * 10^3
end(scope) <- anchor_point + scope_size * 10^3
frags <- GRanges(expRPM(my_obj))
hits <- findOverlaps(frags, scope)
frags <- frags[queryHits(hits)]
x <- as.integer((end(frags) + start(frags)) / 2) - anchor_point
x <- x / 1000
y <- mcols(frags)["RPMs"][, 1]
if (xaxt) {
plot(x, y, type = "l", ylim = ylim, ylab = "RPM", xlab = "")
}
else {
plot(x, y, type = "l", xaxt = "n", ylim = ylim, ylab = "RPM", xlab = "")
}
legend("topright", tag, bty = "n")
}
method_1 <- function(experiment, viewpoint_chromosome, viewpoint_primer_forward, viewpoint_primer_reverse, restrictionEnzyme, distance = 2e5, ylim = c(0, 4000)) {
outfile <- paste(experiment, as.character(distance), sep = "")
outfile <- paste(outfile, ".pdf", sep = "")
pdf(outfile)
tags <- paste(experiment, c("0h", "6h", "12h", "18h", "24h"), sep = "")
# tags <- paste(experiment, c("0h", "6h"), sep="")
par(mfrow = c(length(tags), 1), mai = c(0, 1, 0, 1), omi = c(1, 0, 1, 0))
i <- 0
obj.list <- list()
for (tag in tags) {
expFile <- paste(tag, ".trim.bam", sep = "")
expFile <- paste("bam/", expFile, sep = "")
my_obj <- new("r3Cseq",
organismName = "hg19", alignedReadsBamExpFile = expFile,
isControlInvolved = F, viewpoint_chromosome = viewpoint_chromosome,
viewpoint_primer_forward = viewpoint_primer_forward,
viewpoint_primer_reverse = viewpoint_primer_reverse,
expLabel = tag, restrictionEnzyme = restrictionEnzyme
)
getRawReads(my_obj)
getReadCountPerRestrictionFragment(my_obj)
calculateRPM(my_obj)
getInteractions(my_obj)
expResults <- expInteractionRegions(my_obj)
# export3Cseq2bedGraph(my_obj)
i <- i + 1
obj.list[[i]] <- my_obj
if (i == length(tags)) {
track_plot(my_obj, tag = tag, scope_size = distance, xaxt = TRUE, ylim = ylim)
}
else {
track_plot(my_obj, tag = tag, scope_size = distance, xaxt = FALSE, ylim = ylim)
}
}
dev.off()
return(obj.list)
}
local_interaction <- function(my_obj, scope_size = 1000) {
viewpoint <- getViewpoint(my_obj)
anchor_point <- as.integer((end(viewpoint) + start(viewpoint)) / 2)
scope <- GRanges(viewpoint)
start(scope) <- anchor_point - scope_size * 10^3
end(scope) <- anchor_point + scope_size * 10^3
frags <- GRanges(expRPM(my_obj))
hits <- findOverlaps(frags, scope)
frags <- frags[queryHits(hits)]
x <- as.integer((end(frags) + start(frags)) / 2) - anchor_point
x <- x / 1000
y <- mcols(frags)["RPMs"][, 1]
return(list(x, y))
}
#############################################################################################
timepoints <- c("0h", "6h", "12h", "18h", "24h")
IFNb1_ED <- method_1(
experiment <- "IFNb1-ED_",
viewpoint_chromosome = "chr9",
viewpoint_primer_forward = "CACTCAGTCTATTAGCATTAGATC",
viewpoint_primer_reverse = "AATATGTTTAGGTTTTATGGCTGGCTTTGGGGAAG",
restrictionEnzyme = "EcoRI",
distance = 500
)
save(IFNb1_ED, file = "IFNB1_ED.Robj")
track_plot(IFNb1_ED, "IFNb1-ED_", scope_size = 200, xaxt = FALSE, ylim = c(0, 4500))
IFNb1_HD <- method_1(
experiment <- "IFNb1-HD_",
viewpoint_chromosome = "chr9",
viewpoint_primer_forward = "AACGCAGATAATAAGATC",
viewpoint_primer_reverse = "AAGCTTCAGAGGTCAAAGGGCTAAT",
restrictionEnzyme = "HindIII",
distance = 500
)
save(IFNb1_HD, file = "IFNB1_HD.Robj")
TNFSF10 <- method_1(
experiment <- "TNFSF10-1_",
viewpoint_chromosome = "chr3",
viewpoint_primer_forward = "CTCCATGGGTCAGGTTGTTTTCTTGATGAA",
viewpoint_primer_reverse = "CCAATGCTTCTCCCATGGATC",
restrictionEnzyme = "EcoRI",
distance = 1000,
ylim = c(0, 5000)
)
save(TNFSF10, file = "TNFSF10.Robj")
####################################################################################################
load("IFNB1_ED.Robj")
load("TNFSF10.Robj")
bar.col <- brewer.pal(3, "Blues")[3:1]
pdf("Total_local_ineractions.pdf")
local.sum <- c()
dist.list <- c(10000, 1000, 200)
dist.names <- c("10M", "1M", "200kb")
for (my_obj in IFNb1_ED) {
tp <- local_interaction(my_obj, scope_size = max(dist.list))
for (d in dist.list) {
local.sum <- append(local.sum, sum(tp[[2]][abs(tp[[1]]) < d]))
}
}
local.sum <- matrix(local.sum, nrow = 3, byrow = F, dimnames = list(dist.names, timepoints))
barplot(local.sum,
names.arg = timepoints, beside = T, ylab = "Interactions (RPM)", ylim = c(0, 55000),
col = bar.col, main = "IFNb1"
)
legend("topright", dist.names, col = bar.col, pch = 15)
local.sum <- c()
for (my_obj in TNFSF10) {
tp <- local_interaction(my_obj, scope_size = max(dist.list))
for (d in dist.list) {
local.sum <- append(local.sum, sum(tp[[2]][abs(tp[[1]]) < d]))
}
}
local.sum <- matrix(local.sum, nrow = 3, byrow = F, dimnames = list(dist.names, timepoints))
barplot(local.sum,
names.arg = timepoints, beside = T, ylab = "Interactions (RPM)", ylim = c(0, 55000),
col = bar.col, main = "TNFSF10"
)
legend("topright", dist.names, col = bar.col, pch = 15)
dev.off()
# method_1(
# experiment <- "IFNb1-HD_",
# viewpoint_chromosome='chr9',
# viewpoint_primer_forward='AACGCAGATAATAAGATC',
# viewpoint_primer_reverse='AAGCTTCAGAGGTCAAAGGGCTAAT',
# restrictionEnzyme='HindIII',
# distance=50
# )
# method_1(
# experiment <- "PARP94_",
# viewpoint_chromosome='chr3',
# viewpoint_primer_forward='TTTACCCTATTGTGGCATATATGAATTC',
# viewpoint_primer_reverse='GATCTGCACCTGAGCCAAAGAAAC',
# restrictionEnzyme='EcoRI',
# distance=50
# )