Commit aa28a9d3 authored by Venkat Malladi's avatar Venkat Malladi
Browse files

Chip-seq Figure 4B and C PJ34.

parent a705af10
......@@ -314,7 +314,7 @@ axis(2, at=seq(0,18,2))
dev.off()
# Wilcox Rank Test
# Significance 1.324e-12
# Significance 2.915e-08
wilcox.test(er_alpha[,1], er_alpha[,2],paired=T)
wilcox.test(er_alpha[,1], er_alpha[,3],paired=T)
wilcox.test(er_alpha[,1], er_alpha[,4],paired=T)
......@@ -339,7 +339,7 @@ axis(2, at=seq(0,12,2))
dev.off()
# Wilcox Rank Test
# Significance 1.324e-12
# Significance 2.915e-08
wilcox.test(foxa1[,1], foxa1[,2],paired=T)
wilcox.test(foxa1[,1], foxa1[,3],paired=T)
wilcox.test(foxa1[,1], foxa1[,4],paired=T)
......@@ -364,7 +364,7 @@ axis(2, at=seq(0,12,2))
dev.off()
# Wilcox Rank Test
# Significance 1.324e-12
# Significance 2.915e-08
wilcox.test(h3k27ac[,1], h3k27ac[,2],paired=T)
wilcox.test(h3k27ac[,1], h3k27ac[,3],paired=T)
wilcox.test(h3k27ac[,1], h3k27ac[,4],paired=T)
......@@ -372,7 +372,7 @@ wilcox.test(h3k27ac[,2], h3k27ac[,3],paired=T)
wilcox.test(h3k27ac[,2], h3k27ac[,4],paired=T)
wilcox.test(h3k27ac[,3], h3k27ac[,4],paired=T)
# PAPR-1
# PARP-1
rpkm_l_er_parp <- countOverlaps(closest_er_alpha_luc,l_er_parp)/ (width(closest_er_alpha_luc)/1000 * library_l_er_parp/1000000)
rpkm_p_er_parp <- countOverlaps(closest_er_alpha_luc,p_er_parp)/ (width(closest_er_alpha_luc)/1000 * library_p_er_parp/1000000)
rpkm_l_v_parp <- countOverlaps(closest_er_alpha_luc,l_v_parp)/ (width(closest_er_alpha_luc)/1000 * library_l_v_parp/1000000)
......@@ -389,7 +389,7 @@ axis(2, at=seq(0,1.2,.2))
dev.off()
# Wilcox Rank Test
# Significance 1.324e-12
# Significance 2.915e-08
wilcox.test(parp[,1], parp[,2],paired=T)
wilcox.test(parp[,1], parp[,3],paired=T)
wilcox.test(parp[,1], parp[,4],paired=T)
......@@ -398,21 +398,21 @@ wilcox.test(parp[,2], parp[,4],paired=T)
wilcox.test(parp[,3], parp[,4],paired=T)
# NELF-E
rpkm_l_er_parp <- countOverlaps(closest_er_alpha_luc,l_er_parp)/ (width(closest_er_alpha_luc)/1000 * library_l_er_parp/1000000)
rpkm_p_er_parp <- countOverlaps(closest_er_alpha_luc,p_er_parp)/ (width(closest_er_alpha_luc)/1000 * library_p_er_parp/1000000)
rpkm_l_v_parp <- countOverlaps(closest_er_alpha_luc,l_v_parp)/ (width(closest_er_alpha_luc)/1000 * library_l_v_parp/1000000)
rpkm_p_v_parp <- countOverlaps(closest_er_alpha_luc,p_v_parp)/ (width(closest_er_alpha_luc)/1000 * library_p_v_parp/1000000)
parp <- data.frame(rpkm_l_v_parp,rpkm_l_er_parp,rpkm_p_v_parp,rpkm_p_er_parp)
parp_e <- data.frame(rpkm_l_er_parp,rpkm_p_er_parp)
jpeg('boxplot_parp_1_rpkm.jpg')
boxplot(parp, col=(c("forestgreen","forestgreen","forestgreen","forestgreen")), main="PARP-1", ylab="RPKM",outline=FALSE, notch=TRUE, names=c("L", "L+E", "P", "P+E"), ylim = c(0,1.2),yaxt="n", cex.axis=1,las=2,lwd=4,lty=1)
axis(2, at=seq(0,1.2,.2))
dev.off()
jpeg('boxplot_parp_1_e2_rpkm.jpg')
boxplot(parp_e, col=(c("forestgreen","forestgreen")), main="PARP-1", ylab="RPKM",outline=FALSE, notch=TRUE, names=c("L+E", "P+E"), ylim = c(0,1.2),yaxt="n", cex.axis=1,las=2,lwd=4,lty=1)
axis(2, at=seq(0,1.2,.2))
dev.off()
# NELF-E # To Do
#rpkm_l_er_parp <- countOverlaps(closest_er_alpha_luc,l_er_parp)/ (width(closest_er_alpha_luc)/1000 * library_l_er_parp/1000000)
#rpkm_p_er_parp <- countOverlaps(closest_er_alpha_luc,p_er_parp)/ (width(closest_er_alpha_luc)/1000 * library_p_er_parp/1000000)
#rpkm_l_v_parp <- countOverlaps(closest_er_alpha_luc,l_v_parp)/ (width(closest_er_alpha_luc)/1000 * library_l_v_parp/1000000)
#rpkm_p_v_parp <- countOverlaps(closest_er_alpha_luc,p_v_parp)/ (width(closest_er_alpha_luc)/1000 * library_p_v_parp/1000000)
#parp <- data.frame(rpkm_l_v_parp,rpkm_l_er_parp,rpkm_p_v_parp,rpkm_p_er_parp)
#parp_e <- data.frame(rpkm_l_er_parp,rpkm_p_er_parp)
#jpeg('boxplot_parp_1_rpkm.jpg')
#boxplot(parp, col=(c("forestgreen","forestgreen","forestgreen","forestgreen")), main="PARP-1", ylab="RPKM",outline=FALSE, notch=TRUE, names=c("L", "L+E", "P", "P+E"), ylim = c(0,1.2),yaxt="n", cex.axis=1,las=2,lwd=4,lty=1)
#axis(2, at=seq(0,1.2,.2))
#dev.off()
#jpeg('boxplot_parp_1_e2_rpkm.jpg')
#boxplot(parp_e, col=(c("forestgreen","forestgreen")), main="PARP-1", ylab="RPKM",outline=FALSE, notch=TRUE, names=c("L+E", "P+E"), ylim = c(0,1.2),yaxt="n", cex.axis=1,las=2,lwd=4,lty=1)
#axis(2, at=seq(0,1.2,.2))
#dev.off()
# Wilcox Rank Test
# Significance 1.324e-12
......@@ -535,7 +535,7 @@ axis(2, at=seq(0,12,2))
dev.off()
# Wilcox Rank Test
# Significance 1.324e-12
# Significance 2.915e-08
wilcox.test(ffoxa1[,1], ffoxa1[,2],paired=T)
wilcox.test(ffoxa1[,1], ffoxa1[,3],paired=T)
wilcox.test(ffoxa1[,1], ffoxa1[,4],paired=T)
......@@ -560,7 +560,7 @@ axis(2, at=seq(0,1.4,.2))
dev.off()
# Wilcox Rank Test
# Significance 1.324e-12
# Significance 2.915e-08
wilcox.test(fer_alpha[,1], fer_alpha[,2],paired=T)
wilcox.test(fer_alpha[,1], fer_alpha[,3],paired=T)
wilcox.test(fer_alpha[,1], fer_alpha[,4],paired=T)
......@@ -587,7 +587,7 @@ axis(2, at=seq(0,6,2))
dev.off()
# Wilcox Rank Test
# Significance 1.324e-12
# Significance 2.915e-08
wilcox.test(fh3k27ac[,1], fh3k27ac[,2],paired=T)
wilcox.test(fh3k27ac[,1], fh3k27ac[,3],paired=T)
wilcox.test(fh3k27ac[,1], fh3k27ac[,4],paired=T)
......@@ -612,7 +612,7 @@ axis(2, at=seq(0,1.4,.2))
dev.off()
# Wilcox Rank Test
# Significance 1.324e-12
# Significance 2.915e-08
wilcox.test(fparp[,1], fparp[,2],paired=T)
wilcox.test(fparp[,1], fparp[,3],paired=T)
wilcox.test(fparp[,1], fparp[,4],paired=T)
......@@ -620,8 +620,203 @@ wilcox.test(fparp[,2], fparp[,3],paired=T)
wilcox.test(fparp[,2], fparp[,4],paired=T)
wilcox.test(fparp[,3], fparp[,4],paired=T)
#Figure 4B and C
## Read alignments
e_er_R1 <- as(readGAlignments("/Volumes/project/GCRB/Lee_Lab/s163035/MCF7_ES_PJ34_Shrikanth/E2_ER-alpha/remove-duplicates.sh-1.1.0/E2_Treated_1_ER_S8.filtered.no_dups.bam"), "GRanges")
e_er_R2 <- as(readGAlignments("/Volumes/project/GCRB/Lee_Lab/s163035/MCF7_ES_PJ34_Shrikanth/E2_ER-alpha/remove-duplicates.sh-1.1.0/E2_Treated_2_ER_S11.filtered.no_dups.bam"), "GRanges")
e_foxa1_R1 <- as(readGAlignments("/Volumes/project/GCRB/Lee_Lab/s163035/MCF7_ES_PJ34_Shrikanth/E2_FOXA1/remove-duplicates.sh-1.1.0/E2_Treated_1_FOXA1_S9.filtered.no_dups.bam"), "GRanges")
e_foxa1_R2 <- as(readGAlignments("/Volumes/project/GCRB/Lee_Lab/s163035/MCF7_ES_PJ34_Shrikanth/E2_FOXA1/remove-duplicates.sh-1.1.0/E2_Treated_2_FOXA1_S12.filtered.no_dups.bam"), "GRanges")
e_input_R1 <- as(readGAlignments("/Volumes/project/GCRB/Lee_Lab/s163035/MCF7_ES_PJ34_Shrikanth/E2_Input/remove-duplicates.sh-1.1.0/E2_Treated_1.filtered.no_dups.bam"), "GRanges")
e_input_R2 <- as(readGAlignments("/Volumes/project/GCRB/Lee_Lab/s163035/MCF7_ES_PJ34_Shrikanth/E2_Input/remove-duplicates.sh-1.1.0/E2_Treated_2.filtered.no_dups.bam"), "GRanges")
pj_e_er_R1 <- as(readGAlignments("/Volumes/project/GCRB/Lee_Lab/s163035/MCF7_ES_PJ34_Shrikanth/PJ34_E2_ER-alpha/remove-duplicates.sh-1.1.0/PJ34_E2_Treated_1_ER_S20.filtered.no_dups.bam"), "GRanges")
pj_e_er_R2 <- as(readGAlignments("/Volumes/project/GCRB/Lee_Lab/s163035/MCF7_ES_PJ34_Shrikanth/PJ34_E2_ER-alpha/remove-duplicates.sh-1.1.0/PJ34_E2_Treated_2_ER_S23.filtered.no_dups.bam"), "GRanges")
pj_e_foxa_R1 <- as(readGAlignments("/Volumes/project/GCRB/Lee_Lab/s163035/MCF7_ES_PJ34_Shrikanth/PJ34_E2_FOXA1/remove-duplicates.sh-1.1.0/PJ34_E2_Treated_1_FOXA1_S21.filtered.no_dups.bam"), "GRanges")
pj_e_foxa_R2 <- as(readGAlignments("/Volumes/project/GCRB/Lee_Lab/s163035/MCF7_ES_PJ34_Shrikanth/PJ34_E2_FOXA1/remove-duplicates.sh-1.1.0/PJ34_E2_Treated_2_FOXA1_S24.filtered.no_dups.bam"), "GRanges")
pj_e_input_R1 <- as(readGAlignments("/Volumes/project/GCRB/Lee_Lab/s163035/MCF7_ES_PJ34_Shrikanth/PJ34_E2_Input/remove-duplicates.sh-1.1.0/PJ34_E2_Treated_1.filtered.no_dups.bam"), "GRanges")
pj_e_input_R2 <- as(readGAlignments("/Volumes/project/GCRB/Lee_Lab/s163035/MCF7_ES_PJ34_Shrikanth/PJ34_E2_Input/remove-duplicates.sh-1.1.0/PJ34_E2_Treated_2.filtered.no_dups.bam"), "GRanges")
pj_er_R1 <- as(readGAlignments("/Volumes/project/GCRB/Lee_Lab/s163035/MCF7_ES_PJ34_Shrikanth/PJ34_ER-alpha/remove-duplicates.sh-1.1.0/PJ34_Treated_1_ER_S14.filtered.no_dups.bam"), "GRanges")
pj_er_R2 <- as(readGAlignments("/Volumes/project/GCRB/Lee_Lab/s163035/MCF7_ES_PJ34_Shrikanth/PJ34_ER-alpha/remove-duplicates.sh-1.1.0/PJ34_Treated_2_ER_S17.filtered.no_dups.bam"), "GRanges")
pj_foxa1_R1 <- as(readGAlignments("/Volumes/project/GCRB/Lee_Lab/s163035/MCF7_ES_PJ34_Shrikanth/PJ34_FOXA1/remove-duplicates.sh-1.1.0/PJ34_Treated_1_FOXA1_S15.filtered.no_dups.bam"), "GRanges")
pj_foxa1_R2 <- as(readGAlignments("/Volumes/project/GCRB/Lee_Lab/s163035/MCF7_ES_PJ34_Shrikanth/PJ34_FOXA1/remove-duplicates.sh-1.1.0/PJ34_Treated_2_FOXA1_S18.filtered.no_dups.bam"), "GRanges")
pj_input_R1 <- as(readGAlignments("/Volumes/project/GCRB/Lee_Lab/s163035/MCF7_ES_PJ34_Shrikanth/PJ34_Input/remove-duplicates.sh-1.1.0/PJ34_Treated_1.filtered.no_dups.bam"), "GRanges")
pj_input_R1 <- as(readGAlignments("/Volumes/project/GCRB/Lee_Lab/s163035/MCF7_ES_PJ34_Shrikanth/PJ34_Input/remove-duplicates.sh-1.1.0/PJ34_Treated_2.filtered.no_dups.bam"), "GRanges")
v_er_R1 <- as(readGAlignments("/Volumes/project/GCRB/Lee_Lab/s163035/MCF7_ES_PJ34_Shrikanth/Veh_E2_ER-alpha/remove-duplicates.sh-1.1.0/Vehicle_1_ER_S2.filtered.no_dups.bam"), "GRanges")
v_er_R2 <- as(readGAlignments("/Volumes/project/GCRB/Lee_Lab/s163035/MCF7_ES_PJ34_Shrikanth/Veh_E2_ER-alpha/remove-duplicates.sh-1.1.0/Vehicle_2_ER_S5.filtered.no_dups.bam"), "GRanges")
v_foxa1_R1 <- as(readGAlignments("/Volumes/project/GCRB/Lee_Lab/s163035/MCF7_ES_PJ34_Shrikanth/Veh_E2_FOXA1/remove-duplicates.sh-1.1.0/Vehicle_1_FOXA1_S3.filtered.no_dups.bam"), "GRanges")
v_foxa1_R2 <- as(readGAlignments("/Volumes/project/GCRB/Lee_Lab/s163035/MCF7_ES_PJ34_Shrikanth/Veh_E2_FOXA1/remove-duplicates.sh-1.1.0/Vehicle_2_FOXA1_S6.filtered.no_dups.bam"), "GRanges")
v_input_R1 <- as(readGAlignments("/Volumes/project/GCRB/Lee_Lab/s163035/MCF7_ES_PJ34_Shrikanth/Veh_E2_Input/remove-duplicates.sh-1.1.0/Vehicle_1.filtered.no_dups.bam"), "GRanges")
v_input_R2 <- as(readGAlignments("/Volumes/project/GCRB/Lee_Lab/s163035/MCF7_ES_PJ34_Shrikanth/Veh_E2_Input/remove-duplicates.sh-1.1.0/Vehicle_2.filtered.no_dups.bam"), "GRanges")
## Combine replicates
pj_e_er <- c(pj_e_er_R1, pj_e_er_R2)
pj_e_input <- c(pj_e_input_R1, pj_e_input_R2)
e_er <- c(e_er_R1, e_er_R2)
e_input <- c(e_input_R1, e_input_R2)
pj_e_foxa1 <- c(pj_e_foxa_R1, pj_e_foxa_R2)
e_foxa1 <- c(e_foxa1_R1, e_foxa1_R1)
## Find average
library_pj_e_er <- NROW(pj_e_er)
library_pj_e_input <- NROW(pj_e_input)
library_e_er <- NROW(e_er)
library_e_input <- NROW(e_input)
library_pj_e_foxa1 <- NROW(pj_e_foxa1)
library_e_foxa1 <- NROW(e_foxa1)
## Metagene
#### ER
expReads <- mean(c(library_e_er, library_pj_e_er,library_e_input,library_pj_e_input))
mg_e_er<- runMetaGene(features=closest_er_alpha_luc_summit, reads=e_er, size=100, normCounts=expReads/library_e_er, up = 2000, sampling=FALSE)
mg_e_input<- runMetaGene(features=closest_er_alpha_luc_summit, reads=e_input, size=100, normCounts=expReads/library_e_input, up = 2000, sampling=FALSE)
mg_pj_e_er<- runMetaGene(features=closest_er_alpha_luc_summit, reads=pj_e_er, size=100, normCounts=expReads/library_pj_e_er, up = 2000, sampling=FALSE)
mg_pj_e_input<- runMetaGene(features=closest_er_alpha_luc_summit, reads=pj_e_input, size=100, normCounts=expReads/library_pj_e_input, up = 2000, sampling=FALSE)
### Subtract reads
mg_e_er_enrichment <- mg_e_er$sense / mg_e_input$sense
mg_pj_e_er_enrichment <- mg_pj_e_er$sense / mg_pj_e_input$sense
#### FOXA
expReads <- mean(c(library_e_foxa1, library_pj_e_foxa1,library_e_input,library_pj_e_input))
mg_e_foxa1<- runMetaGene(features=closest_er_alpha_luc_summit, reads=e_foxa1, size=100, normCounts=expReads/library_e_foxa1, up = 2000, sampling=FALSE)
mg_e_input<- runMetaGene(features=closest_er_alpha_luc_summit, reads=e_input, size=100, normCounts=expReads/library_e_input, up = 2000, sampling=FALSE)
mg_pj_e_foxa1<- runMetaGene(features=closest_er_alpha_luc_summit, reads=pj_e_foxa1, size=100, normCounts=expReads/library_pj_e_foxa1, up = 2000, sampling=FALSE)
mg_pj_e_input<- runMetaGene(features=closest_er_alpha_luc_summit, reads=p_er_input, size=100, normCounts=expReads/library_pj_e_input, up = 2000, sampling=FALSE)
### Subtract reads
mg_e_foxa1_enrichment <- mg_e_foxa1$sense / mg_e_input$sense
mg_pj_e_foxa1_enrichment <- mg_pj_e_foxa1$sense / mg_pj_e_input$sense
plotMetaGene <- function(filename, mg_a, mg_b, MIN, MAX, legend){
jpeg(filename)
POS=c(-2000:+1999)
par(font=3, font.lab=2, font.axis=2, mar=c(5.5, 5.5, 2, 2) + 0.3)
plot(POS, mg_a, col=alpha("grey",0.2), type="h", xlim=c(-1000, 1000), cex.axis=1, cex.lab=2, ylim=c(MIN, MAX), ylab="Norm. Read Density", xlab="Position (relative to Peak Summit)")
lines(POS, mg_a, col="black", las=1, lwd=4, type="l")
points(POS, mg_b, col="yellow", type="h",las=1, lwd=4, xlim=c(-1000, 1000))
lines(POS, mg_b, col="yellow2", las=1, lwd=4, type="l")
dev.off()
}
### Plot
legend <- c("E2", "PJ34 + E2")
MAX <- max(c(mg_e_er_enrichment, mg_pj_e_er_enrichment))
MIN <- min(c(mg_e_er_enrichment, mg_pj_e_er_enrichment))
plotMetaGene("closest_er_alpha_peak_pj34.jpg", mg_pj_e_er_enrichment, mg_e_er_enrichment, MIN, MAX, legend)
legend <- c("E2", "PJ34 + E2")
MAX <- max(c(mg_e_foxa1_enrichment, mg_pj_e_foxa1_enrichment))
MIN <- min(c(mg_e_foxa1_enrichment, mg_pj_e_foxa1_enrichment))
plotMetaGene("closest_er_alpha_foxa1_peak_pj34.jpg", mg_pj_e_foxa1_enrichment, mg_e_foxa1_enrichment, MIN, MAX, legend)
## Import files RPKM = numReads / ( geneLength/1000 * totalNumReads/1,000,000 )
# ER-alpha
rpkm_e_er <- countOverlaps(closest_er_alpha_luc,e_er)/ (width(closest_er_alpha_luc)/1000 * library_e_er/1000000)
rpkm_pj_e_er <- countOverlaps(closest_er_alpha_luc,pj_e_er)/ (width(closest_er_alpha_luc)/1000 * library_pj_e_er/1000000)
pj_er_alpha_e <- data.frame(rpkm_e_er,rpkm_pj_e_er)
jpeg('boxplot_er_alpha_e2_rpkm_pj34.jpg')
boxplot(pj_er_alpha_e, col=(c("deepskyblue2","deepskyblue2")), main="ER-alpha", ylab="RPKM",outline=FALSE, notch=TRUE, names=c( "V+E", "PJ+E"), ylim = c(0,35),yaxt="n", cex.axis=1,las=2,lwd=4,lty=1)
axis(2, at=seq(0,35,5))
dev.off()
# Wilcox Rank Test
# Significance 2.915e-08
wilcox.test(pj_er_alpha_e[,1], pj_er_alpha_e[,2],paired=T)
# FOXA1
rpkm_e_er_foxa1 <- countOverlaps(closest_er_alpha_luc,e_foxa1)/ (width(closest_er_alpha_luc)/1000 * library_e_foxa1/1000000)
rpkm_pj_er_foxa1 <- countOverlaps(closest_er_alpha_luc,pj_e_foxa1)/ (width(closest_er_alpha_luc)/1000 * library_pj_e_foxa1/1000000)
pj_foxa1_e <- data.frame(rpkm_e_er_foxa1,rpkm_pj_er_foxa1)
jpeg('boxplot_foxa1_e2_rpkm_pj34.jpg')
boxplot(pj_foxa1_e, col=(c("firebrick3","firebrick3")), main="FOXA1", ylab="RPKM",outline=FALSE, notch=TRUE, names=c("L+E", "P+E"), ylim = c(0,14),yaxt="n", cex.axis=1,las=2,lwd=4,lty=1)
axis(2, at=seq(0,14,2))
dev.off()
# Wilcox Rank Test
# Significance 2.915e-08
wilcox.test(pj_foxa1_e[,1], pj_foxa1_e[,2],paired=T)
## Metagene
#### FOXA1
expReads <- mean(c(library_e_foxa1, library_pj_e_foxa1,library_e_input,library_pj_e_input))
fmg_e_foxa1<- runMetaGene(features=closest_foxa1_luc_summit, reads=e_foxa1, size=100, normCounts=expReads/library_e_foxa1, up = 2000, sampling=FALSE)
fmg_e_input<- runMetaGene(features=closest_foxa1_luc_summit, reads=e_input, size=100, normCounts=expReads/library_e_input, up = 2000, sampling=FALSE)
fmg_pj34_e_foxa1<- runMetaGene(features=closest_foxa1_luc_summit, reads=pj_e_foxa1, size=100, normCounts=expReads/library_pj_e_foxa1, up = 2000, sampling=FALSE)
fmg_pj34_e_input<- runMetaGene(features=closest_foxa1_luc_summit, reads=pj_e_input, size=100, normCounts=expReads/library_pj_e_input, up = 2000, sampling=FALSE)
### Subtract reads
fmg_e_foxa1_enrichment <- fmg_e_foxa1$sense / fmg_e_input$sense
fmg_pj34_e_foxa1_enrichment <- fmg_pj34_e_foxa1$sense / fmg_pj34_e_input$sense
#### ER
expReads <- mean(c(library_e_er, library_pj_e_er,library_e_input,library_pj_e_input))
fmg_e_er<- runMetaGene(features=closest_foxa1_luc_summit, reads=l_er, size=100, normCounts=expReads/library_l_er, up = 2000, sampling=FALSE)
fmg_e_er_input<- runMetaGene(features=closest_foxa1_luc_summit, reads=l_er_input, size=100, normCounts=expReads/library_l_er_input, up = 2000, sampling=FALSE)
fmg_pj34_er<- runMetaGene(features=closest_foxa1_luc_summit, reads=pj_e_er, size=100, normCounts=expReads/library_pj_e_er, up = 2000, sampling=FALSE)
fmg_pj34_input<- runMetaGene(features=closest_foxa1_luc_summit, reads=pj_e_input, size=100, normCounts=expReads/library_pj_e_input, up = 2000, sampling=FALSE)
### Subtract reads
fmg_e_er_enrichment <- fmg_e_er$sense / fmg_e_er_input$sense
fmg_pj34_er_enrichment <- fmg_pj34_er$sense / fmg_pj34_input$sense
plotMetaGene <- function(filename, mg_a, mg_b, MIN, MAX, legend){
jpeg(filename)
POS=c(-2000:+1999)
par(font=3, font.lab=2, font.axis=2, mar=c(5.5, 5.5, 2, 2) + 0.3)
plot(POS, mg_a, col=alpha("cyan",0.2), type="h", xlim=c(-1000, 1000), cex.axis=1, cex.lab=2, ylim=c(MIN, MAX), ylab="Norm. Read Density", xlab="Position (relative to Peak Summit)")
lines(POS, mg_a, col="cyan2", las=1, lwd=4, type="l")
points(POS, mg_b, col="grey", type="h",las=1, lwd=4, xlim=c(-1000, 1000))
lines(POS, mg_b, col="black", las=1, lwd=4, type="l")
dev.off()
}
### Plot
legend <- c("E2", "PJ34 + E2")
MAX <- max(c(fmg_e_foxa1_enrichment, fmg_pj34_e_foxa1_enrichment))
MIN <- min(c(fmg_e_foxa1_enrichment, fmg_pj34_e_foxa1_enrichment))
plotMetaGene("closest_foxa1_peak_pj34.jpg", fmg_e_foxa1_enrichment, fmg_pj34_e_foxa1_enrichment, MIN, MAX, legend)
legend <- c("E2", "PJ34 + E2")
MAX <- max(c(fmg_e_er_enrichment, fmg_pj34_er_enrichment))
MIN <- min(c(fmg_e_er_enrichment, fmg_pj34_er_enrichment))
plotMetaGene("closest_foxa1_er_alpha_peak_pj34.jpg", fmg_e_er_enrichment, fmg_pj34_er_enrichment, MIN, MAX, legend)
## Import files RPKM = numReads / ( geneLength/1000 * totalNumReads/1,000,000 )
# FOXA1
frpkm_e_foxa1 <- countOverlaps(closest_foxa1_luc,e_foxa1)/ (width(closest_foxa1_luc)/1000 * library_e_foxa1/1000000)
frpkm_pj_e_foxa1 <- countOverlaps(closest_foxa1_luc,pj_e_foxa1)/ (width(closest_foxa1_luc)/1000 * library_pj_e_foxa1/1000000)
ffoxa1_e_pj <- data.frame(frpkm_e_foxa1,frpkm_pj_e_foxa1)
jpeg('boxplot_ffoxa1_e2_rpkm_pj34.jpg')
boxplot(ffoxa1_e_pj, col=(c("firebrick3","firebrick3")), main="FOXA1", ylab="RPKM",outline=FALSE, notch=TRUE, names=c("v+E", "PJ+E"), ylim = c(0,14),yaxt="n", cex.axis=1,las=2,lwd=4,lty=1)
axis(2, at=seq(0,14,2))
dev.off()
# Wilcox Rank Test
# Significance 2.915e-08
wilcox.test(ffoxa1_e_pj[,1], ffoxa1_e_pj[,2],paired=T)
# ER-alpha
frpkm_e_er <- countOverlaps(closest_foxa1_luc,e_er)/ (width(closest_foxa1_luc)/1000 * library_pj_e_er/1000000)
frpkm_pj34_er <- countOverlaps(closest_foxa1_luc,pj_e_er)/ (width(closest_foxa1_luc)/1000 * library_pj_e_er/1000000)
fer_alpha_e_pj <- data.frame(frpkm_e_er,frpkm_pj34_er)
jpeg('boxplot_fer_alpha_e2_rpkm_pj34.jpg')
boxplot(fer_alpha_e_pj, col=(c("deepskyblue2","deepskyblue2")), main="ER-alpha", ylab="RPKM",outline=FALSE, notch=TRUE, names=c( "V+E", "PJ+E"), ylim = c(0,1.4),yaxt="n", cex.axis=1,las=2,lwd=4,lty=1)
axis(2, at=seq(0,1.4,.2))
dev.off()
# Wilcox Rank Test
# Significance 2.915e-08
wilcox.test(fer_alpha_e_pj[,1], fer_alpha_e_pj[,2],paired=T)
```
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment