Merge pull request #70 from Mayrlab/dev

Version 0.4.0

Author: mfansler

README.md

@@ -69,7 +69,7 @@ repository directly tests running examples on the GitHub-hosted runners.
 ### Conda/Mamba Mode (MacOS or Linux)
 Snakemake can use Conda to install the needed software. This configuration requires:
 
- - [Snakemake][ref:snakemake] >= 5.11<sup>ª</sup>
+ - [Snakemake][ref:snakemake] >= 6.0, < 8.0<sup>ª</sup>
  - [Conda](https://docs.conda.io/projects/conda/en/latest/)
 
 If Conda is not already installed, we strongly recommend installing 
@@ -84,7 +84,7 @@ conda install -n base -c conda-forge mamba
 Snakemake can use the pre-built scUTRsquant Docker image to provide all additional software.
 This configuration requires installing:
 
- - [Snakemake][ref:snakemake] >= 5.11<sup>ª</sup>
+ - [Snakemake][ref:snakemake] >= 6.0, < 8.0<sup>ª</sup>
  - [Singularity](https://singularity.lbl.gov/index.html)
 
 
@@ -96,7 +96,13 @@ This configuration requires installing:
     git clone git@github.com:Mayrlab/scUTRquant.git
     ```
 
-2. Download the UTRome annotation, kallisto index, and merge file.
+As of scUTRquant v0.4.0, that is all!
+
+<details>
+  
+> *The following steps of prepopulating the annotations and barcode whitelists are now handled directly in the pipeline.*
+  
+2. (**Deprecated**) Download the UTRome annotation, kallisto index, and merge file.
     **Human (hg38)**
     ```
     cd scUTRquant/extdata/targets/utrome_hg38_v1
@@ -112,7 +118,7 @@ This configuration requires installing:
     `utrome_kdx`, and `utrome_merge` to point to the central location. In that
     case, one does not need to redownload the files.
 
-3. (Optional) Download the barcode whitelists.
+3. (**Deprecated** and Optional) Download the barcode whitelists.
     ```
     cd scUTRquant/extdata/bxs
     sh download_10X_whitelists.sh
@@ -121,6 +127,7 @@ This configuration requires installing:
     and referenced by the `bx_whitelist` variable in the `configfile`.
 
 For GitHub runners, it takes ~ 3 mins to clone and download the scUTRquant files.
+</details>
 
 # Running Examples
 Examples are provided in the `scUTRquant/examples` folder. Each includes a script
@@ -219,6 +226,7 @@ pipeline. The following keys are expected:
  - `min_umis`: minimum number of UMIs per cell; cells below this threshold are excluded
  - `cell_annots`: (optional) CSV file with a key column that matches the `<sample_id>_<cell_bx>` format
  - `cell_annots_key`: specifies the name of the key column in the `cell_annots` file; default is `cell_id`
+ - `exclude_unannotated_cells`: boolean indicating whether unannotated cells should be excluded from the final output; default is `False`
  
 ### Default Values
 

Snakefile

@@ -1,11 +1,15 @@
 container: "docker://mfansler/scutr-quant:0.1.5"
 configfile: "config.yaml"
 
-from snakemake.io import load_configfile
-import pandas as pd
 import os
+import pandas as pd
+from snakemake.io import load_configfile
+from snakemake.utils import min_version
 from sys import stderr
 
+# set minimum Snakemake version
+min_version("6.0")
+
 # print to stderr
 def message(*args, **kwargs):
     print(*args, file=stderr, **kwargs)
@@ -66,6 +70,66 @@ def get_outputs():
 rule all:
     input: get_outputs()
 
+################################################################################
+## Downloading and Preprocessing
+################################################################################
+
+## generate downloading rules for target files if a script is provided
+## NB: script should generate files *relative* to the script
+for target_id, target in targets.items():
+    if 'download_script' not in target or target['download_script'] is None:
+        continue
+
+    target_path = target['path']
+    download_script = target_path + target['download_script']
+    
+    FILE_KEYS = ['gtf', 'kdx', 'merge_tsv', 'tx_annots', 'gene_annots']
+    target_files = [target_path + target[k] for k in FILE_KEYS]
+
+    rule:
+        name: f"download_{target_id}"
+        input: f"{download_script}"
+        output: expand("{target_file}", target_file=target_files)
+        conda: "envs/downloading.yaml"
+        shell:
+            """
+            pushd $(dirname {input})
+            sh $(basename {input})
+            popd
+            """
+
+## Import downloading rules for barcodes
+module bxs_workflow:
+    snakefile: "extdata/bxs/download.smk"
+
+use rule * from bxs_workflow
+
+## Convert merge data for tx output
+rule generate_tx_merge:
+    input:
+        tsv=get_target_file('merge_tsv')
+    output:
+        "data/utrs/{target}/tx_merge.tsv"
+    shell:
+        """
+        tail -n+2 {input.tsv} | cut -f1,2 > {output}
+        """
+
+## Convert merge data for gene output
+rule generate_gene_merge:
+    input:
+        tsv=get_target_file('merge_tsv')
+    output:
+        "data/utrs/{target}/gene_merge.tsv"
+    shell:
+        """
+        tail -n+2 {input.tsv} | cut -f1,3 > {output}
+        """
+
+################################################################################
+## kallisto-bustools
+################################################################################
+
 def get_file_type(wildcards):
     return "--bam" if samples.file_type[wildcards.sample_id] == 'bam' else ""
 
@@ -157,26 +221,6 @@ rule bustools_correct_sort:
         bustools sort -t{threads} -T {params.tmpDir} -o {output} {input}
         """
 
-rule generate_tx_merge:
-    input:
-        tsv=get_target_file('merge_tsv')
-    output:
-        "data/utrs/{target}/tx_merge.tsv"
-    shell:
-        """
-        tail -n+2 {input.tsv} | cut -f1,2 > {output}
-        """
-
-rule generate_gene_merge:
-    input:
-        tsv=get_target_file('merge_tsv')
-    output:
-        "data/utrs/{target}/gene_merge.tsv"
-    shell:
-        """
-        tail -n+2 {input.tsv} | cut -f1,3 > {output}
-        """
-
 def get_input_busfile(wildcards):
     if config['correct_bus']:
         return "data/kallisto/%s/%s/output.corrected.sorted.bus" % (wildcards.target, wildcards.sample_id)
@@ -219,18 +263,9 @@ rule bustools_count_genes:
         bustools count -e {input.ec} -t {input.txs} -g {input.merge} -o $prefix --em --genecounts {input.bus}
         """
 
-rule report_umis_per_cell:
-    input:
-        bxs="data/kallisto/{target}/{sample_id}/txs.barcodes.txt",
-        txs="data/kallisto/{target}/{sample_id}/txs.genes.txt",
-        mtx="data/kallisto/{target}/{sample_id}/txs.mtx"
-    params:
-        min_umis=config['min_umis']
-    output:
-        "qc/umi_count/{target}/{sample_id}.umi_count.html"
-    conda: "envs/rmd-reporting.yaml"
-    script:
-        "scripts/report_umi_counts_per_cell.Rmd"
+################################################################################
+## Outputs
+################################################################################
 
 rule mtxs_to_sce_txs:
     input:
@@ -250,6 +285,7 @@ rule mtxs_to_sce_txs:
         sample_ids=samples.index.values,
         min_umis=config['min_umis'],
         cell_annots_key=config['cell_annots_key'],
+        exclude_unannotated_cells=config['exclude_unannotated_cells'],
         tmp_dir=config['tmp_dir'],
         use_hdf5=False
     resources:
@@ -276,6 +312,7 @@ rule mtxs_to_sce_genes:
         sample_ids=samples.index.values,
         min_umis=config['min_umis'],
         cell_annots_key=config['cell_annots_key'],
+        exclude_unannotated_cells=config['exclude_unannotated_cells'],
         tmp_dir=config['tmp_dir'],
         use_hdf5=False
     resources:
@@ -303,6 +340,7 @@ rule mtxs_to_sce_h5_txs:
         sample_ids=samples.index.values,
         min_umis=config['min_umis'],
         cell_annots_key=config['cell_annots_key'],
+        exclude_unannotated_cells=config['exclude_unannotated_cells'],
         tmp_dir=lambda wcs: config['tmp_dir'] + "/sce-txs-" + wcs.target + "-" + config['dataset_name'],
         use_hdf5=True
     resources:
@@ -331,6 +369,7 @@ rule mtxs_to_sce_h5_genes:
         sample_ids=samples.index.values,
         min_umis=config['min_umis'],
         cell_annots_key=config['cell_annots_key'],
+        exclude_unannotated_cells=config['exclude_unannotated_cells'],
         tmp_dir=lambda wcs: config['tmp_dir'] + "/sce-genes-" + wcs.target + "-" + config['dataset_name'],
         use_hdf5=True
     resources:
@@ -339,3 +378,22 @@ rule mtxs_to_sce_h5_genes:
     conda: "envs/bioconductor-sce.yaml"
     script:
         "scripts/mtxs_to_sce_genes.R"
+
+
+################################################################################
+## Reports
+################################################################################
+
+rule report_umis_per_cell:
+    input:
+        bxs="data/kallisto/{target}/{sample_id}/txs.barcodes.txt",
+        txs="data/kallisto/{target}/{sample_id}/txs.genes.txt",
+        mtx="data/kallisto/{target}/{sample_id}/txs.mtx"
+    params:
+        min_umis=config['min_umis']
+    output:
+        "qc/umi_count/{target}/{sample_id}.umi_count.html"
+    conda: "envs/rmd-reporting.yaml"
+    script:
+        "scripts/report_umi_counts_per_cell.Rmd"
+

config.yaml

@@ -2,6 +2,7 @@ dataset_name: "output"
 sample_regex: ".*"
 cell_annots: null
 cell_annots_key: "cell_id"
+exclude_unannotated_cells: False
 target: "utrome_mm10_v2"
 output_type:
   - "txs"

envs/downloading.yaml

@@ -0,0 +1,8 @@
+name: downloading
+channels:
+  - conda-forge
+  - nodefaults
+dependencies:
+  - gzip
+  - tar
+  - wget

extdata/bxs/download.smk

@@ -0,0 +1,22 @@
+## Download Rules
+##
+## Any rules added here will automatically be imported into the main Snakefile.
+## Important notes:
+##   - all `input`/`output` must be relative to main Snakefile
+##   - `conda` is relative to this file
+##   - `shell` will execute relative to main Snakefile
+##   - `download_10X_whitelists` shows how to run script in-place (see shell)
+
+rule download_10X_whitelists:
+    input: "extdata/bxs/download_10X_whitelists.sh"
+    output:
+        "extdata/bxs/737K-april-2014_rc.txt",
+        "extdata/bxs/737K-august-2016.txt",
+        "extdata/bxs/3M-february-2018.txt"
+    conda: "../../envs/downloading.yaml"
+    shell:
+        """
+        pushd $(dirname {input})
+        sh $(basename {input})
+        popd
+        """

extdata/bxs/download_10X_whitelists.sh

@@ -1,12 +1,13 @@
 #!/usr/bin/env bash
 
 # 10xv1
-wget https://github.com/10XGenomics/cellranger/raw/3.0.2/lib/python/cellranger/barcodes/737K-april-2014_rc.txt
+wget https://teichlab.github.io/scg_lib_structs/data/737K-april-2014_rc.txt.gz \
+     && gzip -d 737K-april-2014_rc.txt.gz
 
 # 10xv2
-wget https://github.com/10XGenomics/cellranger/raw/3.0.2/lib/python/cellranger/barcodes/737K-august-2016.txt
+wget https://teichlab.github.io/scg_lib_structs/data/737K-august-2016.txt.gz \
+     && gzip -d 737K-august-2016.txt.gz
 
 # 10xv3
-wget https://github.com/10XGenomics/cellranger/raw/3.0.2/lib/python/cellranger/barcodes/3M-february-2018.txt.gz \
+wget https://teichlab.github.io/scg_lib_structs/data/3M-february-2018.txt.gz \
     && gzip -d 3M-february-2018.txt.gz
-

extdata/targets/targets.yaml

@@ -6,6 +6,7 @@ utrome_mm10_v1:
   merge_tsv: "adult.utrome.e3.t200.f0.999.w500.merge.tsv"
   tx_annots: "utrome_txs_annotation.Rds"
   gene_annots: "utrome_genes_annotation.Rds"
+  download_script: "download_utrome.sh"
 
 utrome_mm10_v2:
   path: "extdata/targets/utrome_mm10_v2/"
@@ -15,6 +16,7 @@ utrome_mm10_v2:
   merge_tsv: "utrome.e30.t5.gc25.pas3.f0.9999.w500.m200.tsv"
   tx_annots: "utrome_txs_annotation.Rds"
   gene_annots: "utrome_genes_annotation.Rds"
+  download_script: "download_utrome.sh"
 
 utrome_hg38_v1:
   path: "extdata/targets/utrome_hg38_v1/"
@@ -24,6 +26,7 @@ utrome_hg38_v1:
   merge_tsv: "utrome.e30.t5.gc39.pas3.f0.9999.w500.m200.tsv"
   tx_annots: "utrome_txs_annotation.Rds"
   gene_annots: "utrome_genes_annotation.Rds"
+  download_script: "download_utrome.sh"
 
 ## example custom target using GENCODE
 ## see https://github.com/Mayrlab/txcutr-db