diffbind works

workflow/main.nf

@@ -4,8 +4,7 @@
 // params.bams="$baseDir/../test/*.bam"
    params.testpath="/project/BICF/BICF_Core/bchen4/chipseq_analysis/test/"
    params.design="/project/BICF/BICF_Core/bchen4/chipseq_analysis/test/samplesheet.csv"
-   params.genebed="/project/BICF/BICF_Core/bchen4/chipseq_analysis/test/genome/hg19/gene.bed"
-
+   params.genomepath="/project/BICF/BICF_Core/bchen4/chipseq_analysis/test/genome/hg19/"
 // design_file = file(params.design)
 // bams=file(params.bams)
 //gtf_file = file("$params.genome/gencode.gtf")
@@ -20,38 +19,37 @@ process processdesign {
 //   file annotation Tdx
    output:
      file "new_design" into deeptools_design
-//     file "new_design" into diffbind_design
+     file "new_design" into diffbind_design
 
      script:
      """
      module load python/2.7.x-anaconda
      source activate /project/shared/bicf_workflow_ref/chipseq_bchen4/
      python $baseDir/scripts/preprocessDesign.py -i ${params.design} 
-"""
+     """
 }
 
-process run_deeptools {
+//process run_deeptools {
 //   publishDir "$baseDir/output", mode: 'copy'
-   input:
-     file deeptools_design_file from deeptools_design
+//   input:
+//     file deeptools_design_file from deeptools_design
 //   file annotation Tdx
-   output:
-     stdout result
-     script:
-     """
-     module load python/2.7.x-anaconda
-     source activate /project/shared/bicf_workflow_ref/chipseq_bchen4/
-     module load deeptools/2.3.5 
-     python $baseDir/scripts/runDeepTools.py -i $deeptools_design_file -g ${params.genebed}
-"""
-}
+//   output:
+//     stdout result
+//     script:
+//     """
+//     module load python/2.7.x-anaconda
+//     source activate /project/shared/bicf_workflow_ref/chipseq_bchen4/
+//     module load deeptools/2.3.5 
+//     python $baseDir/scripts/runDeepTools.py -i $deeptools_design_file -g ${params.genomepath}}
+//"""
+//}
 
 
 process run_diffbind {
 //   publishDir "$baseDir/output", mode: 'copy'
    input:
      file diffbind_design_file from diffbind_design
-//   file annotation Tdx
    output:
      file "*_diffbind.bed" into diffpeaks_chipseeker
      file "*_diffbind.bed" into diffpeaks_meme
@@ -60,15 +58,15 @@ process run_diffbind {
      """
      module load python/2.7.x-anaconda
      source activate /project/shared/bicf_workflow_ref/chipseq_bchen4/
-     Rscript $baseDir/scripts/runDiffBind.R -i $diffbind_design_file
+     Rscript $baseDir/scripts/runDiffBind.R $diffbind_design_file
 """
 }
 
 //process run_chipseeker {
-////   publishDir "$baseDir/output", mode: 'copy'
+//   publishDir "$baseDir/output", mode: 'copy'
 //   input:
 //     file diffbind_design_file from diffbind_design
-////   file annotation Tdx
+//   file annotation Tdx
 //   output:
 //     file "*_diffbind.bed" into diffpeaks_chipseeker
 //     file "*_diffbind.bed" into diffpeaks_meme
@@ -77,8 +75,21 @@ process run_diffbind {
 //     """
 //     module load python/2.7.x-anaconda
 //     source activate /project/shared/bicf_workflow_ref/chipseq_bchen4/
-//     Rscript $baseDir/scripts/runDiffBind.R -i $diffbind_design_file
+//     Rscript $baseDir/scripts/runDiffBind.R $diffbind_design_file
 //"""
 //}
 
-
+//process run_meme {
+//   publishDir "$baseDir/output", mode: 'copy'
+//   input:
+//     file peaks_meme from diffpeaks_meme.flatten()
+//   output:
+//     stdout result
+//     script:
+//     """
+//     module load python/2.7.x-anaconda
+//     source activate /project/shared/bicf_workflow_ref/chipseq_bchen4/
+//     module load meme/4.11.1-gcc-openmpi
+//     python $baseDir/scripts/runMemechip.py -i $peaks_meme
+//"""
+//}

workflow/scripts/runChipseeker.R

@@ -7,11 +7,22 @@ args = commandArgs(trailingOnly=TRUE)
 #}
 
 library(ChIPseeker)
+if args[3]=="hg19"
+{ 
 library(TxDb.Hsapiens.UCSC.hg19.knownGene)
 txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
-#library(clusterProfiler)
+}
+if args[3]=="mm10"
+{ 
+library(TxDb.Hsapiens.UCSC.mm10.knownGene)
+txdb <- TxDb.Hsapiens.UCSC.mm10.knownGene
+}
+if args[3]=="hg38"
+{ 
+library(TxDb.Hsapiens.UCSC.hg38.knownGene)
+txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
+}
 
-#files = list.files(".")
 files<-as.list(unlist(strsplit(args[1],",")))
 names(files)<-as.list(unlist(strsplit(args[2],",")))
 print(files)
@@ -37,10 +48,4 @@ for(index in c(1:length(peakAnnoList)))
 
 
 }
-#promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
-#tagMatrixList <- lapply(files, getTagMatrix, windows=promoter)
-
-#plotAvgProf(tagMatrixList, xlim=c(-3000, 3000), facet="row")
-
-#overlappeakfiles <- as.list(list.files("overlappeaks/"))
 

workflow/scripts/runDeepTools.py

@@ -47,7 +47,7 @@ def bam2bw_wrapper(command):
 
 def run_signal(files, labels, genome):
   #compute matrix and draw profile and heatmap
-  gene_bed = genome#"/project/BICF/BICF_Core/bchen4/chipseq_analysis/test/genome/"+genome+"/gene.bed"
+  gene_bed = genome+"gene.bed"#"/project/BICF/BICF_Core/bchen4/chipseq_analysis/test/genome/"+genome+"/gene.bed"
   bw_commands = []
   for f in files:
     bw_commands.append("bamCoverage -bs 10 -b "+f+" -o "+f.replace("bam","bw"))