Develop_atac

conf/atac/preprocess.config

@@ -1,4 +1,8 @@
 params {
+    atac_preprocess_tools {
+        mark_duplicates_method = 'MarkDuplicates'
+        adapter_trimming_method = 'Trim_Galore'
+    }
     data {
         atac_preprocess {
             metadata = 'metadata.tsv'
@@ -7,10 +11,13 @@ params {
 }
 
 includeConfig './../../src/singlecelltoolkit/singlecelltoolkit.config'
+includeConfig './../../src/singlecelltoolkit/conf/sctk_mapping.config'
 includeConfig './../../src/trimgalore/trimgalore.config'
+includeConfig './../../src/fastp/fastp.config'
 includeConfig './../../src/bwamaptools/bwamaptools.config'
+includeConfig './../../src/gatk/gatk.config'
 includeConfig './../../src/bwamaptools/conf/bwa_mapping.config'
 includeConfig './../../src/sinto/sinto.config'
 includeConfig './../../src/bap/bap.config'
-includeConfig './../../src/bap/conf/bap_biorad_debarcode.config'
+//includeConfig './../../src/bap/conf/bap_biorad_debarcode.config'
 

conf/atac/qc_filtering.config

@@ -1,3 +1,4 @@
-includeConfig './../../src/archr/archr.config'
 includeConfig './../../src/pycistopic/pycistopic.config'
+includeConfig './../../src/singlecelltoolkit/singlecelltoolkit.config'
+includeConfig './../../src/singlecelltoolkit/conf/sctk_saturation.config'
 

docs/index.rst

@@ -20,6 +20,7 @@
    :caption: scATAC-seq
 
    scatac-seq
+   scatac-seq_qc
 
 .. toctree::
    :maxdepth: 2

docs/scatac-seq_qc.rst

@@ -0,0 +1,212 @@
+scATAC-seq QC and Cell Calling
+==============================
+
+This workflow uses the Python implementation of `cisTopic <https://github.com/aertslab/cisTopic>`_ (pycisTopic) to perform quality control and cell calling.
+The inputs here are a fragments and bam file for each sample.
+
+This workflow is currently available in the ``develop_atac`` branch (use ``nextflow pull vib-singlecell-nf/vsn-pipelines -r develop_atac`` to sync this branch).
+
+----
+
+
+Running the workflow
+********************
+
+Technical considerations
+------------------------
+
+1. Direct the Nextflow work directory to an alternate path (e.g. a scratch drive) using the ``NXF_WORK`` environmental variable::
+
+    nwork=/path/to/scratch/example_project
+    mkdir $nwork
+    export NXF_WORK=$nwork
+
+Note that if you start a new shell, ``NXF_WORK`` must be set again, or the pipeline will not resume properly.
+
+
+2. Important for pycisTopic Ray issues: the system default temp location may become full.
+   A workaround is to include a volume mapping to the alternate ``/tmp`` ``-B /alternate/path/to/tmp:/tmp`` using the volume mount options in Docker or Singularity.
+   For example in the container engine options:
+  - Singularity run options: ``runOptions = '--cleanenv -H $PWD -B /data,/alternate/path/to/tmp:/tmp'``
+  - Docker run options: ``runOptions = '-i -v /data:/data -v /alternate/path/to/tmp:/tmp'``
+
+
+3. Use the ``--quiet`` flag with ``nextflow run`` to suppress the printing of each file that is detected by the pipeline.
+
+----
+
+Configuration
+-------------
+
+For each sample, this pipeline take a bam and a fragments file.
+These can be specified separately, or from a Cell Ranger ATAC/ARC ``outs/`` path.
+
+Input with independent bam and fragments files
+______________________________________________
+
+Use the profiles ``bam`` and ``fragments``::
+
+    nextflow config vib-singlecell-nf/vsn-pipelines/main_atac.nf \
+        -profile atac_qc_filtering,bam,fragments,vsc,pycistopic_hg38 \
+        > atac_qc.config
+
+Preset profiles are available for human (``pycistopic_hg38``), mouse (``pycistopic_mm10``), and fly (``pycistopic_dmel``).
+Or, these profiles can be omitted and set manually in the config (biomart, macs2).
+
+
+Input data (bam and fragments files) are specified in the data section::
+
+    data {
+       fragments {
+          file_paths = '/staging/leuven/stg_00002/lcb/cflerin/analysis/asap/20210527_hydrop-atac_asabr/atac_preprocess/out_run1/data/fragments/ASA__*tsv.gz'
+          suffix = '.sinto.fragments.tsv.gz'
+          index_extension = '.tbi'
+       }
+       bam {
+          file_paths = '/staging/leuven/stg_00002/lcb/cflerin/analysis/asap/20210527_hydrop-atac_asabr/atac_preprocess/out_run1/data/bam/ASA*bam'
+          suffix = '.bwa.out.possorted.bam'
+          index_extension = '.bai'
+       }
+    }
+
+
+Multiple files can be specified with ``*`` in ``file_paths`` or by separating the paths with a comma.
+
+.. warning::
+
+    The ``suffix`` for both bam and fragments will be removed from the filename to get sample IDs.
+    The sample names obtained must match between bam and fragments for the files to be paired properly in the workflow.
+
+
+Input with Cell Ranger ATAC data
+________________________________
+
+Use the ``tenx_atac`` profile::
+
+    nextflow config vib-singlecell-nf/vsn-pipelines/main_atac.nf \
+        -profile atac_qc_filtering,tenx_atac,vsc,pycistopic_hg38 \
+        > atac_qc.config
+
+Input data (the Cell Ranger ``outs/`` path) are specified in the data section::
+
+    data {
+        tenx_atac {
+            cellranger_mex = '/data/cellranger_atac_2.0/*/outs,/data/processed/cellranger_arc_2.0.0/*/outs'
+        }
+    }
+
+Multiple files can be specified with ``*`` in ``tenx_atac`` or by separating the paths with a comma.
+
+
+Input directly from the preprocessing pipeline
+______________________________________________
+
+It is also possible to run these QC steps directly after the ``atac_preprocess`` pipeline, with a single command.
+In this case, all the appropriate configuration profiles must be included at the configuration start::
+
+    nextflow config vib-singlecell-nf/vsn-pipelines/main_atac.nf \
+        -profile atac_preprocess,atac_qc_filtering,pycistopic_hg38,vsc \
+        > atac_preprocess_and_qc.config
+
+Note that here, we do not include ``bam`` and ``fragments`` profiles that specify the input data locations to the QC steps since these are piped directly from the preprocessing pipeline.
+One caveat to this is that it could potentially make it harder to run the qc pipeline with ``-resume`` later on, especially if the Nextflow ``work/`` directory is not saved due to disk space concerns.
+
+To execute the preprocessing and mapping pipeline in one step, use the ``atac_preprocess_with_qc`` entry point::
+
+    nextflow -C atac_preprocess_and_qc.config run \
+        vib-singlecell-nf/vsn-pipelines/main_atac.nf \
+        -entry atac_preprocess_with_qc -resume --quiet
+
+
+----
+
+Execution
+---------
+
+After configuring, the workflow can be run with:
+
+.. code:: bash
+
+    nextflow -C atac_qc.config run \
+        vib-singlecell-nf/vsn-pipelines/main_atac.nf \
+        -entry atac_qc_filtering --quiet -resume
+
+After completing, view the report in ``out/notebooks/<project_name>__pycisTopic_QC_report.html``. To change the filtering settings, use the ``params.tools.pycistopic.call_cells`` section.
+
+Adjusting the filter settings
+-----------------------------
+
+In the pycisTopic parameters, filter settings can be applied in this section::
+
+    pycistopic {
+        call_cells {
+            report_ipynb = '/src/pycistopic/bin/pycisTopic_qc_report_template.ipynb'
+            use_density_coloring_on_scatterplot = true
+            use_detailed_title_on_scatterplot = true
+            filter_frags_lower = '1000'
+            filter_frags_upper = ''
+            filter_tss_lower = '8'
+            filter_tss_upper = ''
+            filter_frip_lower = ''
+            filter_frip_upper = ''
+            filter_dup_rate_lower = ''
+            filter_dup_rate_upper = ''
+        }
+    }
+
+If a setting is empty (``''``), this filter will not be applied.
+If set to a single value (i.e. ``filter_frags_lower=1000``), this will apply this filter value to all samples.
+To use sample-specific filters, this can be written as::
+
+    filter_frags_lower = [
+      'default': 1000,
+      'Sample_1': 1500,
+      'Sample_2': 2000,
+    ]
+
+The ``default`` setting (optional) is applied to all samples not listed in array.
+If this default setting is missing, no filter will be applied to samples not listed in the array (all barcodes kept).
+
+After setting the filters, the pipeline can be re-run to apply the new filters (use ``-resume``).
+
+The additional settings control the output of the scatter plots in the report:
+* ``use_density_coloring_on_scatterplot``: Slower when turned on; it can be helpful to set this to ``false`` until the proper thresholds are determined.
+* ``use_detailed_title_on_scatterplot``: Adds the cell count and median values after filtering to the title of each plot.
+
+----
+
+Output
+******
+
+An example output tree is shown here.
+
+.. code:: bash
+
+    out/
+    ├── data
+    │   ├── macs2
+    │   │   ├── sample_1.peaks.narrowPeak
+    │   │   ├── sample_1.summits.bed
+    │   │   ├── sample_2.peaks.narrowPeak
+    │   │   └── sample_2.summits.bed
+    │   └── pycistopic
+    │       └── qc
+    │           ├── benchmark_library_downsampled__metadata.pickle
+    │           ├── benchmark_library_downsampled__profile_data.pickle
+    │           ├── selected_barcodes
+    │           │   ├── sample_1.cell_barcodes.txt
+    │           │   └── sample_2.cell_barcodes.txt
+    │           └── selected_barcodes_nFrag
+    │               ├── sample_1.barcodes_nFrag_thr.txt
+    │               └── sample_2.barcodes_nFrag_thr.txt
+    └── notebooks
+        ├── example_project__pycisTopic_QC_report.html
+        └── example_project__pycisTopic_QC_report.ipynb
+
+
+* ``macs2``: contains the narrowPeak and bed file for each sample.
+* ``pycistopic``:
+  * ``qc``: contains Python objects (in pickle format) for the metadata and profile data computed by pycisTopic.
+    * ``selected_barcodes``: contains a text file with selected cell barcodes (one per line) based on the thresholds set in the config file.
+    * ``selected_barcodes_nFrag``: contains a text file with barcodes (one per line) that have unique fragment counts greater than the ``params.tools.pycistopic.compute_qc_stats.n_frag`` setting in the pycisTopic parameters.
+

main_atac.nf

@@ -37,53 +37,42 @@ workflow cellranger_atac {
 }
 
 
-workflow cistopic {
-
-    include {
-        cistopic as CISTOPIC
-    } from './src/cistopic/main' params(params)
-
-    getDataChannel | CISTOPIC
-
-}
-
-
 workflow atac_preprocess {
 
     // generic ATAC-seq preprocessing pipeline: adapter trimming, mapping, fragments file generation
     include {
-        ATAC_PREPROCESS_WITH_METADATA;
+        ATAC_PREPROCESS;
     } from './workflows/atac/preprocess.nf' params(params)
 
-    ATAC_PREPROCESS_WITH_METADATA(file(params.data.atac_preprocess.metadata))
+    ATAC_PREPROCESS(file(params.data.atac_preprocess.metadata))
 
 }
 
 
 workflow atac_preprocess_bap {
 
     include {
-        ATAC_PREPROCESS_WITH_METADATA;
+        ATAC_PREPROCESS;
     } from './workflows/atac/preprocess.nf' params(params)
     include {
         get_bam;
-        BAP__BARCODE_MULTIPLET_PIPELINE;
+        BAP__BARCODE_MULTIPLET_WF;
     } from './src/bap/main.nf' params(params)
 
-    ATAC_PREPROCESS_WITH_METADATA(file(params.data.atac_preprocess.metadata)) |
+    ATAC_PREPROCESS(file(params.data.atac_preprocess.metadata)) |
         get_bam |
-        BAP__BARCODE_MULTIPLET_PIPELINE
+        BAP__BARCODE_MULTIPLET_WF
 
 }
 
 
 workflow bap {
 
     include {
-        BAP__BARCODE_MULTIPLET_PIPELINE;
+        BAP__BARCODE_MULTIPLET_WF;
     } from './src/bap/main.nf' params(params)
 
-    getDataChannel | BAP__BARCODE_MULTIPLET_PIPELINE
+    getDataChannel | BAP__BARCODE_MULTIPLET_WF
 
 }
 
@@ -101,6 +90,21 @@ workflow atac_qc_filtering {
 
 }
 
+workflow atac_preprocess_with_qc {
+
+    // generic ATAC-seq preprocessing pipeline: adapter trimming, mapping, fragments file generation
+    include {
+        ATAC_PREPROCESS;
+    } from './workflows/atac/preprocess.nf' params(params)
+    include {
+        ATAC_QC_PREFILTER;
+    } from './workflows/atac/qc_filtering.nf' params(params)
+
+    pp = ATAC_PREPROCESS(file(params.data.atac_preprocess.metadata))
+    ATAC_QC_PREFILTER(pp.bam.mix(pp.fragments))
+
+}
+
 workflow atac_preprocess_freemuxlet {
 
     // generic ATAC-seq preprocessing pipeline: adapter trimming, mapping, fragments file generation

nextflow.config

@@ -186,6 +186,7 @@ profiles {
         includeConfig 'conf/nemesh.config'
         includeConfig 'src/fastp/fastp.config'
         includeConfig 'src/picard/picard.config'
+        includeConfig 'src/picard/conf/nemesh.config'
         includeConfig 'src/dropseqtools/dropseqtools.config'
         includeConfig 'src/star/star.config'
         includeConfig 'src/dropletutils/dropletutils.config'
@@ -465,9 +466,6 @@ profiles {
         includeConfig 'src/cellranger-atac/conf/count.config'
         includeConfig 'src/cellranger-atac/conf/count_metadata.config'
     }
-    cistopic {
-        includeConfig 'src/cistopic/cistopic.config'
-    }
     atac_preprocess {
         includeConfig 'conf/atac/preprocess.config'
     }
@@ -482,6 +480,10 @@ profiles {
         includeConfig 'conf/atac/preprocess.config'
         includeConfig 'src/popscle/popscle.config'
     }
+    /* pycistopic (biomart/macs2) profiles (load after pycistopic config) */
+    pycistopic_hg38 { includeConfig 'src/pycistopic/conf/pycistopic_hg38.config' }
+    pycistopic_mm10 { includeConfig 'src/pycistopic/conf/pycistopic_mm10.config' }
+    pycistopic_dmel { includeConfig 'src/pycistopic/conf/pycistopic_dmel.config' }
 
 
     /*

src/bap/conf/bap_barcode_multiplet.config

@@ -2,7 +2,7 @@ params {
     tools {
         bap {
             barcode_multiplet {
-                report_ipynb = './src/bap/bin/bap_barcode_multiplet_report.ipynb'
+                report_ipynb = '/src/bap/bin/bap_barcode_multiplet_report.ipynb'
                 bead_tag = 'CB' // existing tag in bam file
                 drop_tag = 'DB' // new tag with merged barcodes
                 // filters:

src/bap/main.nf

@@ -4,7 +4,7 @@ nextflow.enable.dsl=2
 //  Import sub-workflows from the modules:
 
 include {
-    SC__BAP__BARCODE_MULTIPLET_PIPELINE;
+    BAP__BARCODE_MULTIPLET_PIPELINE as BARCODE_MULTIPLET_PIPELINE;
 } from './processes/barcode_multiplet.nf' params(params)
 
 include {
@@ -25,15 +25,14 @@ workflow get_bam {
         bam
 }
 
-
-workflow BAP__BARCODE_MULTIPLET_PIPELINE {
+workflow BAP__BARCODE_MULTIPLET_WF {
 
     take:
         bam
 
     main:
 
-        bap = SC__BAP__BARCODE_MULTIPLET_PIPELINE(bam.map { it -> tuple(it[0], it[1], it[2]) })
+        bap = BARCODE_MULTIPLET_PIPELINE(bam.map { it -> tuple(it[0], it[1][0], it[1][1]) })
 
         GENERATE_REPORT(
             file(workflow.projectDir + params.tools.bap.barcode_multiplet.report_ipynb),