Homec4science
Diffusion PhD_thesis (master) Browse

R9524/scripts/ch_encode_peaks/figure_motif_association.R
da4a342448bamaster

Last Change
Raw File

figure_motif_association.R
View Options

# libraries
library(gplots)
library(RColorBrewer)
setwd("/local/groux/encode_uniform_peaks/")
get.sample.labels = function(DATA_DF, SAMPLE_DF)
{
LABELS = vector(length=nrow(DATA_DF), mode="character")
for(i in 1:nrow(DATA_DF))
{ index = grep(DATA_DF[i,1], SAMPLE_DF[,2])
LABELS[i] = SAMPLE_DF[index,7]
}
return(LABELS)
}
density.mean = function(P.X, X)
{ # Given a vector of values and their associates probabilities,
# this function returns the weighted mean.
#
# P.X a vector of numerical containing the prob. of X
# X a vector of numerical values
#
# Returns a numerical, the weighted mean of X
# if the density does not sum up to 1
P.X = P.X / sum(P.X)
m = P.X %*% X
return(m)
}
density.var = function(P.X, X)
{ # Given a vector of values and their associates probabilities,
# this function returns the weighted mean.
#
# P.X a vector of numerical containing the prob. of X
# X a vector of numerical values
#
# Returns a numerical, the weighted variance of X
# if the density does not sum up to 1
# if the density does not sum up to 1
P.X = P.X / sum(P.X)
m = density.mean(X, P.X)
# v = ((X - m)**2) %*% Y
v = ((X**2) %*% P.X) - (m**2)
return(v)
}
process_data = function(file_data_all, file_data_with, file_data_wo, file_table, file_sample, key)
{
# density reading and row labelling
data_all = read.table(file_data_all, header=F, stringsAsFactors=F)
data_with = read.table(file_data_with, header=F, stringsAsFactors=F)
data_wo = read.table(file_data_wo, header=F, stringsAsFactors=F)
sample_df = read.table(file_sample, header=F, stringsAsFactors=F)
labels = get.sample.labels(data_all, sample_df)
rownames(data_all) = labels
rownames(data_with) = labels
rownames(data_wo) = labels
data_all = as.matrix(data_all[,-1])
data_with = as.matrix(data_with[,-1])
data_wo = as.matrix(data_wo[,-1])
n.row = nrow(data_all)
n.col = ncol(data_all)
# association table reading and row labelling
table = read.table(file_table, header=T, stringsAsFactors=F)
labels = get.sample.labels(table, sample_df)
table$Label = labels
labels = sample_df = NULL
# select only the rows of table which have "MA0139.1" as value in MotifB
index = which(table$MotifB == key)
table = table[index,]
index = NULL
# force data and table rows to have the same order
# normally, data rowname values should also be table$Peaklist values
ord1 = order(rownames(data_all))
ord2 = order(table$Label)
data_all = data_all[ord1,]
data_with = data_with[ord1,]
data_wo = data_wo[ord1,]
table = table[ord2,]
# compute mean distance to peak and its variance
x_lab = c(0:(ncol(data_all[,1:501])-1))
mean_all = apply(data_all[,1:501], 1, density.mean, x_lab)
mean_with = apply(data_with[,1:501], 1, density.mean, x_lab)
mean_wo = apply(data_wo[,1:501], 1, density.mean, x_lab)
sd_all = sqrt(apply(data_all[,1:501], 1, density.var, x_lab))
sd_with = sqrt(apply(data_with[,1:501], 1, density.var, x_lab))
sd_wo = sqrt(apply(data_wo[,1:501], 1, density.var, x_lab))
# determine significant association/repulsion
index.ass.pos = which((table$Odd.ratio > 1) & (table$Odd.ratio.CI.lower > 1))
index.ass.neg = which((table$Odd.ratio < 1) & (table$Odd.ratio.CI.upper < 1))
# set colors for odd-ratio (<1, >1 and ==1)
color.neg = brewer.pal(11, "RdYlGn")[2] # red
color.pos = brewer.pal(11, "RdYlGn")[10] # green
color.neutral = "black"
colors = rep(color.neutral, nrow(data_all))
colors[index.ass.pos] = color.pos
colors[index.ass.neg] = color.neg
# point radius
# proportional to the propotion of peak having only the other motif
point.cex = table$With.motifB / (table$With.motifA+table$With.motifB+table$With.both.motif+table$With.no.motif)
constant = 2/max(point.cex)
point.cex = constant * point.cex
# number each point for legend purposes
x.lim = c(0.5*min(table$Odd.ratio), max(table$Odd.ratio[which(!is.infinite(table$Odd.ratio))]))
y.lim = c(0, max(mean_all))
numb = 1:nrow(table)
numb.dy = 0 # (y.lim[2] - y.lim[1]) /30
numb.dx = (x.lim[2] - x.lim[1]) /30
return(list(table=table,
data_all=data_all,
data_with=data_with,
data_wo=data_wo,
mean_all=mean_all,
mean_with=mean_with,
mean_wo=mean_wo,
sd_all=sd_all,
sd_with=sd_with,
sd_wo=sd_wo,
x.lim=x.lim,
y.lim=y.lim,
colors=colors,
point.cex=point.cex,
numb=numb))
}
# list of TF with their respective PWM ID which could be given as key
# 1st column is for labelling and naming,
# 2nd column is an ID that should be found in DATA$MotifB using the given key
assoc_table = matrix(data=c("CTCF", "MA0139.1",
"Jund", "JUND_HUMAN.H10MO.A"),
ncol=2, byrow=T)
file_table = "results/analysis_6_motif_colocalization/with_rmsk_motif/GM12878_association_scores_500bp.tsv"
file_sample = "ENCODE_Uniform_GM12878_data.tsv"
ctcf_data = process_data("results/analysis_6_motif_colocalization/with_rmsk_motif/motif_density/GM12878_CTCF_motif_density_allpeak.tsv",
"results/analysis_6_motif_colocalization/with_rmsk_motif/motif_density/GM12878_CTCF_motif_density_peakwithmotif.tsv",
"results/analysis_6_motif_colocalization/with_rmsk_motif/motif_density/GM12878_CTCF_motif_density_peaknomotif.tsv",
file_table,
file_sample,
"MA0139.1")
jund_data = process_data("results/analysis_6_motif_colocalization/with_rmsk_motif/motif_density/GM12878_JunD_motif_density_allpeak.tsv",
"results/analysis_6_motif_colocalization/with_rmsk_motif/motif_density/GM12878_JunD_motif_density_peakwithmotif.tsv",
"results/analysis_6_motif_colocalization/with_rmsk_motif/motif_density/GM12878_JunD_motif_density_peaknomotif.tsv",
file_table,
file_sample,
"JUND_HUMAN.H10MO.A")
#####################################################################
# Figure CTCF motif odd radio vs mean distance (peak to CTCF motif) #
#####################################################################
# X11(height=5, width=18)
png("/local/groux/phd_thesis/images/ch_encode_peaks/ctcf_motif_odd_ratio_vs_meandist.png",
height=5, width=18, res=720, units="in")
l = layout(matrix(nrow=1, ncol=4, data=c(1,2,3,4)), widths=c(3,3,3,1.5))
layout.show(l)
p = par(mar=c(5.1, 4.1, 4.1, 2.1))
# all
plot(x=ctcf_data$table$Odd.ratio, y=ctcf_data$mean_all, col="white", xlim=ctcf_data$x.lim, ylim=ctcf_data$y.lim,
main="All peaks", xlab="Odd ratio", ylab="mean dist. to motif",
cex.main=2, cex.axis=1.5, cex.lab=1.5)
# add Odd.ratio 95% CI
segments(x0=ctcf_data$table$Odd.ratio, y0=ctcf_data$mean_all, x1=ctcf_data$table$Odd.ratio.CI.lower, y1=ctcf_data$mean_all, col=ctcf_data$colors, lty=3, lwd=0.5)
segments(x0=ctcf_data$table$Odd.ratio, y0=ctcf_data$mean_all, x1=ctcf_data$table$Odd.ratio.CI.upper, y1=ctcf_data$mean_all, col=ctcf_data$colors, lty=3, lwd=0.5)
# add standard deviation around mean
segments(x0=ctcf_data$table$Odd.ratio, y0=ctcf_data$mean_all, x1=ctcf_data$table$Odd.ratio, y1=ctcf_data$mean_all - ctcf_data$sd_all, col=ctcf_data$colors, lty=3, lwd=0.5)
segments(x0=ctcf_data$table$Odd.ratio, y0=ctcf_data$mean_all, x1=ctcf_data$table$Odd.ratio, y1=ctcf_data$mean_all + ctcf_data$sd_all, col=ctcf_data$colors, lty=3, lwd=0.5)
# plot the points on top of the segments
points(x=ctcf_data$table$Odd.ratio, y=ctcf_data$mean_all, col=ctcf_data$colors, cex=ctcf_data$point.cex, pch=19)
# plot a number above each point for legend
text(x=ctcf_data$table$Odd.ratio-0.05, y=ctcf_data$mean_all, label=ctcf_data$numb, col=ctcf_data$colors, cex=1)
p = par(xpd=NA)
text(0.5*min(ctcf_data$x.lim), 1.12*max(ctcf_data$y.lim), "A", cex=6)
par(p)
# with motif
plot(x=ctcf_data$table$Odd.ratio, y=ctcf_data$mean_with, col="white", xlim=ctcf_data$x.lim, ylim=ctcf_data$y.lim,
main="Peaks with motif", xlab="Odd ratio", ylab="mean dist. to motif",
cex.main=2, cex.axis=1.5, cex.lab=1.5)
# add Odd.ratio 95% CI
segments(x0=ctcf_data$table$Odd.ratio, y0=ctcf_data$mean_with, x1=ctcf_data$table$Odd.ratio.CI.lower, y1=ctcf_data$mean_with, col=ctcf_data$colors, lty=3, lwd=0.5)
segments(x0=ctcf_data$table$Odd.ratio, y0=ctcf_data$mean_with, x1=ctcf_data$table$Odd.ratio.CI.upper, y1=ctcf_data$mean_with, col=ctcf_data$colors, lty=3, lwd=0.5)
# add standard deviation around mean
segments(x0=ctcf_data$table$Odd.ratio, y0=ctcf_data$mean_with, x1=ctcf_data$table$Odd.ratio, y1=ctcf_data$mean_with - ctcf_data$sd_with, col=ctcf_data$colors, lty=3, lwd=0.5)
segments(x0=ctcf_data$table$Odd.ratio, y0=ctcf_data$mean_with, x1=ctcf_data$table$Odd.ratio, y1=ctcf_data$mean_with + ctcf_data$sd_with, col=ctcf_data$colors, lty=3, lwd=0.5)
# plot the points on top of the segments
points(x=ctcf_data$table$Odd.ratio, y=ctcf_data$mean_with, col=ctcf_data$colors, cex=ctcf_data$point.cex, pch=19)
# plot a number above each point for legend
text(x=ctcf_data$table$Odd.ratio-0.05, y=ctcf_data$mean_with, label=ctcf_data$numb, col=ctcf_data$colors, cex=1)
p = par(xpd=NA)
text(0.5*min(ctcf_data$x.lim), 1.12*max(ctcf_data$y.lim), "B", cex=6)
par(p)
# without motif
plot(x=ctcf_data$table$Odd.ratio, y=ctcf_data$mean_wo, col="white", xlim=ctcf_data$x.lim, ylim=ctcf_data$y.lim,
main="Peaks without motif", xlab="Odd ratio", ylab="mean dist. to motif",
cex.main=2, cex.axis=1.5, cex.lab=1.5)
# add Odd.ratio 95% CI
segments(x0=ctcf_data$table$Odd.ratio, y0=ctcf_data$mean_wo, x1=ctcf_data$table$Odd.ratio.CI.lower, y1=ctcf_data$mean_wo, col=ctcf_data$colors, lty=3, lwd=0.5)
segments(x0=ctcf_data$table$Odd.ratio, y0=ctcf_data$mean_wo, x1=ctcf_data$table$Odd.ratio.CI.upper, y1=ctcf_data$mean_wo, col=ctcf_data$colors, lty=3, lwd=0.5)
# add standard deviation around mean
segments(x0=ctcf_data$table$Odd.ratio, y0=ctcf_data$mean_wo, x1=ctcf_data$table$Odd.ratio, y1=ctcf_data$mean_wo - ctcf_data$sd_wo, col=ctcf_data$colors, lty=3, lwd=0.5)
segments(x0=ctcf_data$table$Odd.ratio, y0=ctcf_data$mean_wo, x1=ctcf_data$table$Odd.ratio, y1=ctcf_data$mean_wo + ctcf_data$sd_wo, col=ctcf_data$colors, lty=3, lwd=0.5)
# plot the points on top of the segments
points(x=ctcf_data$table$Odd.ratio, y=ctcf_data$mean_wo, col=ctcf_data$colors, cex=ctcf_data$point.cex, pch=19)
# plot a number above each point for legend
text(x=ctcf_data$table$Odd.ratio-0.05, y=ctcf_data$mean_wo, label=ctcf_data$numb, col=ctcf_data$colors, cex=1)
p = par(xpd=NA)
text(0.5*min(ctcf_data$x.lim), 1.12*max(ctcf_data$y.lim), "C", cex=6)
par(p)
p = par(mar=c(5.1,0,4.1,0))
plot(0, xlim=c(0,1), ylim=ctcf_data$y.lim, col="white", bty='n', main="", xlab="", ylab="", xaxt='n', yaxt='n')
# legend in the margin, each point number is written with the dataset name on its right
x.text = c(rep(0.0, floor(nrow(ctcf_data$table)/2)), rep(0.5, ceiling(nrow(ctcf_data$table)/2)))
y.text = c(seq(max(ctcf_data$y.lim), max(ctcf_data$y.lim)/20, length.out=floor(nrow(ctcf_data$table)/2)),
seq(max(ctcf_data$y.lim), max(ctcf_data$y.lim)/20, length.out=ceiling(nrow(ctcf_data$table)/2)))
text(x=x.text, y=y.text, labels=paste(ctcf_data$numb, ctcf_data$table$Label, sep=" "),
cex=1.45, pos=4, xpd=NA, col=ctcf_data$colors)
par(p)
dev.off()
#####################################################################
# Figure jund motif odd radio vs mean distance (peak to jund motif) #
#####################################################################
jund_data$x.lim[2] = ctcf_data$x.lim[2]
# X11(height=5, width=18)
png("/local/groux/phd_thesis/images/ch_encode_peaks/jund_motif_odd_ratio_vs_meandist.png",
height=5, width=18, res=720, units="in")
l = layout(matrix(nrow=1, ncol=4, data=c(1,2,3,4)), widths=c(3,3,3,1.5))
layout.show(l)
p = par(mar=c(5.1, 4.1, 4.1, 2.1))
# all
plot(x=jund_data$table$Odd.ratio, y=jund_data$mean_all, col="white", xlim=jund_data$x.lim, ylim=jund_data$y.lim,
main="All peaks", xlab="Odd ratio", ylab="mean dist. to motif",
cex.main=2, cex.axis=1.5, cex.lab=1.5)
# add Odd.ratio 95% CI
segments(x0=jund_data$table$Odd.ratio, y0=jund_data$mean_all, x1=jund_data$table$Odd.ratio.CI.lower, y1=jund_data$mean_all, col=jund_data$colors, lty=3, lwd=0.5)
segments(x0=jund_data$table$Odd.ratio, y0=jund_data$mean_all, x1=jund_data$table$Odd.ratio.CI.upper, y1=jund_data$mean_all, col=jund_data$colors, lty=3, lwd=0.5)
# add standard deviation around mean
segments(x0=jund_data$table$Odd.ratio, y0=jund_data$mean_all, x1=jund_data$table$Odd.ratio, y1=jund_data$mean_all - jund_data$sd_all, col=jund_data$colors, lty=3, lwd=0.5)
segments(x0=jund_data$table$Odd.ratio, y0=jund_data$mean_all, x1=jund_data$table$Odd.ratio, y1=jund_data$mean_all + jund_data$sd_all, col=jund_data$colors, lty=3, lwd=0.5)
# plot the points on top of the segments
points(x=jund_data$table$Odd.ratio, y=jund_data$mean_all, col=jund_data$colors, cex=jund_data$point.cex, pch=19)
# plot a number above each point for legend
text(x=jund_data$table$Odd.ratio-0.05, y=jund_data$mean_all, label=jund_data$numb, col=jund_data$colors, cex=1)
p = par(xpd=NA)
text(0.5*min(jund_data$x.lim), 1.12*max(jund_data$y.lim), "A", cex=6)
par(p)
# with motif
plot(x=jund_data$table$Odd.ratio, y=jund_data$mean_with, col="white", xlim=jund_data$x.lim, ylim=jund_data$y.lim,
main="Peaks with motif", xlab="Odd ratio", ylab="mean dist. to motif",
cex.main=2, cex.axis=1.5, cex.lab=1.5)
# add Odd.ratio 95% CI
segments(x0=jund_data$table$Odd.ratio, y0=jund_data$mean_with, x1=jund_data$table$Odd.ratio.CI.lower, y1=jund_data$mean_with, col=jund_data$colors, lty=3, lwd=0.5)
segments(x0=jund_data$table$Odd.ratio, y0=jund_data$mean_with, x1=jund_data$table$Odd.ratio.CI.upper, y1=jund_data$mean_with, col=jund_data$colors, lty=3, lwd=0.5)
# add standard deviation around mean
segments(x0=jund_data$table$Odd.ratio, y0=jund_data$mean_with, x1=jund_data$table$Odd.ratio, y1=jund_data$mean_with - jund_data$sd_with, col=jund_data$colors, lty=3, lwd=0.5)
segments(x0=jund_data$table$Odd.ratio, y0=jund_data$mean_with, x1=jund_data$table$Odd.ratio, y1=jund_data$mean_with + jund_data$sd_with, col=jund_data$colors, lty=3, lwd=0.5)
# plot the points on top of the segments
points(x=jund_data$table$Odd.ratio, y=jund_data$mean_with, col=jund_data$colors, cex=jund_data$point.cex, pch=19)
# plot a number above each point for legend
text(x=jund_data$table$Odd.ratio-0.05, y=jund_data$mean_with, label=jund_data$numb, col=jund_data$colors, cex=1)
p = par(xpd=NA)
text(0.5*min(jund_data$x.lim), 1.12*max(jund_data$y.lim), "B", cex=6)
par(p)
# without motif
plot(x=jund_data$table$Odd.ratio, y=jund_data$mean_wo, col="white", xlim=jund_data$x.lim, ylim=jund_data$y.lim,
main="Peaks without motif", xlab="Odd ratio", ylab="mean dist. to motif",
cex.main=2, cex.axis=1.5, cex.lab=1.5)
# add Odd.ratio 95% CI
segments(x0=jund_data$table$Odd.ratio, y0=jund_data$mean_wo, x1=jund_data$table$Odd.ratio.CI.lower, y1=jund_data$mean_wo, col=jund_data$colors, lty=3, lwd=0.5)
segments(x0=jund_data$table$Odd.ratio, y0=jund_data$mean_wo, x1=jund_data$table$Odd.ratio.CI.upper, y1=jund_data$mean_wo, col=jund_data$colors, lty=3, lwd=0.5)
# add standard deviation around mean
segments(x0=jund_data$table$Odd.ratio, y0=jund_data$mean_wo, x1=jund_data$table$Odd.ratio, y1=jund_data$mean_wo - jund_data$sd_wo, col=jund_data$colors, lty=3, lwd=0.5)
segments(x0=jund_data$table$Odd.ratio, y0=jund_data$mean_wo, x1=jund_data$table$Odd.ratio, y1=jund_data$mean_wo + jund_data$sd_wo, col=jund_data$colors, lty=3, lwd=0.5)
# plot the points on top of the segments
points(x=jund_data$table$Odd.ratio, y=jund_data$mean_wo, col=jund_data$colors, cex=jund_data$point.cex, pch=19)
# plot a number above each point for legend
text(x=jund_data$table$Odd.ratio-0.05, y=jund_data$mean_wo, label=jund_data$numb, col=jund_data$colors, cex=1)
p = par(xpd=NA)
text(0.5*min(jund_data$x.lim), 1.12*max(jund_data$y.lim), "C", cex=6)
par(p)
p = par(mar=c(5.1,0,4.1,0))
plot(0, xlim=c(0,1), ylim=jund_data$y.lim, col="white", bty='n', main="", xlab="", ylab="", xaxt='n', yaxt='n')
# legend in the margin, each point number is written with the dataset name on its right
x.text = c(rep(0.0, floor(nrow(jund_data$table)/2)), rep(0.5, ceiling(nrow(jund_data$table)/2)))
y.text = c(seq(max(jund_data$y.lim), max(jund_data$y.lim)/20, length.out=floor(nrow(jund_data$table)/2)),
seq(max(jund_data$y.lim), max(jund_data$y.lim)/20, length.out=ceiling(nrow(jund_data$table)/2)))
text(x=x.text, y=y.text, labels=paste(jund_data$numb, jund_data$table$Label, sep=" "),
cex=1.45, pos=4, xpd=NA, col=jund_data$colors)
par(p)
dev.off()
######################
# CTCF motif density #
######################
# peaks with motif
# X11(height=7, width=6)
png("/local/groux/phd_thesis/images/ch_encode_peaks/ctcf_withmotif_density.png",
height=7, width=6, res=720, units="in")
colors = colorRampPalette(c("blue", "white", "red"))(n = 100)
heatmap.2(ctcf_data$data_with[,1:200],
colRow=ctcf_data$colors,
Rowv=T,
Colv=F,
dendrogram="row",
scale="row",
trace="none",
margins=c(5,10),
col=colors,
cexCol=0.3,
labRow=ctcf_data$table$Label,
main="Peak with motif",
xlab="position to peak [bp]")
dev.off()
# peaks without motif
# X11(height=7, width=6)
png("/local/groux/phd_thesis/images/ch_encode_peaks/ctcf_nomotif_density.png",
height=7, width=6, res=720, units="in")
colors = colorRampPalette(c("blue", "white", "red"))(n = 100)
heatmap.2(ctcf_data$data_wo[,1:200],
colRow=ctcf_data$colors,
Rowv=T,
Colv=F,
dendrogram="row",
scale="row",
trace="none",
margins=c(5,10),
col=colors,
cexCol=0.3,
labRow=ctcf_data$table$Label,
main="Peak with motif",
xlab="position to peak [bp]")
dev.off()
######################
# jund motif density #
######################
# peaks with motif
# X11(height=7, width=6)
png("/local/groux/phd_thesis/images/ch_encode_peaks/jund_withmotif_density.png",
height=7, width=6, res=720, units="in")
colors = colorRampPalette(c("blue", "white", "red"))(n = 100)
heatmap.2(ctcf_data$data_with[,1:200],
colRow=ctcf_data$colors,
Rowv=T,
Colv=F,
dendrogram="row",
scale="row",
trace="none",
margins=c(5,10),
col=colors,
cexCol=0.3,
labRow=ctcf_data$table$Label,
main="Peak with motif",
xlab="position to peak [bp]")
dev.off()
# peaks without motif
# X11(height=7, width=6)
png("/local/groux/phd_thesis/images/ch_encode_peaks/jund_nomotif_density.png",
height=7, width=6, res=720, units="in")
colors = colorRampPalette(c("blue", "white", "red"))(n = 100)
heatmap.2(ctcf_data$data_wo[,1:200],
colRow=ctcf_data$colors,
Rowv=T,
Colv=F,
dendrogram="row",
scale="row",
trace="none",
margins=c(5,10),
col=colors,
cexCol=0.3,
labRow=ctcf_data$table$Label,
main="Peak with motif",
xlab="position to peak [bp]")
dev.off()
#########################
# CTCF motif odd-ratios #
#########################
# X11(height=4, width=16)
png("/local/groux/phd_thesis/images/ch_encode_peaks/ctcf_odd_ratio.png",
height=4, width=16, res=720, units="in")
ord = order(ctcf_data$table$Odd.ratio, decreasing=T)
y = ctcf_data$table$Odd.ratio[ord]
y0 = ctcf_data$table$Odd.ratio.CI.lower[ord]
y1 = ctcf_data$table$Odd.ratio.CI.upper[ord]
x = 1:length(y)
sd = ctcf_data$sd_all[ord]
p = par(cex.axis=0.5, mar=c(10.1, 4.1, 1.1, 4.1), cex.axis=1.1)
plot(x, y, xaxt='n', col=ctcf_data$colors[ord],
ylim=c(0,3),
main="CTCF motif odd-ratios", xlab="", ylab="Odd ratio")
# confidence interval
segments(x0=x, y0=y0, x1=x, y1=y1, col=ctcf_data$colors[ord], lwd=1.5)
# make X-axis
index.ass.pos = ((ctcf_data$table$Odd.ratio > 1) & (ctcf_data$table$Odd.ratio.CI.lower > 1))[ord]
index.ass.neg = ((ctcf_data$table$Odd.ratio < 1) & (ctcf_data$table$Odd.ratio.CI.upper < 1))[ord]
index.ass.no = !(index.ass.pos | index.ass.neg)
axis(side=1, at=x[index.ass.no], labels=ctcf_data$table$Label[ord][index.ass.no], las=2, hadj=1, col.axis=ctcf_data$colors[ord][index.ass.no][1])
axis(side=1, at=x[index.ass.pos], labels=ctcf_data$table$Label[ord][index.ass.pos], las=2, hadj=1, col.axis=ctcf_data$colors[ord][index.ass.pos][1])
axis(side=1, at=x[index.ass.neg], labels=ctcf_data$table$Label[ord][index.ass.neg], las=2, hadj=1, col.axis=ctcf_data$colors[ord][index.ass.neg][1])
# label 1
abline(h=1, lty=2)
dev.off()
#########################
# JunD motif odd-ratios #
#########################
# X11(height=4, width=16)
png("/local/groux/phd_thesis/images/ch_encode_peaks/jund_odd_ratio.png",
height=4, width=16, res=720, units="in")
ord = order(jund_data$table$Odd.ratio, decreasing=T)
y = jund_data$table$Odd.ratio[ord]
y0 = jund_data$table$Odd.ratio.CI.lower[ord]
y1 = jund_data$table$Odd.ratio.CI.upper[ord]
x = 1:length(y)
sd = jund_data$sd_all[ord]
p = par(cex.axis=0.5, mar=c(10.1, 4.1, 1.1, 4.1), cex.axis=1.1)
plot(x, y, xaxt='n', col=jund_data$colors[ord],
ylim=c(0,3),
main="JunD motif odd-ratios", xlab="", ylab="Odd ratio")
# confidence interval
segments(x0=x, y0=y0, x1=x, y1=y1, col=jund_data$colors[ord], lwd=1.5)
# make X-axis
index.ass.pos = ((jund_data$table$Odd.ratio > 1) & (jund_data$table$Odd.ratio.CI.lower > 1))[ord]
index.ass.neg = ((jund_data$table$Odd.ratio < 1) & (jund_data$table$Odd.ratio.CI.upper < 1))[ord]
index.ass.no = !(index.ass.pos | index.ass.neg)
axis(side=1, at=x[index.ass.no], labels=jund_data$table$Label[ord][index.ass.no], las=2, hadj=1, col.axis=jund_data$colors[ord][index.ass.no][1])
axis(side=1, at=x[index.ass.pos], labels=jund_data$table$Label[ord][index.ass.pos], las=2, hadj=1, col.axis=jund_data$colors[ord][index.ass.pos][1])
axis(side=1, at=x[index.ass.neg], labels=jund_data$table$Label[ord][index.ass.neg], las=2, hadj=1, col.axis=jund_data$colors[ord][index.ass.neg][1])
# label 1
abline(h=1, lty=2)
dev.off()
c4science · Help