add facets process to pipeline

bin/facets.R

@@ -0,0 +1,34 @@
+#!/usr/bin/env Rscript
+
+library(facets)
+
+args <- commandArgs(TRUE)
+snp_pileup_file <- args[1] # standard snp-pilup output
+output_pdf_file <- args[2] # file to save Plotly PDF to
+output_seg_file <- args[3] # file to save the csv CNV segment
+
+datafile = snp_pileup_file # path to the snp-pileup file
+rcmat = readSnpMatrix(datafile) #read in the matrix
+preproc_sample = preProcSample(rcmat)
+
+## A bivariate genome segmentation is performed on logR and logOR by extending the CBS algotithm (Olshen et al., 2004; Venkatraman and Olshen, 2007) to the bivariate scenario using a Tsquared statistic for identifiying change points ##
+
+oo=procSample(preproc_sample,cval=150) ##  Lower cval lead to higher sensitivity for small changes
+
+## Call allele-specific copy number and associated cellular fraction, estimate tumor purity and ploidy
+
+fit=emcncf(oo)
+
+#head(fit$cncf) ## The segmentation result and the EM fit output looks like this
+#fit$purity ## estimated sample purity
+#fit$ploidy ## estimated sample ploidy
+
+save.image("loaded.Rdata")
+write.csv(fit$cncf,file=output_seg_file,quote = F,row.names = F)
+
+pdf(output_pdf_file,width = 10, height = 10)
+sname<- sprintf('ploidy= %.2f; purity= %.2f', fit$ploidy, fit$purity)
+plotSample(x=oo,emfit=fit,sname=sname)
+dev.off()
+
+## The top panel of the figure displays logR with chromosomes alternating in blue and gray. The green line indicates the median logR in the sample. The purple line indicates the logR of the diploid state. The second panel displays logOR. Segment means are ploted in red lines. The third panel plots the total (black) and minor (red) copy number for each segment. The bottom bar shows the associated cellular fraction (cf). Dark blue indicates high cf. Light blue indicates low cf. Beige indicates a normal segment (total=2,minor=1) ##

main.nf

@@ -3927,12 +3927,13 @@ process cnvkit_plotly {
 }
 
 process snp_pileup {
+    publishDir "${params.outputDir}/cnv", mode: 'copy'
 
     input:
     set val(comparisonID), val(tumorID), file(tumorBam), file(tumorBai), val(normalID), file(normalBam), file(normalBai), file(snp_vcf), file(snp_vcf_tbi) from samples_dd_bam_noHapMap_pairs2.combine(common_snp_vcf).combine(common_snp_vcf_tbi)
 
     output:
-    set val(tumorBamID), file(output_cnvsnp) into snp_pileup
+    set val(prefix), file(output_cnvsnp) into snp_pileup
 
     script:
     caller = "FACETS"
@@ -3954,6 +3955,24 @@ process snp_pileup {
     """
 }
 
+process facets {
+    publishDir "${params.outputDir}/cnv", mode: 'copy'
+
+    input:
+    set val(prefix), file(output_cnvsnp) from snp_pileup
+
+    output:
+    file("${output_segment}")
+    file("${output_pdf}")
+
+    script:
+    output_segment = "${prefix}.segment.csv"
+    output_pdf = "${prefix}.plot.pdf"
+    """
+    facets.R "${output_cnvsnp}" "${output_pdf}" "${output_segment}"
+    """
+}
+
 
 
 

nextflow.config

@@ -549,6 +549,9 @@ profiles {
             withName: snp_pileup {
                 container = "${params.containerDir}/cnv_facets-0.14.0.simg"
             }
+            withName: facets {
+                container = "${params.containerDir}/cnv_facets-0.14.0.simg"
+            }
         }
     }