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_Nov_2017.R
95 lines (80 loc) · 3.5 KB
/
4C_analysis_Nov_2017.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
setwd("~/Documents/IFN_enhancer/R/Data/4C_Nov_2017/")
library(r3Cseq)
library(BSgenome.Hsapiens.UCSC.hg19)
#######################################################
# Method_1
#######################################################
tags <- paste("L2_", c("0h", "6h", "12h", "18h", "24h"), sep = "")
for (tag in tags) {
expFile <- paste(tag, ".sorted.bam", sep = "")
my_obj <- new("r3Cseq",
organismName = "hg19", alignedReadsBamExpFile = expFile,
isControlInvolved = F, viewpoint_chromosome = "chr9",
viewpoint_primer_forward = "gccccattcccagagtaaggaatgttgcacaaacaaacccagaagaatctag",
viewpoint_primer_reverse = "ggaggaggttgaagtaatatgtttaggttttatggctggctttggggaag",
expLabel = tag, restrictionEnzyme = "EcoRI"
)
getRawReads(my_obj)
getReadCountPerRestrictionFragment(my_obj)
calculateRPM(my_obj)
getInteractions(my_obj)
expResults <- expInteractionRegions(my_obj)
export3Cseq2bedGraph(my_obj)
pdf(paste(tag, ".pdf", sep = ""))
par(mfrow = c(1, 1))
plotOverviewInteractions(my_obj)
plotInteractionsNearViewpoint(my_obj, distance = 2e5)
dev.off()
tp <- as.data.frame(expResults)
tp <- tp[tp["space"] == "chr9", ]
bed_data <- data.frame(chr = tp[, 1], start = tp[, 2], end = tp[, 3], name = row.names(tp), score = tp[, 6] / max(tp[, 6]) * 1000)
write.table(bed_data, file = paste(tag, ".chr9_interactions.bed", sep = ""), quote = F, row.names = F, col.names = F)
}
#######################################################
# Method_1: time points as replicates
#######################################################
# 'r3CSeqInBatch' kept asking control bam files, which doesn't make sense for our purpose. So I decided to go with a merged bam file
my_obj <- new("r3Cseq",
organismName = "hg19", alignedReadsBamExpFile = "merged.bam",
isControlInvolved = F, viewpoint_chromosome = "chr9",
viewpoint_primer_forward = "gccccattcccagagtaaggaatgttgcacaaacaaacccagaagaatctag",
viewpoint_primer_reverse = "ggaggaggttgaagtaatatgtttaggttttatggctggctttggggaag",
expLabel = tag, restrictionEnzyme = "EcoRI"
)
getRawReads(my_obj)
getReadCountPerRestrictionFragment(my_obj)
calculateRPM(my_obj)
getInteractions(my_obj)
expResults <- expInteractionRegions(my_obj)
export3Cseq2bedGraph(my_obj)
pdf(paste("merged.pdf", sep = ""))
par(mfrow = c(1, 1))
plotOverviewInteractions(my_obj)
plotInteractionsNearViewpoint(my_obj, distance = 2e5)
dev.off()
tp <- as.data.frame(expResults)
tp <- tp[tp["space"] == "chr9", ]
bed_data <- data.frame(chr = tp[, 1], start = tp[, 2], end = tp[, 3], name = row.names(tp), score = tp[, 6] / max(tp[, 6]) * 1000)
write.table(bed_data, file = "merged.chr9_interactions.bed", quote = F, row.names = F, col.names = F)
#######################################################
# Method_2
#######################################################
expFile <- "4C12h.bam"
contrFile <- "4C0h.bam"
my_obj <- new("r3Cseq",
organismName = "hg19", alignedReadsBamExpFile = expFile,
alignedReadsBamContrFile = contrFile, isControlInvolved = TRUE, viewpoint_chromosome = "chr9",
viewpoint_primer_forward = "tgttgcacaaacaaacccaga",
viewpoint_primer_reverse = "TATGGCTGGCTTTGGGGAAG",
expLabel = "GM12h", contrLabel = "GM0h", restrictionEnzyme = "EcoRI"
)
getRawReads(my_obj)
getReadCountPerRestrictionFragment(my_obj)
calculateRPM(my_obj)
getInteractions(my_obj)
expResults <- expInteractionRegions(my_obj)
contrResults <- contrInteractionRegions(my_obj)
export3CseqRawReads2bedGraph(my_obj)
export3Cseq2bedGraph(my_obj)
plotOverviewInteractions(my_obj)
plotInteractionsNearViewpoint(my_obj, distance = 2e5)