Exercise part4 – Alternative approach in R to plot and visualize the data

As an alternative to deeptools2, we can always use R starting from the BigWig files. Below is the R code that achieves this.

You will need to use the library CoverageView (available through bioconductor).

The code has “hard-coded” paths in it, but you should, by now, be able to interact with R enough to change these accordingly.

(Note: You will obviously need to download the BigWig files to your local machine for this to work).

options(width=40)
library(CoverageView)


# Setting the working directory in R to where the BED Files are

setwd("/Users/nd48/Documents/Analysis/NCB-216/bedFiles")


# Loading the BigWig files for each sample and then creating a 
# the object for each sample using the "CreateBigWigFile"

kobrg1BigWIGfile<-("Sample_KO_Brg1_Chip.bw")
kobrg1<-CoverageBigWigFile(kobrg1BigWIGfile,reads_mapped=25626121)

koh3k9BigWIGfile<-("Sample_KO_H3K9Me3_Chip.bw")
koh3k9<-CoverageBigWigFile(koh3k9BigWIGfile,reads_mapped=23568151)

koinputBigWIGfile<-("Sample_KO_Input.bw")
koinput<-CoverageBigWigFile(koinputBigWIGfile,reads_mapped=42594560)

wtbrg1BigWIGfile<-("Sample_WT_Brg1_Chip.bw")
wtbrg1<-CoverageBigWigFile(wtbrg1BigWIGfile,reads_mapped=14848901)

wth3k9BigWIGfile<-("Sample_WT_H3K9Me3_Chip.bw")
wth3k9<-CoverageBigWigFile(wth3k9BigWIGfile,reads_mapped=41502848)

wtinputBigWIGfile<-("Sample_WT_Input.bw")
wtinput<-CoverageBigWigFile(wtinputBigWIGfile,reads_mapped=78286472)


# Loading the BED files for the UP, Down, and Control genes (no DE)

upGenes<-("/Users/nd48/Documents/Analysis/NCB-216/Gene_list/TSSbeds/upRegTSS.bed")
downGenes<-("/Users/nd48/Documents/Analysis/NCB-216/Gene_list/TSSbeds/downRegTSS.bed")
noGenes<-("/Users/nd48/Documents/Analysis/NCB-216/Gene_list/TSSbeds/noRegTSS.bed")

##################################################################
########################### Up Genes #############################
##################################################################

# Calculating the coverage matrix for each BigWig sample file around
# the TSS of UP regulated genes. We are scanning 5kb upstream and downstream
# of the defined TSS at a window size (bin) of every 10bp. We are then printing
# the coverage prfile (write.profile) to a text file

kobrg1upgenes<-cov.matrix(kobrg1,coordfile=upGenes,extend=5000,bin_width=10)
write.profile(kobrg1upgenes,outfile = "/Users/nd48/Documents/Analysis/NCB-216/Gene_list/Coview/kobrg1upgenes5kbTSS.txt")

koh3k9upgenes<-cov.matrix(koh3k9,coordfile=upGenes,extend=5000,bin_width=10)
write.profile(koh3k9upgenes,outfile = "/Users/nd48/Documents/Analysis/NCB-216/Gene_list/Coview/koh3k9upgenes5kbTSS.txt")

koinputupgenes<-cov.matrix(koinput,coordfile=upGenes,extend=5000,bin_width=10)
write.profile(koinputupgenes,outfile = "/Users/nd48/Documents/Analysis/NCB-216/Gene_list/Coview/koinputupgenes5kbTSS.txt")

wtbrg1upgenes<-cov.matrix(wtbrg1,coordfile=upGenes,extend=5000,bin_width=10)
write.profile(wtbrg1upgenes,outfile = "/Users/nd48/Documents/Analysis/NCB-216/Gene_list/Coview/wtbrg1upgenes5kbTSS.txt")

wth3k9upgenes<-cov.matrix(wth3k9,coordfile=upGenes,extend=5000,bin_width=10)
write.profile(wth3k9upgenes,outfile = "/Users/nd48/Documents/Analysis/NCB-216/Gene_list/Coview/wth3k9upgenes5kbTSS.txt")

wtinputupgenes<-cov.matrix(wtinput,coordfile=upGenes,extend=5000,bin_width=10)
write.profile(wtinputupgenes,outfile = "/Users/nd48/Documents/Analysis/NCB-216/Gene_list/Coview/wtinputupgenes5kbTSS.txt")


##################################################################
########################### Down Genes #############################
##################################################################

# Same as above but for DOWN regulated genes

kobrg1downgenes<-cov.matrix(kobrg1,coordfile=downGenes,extend=5000,bin_width=10)
write.profile(kobrg1downgenes,outfile = "/Users/nd48/Documents/Analysis/NCB-216/Gene_list/Coview/kobrg1downgenes5kbTSS.txt")

koh3k9downgenes<-cov.matrix(koh3k9,coordfile=downGenes,extend=5000,bin_width=10)
write.profile(koh3k9downgenes,outfile = "/Users/nd48/Documents/Analysis/NCB-216/Gene_list/Coview/koh3k9downgenes5kbTSS.txt")

koinputdowngenes<-cov.matrix(koinput,coordfile=downGenes,extend=5000,bin_width=10)
write.profile(koinputdowngenes,outfile = "/Users/nd48/Documents/Analysis/NCB-216/Gene_list/Coview/koinputdowngenes5kbTSS.txt")

wtbrg1downgenes<-cov.matrix(wtbrg1,coordfile=downGenes,extend=5000,bin_width=10)
write.profile(wtbrg1downgenes,outfile = "/Users/nd48/Documents/Analysis/NCB-216/Gene_list/Coview/wtbrg1downgenes5kbTSS.txt")

wth3k9downgenes<-cov.matrix(wth3k9,coordfile=downGenes,extend=5000,bin_width=10)
write.profile(wth3k9downgenes,outfile = "/Users/nd48/Documents/Analysis/NCB-216/Gene_list/Coview/wth3k9downgenes5kbTSS.txt")

wtinputdowngenes<-cov.matrix(wtinput,coordfile=downGenes,extend=5000,bin_width=10)
write.profile(wtinputdowngenes,outfile = "/Users/nd48/Documents/Analysis/NCB-216/Gene_list/Coview/wtinputdowngenes5kbTSS.txt")


##################################################################
########################### No Genes #############################
##################################################################


# And again for our CONTROL (noDE) genes.

kobrg1nogenes<-cov.matrix(kobrg1,coordfile=noGenes,extend=5000,bin_width=10)
write.profile(kobrg1nogenes,outfile = "/Users/nd48/Documents/Analysis/NCB-216/Gene_list/Coview/kobrg1nogenes5kbTSS.txt")

koh3k9nogenes<-cov.matrix(koh3k9,coordfile=noGenes,extend=5000,bin_width=10)
write.profile(koh3k9nogenes,outfile = "/Users/nd48/Documents/Analysis/NCB-216/Gene_list/Coview/koh3k9nogenes5kbTSS.txt")

koinputnogenes<-cov.matrix(koinput,coordfile=noGenes,extend=5000,bin_width=10)
write.profile(koinputnogenes,outfile = "/Users/nd48/Documents/Analysis/NCB-216/Gene_list/Coview/koinputnogenes5kbTSS.txt")

wtbrg1nogenes<-cov.matrix(wtbrg1,coordfile=noGenes,extend=5000,bin_width=10)
write.profile(wtbrg1nogenes,outfile = "/Users/nd48/Documents/Analysis/NCB-216/Gene_list/Coview/wtbrg1nogenes5kbTSS.txt")

wth3k9nogenes<-cov.matrix(wth3k9,coordfile=noGenes,extend=5000,bin_width=10)
write.profile(wth3k9nogenes,outfile = "/Users/nd48/Documents/Analysis/NCB-216/Gene_list/Coview/wth3k9nogenes5kbTSS.txt")

wtinputnogenes<-cov.matrix(wtinput,coordfile=noGenes,extend=5000,bin_width=10)
write.profile(wtinputnogenes,outfile = "/Users/nd48/Documents/Analysis/NCB-216/Gene_list/Coview/wtinputnogenes5kbTSS.txt")



# Finally, we are plotting the coverage profiles and heatmaps 
# (in this case for the DOWN regulated genes), from the objects that
# we calculated earlier. Repeat this process for the rest of the objects
# (UP and noDE gene objects).


# Draw profiles

draw.profile(kobrg1nogenes,ylab="avg coverage",outfile="/Users/nd48/Documents/Analysis/NCB-216/Gene_list/BEDs/kobrg1nogenes.png",main="KO Brg1 NoDE Genes")

draw.profile(koh3k9nogenes,ylab="avg coverage",outfile="/Users/nd48/Documents/Analysis/NCB-216/Gene_list/BEDs/koh3k9nogenes.png",main="KO H3K9 NoDE Genes")

draw.profile(koinputnogenes,ylab="avg coverage",outfile="/Users/nd48/Documents/Analysis/NCB-216/Gene_list/BEDs/koinputnogenes.png",main="KO INPUT NoDE Genes")

draw.profile(wtbrg1nogenes,ylab="avg coverage",outfile="/Users/nd48/Documents/Analysis/NCB-216/Gene_list/BEDs/wtbrg1nogenes.png",main="WT Brg1 NoDE Genes")

draw.profile(wth3k9nogenes,ylab="avg coverage",outfile="/Users/nd48/Documents/Analysis/NCB-216/Gene_list/BEDs/wth3k9nogenes.png",main="WT H3K9 NoDE Genes")

draw.profile(wtinputnogenes,ylab="avg coverage",outfile="/Users/nd48/Documents/Analysis/NCB-216/Gene_list/BEDs/wtinputnogenes.png",main="WT INPUT NoDE Genes")


# Draw heatmaps

draw.heatmap(kobrg1nogenes,outfile="/Users/nd48/Documents/Analysis/NCB-216/Gene_list/BEDs/kobrg1nogenes_heatmap.png",main="KO Brg1 NoDE Genes")

draw.heatmap(koh3k9nogenes,outfile="/Users/nd48/Documents/Analysis/NCB-216/Gene_list/BEDs/koh3k9nogenes_heatmap.png",main="KO H3K9 NoDE Genes")

draw.heatmap(koinputnogenes,outfile="/Users/nd48/Documents/Analysis/NCB-216/Gene_list/BEDs/koinputnogenes_heatmap.png",main="KO INPUT NoDE Genes")

draw.heatmap(wtbrg1nogenes,outfile="/Users/nd48/Documents/Analysis/NCB-216/Gene_list/BEDs/wtbrg1nogenes_heatmap.png",main="WT Brg1 NoDE Genes")

draw.heatmap(wth3k9nogenes,outfile="/Users/nd48/Documents/Analysis/NCB-216/Gene_list/BEDs/wth3k9nogenes_heatmap.png",main="WT H3K9 NoDE Genes")

draw.heatmap(wtinputnogenes,outfile="/Users/nd48/Documents/Analysis/NCB-216/Gene_list/BEDs/wtinputnogenes_heatmap.png",main="WT INPUT NoDE Genes")