-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathabsenrich.R
40 lines (35 loc) · 1.56 KB
/
absenrich.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
batchMode = TRUE
if (batchMode) {
cmdArgs = commandArgs(TRUE)
print(cmdArgs)
}
GENOME.STATS <- "/gpfs/home/wallen/data/genome_data/mm9_window_1kb_stats.txt"
compareCpGAndCpAMeth <- function() {
# setwd("/gpfs/home/wallen/data/genome_data")
# window.1kb.stats <- read.table("mm9_window_1kb_stats.txt",
# header=TRUE, sep="\t")
setwd("~/experiment/experiment/stavros_data")
load("/gpfs/home/wallen/data/genome_data/mm9_window_1kb_stats.Rdata")
omp.medip <- load.DipData("omp_medip_1kb")
omp.hmedip <- load.DipData("omp_hmedip_1kb")
cat("Sampling...\n")
idx <- sample(1:length(omp.medip$genome.data[,'raw']), 1000000)
png("cpn_corr.png")
par(mfrow=c(2,2))
cat("Drawing plot...\n")
plot(mm9.window.1kb.stats$CpA[idx],omp.medip$genome.data[idx, 'raw'],
pch=20, cex=.1, col=rgb(0, 0, 100, 50, maxColorValue=255),
ylim=c(0, 100), main="CpA meC")
plot(mm9.window.1kb.stats$CpG[idx],omp.medip$genome.data[idx, 'raw'],
pch=20, cex=.1, col=rgb(0, 0, 100, 50, maxColorValue=255),
ylim=c(0, 100), main="CpG meC")
plot(mm9.window.1kb.stats$CpA[idx],omp.hmedip$genome.data[idx, 'raw'],
pch=20, cex=.1, col=rgb(0,0,100, 50, maxColorValue=255),
ylim=c(0, 100), main="CpA hmeC")
plot(mm9.window.1kb.stats$CpG[idx],omp.hmedip$genome.data[idx, 'raw'],
pch=20, cex=.1, col=rgb(0,0,100,50, maxColorValue=255),
ylim=c(0, 100), main="CpG hmeC")
dev.off()
}
load("/gpfs/home/wallen/data/genome_data/mm9_window_1kb_stats.Rdata")
omp.medip <- load.DipData("omp_medip_1kb")$genome.data[,'raw']