Mercurial > repos > bgruening > diffbind

comparison diffbind.R @ 11:4c7ab9995f9e draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/diffbind commit cc4c1c4131518b9cbf986a1f252767ff73ca938e
author iuc
date Sat, 07 Apr 2018 15:45:41 -0400
parents d7725c5596ab
children 1de83981d43c
comparison
equal deleted inserted replaced
10:d7725c5596ab 11:4c7ab9995f9e
2 options( show.error.messages=F, error = function () { cat( geterrmessage(), file=stderr() ); q( "no", 1, F ) } ) 2 options( show.error.messages=F, error = function () { cat( geterrmessage(), file=stderr() ); q( "no", 1, F ) } )
3 # we need that to not crash galaxy with an UTF8 error on German LC settings. 3 # we need that to not crash galaxy with an UTF8 error on German LC settings.
4 Sys.setlocale("LC_MESSAGES", "en_US.UTF-8") 4 Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")
5 5
6 suppressPackageStartupMessages({ 6 suppressPackageStartupMessages({
7 library('getopt') 7 library('getopt')
8 library('DiffBind') 8 library('DiffBind')
9 library('rjson')
9 }) 10 })
10 11
11 options(stringAsfactors = FALSE, useFancyQuotes = FALSE) 12 options(stringAsfactors = FALSE, useFancyQuotes = FALSE)
12 args <- commandArgs(trailingOnly = TRUE) 13 args <- commandArgs(trailingOnly = TRUE)
13 14
14 #get options, using the spec as defined by the enclosed list. 15 #get options, using the spec as defined by the enclosed list.
15 #we read the options from the default: commandArgs(TRUE). 16 #we read the options from the default: commandArgs(TRUE).
16 spec = matrix(c( 17 spec = matrix(c(
18 'infile' , 'i', 1, "character",
19 'outfile' , 'o', 1, "character",
20 'scorecol', 'n', 1, "integer",
21 'lowerbetter', 'l', 1, "logical",
22 'summits', 's', 1, "integer",
23 'th', 't', 1, "double",
24 'format', 'f', 1, "character",
25 'plots' , 'p', 2, "character",
26 'bmatrix', 'b', 0, "logical",
27 "rdaOpt", "r", 0, "logical",
28 'infoOpt' , 'a', 0, "logical",
17 'verbose', 'v', 2, "integer", 29 'verbose', 'v', 2, "integer",
18 'help' , 'h', 0, "logical", 30 'help' , 'h', 0, "logical"
19 'outfile' , 'o', 1, "character",
20 'plots' , 'p', 2, "character",
21 'infile' , 'i', 1, "character",
22 'format', 'f', 1, "character",
23 'th', 't', 1, "double",
24 'bmatrix', 'b', 0, "logical",
25 "rdaOpt", "r", 0, "logical"
26 ), byrow=TRUE, ncol=4); 31 ), byrow=TRUE, ncol=4);
27 32
28 opt = getopt(spec); 33 opt = getopt(spec);
29 34
30 # if help was asked for print a friendly message 35 # if help was asked for print a friendly message
32 if ( !is.null(opt$help) ) { 37 if ( !is.null(opt$help) ) {
33 cat(getopt(spec, usage=TRUE)); 38 cat(getopt(spec, usage=TRUE));
34 q(status=1); 39 q(status=1);
35 } 40 }
36 41
42 parser <- newJSONParser()
43 parser$addData(opt$infile)
44 factorList <- parser$getObject()
45 filenamesIn <- unname(unlist(factorList[[1]][[2]]))
46 peaks <- filenamesIn[grepl("peaks.bed", filenamesIn)]
47 bams <- filenamesIn[grepl("bamreads.bam", filenamesIn)]
48 ctrls <- filenamesIn[grepl("bamcontrol.bam", filenamesIn)]
49
50 # get the group and sample id from the peaks filenames
51 groups <- sapply(strsplit(peaks,"-"), `[`, 1)
52 samples <- sapply(strsplit(peaks,"-"), `[`, 2)
53
54 if ( length(ctrls) != 0 ) {
55 sampleTable <- data.frame(SampleID=samples,
56 Condition=groups,
57 bamReads=bams,
58 bamControl=ctrls,
59 Peaks=peaks,
60 Tissue=samples, # using "Tissue" column to display ids as labels in PCA plot
61 stringsAsFactors=FALSE)
62 } else {
63 sampleTable <- data.frame(SampleID=samples,
64 Replicate=samples,
65 Condition=groups,
66 bamReads=bams,
67 Peaks=peaks,
68 Tissue=samples,
69 stringsAsFactors=FALSE)
70 }
71
72 sample = dba(sampleSheet=sampleTable, peakFormat='bed', scoreCol=opt$scorecol, bLowerScoreBetter=opt$lowerbetter)
73
74 if ( !is.null(opt$summits) ) {
75 sample_count = dba.count(sample, summits=opt$summits)
76 } else {
77 sample_count = dba.count(sample)
78 }
79
80 sample_contrast = dba.contrast(sample_count, categories=DBA_CONDITION, minMembers=2)
81 sample_analyze = dba.analyze(sample_contrast)
82 diff_bind = dba.report(sample_analyze, th=opt$th)
83
84 # Generate plots
37 if ( !is.null(opt$plots) ) { 85 if ( !is.null(opt$plots) ) {
38 pdf(opt$plots) 86 pdf(opt$plots)
87 orvals = dba.plotHeatmap(sample_analyze, contrast=1, correlations=FALSE, cexCol=0.8, th=opt$th)
88 dba.plotPCA(sample_analyze, contrast=1, th=opt$th, label=DBA_TISSUE, labelSize=0.3)
89 dba.plotMA(sample_analyze, th=opt$th)
90 dba.plotVolcano(sample_analyze, th=opt$th)
91 dba.plotBox(sample_analyze, th=opt$th)
92 dev.off()
39 } 93 }
40 94
41 sample = dba(sampleSheet=opt$infile, peakFormat='bed') 95 # Output differential binding sites
42 sample_count = dba.count(sample)
43 sample_contrast = dba.contrast(sample_count, categories=DBA_CONDITION, minMembers=2)
44 sample_analyze = dba.analyze(sample_contrast)
45 diff_bind = dba.report(sample_analyze)
46 orvals = dba.plotHeatmap(sample_analyze, contrast=1, correlations=FALSE)
47 dev.off()
48
49 resSorted <- diff_bind[order(diff_bind$FDR),] 96 resSorted <- diff_bind[order(diff_bind$FDR),]
50 write.table(as.data.frame(resSorted), file = opt$outfile, sep="\t", quote = FALSE, append=TRUE, row.names = FALSE, col.names = FALSE) 97 write.table(as.data.frame(resSorted), file = opt$outfile, sep="\t", quote = FALSE, append=TRUE, row.names = FALSE)
51 98
52 # Output binding affinity scores 99 # Output binding affinity scores
53 if (!is.null(opt$bmatrix)) { 100 if (!is.null(opt$bmatrix)) {
54 bmat <- dba.peakset(sample_count, bRetrieve=TRUE, DataType=DBA_DATA_FRAME) 101 bmat <- dba.peakset(sample_count, bRetrieve=TRUE, DataType=DBA_DATA_FRAME)
55 write.table(as.data.frame(bmat), file="bmatrix.tab", sep="\t", quote=FALSE, row.names=FALSE, col.names=FALSE) 102 write.table(as.data.frame(bmat), file="bmatrix.tab", sep="\t", quote=FALSE, row.names=FALSE)
56 } 103 }
57 104
58 ## Output RData file 105 # Output RData file
59
60 if (!is.null(opt$rdaOpt)) { 106 if (!is.null(opt$rdaOpt)) {
61 save.image(file = "DiffBind_analysis.RData") 107 save.image(file = "DiffBind_analysis.RData")
62 } 108 }
63 109
64 sessionInfo() 110 # Output analysis info
111 if (!is.null(opt$infoOpt)) {
112 info <- "DiffBind_analysis_info.txt"
113 cat("dba.count Info\n\n", file=info, append = TRUE)
114 capture.output(sample, file=info, append=TRUE)
115 cat("\ndba.analyze Info\n\n", file=info, append = TRUE)
116 capture.output(sample_analyze, file=info, append=TRUE)
117 cat("\nSessionInfo\n\n", file=info, append = TRUE)
118 capture.output(sessionInfo(), file=info, append=TRUE)
119 }