draw coverage using nextflow

coverage.nf

bedIn = file( params.bed )
reference = file( params.reference )

bamChannel = Channel.fromPath( params.bams ).
		splitCsv(header: false,sep:',',strip:true).
		map{T->T[0]}

String jvarkit(String tool) {
	return "/path/to/jvarkit-git/dist/"+tool+".jar";
	}

process extendBed {
	tag "extend bed"
	input:
		file bed from bedIn
		file ref from reference
	output:
		file("norm.bed") into normalized_bed
	script:
	"""
	module load bedtools

	cut -f1,2 "${ref.toRealPath()}.fai" > jeter.genome
	grep -v '^(browser|track|#)' "${bed}" | cut -f1,2,3 | sed 's/^chr//;s/^/chr/'  | sort -T . -t '\t' -k1,1V -k2,2n > tmp1.bed

	# just get extended start and end
	bedtools slop -i tmp1.bed -g jeter.genome -pct -b "${params.extend}" | cut -f 2,3 > tmp2.bed

	# paste to get chrom/original-start/original-end/x- start/x-end
	paste tmp1.bed tmp2.bed > norm.bed

	rm jeter.genome tmp1.bed tmp2.bed
	"""

	}

process genesInBed {
	tag "bed ${chrom}:${start}-${end}"
	input:
		set chrom,delstart,delend,start,end from normalized_bed.splitCsv(header: false,sep:'\t',strip:true)
	output:
		set chrom,delstart,delend,start,end,file("${chrom}_${start}_${end}.R") into region_R
		
	"""
	set -o pipefail

	java -jar ${jvarkit("kg2bed")} /path/to/wgEncodeGencodeCompV19.txt.gz |\
	awk -F '	' '(\$6=="EXON" && \$1=="${chrom}" && !(int(\$2)>${end} || int(\$3) < ${start}))' |\
	cut -f 1,2,3 |\
	LC_ALL=C sort -T . -t '\t' -k1,1 -k2,2n |\
	bedtools merge |\
	awk -F '\t' '{printf("rect(%d, 0,%d, -1, col=\\\"green\\\")\\n",\$2,\$3);}' > "${chrom}_${start}_${end}.R"
	
	"""
	}

process drawCoverage {
	tag "call ${chrom},${start}-${end} : ${bam}"
	input:
		file ref from reference
		set chrom,delstart,delend,start,end,r_file,bam from region_R.combine(bamChannel)
	output:
		set val("${chrom}:${start}-${end}"),file("out.pdf") into region_pdfs

	script:
	"""
	set -o pipefail
	
	module load samtools r

	samtools view --reference "${ref.toRealPath()}" -H "${bam}" | grep "^@RG" | tr "\t" "\\n" | grep ^SM: -m1 | awk -F ':' '{printf("SM <- \\\"%s\\\"\\n",\$2);}' > jeter.R


samtools depth --reference "${ref.toRealPath()}" -a -r "${chrom}:${start}-${end}" "${bam}" | cut -f 2,3 > jeter.tsv


cat << __EOF__ >> jeter.R

T<-read.table("jeter.tsv", header = FALSE, sep = "\t",colClasses = c("integer","integer"),col.names=c("pos","cov"))

pdf(paste("out.pdf",sep=""), 15, 5)
plot(x=T\\\$pos,y=T\\\$cov,
   main=paste(SM," Coverage ${chrom}:${start}-${end}.",sep=""),
  sub=basename("${bam}"),
   xlab="Position", ylab="Depth ",
  xlim = c(min(T[1]),max(T[1])),
  ylim = c(0,max(T\\\$cov)),
  col=rgb(0.2,0.1,0.5,0.9) , 
  type="l", 
  pch=19)
abline(h=median(T\\\$cov), col="blue")
abline(h=median(T\\\$cov)*0.5, col="orange")
abline(h=median(T\\\$cov)*1.5, col="orange")
abline(v=${delstart}, col="green")
abline(v=${delend}, col="green")
__EOF__

cat ${r_file} >> jeter.R

echo 'dev.off()' >> jeter.R

R --vanilla < jeter.R

	rm -f jeter.R jeter.tsv
	"""
	}

process mergeRegionPdf {
	tag "merge ${rgn} N=${pdfs.size()}"
	input:
		set rgn,pdfs from region_pdfs.groupTuple()
	output:
		file("rgn.pdf") into merged1_pdf
	script:
	"""
	echo -e ".ad c\\n${rgn}" | groff -Tps | ps2pdf - cover.pdf
	gs -dBATCH -dNOPAUSE -q -sDEVICE=pdfwrite -dPDFSETTINGS=/prepress -sOutputFile=rgn.pdf cover.pdf ${pdfs.join(" ")}
	rm cover.pdf	
	"""

	}

process mergePDFs {
	tag "merge pdfs N=${vcfs.size()}"

	input:
		val vcfs from merged1_pdf.collect()
	output:
		file("${params.prefix}pdf") into merged_pdf
	script:
	"""
	gs -dBATCH -dNOPAUSE -q -sDEVICE=pdfwrite -dPDFSETTINGS=/prepress -sOutputFile=${params.prefix}pdf ${vcfs.join(" ")}
	"""	
	}

run.sh

nextflow  run -with-dag workflow.dot -with-trace -resume -work-dir $(OUTDIR) \
		--prefix "biostar" \
		--reference "/path/to/reference.fasta" \
		--extend 0.3 \
		--bams input/bam.list \
		--bed input.bed \
		 coverage.nf