library(DiffBind)#set working directory and place sample.csv file in directory
setwd("/Volumes/Terranova2/PROSTATE_12-19-16/PRAD_FINAL_Analysis_June2018/PRAD_Tumors_Normals/04aln_downsample/")
#create a dba.object#
H3K4me1 <- dba(sampleSheet = "PRAD_H3K4me1_DiffBind.csv", minOverlap = 2)
H3K27ac <- dba(sampleSheet = "PRAD_H3K27ac_DiffBind.csv", minOverlap = 2)
#count reads in binding site intervals: Scoring metric can be changed and does not influence differential analysis#
H3K4me1_counts_READS_MINUS <- dba.count(H3K4me1, minOverlap = 2, score = DBA_SCORE_READS_MINUS, fragmentSize = 200)
H3K27ac_counts_READS_MINUS <- dba.count(H3K27ac, minOverlap = 2, score = DBA_SCORE_READS_MINUS, fragmentSize = 200)
#establish contrast for differential analysis#
H3K4me1_contrast <- dba.contrast(H3K4me1_counts_READS_MINUS, categories = DBA_TISSUE, minMembers = 2)
H3K27ac_contrast <- dba.contrast(H3K27ac_counts_READS_MINUS, categories = DBA_TISSUE, minMembers = 2)
#perform differential analysis#
H3K4me1_READS_MINUS_analysis <- dba.analyze(H3K4me1_contrast, method = DBA_DESEQ2, bSubControl = TRUE, bFullLibrarySize = TRUE)
H3K27ac_READS_MINUS_analysis <- dba.analyze(H3K27ac_contrast, method = DBA_DESEQ2, bSubControl = TRUE, bFullLibrarySize = TRUE)#MAPlot
dba.plotMA(H3K4me1_READS_MINUS_analysis, contrast = 1, method = DBA_DESEQ2, th = 0.005, bUsePval = TRUE, bNormalized = TRUE)dba.plotMA(H3K27ac_READS_MINUS_analysis, contrast = 1, method = DBA_DESEQ2, th = 0.005, bUsePval = TRUE, bNormalized = TRUE)Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.