From 496066b20dc67e6d13f35889b49ef93d2550fa81 Mon Sep 17 00:00:00 2001
From: Beibei Chen <beibei.chen@utsouthwestern.edu>
Date: Mon, 21 Nov 2016 14:41:39 -0600
Subject: [PATCH] add diffbind running script
---
workflow/scripts/runDiffBind.R | 30 ++++++++++++++++++++++++++++++
1 file changed, 30 insertions(+)
create mode 100644 workflow/scripts/runDiffBind.R
diff --git a/workflow/scripts/runDiffBind.R b/workflow/scripts/runDiffBind.R
new file mode 100644
index 0000000..e95d454
--- /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)
+}
+
--
GitLab