🛠️ Modification to bin cobalt outputs and add purity forcing

.github/workflows/ci.yml

@@ -19,11 +19,12 @@ jobs:
  sudo mv nf-test /usr/local/bin/
  - name: Pull Docker image and cache
  run: |
- docker pull papaemmelab/purple:v0.1.0
+ docker pull papaemmelab/purple:v0.1.1
  - name: Run unit tests of each process for Amber, Cobalt, Purple
  run: |
  nf-test test tests/main.runamber.nf.test
  nf-test test tests/main.runcobalt.nf.test
+ nf-test test tests/main.bincobalt.nf.test
  nf-test test tests/main.runpurple.nf.test
  - name: Run pipeline end-to-end test
  run: |

Dockerfile

@@ -1,7 +1,15 @@
 FROM papaemmelab/docker-hmftools:v1.0.0
 
+# Clean up to free space
+RUN apt-get clean
+
 # Install dependencies
-RUN apt-get update && apt-get install -y tar curl python3
+RUN apt-get update && apt-get install -y \
+ tar \
+ curl \
+ python3 \
+ python3-pip \
+ python3-venv
 
 # Install Google Cloud SDK
 RUN curl -sSL https://sdk.cloud.google.com | bash
@@ -20,6 +28,10 @@ RUN \
  tar -xzf /data/${REF_DIR}.tar.gz -C /data && \
  rm /data/${REF_DIR}.tar.gz
 
+# Link python3 to python and install packages
+RUN ln -s /usr/bin/python3 /usr/bin/python && \
+ pip install --no-cache-dir numpy pandas
+
 COPY . /app
 WORKDIR /app
 

README.md

@@ -44,6 +44,11 @@ More information on their [docs](https://github.com/hartwigmedical/hmftools/blob
 sh tests/run_test.sh
 ```
 
+To check tests and update snapshots run:
+```bash
+nf-test test --update-snapshot
+```
+
 ### Docker
 
 [Purple Docker image](https://hub.docker.com/r/papaemmelab/purple) was built for several platforms using [docker buildx](https://docs.docker.com/buildx/working-with-buildx/).

main.nf

@@ -1,4 +1,5 @@
-params.cores = 4
+params.cores = 1 
+params.memory = '4 GB'
 
 // Params Defaults in juno
 params.refGenome = "/work/isabl/ref/homo_sapiens/GRCh37d5/gr37.fasta"
@@ -8,17 +9,26 @@ params.loci = "/data/copy_number/GermlineHetPon.37.vcf.gz"
 params.gcProfile = "/data/copy_number/GC_profile.1000bp.37.cnp"
 params.ensemblDataDir = "/data/common/ensembl_data"
 params.diploidRegions = "/data/copy_number/DiploidRegions.37.bed.gz"
+params.binProbes = 0
+params.binLogR = 0
+params.minPurity = 0.08
+params.maxPurity = 1.0
 
 
 log.info """\
  HMFTOOLS - PURPLE
  ========================================
  Params:
  ----------------------------------------
- tumor : ${params.tumor}
- tumorBam : ${params.tumorBam}
- outdir : ${params.outdir}
- cores : ${params.cores}
+ tumor : ${params.tumor}
+ tumorBam : ${params.tumorBam}
+ outdir : ${params.outdir}
+ cores : ${params.cores}
+ memory : ${params.memory}
+ binProbes : ${params.binProbes}
+ binLogR : ${params.binLogR}
+ minPurity : ${params.minPurity}
+ maxPurity : ${params.maxPurity}
  ========================================
  Workflow:
  ----------------------------------------
@@ -32,7 +42,7 @@ process runAmber {
  tag "AMBER on ${params.tumor}"
  publishDir "${params.outdir}/amber", mode: 'copy'
  cpus params.cores
- memory '32 GB'
+ memory params.memory
  time '1h'
 
  input:
@@ -46,21 +56,30 @@ process runAmber {
 
  script:
  """
- amber \
- -tumor ${tumor} \
- -tumor_bam ${tumorBam} \
- -output_dir \$PWD \
- -threads ${params.cores} \
- -loci ${params.loci} \
- -ref_genome_version ${params.genomeVersion}
+ if [ -f "${params.outdir}/amber/${tumor}.amber.baf.tsv.gz" ] && \
+ [ -f "${params.outdir}/amber/${tumor}.amber.baf.pcf" ] && \
+ [ -f "${params.outdir}/amber/${tumor}.amber.qc" ]; then
+ echo "Output files already exist. Skipping amber execution."
+ ln -s ${params.outdir}/amber/${tumor}.amber.baf.tsv.gz ${tumor}.amber.baf.tsv.gz
+ ln -s ${params.outdir}/amber/${tumor}.amber.baf.pcf ${tumor}.amber.baf.pcf
+ ln -s ${params.outdir}/amber/${tumor}.amber.qc ${tumor}.amber.qc
+ else
+ amber \
+ -tumor ${tumor} \
+ -tumor_bam ${tumorBam} \
+ -output_dir \$PWD \
+ -threads ${params.cores} \
+ -loci ${params.loci} \
+ -ref_genome_version ${params.genomeVersion}
+ fi
  """.stripIndent()
 }
 
 process runCobalt {
  tag "COBALT on ${params.tumor}"
  publishDir "${params.outdir}/cobalt", mode: 'copy'
  cpus params.cores
- memory '32 GB'
+ memory params.memory
  time '1h'
 
  input:
@@ -73,21 +92,102 @@ process runCobalt {
 
  script:
  """
- cobalt \
+ if [ -f "${params.outdir}/cobalt/${tumor}.cobalt.ratio.tsv.gz" ] && \
+ [ -f "${params.outdir}/cobalt/${tumor}.cobalt.ratio.pcf" ]; then
+ echo "Output files already exist. Skipping cobalt execution."
+ ln -s ${params.outdir}/cobalt/${tumor}.cobalt.ratio.tsv.gz ${tumor}.cobalt.ratio.tsv.gz
+ ln -s ${params.outdir}/cobalt/${tumor}.cobalt.ratio.pcf ${tumor}.cobalt.ratio.pcf
+ else
+ cobalt \
  -tumor ${tumor} \
  -tumor_bam ${tumorBam} \
  -output_dir \$PWD \
  -threads ${params.cores} \
  -gc_profile ${params.gcProfile} \
  -tumor_only_diploid_bed ${params.diploidRegions}
+ fi
+ """.stripIndent()
+}
+
+process binCobalt {
+ tag "COBALT BIN on ${params.tumor}"
+ publishDir "${params.outdir}/cobalt/binned_${params.binProbes}_probes_${params.binLogR}_LogR", mode: 'copy'
+ cpus params.cores
+ memory params.memory
+ time '1h'
+
+ input:
+ val tumor
+ val binProbes
+ val binLogR
+ path cobalt_ratio_pcf
+ path cobalt_ratio_tsv
+
+ output:
+ path "${tumor}.cobalt.ratio.tsv.gz", emit: cobalt_ratio_tsv
+ path "${tumor}.cobalt.ratio.pcf", emit: cobalt_ratio_pcf
+
+ script:
+ """
+ #!/usr/bin/env python
+ import pandas as pd
+ import numpy as np
+
+ cobalt_ratio_pcf = pd.read_csv('${cobalt_ratio_pcf}', sep='\\t')
+ cobalt_ratio_pcf_probes = pd.DataFrame(columns=cobalt_ratio_pcf.columns)
+
+ # First bin by probes
+ chrom_arm = None
+ last_idx = None
+ for idx, seg in cobalt_ratio_pcf.iterrows():
+ if chrom_arm != '_'.join(seg[['chrom','arm']].astype(str)):
+ chrom_arm = '_'.join(seg[['chrom','arm']].astype(str))
+ cobalt_ratio_pcf_probes = pd.concat([cobalt_ratio_pcf_probes, seg.to_frame().T], ignore_index=True)
+ last_idx = cobalt_ratio_pcf_probes.index[-1]
+ continue
+ if (
+ cobalt_ratio_pcf_probes.loc[last_idx, 'n.probes'] <= ${binProbes}
+ ) or (
+ seg['n.probes'] <= ${binProbes}
+ ):
+ means = [cobalt_ratio_pcf_probes.loc[last_idx, 'mean']] * cobalt_ratio_pcf_probes.loc[last_idx, 'n.probes']
+ means.extend([seg['mean']] * seg['n.probes'])
+ cobalt_ratio_pcf_probes.loc[last_idx, 'mean'] = np.mean(means)
+ cobalt_ratio_pcf_probes.loc[last_idx, 'n.probes'] += seg['n.probes']
+ cobalt_ratio_pcf_probes.loc[last_idx, 'end.pos'] = seg['end.pos']
+ else:
+ cobalt_ratio_pcf_probes = pd.concat([cobalt_ratio_pcf_probes, seg.to_frame().T], ignore_index=True)
+ last_idx = cobalt_ratio_pcf_probes.index[-1]
+
+ # Then bin by logR mean
+ cobalt_ratio_pcf_probes = cobalt_ratio_pcf_probes.reset_index().drop(columns="index")
+ cobalt_ratio_pcf_probes_logR = pd.DataFrame(columns=cobalt_ratio_pcf_probes.columns)
+ chrom_arm = None
+ for idx, seg in cobalt_ratio_pcf_probes.iterrows():
+ if chrom_arm != '_'.join(seg[['chrom','arm']].astype(str)):
+ chrom_arm = '_'.join(seg[['chrom','arm']].astype(str))
+ cobalt_ratio_pcf_probes_logR = pd.concat([cobalt_ratio_pcf_probes_logR, seg.to_frame().T], ignore_index=True)
+ last_idx = cobalt_ratio_pcf_probes_logR.index[-1]
+ continue
+ if abs(cobalt_ratio_pcf_probes.loc[last_idx, 'mean'] - seg['mean']) <= ${binLogR}:
+ means = [cobalt_ratio_pcf_probes_logR.loc[last_idx, 'mean']] * cobalt_ratio_pcf_probes_logR.loc[last_idx, 'n.probes']
+ means.extend([seg['mean']] * seg['n.probes'])
+ cobalt_ratio_pcf_probes_logR.loc[last_idx, 'mean'] = np.mean(means)
+ cobalt_ratio_pcf_probes_logR.loc[last_idx, 'n.probes'] += seg['n.probes']
+ cobalt_ratio_pcf_probes_logR.loc[last_idx, 'end.pos'] = seg['end.pos']
+ else:
+ cobalt_ratio_pcf_probes_logR = pd.concat([cobalt_ratio_pcf_probes_logR, seg.to_frame().T], ignore_index=True)
+ last_idx = cobalt_ratio_pcf_probes_logR.index[-1]
+
+ cobalt_ratio_pcf_probes_logR.to_csv("${tumor}.cobalt.ratio.pcf", sep='\\t', index=False)
  """.stripIndent()
 }
 
 process runPurple {
  tag "PURPLE on ${params.tumor}"
- publishDir "${params.outdir}/purple", mode: 'copy'
+ publishDir "${params.outdir}/purple/purple_${params.minPurity}_${params.maxPurity}", mode: 'copy'
  cpus params.cores
- memory '32 GB'
+ memory params.memory
  time '1h'
 
  input:
@@ -97,6 +197,7 @@ process runPurple {
  path amber_qc
  path cobalt_ratio_tsv
  path cobalt_ratio_pcf
+ path cobalt_path
 
  output:
  path "${tumor}.purple.purity.tsv", emit: purple_purity_tsv
@@ -118,21 +219,34 @@ process runPurple {
  purple \
  -tumor ${tumor} \
  -amber ${params.outdir}/amber \
- -cobalt ${params.outdir}/cobalt \
+ -cobalt ${cobalt_path} \
  -output_dir \$PWD \
  -gc_profile ${params.gcProfile} \
  -ref_genome ${params.refGenome} \
  -ref_genome_version ${params.genomeVersion} \
  -ensembl_data_dir ${params.ensemblDataDir} \
- -circos ${params.circos}
+ -circos ${params.circos} \
+ -min_purity ${params.minPurity} \
+ -max_purity ${params.maxPurity}
+
+ rsync -a --no-links \$PWD/ ${params.outdir}/purple/
  """.stripIndent()
 }
 
 workflow {
  tumor = Channel.value(params.tumor)
  tumorBam = Channel.fromPath(params.tumorBam)
-
- runAmber(tumor, tumorBam)
- runCobalt(tumor, tumorBam)
- runPurple(tumor, runAmber.out, runCobalt.out)
+ binProbes = Channel.value(params.binProbes)
+ binLogR = Channel.value(params.binLogR)
+
+ amberOutput = runAmber(tumor, tumorBam)
+ cobaltOutput = runCobalt(tumor, tumorBam)
+
+ if (binProbes != 0 || binLogR != 0) {
+ binCobaltOutput = binCobalt(tumor, binProbes, binLogR, cobaltOutput.cobalt_ratio_pcf, cobaltOutput.cobalt_ratio_tsv)
+ runPurple(tumor, amberOutput, binCobaltOutput, "${params.outdir}/cobalt/binned_${params.binProbes}_probes_${params.binLogR}_LogR")
+ } else {
+ runPurple(tumor, amberOutput, cobaltOutput, "${params.outdir}/cobalt")
+ }
 }
+

nextflow.config

@@ -48,7 +48,7 @@ profiles {
  }
  process {
  executor = 'local'
- container = 'papaemmelab/purple:v0.1.0'
+ container = 'papaemmelab/purple:v0.1.1'
  }
  }
 }

tests/data/TEST.cobalt.ratio.pcf

@@ -0,0 +1,3 @@
+sampleID chrom arm start.pos end.pos n.probes mean
+S1 2 p 14001 397001 340 -7.24
+S1 17 p 1 397001 341 -8.813389396729216

tests/main.bincobalt.nf.test

@@ -0,0 +1,41 @@
+nextflow_process {
+
+ name "Test Process binCobalt"
+ script "main.nf"
+ process "binCobalt"
+
+ test("Should bin Cobalt without failures") {
+
+ when {
+ params {
+ tumor = "TEST"
+ binProbes = 100
+ binLogR = 0.5
+ cobalt_ratio_pcf = "${projectDir}/tests/data/TEST.cobalt.ratio.pcf"
+ cobalt_ratio_tsv = "${projectDir}/tests/data/TEST.cobalt.ratio.tsv.gz"
+ }
+ process {
+ """
+ input[0] = Channel.value(params.tumor)
+ input[1] = Channel.value(params.binProbes)
+ input[2] = Channel.value(params.binLogR)
+ input[3] = Channel.fromPath(params.cobalt_ratio_pcf)
+ input[4] = Channel.fromPath(params.cobalt_ratio_tsv)
+ """
+ }
+ }
+
+ then {
+ assert process.success
+ assert snapshot(process.out).match()
+
+ // check expected files
+ with(process.out) {
+ assert cobalt_ratio_tsv.size() == 1
+ assert cobalt_ratio_pcf.size() == 1
+ }
+ }
+
+ }
+
+}

tests/main.bincobalt.nf.test.snap

@@ -0,0 +1,25 @@
+{
+ "Should bin Cobalt without failures": {
+ "content": [
+ {
+ "0": [
+ "TEST.cobalt.ratio.tsv.gz:md5,87cf8451b04e4373fc4dd7ccda3e6afe"
+ ],
+ "1": [
+ "TEST.cobalt.ratio.pcf:md5,cd95f6b56cecfd80461bf1bd13c7a9fd"
+ ],
+ "cobalt_ratio_pcf": [
+ "TEST.cobalt.ratio.pcf:md5,cd95f6b56cecfd80461bf1bd13c7a9fd"
+ ],
+ "cobalt_ratio_tsv": [
+ "TEST.cobalt.ratio.tsv.gz:md5,87cf8451b04e4373fc4dd7ccda3e6afe"
+ ]
+ }
+ ],
+ "meta": {
+ "nf-test": "0.9.0",
+ "nextflow": "24.04.4"
+ },
+ "timestamp": "2024-08-07T11:15:21.745138"
+ }
+}

tests/main.nf.test

@@ -16,8 +16,8 @@ nextflow_pipeline {
  with(workflow) {
  assert failed
  assert exitStatus == 1
- assert trace.tasks().size() == 3
- assert workflow.trace.succeeded().size() == 2
+ assert trace.tasks().size() == 4
+ assert workflow.trace.succeeded().size() == 3
  assert workflow.trace.failed().size() == 1
  }
  }