10x scRNA-seq: Cell Ranger → GeneExt → CellBender

scRNA-seq

Genome annotation

HPC & computing

From FASTQs to ambient-RNA-corrected count matrices for non-model invertebrates, with GeneExt-extended gene models. Runs as an LSF job chain on an HPC cluster.

Published

May 7, 2026

Modified

October 16, 2026

Pipeline·Tools Cell Ranger 8.0.0 samtools GeneExt CellBender

Overview

10x Genomics 3’ chemistries capture the 3’ end of transcripts. In non-model genomes the gene models usually have short or missing 3’ UTRs, so many reads land just downstream of the annotated gene and are thrown away. GeneExt uses the reads themselves to extend gene models at the 3’ end. More reads are then assigned to genes, giving more UMIs and genes per cell. (Background and first tests: GeneExt on scRNA-seq data.)

The pipeline runs Cell Ranger twice: once with the original annotation (only to get BAMs for GeneExt), and once with the extended annotation (the counts you analyze). CellBender then removes ambient RNA.

flowchart LR
  F[FASTQs<br/>per library] --> C1[1. cellranger count<br/>original ref]
  C1 --> M[2. samtools merge<br/>all BAMs]
  M --> G[3. GeneExt<br/>extend 3' ends]
  G --> R[4. cellranger mkref<br/>extended GTF]
  R --> C2[5. cellranger count<br/>extended ref]
  F --> C2
  C2 --> B[6. CellBender<br/>ambient RNA]
  B --> O[(filtered .h5<br/>for Seurat / scanpy)]

Every step is submitted as an LSF job that waits (#BSUB -w "done(...)") on the step before it. You submit the whole chain once and get an email as each job finishes.

Inputs

All of these go in your config file (see Run it).

Input	Config setting	Notes
Genome FASTA	`FASTA`
Gene models (GTF)	`GTF`	Exon-only, validated GTF. Cell Ranger and GeneExt both need consistent `gene_id`/`transcript_id` attributes.
Cell Ranger reference built from the original GTF	`REF_ORIG`	Reuse one if it exists; otherwise set `RUN_MKREF_ORIG=true` and the pipeline builds it.
FASTQ folder per library	`LIB_FASTQS`	One folder per library
FASTQ sample prefix per library	`LIB_SAMPLES`	The start of the FASTQ file names, passed to `cellranger count --sample`
Library IDs	`LIB_IDS`	Short names used for output folders, e.g. `lib-001`

Software

Tool	How it’s loaded (defaults in the config; change them for your cluster)
Cell Ranger 8.0.0	`module load cellranger/8.0.0`
samtools	conda env `samtools_env` (keep it up to date, see Known issues)
GeneExt	git clone of GeneExt (`GENEEXT_INSTALL`) + conda env `geneext` (install notes)
CellBender	conda env `CellBender`

How the conda environments were set up on the University of Miami Pegasus cluster: Pegasus conda environments.

Default resources for every step: queue bigmem, 16 cores, rusage[mem=15000], 24–120 h walltime. Change QUEUE, CORES and MEM_MB in the config for your cluster.

Run it

The script never needs editing. Everything specific to your run (paths, library names, email, LSF project) goes in a small config file.

Make a config file on the cluster login node:

bash scrnaseq_10x_geneext_cellbender.sh --init my_run.config

Fill it in. Open my_run.config (nano my_run.config) and replace every <fill in: ...> value. Use absolute paths. Software settings (module names, conda env names) have defaults you can change.

RUN="<fill in: short run name, used in job and reference names, e.g. species_runID>"
EMAIL="<fill in: email for LSF job notifications>"
FASTA="<fill in: /path/to/genome.fa>"
LIB_IDS=(     "<fill in: library ID 1, e.g. lib-001>"     "<fill in: library ID 2>" )
# ...

Submit the chain:
```
bash scrnaseq_10x_geneext_cellbender.sh my_run.config
bjobs -w | grep <your RUN name>        # check progress
```
Before submitting anything, the script stops if any <fill in> value is left, if the library lists have different lengths, or if the FASTA, GTF, reference, FASTQ folders or GeneExt install can’t be found.

If a step fails, fix the problem, bkill the jobs still waiting, and resume from that step. Earlier outputs are reused:

START_AT=3 bash scrnaseq_10x_geneext_cellbender.sh my_run.config   # 2 = merge, 3 = GeneExt, 5 = second count, 6 = CellBender

The exact job file for every step is saved in <OUTDIR>/logs/ next to its .out/.err logs, so each run keeps a record of what was submitted. Keep the config file with your project notes: it is the full record of a run’s settings.

Full script

Download scrnaseq_10x_geneext_cellbender.sh

scrnaseq_10x_geneext_cellbender.sh

#!/bin/bash
# =============================================================================
# 10x scRNA-seq: Cell Ranger -> GeneExt -> Cell Ranger -> CellBender  (LSF clusters, e.g. Pegasus)
# Pipeline page: https://kevinhwong1.github.io/KevinHWong_Notebook/pipelines/scrnaseq-10x-geneext-cellbender.html
#
# Steps (each is an LSF job that waits on the one before it):
#   0.  cellranger mkref   original GTF            (optional: RUN_MKREF_ORIG=true)
#   1.  cellranger count   every library, original ref   (BAMs for GeneExt)
#   2.  samtools merge     all library BAMs -> merge.bam
#   3.  GeneExt            extend 3' ends of gene models using merge.bam
#   4.  cellranger mkref   GeneExt GTF
#   5.  cellranger count   every library, GeneExt ref
#   6.  CellBender         ambient RNA removal on each library
#
# Usage:
#   bash scrnaseq_10x_geneext_cellbender.sh --init my_run.config    # 1. write a config template
#   nano my_run.config                                              # 2. replace every placeholder
#   bash scrnaseq_10x_geneext_cellbender.sh my_run.config           # 3. submit the whole chain
#   START_AT=3 bash scrnaseq_10x_geneext_cellbender.sh my_run.config  # resume from a step (0-6)
#
# Nothing in this file needs editing; all run-specific settings live in the config file.
# =============================================================================
set -euo pipefail

usage() { sed -n '15,19p' "$0" | sed 's/^# \{0,1\}//'; }

write_template() {
cat > "$1" <<'TEMPLATE'
# Config for scrnaseq_10x_geneext_cellbender.sh
# Replace every "<fill in: ...>" value. The script refuses to run while any are left.
# This file is sourced by bash: no spaces around "=", and quote values.

# --- Run and cluster ----------------------------------------------------------
RUN="<fill in: short run name, used in job and reference names, e.g. species_runID>"
LSF_PROJECT="<fill in: LSF project/allocation passed to bsub -P>"
EMAIL="<fill in: email for LSF job notifications>"
QUEUE="bigmem"          # LSF queue
CORES=16                # cores per job (bsub -n)
MEM_MB=15000            # bsub -R "rusage[mem=...]"

# --- Genome -------------------------------------------------------------------
FASTA="<fill in: /path/to/genome.fa>"
GTF="<fill in: /path/to/genes.gtf (exons only, consistent gene_id/transcript_id)>"
REF_ORIG="<fill in: /path/to/existing/cellranger/reference (or where step 0 should build it)>"
RUN_MKREF_ORIG=false    # true = build REF_ORIG from FASTA + GTF first (step 0)

# --- Libraries: one entry per library, same order in all three lists -----------
LIB_IDS=(     "<fill in: library ID 1, e.g. lib-001>"     "<fill in: library ID 2>" )
LIB_SAMPLES=( "<fill in: FASTQ sample prefix 1 (cellranger --sample)>" "<fill in: FASTQ sample prefix 2>" )
LIB_FASTQS=(  "<fill in: /path/to/fastq/folder/1>"         "<fill in: /path/to/fastq/folder/2>" )

# --- Output and software --------------------------------------------------------
OUTDIR="<fill in: /path/to/scratch/output/folder>"
GENEEXT_INSTALL="<fill in: /path/to/GeneExt (git clone containing geneext.py)>"

CELLRANGER_MODULE="cellranger/8.0.0"
CONDA_SH="${HOME}/anaconda3/bin/activate"   # conda activation script
SAMTOOLS_ENV="samtools_env"                 # conda env names
GENEEXT_ENV="geneext"
CELLBENDER_ENV="CellBender"

# --- Tool settings ----------------------------------------------------------------
GENEEXT_ARGS="-v 3 -j 8 --clip_strand both"        # --clip_strand both: no same-strand overlaps
CELLBENDER_ARGS="--fpr 0.01 --epochs 150 --learning-rate 0.00005"
TEMPLATE
}

# ----------------------------------------------------------------------------- read config
if [[ "${1:-}" == "--init" ]]; then
  cfg="${2:-run.config}"
  if [[ -e "${cfg}" ]]; then echo "${cfg} already exists; not overwriting." >&2; exit 1; fi
  write_template "${cfg}"
  echo "Wrote ${cfg}. Fill in every <fill in: ...> value, then run:"
  echo "  bash $0 ${cfg}"
  exit 0
fi

CONFIG="${1:-}"
if [[ -z "${CONFIG}" || ! -f "${CONFIG}" ]]; then usage; exit 1; fi

if grep -nE '^[^#]*<fill in' "${CONFIG}"; then
  echo "" >&2
  echo "Replace the placeholder values above in ${CONFIG} first." >&2
  exit 1
fi

# shellcheck source=/dev/null
source "${CONFIG}"
START_AT="${START_AT:-0}"

# ----------------------------------------------------------------------------- sanity checks
problems=0
need_file() { [[ -e "$1" ]] || { echo "Not found: $1  ($2)" >&2; problems=$((problems+1)); }; }
if [[ ${#LIB_IDS[@]} -ne ${#LIB_SAMPLES[@]} || ${#LIB_IDS[@]} -ne ${#LIB_FASTQS[@]} ]]; then
  echo "LIB_IDS, LIB_SAMPLES and LIB_FASTQS must have the same number of entries." >&2; problems=$((problems+1))
fi
need_file "${FASTA}" FASTA
need_file "${GTF}" GTF
need_file "${GENEEXT_INSTALL}/geneext.py" "GeneExt install"
if (( START_AT <= 1 )) && [[ "${RUN_MKREF_ORIG}" != true ]]; then need_file "${REF_ORIG}" "REF_ORIG (or set RUN_MKREF_ORIG=true)"; fi
if (( START_AT <= 5 )); then
  for d in "${LIB_FASTQS[@]}"; do need_file "${d}" "FASTQ folder"; done
fi
if (( problems > 0 )); then echo "Fix the ${problems} problem(s) above in ${CONFIG}." >&2; exit 1; fi

GENEEXT_DIR="${OUTDIR}/GeneExt"
GENEEXT_GTF="${OUTDIR}/${RUN}_GeneExt.gtf"
REF_GENEEXT="${OUTDIR}/${RUN}_geneext"
LOGDIR="${OUTDIR}/logs"
mkdir -p "${OUTDIR}" "${GENEEXT_DIR}" "${LOGDIR}"

# deps JID [JID ...] -> 'done(1) && done(2)', skipping empty IDs (steps not run this time)
deps() {
  local out="" j
  for j in "$@"; do
    [[ -z "${j}" ]] && continue
    out+="${out:+ && }done(${j})"
  done
  echo "${out}"
}

# submit NAME WALLTIME "DEPS" <<EOF ... EOF   -> prints the LSF job ID
# The job file is kept in ${LOGDIR} (not /tmp) so every run has a record of what was submitted.
submit() {
  local name=$1 walltime=$2 wait=$3 jobfile="${LOGDIR}/${1}.lsf"
  {
    echo "#!/bin/bash"
    echo "#BSUB -J ${RUN}_${name}"
    echo "#BSUB -q ${QUEUE}"
    echo "#BSUB -P ${LSF_PROJECT}"
    echo "#BSUB -n ${CORES}"
    echo "#BSUB -W ${walltime}"
    echo "#BSUB -R \"rusage[mem=${MEM_MB}]\""
    [[ -n "${wait}" ]] && echo "#BSUB -w \"${wait}\""
    echo "#BSUB -u ${EMAIL}"
    echo "#BSUB -o ${LOGDIR}/${name}_%J.out"
    echo "#BSUB -e ${LOGDIR}/${name}_%J.err"
    echo "#BSUB -N"
    echo "set -eo pipefail   # no -u: conda activate scripts reference unset variables"
    cat                                   # job body from the heredoc
  } > "${jobfile}"
  bsub < "${jobfile}" | grep -oP '(?<=Job <)\d+'
}

echo "== ${RUN}: submitting LSF job chain from step ${START_AT} ($(date)) =="

# ----------------------------------------------------------------------------- STEP 0
JID_MKREF_ORIG=""
if (( START_AT <= 0 )) && [[ "${RUN_MKREF_ORIG}" == true ]]; then
  JID_MKREF_ORIG=$(submit step0_mkref_orig 24:00 "" <<EOF
module load ${CELLRANGER_MODULE}
cd $(dirname "${REF_ORIG}")
cellranger mkref --genome=$(basename "${REF_ORIG}") --fasta=${FASTA} --genes=${GTF}
EOF
)
  echo "step 0  mkref (original)      ${JID_MKREF_ORIG}"
fi

# ----------------------------------------------------------------------------- STEP 1
JIDS_COUNT_ORIG=()
if (( START_AT <= 1 )); then
  for i in "${!LIB_IDS[@]}"; do
    jid=$(submit "step1_count_${LIB_IDS[$i]}" 120:00 "$(deps "${JID_MKREF_ORIG}")" <<EOF
module load ${CELLRANGER_MODULE}
cd ${OUTDIR}
cellranger count \\
  --id=${LIB_IDS[$i]} \\
  --transcriptome=${REF_ORIG} \\
  --fastqs=${LIB_FASTQS[$i]} \\
  --sample=${LIB_SAMPLES[$i]} \\
  --create-bam=true
EOF
)
    JIDS_COUNT_ORIG+=("${jid}")
    echo "step 1  count ${LIB_IDS[$i]} (original) ${jid}"
  done
fi

# ----------------------------------------------------------------------------- STEP 2
JID_MERGE=""
if (( START_AT <= 2 )); then
  BAMS=""
  for id in "${LIB_IDS[@]}"; do BAMS+=" ${OUTDIR}/${id}/outs/possorted_genome_bam.bam"; done
  JID_MERGE=$(submit step2_bam_merge 24:00 "$(deps "${JIDS_COUNT_ORIG[@]:-}")" <<EOF
source ${CONDA_SH}
conda activate ${SAMTOOLS_ENV}
samtools --version | head -1
cd ${GENEEXT_DIR}
samtools merge -f -@ ${CORES} merge.bam ${BAMS}
samtools quickcheck merge.bam && echo "merge.bam OK"
ls -lh ${BAMS} merge.bam
EOF
)
  echo "step 2  samtools merge        ${JID_MERGE}"
fi

# ----------------------------------------------------------------------------- STEP 3
JID_GENEEXT=""
if (( START_AT <= 3 )); then
  JID_GENEEXT=$(submit step3_geneext 120:00 "$(deps "${JID_MERGE}")" <<EOF
source ${CONDA_SH}
conda activate ${GENEEXT_ENV}
# GeneExt writes tmp/ and calls its R scripts with paths relative to the install dir,
# so it has to run from there, and a stale tmp/ from an earlier run must be removed first.
cd ${GENEEXT_INSTALL}
rm -rf tmp
python geneext.py -g ${GTF} -b ${GENEEXT_DIR}/merge.bam -o ${GENEEXT_GTF} ${GENEEXT_ARGS}
tail -20 ${GENEEXT_GTF}.GeneExt.log
EOF
)
  echo "step 3  GeneExt               ${JID_GENEEXT}"
fi

# ----------------------------------------------------------------------------- STEP 4
JID_MKREF_GE=""
if (( START_AT <= 4 )); then
  JID_MKREF_GE=$(submit step4_mkref_geneext 120:00 "$(deps "${JID_GENEEXT}")" <<EOF
module load ${CELLRANGER_MODULE}
cd ${OUTDIR}
cellranger mkref --genome=$(basename "${REF_GENEEXT}") --fasta=${FASTA} --genes=${GENEEXT_GTF}
EOF
)
  echo "step 4  mkref (GeneExt)       ${JID_MKREF_GE}"
fi

# ----------------------------------------------------------------------------- STEPS 5 + 6
for i in "${!LIB_IDS[@]}"; do
  id="${LIB_IDS[$i]}"
  JID_COUNT_GE=""
  if (( START_AT <= 5 )); then
    JID_COUNT_GE=$(submit "step5_count_${id}_geneext" 120:00 "$(deps "${JID_MKREF_GE}")" <<EOF
module load ${CELLRANGER_MODULE}
cd ${OUTDIR}
cellranger count \\
  --id=${id}_geneext \\
  --transcriptome=${REF_GENEEXT} \\
  --fastqs=${LIB_FASTQS[$i]} \\
  --sample=${LIB_SAMPLES[$i]} \\
  --create-bam=true
EOF
)
    echo "step 5  count ${id} (GeneExt)  ${JID_COUNT_GE}"
  fi

  JID_CB=$(submit "step6_cellbender_${id}" 120:00 "$(deps "${JID_COUNT_GE}")" <<EOF
source ${CONDA_SH}
conda activate ${CELLBENDER_ENV}
cd ${OUTDIR}/${id}_geneext
cellbender remove-background \\
  --input  ${OUTDIR}/${id}_geneext/outs/raw_feature_bc_matrix.h5 \\
  --output ${OUTDIR}/${id}_geneext/${id}_geneext_cellbender_output.h5 \\
  ${CELLBENDER_ARGS}
EOF
)
  echo "step 6  CellBender ${id}      ${JID_CB}"
done

cat <<EOF

All jobs submitted. Monitor with:  bjobs -w | grep ${RUN}
Logs and the exact job scripts:     ${LOGDIR}/
Outputs to copy back, per library:  ${OUTDIR}/<lib>_geneext/<lib>_geneext_cellbender_output{_filtered.h5,_metrics.csv,_report.html,.pdf}
                                    ${OUTDIR}/<lib>_geneext/outs/web_summary.html
EOF

Outputs

For each library <lib>, in ${OUTDIR}/<lib>_geneext/:

File	What it is
`outs/web_summary.html`	Cell Ranger QC report (compare with the step-1 report in `${OUTDIR}/<lib>/outs/`)
`<lib>_geneext_cellbender_output_filtered.h5`	The matrix to analyze: CellBender-corrected, cell barcodes only
`<lib>_geneext_cellbender_output_metrics.csv`	CellBender summary metrics
`<lib>_geneext_cellbender_output_report.html` / `.pdf`	CellBender training and QC report

Plus ${OUTDIR}/<run>_GeneExt.gtf (the extended annotation), with .extension_length.pdf and .peak_coverage.pdf plots from GeneExt.

Copy them back to your laptop:

H=<fill in: user@cluster.address>
D=<fill in: OUTDIR from your config>
for lib in <fill in: your LIB_IDS, space-separated>; do
  scp $H:$D/${lib}/outs/web_summary.html            ./${lib}_web_summary.html
  scp $H:$D/${lib}_geneext/outs/web_summary.html    ./${lib}_web_summary_geneext.html
  for f in .pdf _filtered.h5 _metrics.csv _report.html; do
    scp $H:$D/${lib}_geneext/${lib}_geneext_cellbender_output${f} ./
  done
done

Known issues

GeneExt must run from its install directory, with no old tmp/

GeneExt writes intermediate files to ./tmp and calls its plotting scripts with relative paths (Rscript geneext/plot_extensions.R …). Run it from anywhere else and those steps fail. A tmp/ left over from a crashed run also breaks the next run. The script cds into GENEEXT_INSTALL and runs rm -rf tmp first. Because of the shared tmp/, don’t run two GeneExt jobs at the same time from the same install.

samtools merge failed with the old samtools in samtools_env

The first full run (Pdam 31298) failed at the merge step. Updating samtools in the environment fixed it. The job now prints samtools --version and runs samtools quickcheck on the merged BAM, so a bad merge stops the chain before GeneExt.

A failed step leaves the rest of the chain stuck in PEND

LSF never satisfies done(<jobid>) for a job that exited with an error, so the downstream jobs wait forever. Check with bjobs -l <jobid> (“dependency condition invalid or never satisfied”), bkill them, fix the problem, and resubmit with START_AT=<step>.

Checking BAMs: use samtools, not wc -l

BAMs are compressed binary files, so wc -l doesn’t count reads. Use samtools quickcheck file.bam (is the file intact?) or samtools view -c file.bam (how many reads?).

Example run

Pdam 31298 (all cells + ALDH+), May 2026. GeneExt extended 13,283 / 25,422 genes, median extension 1,297 bp.

Changelog

Date	Change
2026-10-05	Run settings moved to a config file (`--init` writes a template with placeholders; the script checks inputs before submitting).
2026-10-05	Turned the Pdam 31298 scripts into this pipeline page: one script for any number of libraries, `START_AT` to resume, job files kept in `logs/`, merge checks, GeneExt `tmp/` fix.
2026-05-07	First full run on Pdam 31298 (three separate scripts: full, post-count, post-merge).