33from defaults import *
44from utils import *
55import csv
6+ import glob
67
78FORMAT = '%(levelname)s %(asctime)-15s %(name)-20s %(message)s'
89logFormatter = logging .Formatter (FORMAT )
@@ -182,6 +183,7 @@ def run_deseq2(quant_files="", alignments="",
182183 if start <= step :
183184 logger .info ("--------------------------STEP %s--------------------------" % step )
184185 command = "%s -e \" library('DESeq2'); load('%s/txi.rda'); \
186+ samples <- c(%s); \
185187 condition <- factor(c(%s)); \
186188 (colData <- data.frame(row.names=colnames(txi$count), condition));\
187189 counts <- round(txi$counts); mode(counts) <- 'integer'; \
@@ -193,12 +195,24 @@ def run_deseq2(quant_files="", alignments="",
193195 lengths <- txi$length; dimnames(lengths) <- dimnames(dds);\
194196 assays(dds)[['avgTxLength']] <- lengths; }; \
195197 dds <- dds[ rowSums(counts(dds)) >= %d, ]; \
196- dds$condition <- relevel(dds$condition, ref='C1'); \
197- dds <- DESeq(dds); res <- results(dds, alpha=%f); \
198- (summary(res)); save(txi,colData,condition,dds,res, \
199- file='%s/deseq2.rda'); write.table(res, file = '%s/deseq2_res.tab', \
200- quote = FALSE, sep='\\ t');\" "% (
201- R , work_deseq2 , "," .join (map (lambda i :"rep('C%d', %d)" % (i ,n_replicates [i ]),range (len (samples )))),
198+ dds <- DESeq(dds); \
199+ for (i in seq_along(samples)){ \
200+ for (j in seq_along(samples)){ \
201+ if (i < j){\
202+ sample1 <- samples[i]; \
203+ sample2 <- samples[j]; \
204+ res <- results(dds, contrast=c('condition',sample1,sample2), alpha=%f); \
205+ (summary(res)); \
206+ res_file= sprintf('%s/deseq2_res_%%s_vs_%%s.tab',sample1,sample2);\
207+ write.table(res, file = res_file, \
208+ quote = FALSE, sep='\\ t'); \
209+ } \
210+ } \
211+ }; \
212+ save(txi,colData,condition,dds,res, \
213+ file='%s/deseq2.rda');\" "% (
214+ R , work_deseq2 , "," .join (map (lambda i :"'sample%d'" % (i ),range (len (samples )))),
215+ "," .join (map (lambda i :"rep('sample%d', %d)" % (i ,n_replicates [i ]),range (len (samples )))),
202216 mincount , alpha , work_deseq2 , work_deseq2 )
203217 cmd = TimedExternalCmd (command , logger , raise_exception = True )
204218 retcode = cmd .run (cmd_log_fd_out = deseq2_log_fd , cmd_log = deseq2_log , msg = msg , timeout = timeout )
@@ -233,17 +247,29 @@ def run_deseq2(quant_files="", alignments="",
233247 command = "%s -e \" library('DESeq2'); countData <- read.table('%s/featureCounts.txt', \
234248 header=TRUE, row.names=1); countData <- countData[ ,6:ncol(countData)]; \
235249 countData <- as.matrix(countData); \
250+ samples <- c(%s); \
236251 condition <- factor(c(%s)); \
237252 (colData <- data.frame(row.names=colnames(countData), condition));\
238253 dds <- DESeqDataSetFromMatrix(countData=countData, colData=colData, design=~ condition);\
239254 dds <- dds[ rowSums(counts(dds)) >= %d, ]; \
240- dds$condition <- relevel(dds$condition, ref='C1'); \
241- dds <- DESeq(dds); res <- results(dds, alpha=%f); \
242- (summary(res)); save(countData,colData,condition,dds,res, \
243- file='%s/deseq2.rda'); write.table(res, file = '%s/deseq2_res.tab', \
244- quote = FALSE, sep='\\ t');\" "% (
245- R , work_deseq2 , "," .join (map (lambda i :"rep('C%d', %d)" % (i ,n_replicates [i ]),
246- range (len (samples )))),
255+ dds <- DESeq(dds); \
256+ for (i in seq_along(samples)){ \
257+ for (j in seq_along(samples)){ \
258+ if (i < j){\
259+ sample1 <- samples[i]; \
260+ sample2 <- samples[j]; \
261+ res <- results(dds, contrast=c('condition',sample1,sample2), alpha=%f); \
262+ (summary(res)); \
263+ res_file= sprintf('%s/deseq2_res_%%s_vs_%%s.tab',sample1,sample2);\
264+ write.table(res, file = res_file, \
265+ quote = FALSE, sep='\\ t'); \
266+ } \
267+ } \
268+ }; \
269+ save(countData,colData,condition,dds,res, \
270+ file='%s/deseq2.rda');\" "% (
271+ R , work_deseq2 , "," .join (map (lambda i :"'sample%d'" % (i ),range (len (samples )))),
272+ "," .join (map (lambda i :"rep('sample%d', %d)" % (i ,n_replicates [i ]),range (len (samples )))),
247273 mincount , alpha , work_deseq2 , work_deseq2 )
248274 cmd = TimedExternalCmd (command , logger , raise_exception = True )
249275 retcode = cmd .run (cmd_log_fd_out = deseq2_log_fd , cmd_log = deseq2_log , msg = msg , timeout = timeout )
@@ -299,17 +325,29 @@ def run_deseq2(quant_files="", alignments="",
299325 command = "%s -e \" library('DESeq2'); countData <- read.table('%s/featureCounts.txt', \
300326 header=TRUE, row.names=1); countData <- countData[ ,6:ncol(countData)]; \
301327 countData <- as.matrix(countData); \
328+ samples <- c(%s); \
302329 condition <- factor(c(%s)); \
303330 (colData <- data.frame(row.names=colnames(countData), condition));\
304331 dds <- DESeqDataSetFromMatrix(countData=countData, colData=colData, design=~ condition);\
305332 dds <- dds[ rowSums(counts(dds)) >= %d, ]; \
306- dds$condition <- relevel(dds$condition, ref='C1'); \
307- dds <- DESeq(dds); res <- results(dds, alpha=%f); \
308- (summary(res)); save(countData,colData,condition,dds,res, \
309- file='%s/deseq2.rda'); write.table(res, file = '%s/deseq2_res.tab', \
310- quote = FALSE, sep='\\ t');\" "% (
311- R , work_deseq2 , "," .join (map (lambda i :"rep('C%d', %d)" % (i ,n_replicates [i ]),
312- range (len (samples )))),
333+ dds <- DESeq(dds); \
334+ for (i in seq_along(samples)){ \
335+ for (j in seq_along(samples)){ \
336+ if (i < j){\
337+ sample1 <- samples[i]; \
338+ sample2 <- samples[j]; \
339+ res <- results(dds, contrast=c('condition',sample1,sample2), alpha=%f); \
340+ (summary(res)); \
341+ res_file= sprintf('%s/deseq2_res_%%s_vs_%%s.tab',sample1,sample2);\
342+ write.table(res, file = res_file, \
343+ quote = FALSE, sep='\\ t'); \
344+ } \
345+ } \
346+ }; \
347+ save(countData,colData,condition,dds, \
348+ file='%s/deseq2.rda');\" "% (
349+ R , work_deseq2 , "," .join (map (lambda i :"'sample%d'" % (i ),range (len (samples )))),
350+ "," .join (map (lambda i :"rep('sample%d', %d)" % (i ,n_replicates [i ]),range (len (samples )))),
313351 mincount , alpha , work_deseq2 , work_deseq2 )
314352 cmd = TimedExternalCmd (command , logger , raise_exception = True )
315353 retcode = cmd .run (cmd_log_fd_out = deseq2_log_fd , cmd_log = deseq2_log , msg = msg , timeout = timeout )
@@ -322,20 +360,21 @@ def run_deseq2(quant_files="", alignments="",
322360 msg = "Copy predictions to output directory for %s." % samples_txt
323361 if start <= step :
324362 logger .info ("--------------------------STEP %s--------------------------" % step )
325- if os .path .exists ("%s/deseq2_res.tab" % work_deseq2 ):
326- command = "cp %s/deseq2_res.tab %s/deseq2_res.tab" % (
327- work_deseq2 , out_deseq2 )
328- cmd = TimedExternalCmd (command , logger , raise_exception = True )
329- retcode = cmd .run (cmd_log_fd_out = deseq2_log_fd , cmd_log = deseq2_log , msg = msg , timeout = timeout )
363+ if len (glob .glob ("%s/deseq2_res*.tab" % work_deseq2 ))> 0 :
364+ for out_file in glob .glob ("%s/deseq2_res*.tab" % work_deseq2 ):
365+ command = "cp %s %s/" % (
366+ out_file , out_deseq2 )
367+ cmd = TimedExternalCmd (command , logger , raise_exception = True )
368+ retcode = cmd .run (cmd_log_fd_out = deseq2_log_fd , cmd_log = deseq2_log , msg = msg , timeout = timeout )
330369 else :
331370 logger .info ("Skipping step %d: %s" % (step ,msg ))
332371 step += 1
333372
334373 diff = ""
335- if os . path . exists ("%s/deseq2_res.tab" % out_deseq2 ):
374+ if len ( glob . glob ("%s/deseq2_res* .tab" % out_deseq2 )) > 0 :
336375 logger .info ("DESeq2 was successfull!" )
337- logger .info ("Output differential expressions: %s/deseq2_res.tab" % out_deseq2 )
338- diff = "%s/deseq2_res.tab" % out_deseq2
376+ logger .info ("Output differential expressions: %s" % ( glob . glob ( "%s /deseq2_res* .tab"% out_deseq2 )) )
377+ diff = glob . glob ( "%s/deseq2_res* .tab" % out_deseq2 )
339378 else :
340379 logger .info ("DESeq2 failed!" )
341380 return diff
0 commit comments