library(DiffBind)setwd("/Volumes/Terranova2/KS_SMINT_FINAL_Analysis_May2018/KS-SmINT_FED_FASTED_REFED_Fastq/KS_Intestine_ChIP-seq_FINAL_BigWig_20Million_Reads_e5__02-01-18/04aln_downsample/")
#create a dba.object#
H3K9bhb <- dba(sampleSheet = "H3K9bhb_singlemark.csv", minOverlap = 2)
H4K12bhb <- dba(sampleSheet = "H4K12bhb_singlemark.csv", minOverlap = 2)
#count reads in binding site intervals: Scoring metric can be changed#
H3K9bhb_RPKM <- dba.count(H3K9bhb, minOverlap = 2, score = DBA_SCORE_RPKM, fragmentSize = 200)
H4K12bhb_RPKM <- dba.count(H4K12bhb, minOverlap = 2, score = DBA_SCORE_RPKM, fragmentSize = 200)
#establish contrast for differential analysis#
H3K9bhb_RPKM_contrast <- dba.contrast(H3K9bhb_RPKM, categories = DBA_TISSUE, minMembers = 2)
H4K12bhb_RPKM_contrast <- dba.contrast(H4K12bhb_RPKM, categories = DBA_TISSUE, minMembers = 2)
#perform differential analysis#
H3K9bhb_RPKM_analysis <- dba.analyze(H3K9bhb_RPKM_contrast, method = DBA_DESEQ2, bSubControl = TRUE, bFullLibrarySize = TRUE)
H4K12bhb_RPKM_analysis <- dba.analyze(H4K12bhb_RPKM_contrast, method = DBA_DESEQ2, bSubControl = TRUE, bFullLibrarySize = TRUE)#MAPlot
dba.plotMA(H3K9bhb_RPKM_analysis, method = DBA_DESEQ2, th = 0.01)dba.plotMA(H4K12bhb_RPKM_analysis, method = DBA_DESEQ2, th = 0.01)Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.