snRNA-seq Preprocessing - FunGen-xQTL Consortium

Takes CellRanger counts through quality control and cell-type annotation to a Seurat object ready for downstream analysis.

Overview¶

Single-nuclei RNA-seq counts arrive with artefacts that would otherwise be read as biology: dying cells with high mitochondrial content, ambient RNA carried over from the suspension, and droplets holding two nuclei rather than one. Filtering those out, and then deciding what cell type each surviving nucleus is, has to happen before any per-cell-type analysis can start.

Quality control runs through SCTK and Seurat: cells are dropped on mitochondrial percent, total counts (nUMI) and detected genes (nFeature), ambient RNA is removed with decontX, and doublets are removed with a user-selected method, scds by default. Cell types are then assigned by transferring labels from an annotated reference dataset onto the filtered object.

When to run it. On CellRanger output, before pseudobulk aggregation or any per-cell-type phenotype work.

Input¶

--input-dir: the CellRanger output directory, holding the per-sample counts folder with filtered_feature_bc_matrix. Example input/snrnaseq/protocol_example.snrnaseq.cellranger/.
--sample-meta: a CSV mapping cell barcodes to individuals, with columns libraryBatch, cellBarcode and individualID. Example tests/fixtures/snrnaseq_preprocessing/protocol_example.snrnaseq.id_mapping.csv:
```
libraryBatch,cellBarcode,individualID
190403-B4-A,ATTATCCCAGGACATG-1,R2670295
190403-B4-A,GAGCTGCCACACGGAA-1,R1743384
```
--sctk-rds: for cell_annotation, the QC-filtered Seurat object written by sctk_qc.
--seurat-ref: the annotated reference dataset the cell types are transferred from. Example tests/fixtures/snrnaseq_preprocessing/protocol_example.snrnaseq.seurat_ref_SE.rds.
--output-dir: the directory outputs are written to.
--name: the prefix given to every output filename, protocol_example.snrnaseq by default, so results stay traceable to the dataset.

Output¶

{output_dir}/SCTK_results/{name}.filtered_seuratobj.rds - the QC-filtered, clustered Seurat object from sctk_qc, and the input to cell_annotation. No example ships: output/snrna_seq/ has no SCTK_results/ directory.
{output_dir}/QC_table/{name}.SCTK_QC_table.csv - per-sample SCTK QC metrics. Example output/snrna_seq/QC_table/protocol_example.snrnaseq.SCTK_QC_table.csv, 10 columns and 10 rows:
```
"","R6249698","R7840275","R9818080","R8054504"
"Number of Cells,361,269,246,151
"Mean counts",10237,10117,10514,8989.2
```
{output_dir}/SCTK_results/{name}.QC_summary.csv - cell and gene counts at each QC step. No example ships.
{output_dir}/{name}.celltyped_seuratobj.rds - the cell-type annotated Seurat object from cell_annotation. Example output/snrna_seq/protocol_example.snrnaseq.celltyped_seuratobj.rds.

Minimal Working Example¶

Two workflows run in order: sctk_qc filters the CellRanger counts, then cell_annotation labels the surviving cells. They are chained by hand - --sctk-rds in the second command is the object the first one wrote.

QC filtering¶

sctk_qc drops cells on the mitochondrial-percent, minimum-count and minimum-gene thresholds, removes ambient RNA with decontX, and removes doublets with the method --doublet-method selects.

Timing: TBD (on toy dataset)

sos run pipeline/snRNAseq_preprocessing.ipynb sctk_qc \
    --input-dir <path/to/protocol_example.snrnaseq.cellranger> \
    --output-dir output/snrna_seq \
    --sample-meta tests/fixtures/snrnaseq_preprocessing/protocol_example.snrnaseq.id_mapping.csv \
    --name protocol_example.snrnaseq

Cell type annotation¶

Timing: TBD (on toy dataset)

sos run pipeline/snRNAseq_preprocessing.ipynb cell_annotation \
    --sctk-rds output/snrna_seq/SCTK_results/protocol_example.snrnaseq.filtered_seuratobj.rds \
    --output-dir output/snrna_seq \
    --seurat-ref tests/fixtures/snrnaseq_preprocessing/protocol_example.snrnaseq.seurat_ref_SE.rds \
    --name protocol_example.snrnaseq

Command Interface¶

sos run pipeline/snRNAseq_preprocessing.ipynb -h

usage: sos run code/SoS/molecular_phenotypes/snRNAseq_preprocessing.ipynb
               [workflow_name | -t targets] [options] [workflow_options]
  workflow_name:        Single or combined workflows defined in this script
  targets:              One or more targets to generate
  options:              Single-hyphen sos parameters (see "sos run -h" for details)
  workflow_options:     Double-hyphen workflow-specific parameters

Workflows:
  sctk_qc
  cell_annotation

Global Workflow Options:
  --modular-script-dir code/script (as path)
  --cwd output (as path)
  --name 'protocol_example.snrnaseq'
                        Prefix for all output filenames; keeps results traceable
                        to the dataset
  --job-size 1 (as int)
  --walltime 5h
  --mem 16G
  --numThreads 8 (as int)
  --container ''

Sections
  sctk_qc:
    Workflow Options:
      --input-dir VAL (as str, required)
      --output-dir VAL (as str, required)
      --sample-meta VAL (as str, required)
      --max-mito-percent 5 (as int)
      --min-ncounts 250 (as int)
      --min-ngenes 500 (as int)
      --doublet-method scds
      --scds-call Singlet
      --scrublet-threshold 0.25 (as float)
      --ambient-method decontX
      --max-decontXsc 0.5 (as float)
  cell_annotation:
    Workflow Options:
      --sctk-rds VAL (as str, required)
      --output-dir VAL (as str, required)
      --seurat-ref VAL (as str, required)
      --downsample-n 1000 (as int)
      --cluster-threshold 0.9 (as float)

Workflow implementation¶

[global]
parameter: modular_script_dir = path('code/script')  # override with --modular-script-dir
parameter: cwd = path("output")
# Prefix for all output filenames; keeps results traceable to the dataset
parameter: name = "protocol_example.snrnaseq"
parameter: job_size = 1
parameter: walltime = "5h"
parameter: mem = "16G"
parameter: numThreads = 8
parameter: container = ""

cwd = path(f'{cwd:a}')

[sctk_qc]
parameter: input_dir          = str
parameter: output_dir         = str
parameter: sample_meta        = str
parameter: max_mito_percent   = 5
parameter: min_ncounts        = 250
parameter: min_ngenes         = 500
parameter: doublet_method     = 'scds'      # options: 'scds', 'scrublet', 'doubletFinder', 'doubletCells', 'none'
parameter: scds_call          = 'Singlet'
parameter: scrublet_threshold = 0.25
parameter: ambient_method     = 'decontX'   # options: 'decontX', 'none'
parameter: max_decontXsc      = 0.5

input:  []
output: f'{output_dir}/SCTK_results/{name}.filtered_seuratobj.rds',
        f'{output_dir}/QC_table/{name}.SCTK_QC_table.csv'
task: trunk_workers = 1, trunk_size = 1, walltime = '36:00:00', mem = '256G', cores = 16
bash: expand = "${ }", stdout = f'{_output[0]:n}.stdout', stderr = f'{_output[0]:n}.stderr'
    Rscript ${modular_script_dir}/molecular_phenotypes/snRNAseq_preprocessing.R \
        --step sctk_qc \
        --input-dir "${input_dir}" \
        --output-dir "${output_dir}" \
        --name "${name}" \
        --sample-meta "${sample_meta}" \
        --max-mito-percent ${max_mito_percent} \
        --min-ncounts ${min_ncounts} \
        --min-ngenes ${min_ngenes} \
        --doublet-method "${doublet_method}" \
        --scds-call "${scds_call}" \
        --scrublet-threshold ${scrublet_threshold} \
        --ambient-method "${ambient_method}" \
        --max-decontxsc ${max_decontXsc} \
        --num-threads ${numThreads}

[cell_annotation]
parameter: sctk_rds          = str
parameter: output_dir        = str
parameter: seurat_ref        = str
parameter: downsample_n      = 1000
parameter: cluster_threshold = 0.90

input:  []
output: f'{output_dir}/{name}.celltyped_seuratobj.rds'
task: trunk_workers = 1, trunk_size = 1, walltime = '36:00:00', mem = '256G', cores = 16
bash: expand = "${ }", stdout = f'{_output:n}.stdout', stderr = f'{_output:n}.stderr'
    Rscript ${modular_script_dir}/molecular_phenotypes/snRNAseq_preprocessing.R \
        --step cell_annotation \
        --sctk-rds "${sctk_rds}" \
        --output-dir "${output_dir}" \
        --name "${name}" \
        --seurat-ref "${seurat_ref}" \
        --downsample-n ${downsample_n} \
        --cluster-threshold ${cluster_threshold} \
        --num-threads ${numThreads}