Single-nuclei Pseudobulk Preprocessing (RNA-seq and ATAC-seq)

Single-nuclei Pseudobulk Preprocessing (RNA-seq and ATAC-seq)#

Aggregates per-nucleus single-cell counts into per-sample pseudobulk matrices, for RNA-seq and ATAC-seq.

Overview#

Single-nuclei pseudobulk preprocessing for RNA-seq and ATAC-seq. Aggregates per-nucleus counts into per-sample pseudobulk matrices, harmonizes sample IDs, and regresses out technical covariates to produce QTL-ready phenotype values.

The pipeline has three stages:

  • pseudobulk_counts – aggregate a Seurat object into a raw pseudobulk count matrix for one cell type.

  • sampleid_mapping – remap individualID headers to standardized sampleid across metadata and count matrices.

  • pseudobulk_qc – filter, TMM-normalize, fit a technical-covariate model, and write residuals.

(phenotype_formatting reformats residuals into a BED for snATAC-seq caQTL mapping.)

When to run it. Run this module after cell-level counts and sample metadata are available and before phenotype formatting or pseudobulk QTL analysis.

When to run it. After single-nuclei data have been clustered and cell types assigned, and before any QTL scan – the association modules read the matrices this module writes.

Input#

pseudobulk_counts – aggregate single-nuclei counts into a pseudobulk matrix:

  • --seurat-files input/snrnaseq/protocol_example.snrnaseq.seurat_MIC.rds: the Seurat objects to aggregate, each carrying celltype and sample in meta.data.

  • --output-dir (required): where the pseudobulk matrix is written.

  • --celltype (default MIC): the cell type to subset to. Objects lacking it are skipped.

  • --min-cells (default 10): drop samples contributing fewer nuclei than this.

sampleid_mapping – rename matrix columns to study sample IDs:

  • --map-file input/snrnaseq/protocol_example.snrnaseq.id_map.csv (required): two-column map from nuclei-level IDs to study sample IDs.

    individualID,sampleid
DONOR_001,SAMPLE_001
DONOR_002,SAMPLE_002

  • --meta-files input/snrnaseq/protocol_example.snrnaseq.metadata_MIC.csv and --count-files input/snrnaseq/protocol_example.snrnaseq.pseudobulk_counts_MIC.csv.gz: the metadata and count matrices to relabel. The count matrix is genes by samples:

    ,SAMPLE_001,SAMPLE_002,SAMPLE_003,SAMPLE_004,SAMPLE_005,SAMPLE_006,SAMPLE_007,SAMPLE_008
ENSG00000000001,299,172,139,269,742,166,696,820,428,811,188,593,131,270,873,49,1244,259,
ENSG00000000002,533,624,49,2,1911,1822,679,333,147,160,664,92,382,142,160,1323,214,468,1

  • --output-dir (required): where the relabelled files are written.

pseudobulk_qc – filter, optionally batch-correct and normalise:

  • --meta-files and --count-files: the pseudobulk matrices and their metadata.

  • --output-dir (required) and --tech-vars-file input/snrnaseq/protocol_example.snrnaseq.tech_vars_MIC.csv (required): output location, and the technical covariates to regress out.

  • --min-count (default 5), --min-total-count (default 15), --min-prop (default 0.1): edgeR-style expression filters.

  • --min-nuclei (default 20): drop samples built from fewer nuclei than this.

  • --batch-correction (default FALSE) and --batch-method (default limma): whether to correct for batch, and with what.

  • --quant-norm (default FALSE): quantile-normalise the filtered matrix.

  • --blacklist-file, --regions, --gene-list (all empty by default): optional restrictions on the features carried forward.

phenotype_formatting – reshape residuals into the phenotype format (snATAC-seq only):

  • --residual-files input/snrnaseq/atac_residuals/MIC/protocol_example.snrnaseq.MIC_residuals.txt: the residualised matrices to reformat.

  • --gtf-file input/reference_data/Homo_sapiens.GRCh38.103.chr.reformatted.collapse_only.gene.ERCC.gtf (required): gene model supplying the coordinates.

  • --output-dir (required): where the phenotype files are written.

Shared: --cwd (default output), --modular-script-dir (default code/script), and the cluster options --numThreads (default 8), --job-size (default 1), --walltime (default 5h), --mem (default 16G), --container.

Output#

  • pseudobulk_counts_<celltype>.csv.gz – from pseudobulk_counts: raw aggregated counts, one column per sample. Normalisation is left to pseudobulk_qc.

    gene_id,SAMPLE_001,SAMPLE_002,SAMPLE_003,SAMPLE_004,SAMPLE_005,SAMPLE_006,SAMPLE_007,SAM
NOC2L,29,7,4,16,24,16,7,13,3,3,12,3,9,0,12,0,33,37,30,9
KLHL17,8,1,0,6,10,16,5,5,0,1,6,0,5,1,3,2,15,14,19,6

  • relabelled count and metadata files – from sampleid_mapping: the same matrices with columns renamed to study sample IDs via --map-file.

  • filtered, optionally batch-corrected and normalised matrices plus their residuals – from pseudobulk_qc, as GCT:

    ID	SAMPLE_011	SAMPLE_012	SAMPLE_013	SAMPLE_014	SAMPLE_015	SAMPLE_016	SAMPLE_017	SAMPLE_0
12	7.218285529299902	7.765796037583676	6.875518837530288	4.3191760256225535	7.1895314838
18	6.90933528841075	7.765796037583676	6.828976251593258	6.641104120509915	6.252725310296

  • phenotype-format files – from phenotype_formatting, for snATAC-seq only.

  • <name>.<celltype>.residuals.txt – from pseudobulk_qc: the residualised matrix per cell type, which phenotype_formatting reads.

    	SAMPLE_001	SAMPLE_002	SAMPLE_003	SAMPLE_004	SAMPLE_005	SAMPLE_006	SAMPLE_007	SAMPLE_008
ENSG00000000001	10.5316181036712	10.0782299856604	9.69847909370219	10.0321840665649	11.3
ENSG00000000002	9.94704712753556	10.5923588918617	8.97464386590436	4.12683911798198	12.9

  • a .stdout / .stderr log beside each target.

All paths are relative to the --output-dir given to each workflow.

Example output for this module lives under output/snrna_seq/ (counts and residuals) and output/pseudobulk_qc/ (filtered and normalised matrices).

Minimal Working Example#

Step 1. Pseudobulk count matrix generation#

Aggregates single-nuclei counts into pseudobulk count matrices per cell type from Seurat objects.

Timing: TBD (on toy dataset)

sos run pipeline/pseudobulk_preprocessing.ipynb pseudobulk_counts \
    --seurat-files input/snrnaseq/protocol_example.snrnaseq.seurat_MIC.rds \
    --celltype MIC \
    --output-dir output/snrna_seq

Step 2. Sample ID mapping#

Maps original sample identifiers (individualID) to standardized sample IDs (sampleid) across metadata and count matrix files.

Timing: TBD (on toy dataset)

sos run pipeline/pseudobulk_preprocessing.ipynb sampleid_mapping \
    --map-file input/snrnaseq/protocol_example.snrnaseq.id_map.csv \
    --meta-files input/snrnaseq/protocol_example.snrnaseq.metadata_MIC.csv \
    --output-dir output/snrna_seq

Step 3. Pseudobulk QC#

Regresses out technical covariates for downstream QTL analysis. Works for both snATAC-seq and snRNA-seq.

Timing: TBD (on toy dataset)

sos run pipeline/pseudobulk_preprocessing.ipynb pseudobulk_qc \
    --meta-files input/snrnaseq/protocol_example.snrnaseq.metadata_MIC.csv \
    --count-files input/snrnaseq/protocol_example.snrnaseq.pseudobulk_counts_MIC.csv.gz \
    --tech-vars-file input/snrnaseq/protocol_example.snrnaseq.tech_vars_MIC.csv \
    --output-dir output/snrna_seq

Step 4. Phenotype reformatting (snATAC-seq only)#

Converts residuals into a QTL-ready BED format for genome-wide caQTL mapping.

Timing: TBD (on toy dataset)

sos run pipeline/pseudobulk_preprocessing.ipynb phenotype_formatting \
    --residual-files input/snrnaseq/atac_residuals/MIC/protocol_example.snrnaseq.MIC_residuals.txt \
    --output-dir output/snrna_seq \
    --gtf-file input/reference_data/Homo_sapiens.GRCh38.103.chr.reformatted.collapse_only.gene.ERCC.gtf

Command Interface#

sos run pipeline/pseudobulk_preprocessing.ipynb -h

usage: sos run pipeline/pseudobulk_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:
  pseudobulk_counts
  sampleid_mapping
  pseudobulk_qc
  phenotype_formatting

Global Workflow Options:
  --modular-script-dir code/script (as path)
  --cwd output (as path)
  --job-size 1 (as int)
  --walltime 5h
  --mem 16G
  --numThreads 8 (as int)
  --container ''

Sections
  pseudobulk_counts:
    Workflow Options:
      --seurat-files  (as list)
      --output-dir VAL (as str, required)
      --celltype MIC
      --min-cells 10 (as int)
  sampleid_mapping:
    Workflow Options:
      --map-file VAL (as str, required)
      --output-dir VAL (as str, required)
      --meta-files  (as list)
      --count-files  (as list)
  pseudobulk_qc:
    Workflow Options:
      --meta-files  (as list)
      --count-files  (as list)
      --output-dir VAL (as str, required)
      --tech-vars-file VAL (as str, required)
      --blacklist-file ''
      --batch-correction FALSE
      --batch-method limma
      --quant-norm FALSE
      --min-count 5 (as int)
      --min-total-count 15 (as int)
      --min-prop 0.1 (as float)
      --min-nuclei 20 (as int)
      --regions ''
      --gene-list ''
  phenotype_formatting:
    Workflow Options:
      --residual-files  (as list)
      --output-dir VAL (as str, required)
      --gtf-file VAL (as path, required)

Workflow implementation#

Default Parameters#

[global]
parameter: modular_script_dir = path('code/script')  # override with --modular-script-dir
parameter: cwd = path("output")
parameter: name = "protocol_example"
parameter: job_size = 1
parameter: walltime = "5h"
parameter: mem = "16G"
parameter: numThreads = 8
parameter: container = ""


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

usage: sos run pipeline/pseudobulk_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:
  sampleid_mapping
  pseudobulk_qc
  phenotype_formatting
Global Workflow Options:
  --cwd output (as path)
  --job-size 1 (as int)
  --walltime 5h
  --mem 16G
  --numThreads 8 (as int)
  --container ''
Sections
  sampleid_mapping:
    Workflow Options:
      --map-file VAL (as str, required)
      --output-dir VAL (as str, required)
      --meta-files  (as list)
      --count-files  (as list)
  pseudobulk_qc:
    Workflow Options:
      --meta-files  (as list)
      --count-files  (as list)
      --output-dir VAL (as str, required)
      --tech-vars-file VAL (as str, required)
      --blacklist-file ''
      --batch-correction FALSE
      --batch-method limma
      --quant-norm FALSE
      --min-count 5 (as int)
      --min-total-count 15 (as int)
      --min-prop 0.1 (as float)
      --min-nuclei 20 (as int)
      --regions ''
      --gene-list ''
  phenotype_formatting:
    Workflow Options:
      --residual-files  (as list)
      --output-dir VAL (as str, required)
      --gtf-file VAL (as path, required)

pseudobulk_counts#

[pseudobulk_counts]
parameter: seurat_files = []
parameter: output_dir   = str
parameter: celltype     = 'MIC'
parameter: min_cells    = 10

import os

input:  seurat_files
output: f'{output_dir}/{name}.pseudobulk_counts_{celltype}.csv.gz'

task: trunk_workers = 1, trunk_size = 1, walltime = '4:00:00', mem = '64G', cores = 4

bash: expand = "${ }", stdout = f'{_output:n}.stdout', stderr = f'{_output:n}.stderr'
    Rscript ${modular_script_dir}/molecular_phenotypes/QC/pseudobulk_preprocessing.R \
        --name "${name}" \
        --step pseudobulk_counts \
        --seurat-files ${' '.join([f'"{f}"' for f in seurat_files])} \
        --celltype "${celltype}" \
        --min-cells ${min_cells} \
        --output-dir "${output_dir}"

sampleid_mapping#

[sampleid_mapping]
parameter: map_file    = str
parameter: output_dir  = str
parameter: meta_files  = []
parameter: count_files = []

import os

input:  meta_files + count_files
output: [f'{output_dir}/{os.path.basename(f)}' for f in meta_files + count_files]
         
bash: expand = "${ }"
    Rscript ${modular_script_dir}/molecular_phenotypes/QC/pseudobulk_preprocessing.R \
        --name "${name}" \
        --step sampleid_mapping \
        --map-file "${map_file}" \
        ${('--meta-files ' + ' '.join([f'"{f}"' for f in meta_files])) if meta_files else ''} \
        ${('--count-files ' + ' '.join([f'"{f}"' for f in count_files])) if count_files else ''} \
        --output-dir "${output_dir}"

pseudobulk_qc#

[pseudobulk_qc]
parameter: meta_files       = []
parameter: count_files      = []
parameter: output_dir       = str
parameter: tech_vars_file   = str
parameter: blacklist_file   = ''
parameter: batch_correction = "FALSE"
parameter: batch_method     = "limma"
parameter: quant_norm       = "FALSE"
parameter: min_count        = 5
parameter: min_total_count  = 15
parameter: min_prop         = 0.1
parameter: min_nuclei       = 20
parameter: regions          = ''
parameter: gene_list        = ''

import os

_cts = [os.path.basename(f).split('metadata_')[-1].replace('.csv','') for f in meta_files]

input:  meta_files + count_files
output: [f'{output_dir}/{name}.{ct}.residuals.txt' for ct in _cts]

task: trunk_workers = 1, trunk_size = 1, walltime = '6:00:00', mem = '64G', cores = 4

bash: expand = "${ }", stdout = f'{_output[0]:n}.stdout', stderr = f'{_output[0]:n}.stderr'
    Rscript ${modular_script_dir}/molecular_phenotypes/QC/pseudobulk_preprocessing.R \
        --name "${name}" \
        --step pseudobulk_qc \
        --meta-files ${' '.join([f'"{f}"' for f in meta_files])} \
        --count-files ${' '.join([f'"{f}"' for f in count_files])} \
        --tech-vars-file "${tech_vars_file}" \
        --blacklist-file "${blacklist_file}" \
        --batch-correction "${batch_correction}" \
        --batch-method "${batch_method}" \
        --quant-norm "${quant_norm}" \
        --min-count ${min_count} \
        --min-total-count ${min_total_count} \
        --min-prop ${min_prop} \
        --min-nuclei ${min_nuclei} \
        --regions "${regions}" \
        --gene-list "${gene_list}" \
        --output-dir "${output_dir}"

phenotype_reformatting#

[phenotype_formatting]
parameter: residual_files = []
parameter: output_dir     = str
parameter: gtf_file       = path

import os

_cts = [os.path.basename(f).replace(f'{name}.','').replace('.residuals.txt','').replace('.residuals_qn.txt','') for f in residual_files]

input:  residual_files
output: [f'{output_dir}/{name}.{ct}.phenotype.bed.gz' for ct in _cts]

task: trunk_workers = 1, trunk_size = 1, walltime = '2:00:00', mem = '16G', cores = 2

bash: expand = "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout'
    Rscript ${modular_script_dir}/molecular_phenotypes/QC/pseudobulk_preprocessing.R \
        --name "${name}" \
        --step phenotype_formatting \
        --residual-files ${' '.join([f'"{f}"' for f in residual_files])} \
        --gtf-file "${gtf_file}" \
        --output-dir "${output_dir}"