diff --git a/workflow/scripts/runDiffBind.R b/workflow/scripts/runDiffBind.R
new file mode 100644
index 0000000000000000000000000000000000000000..e95d45421cead97070c3c93cf60d29725c9763fc
--- /dev/null
+++ b/workflow/scripts/runDiffBind.R
@@ -0,0 +1,30 @@
+library("DiffBind")
+
+#build dba object from sample sheet and do analysis
+data <- dba(sampleSheet="samplesheet.csv")
+data <- dba.count(data, summits=250)
+data <- dba.contrast(data, minMembers = 2, categories=DBA_CONDITION)
+data <- dba.analyze(data)
+
+#Draw figures
+pdf("diffbind.samples.heatmap.pdf")
+plot(data)
+dev.off()
+
+pdf("diffbind.samples.pca.pdf")
+dba.plotPCA(data, DBA_TISSUE, label=DBA_CONDITION)
+dev.off()
+
+#Save peak reads count
+normcount <- dba.peakset(dta, bRetrieve=T)
+write.table(as.data.frame(normcount),"diffbind.normcount.txt",sep="\t",quote=F,row.names=F)
+
+#Retriving the differentially bound sites
+for (i in c(1:length(data$contrasts)))
+{
+ contrast_name = paste(data$contrasts[[i]]$name1,"vs",
+                      data$contrasts[[i]]$name2,"diffbind.xls",sep="_")
+ report <- dba.report(data, contrast=i, th=1, bCount=TRUE)
+ write.table(as.data.frame(report),contrast_name,sep="\t",quote=F,row.names=F)
+}
+