Preparing spliced/unspliced data

Preparing spliced/unspliced data#

partition_de_by_mechanism needs spliced and unspliced layers (or mature / nascent). This page is the recipes: velocyto, STARsolo, kb-python, alevin-fry, merge into an existing AnnData, then regime_diagnosis.

Sections 1–2 are commands, not executed. Sections 3–5 run on a tiny synthetic object so you can see barcode mismatch and the capture check.

Already have the layers? Quickstart.

import anndata as ad
import numpy as np
import pandas as pd

import scatrans as scat

print("scatrans", scat.__version__)

scatrans 0.10.12

Pick a quantifier#

Cell Ranger count (and STARsolo with only --soloFeatures Gene) drops the intron/exon split. You need a velocity-aware mode. Versions below are current as of August 2026; check each tool before a new project.

Tool

Version

When

velocyto

0.17.17

You already have a Cell Ranger BAM

STARsolo

2.7.11b

FASTQ → Gene + Velocyto in one pass

kb-python

0.30.2, --workflow nac

Fast pseudoalignment

alevin-fry + pyroe

0.11.2 / 0.9.3

Large cohorts, USA mode

Those four score intronic vs exonic reads. Metabolic labeling (scNT-seq, sci-fate) already has new/old — skip to the last section. Do not run velocyto on labeling data.

velocyto (0.17.17)#

Runs directly on a completed cellranger count output folder — no separate alignment step.

pip install velocyto==0.17.17   # last velocyto.py release (2019); still the
                                 # standard wrapper around a Cell Ranger BAM

# SAMPLE_DIR is the cellranger count output folder (contains outs/possorted_genome_bam.bam,
# outs/filtered_feature_bc_matrix/, outs/raw_feature_bc_matrix/)
velocyto run10x SAMPLE_DIR /path/to/refdata-gex-GRCh38-2024-A/genes/genes.gtf

# recommended: mask repeat regions to cut spurious intronic reads
# (download a repeat-masker GTF for your genome build first)
# velocyto run10x -m repeat_msk.gtf SAMPLE_DIR genes.gtf

Output: SAMPLE_DIR/velocyto/<sample>.loom, with spliced / unspliced / ambiguous layers already inside.

STARsolo (STAR 2.7.11b)#

One alignment pass produces both the standard gene count matrix and the velocity split.

# 1) build a genome index once per reference (skip if you already have one)
STAR --runMode genomeGenerate \
     --genomeDir STAR_index \
     --genomeFastaFiles genome.fa \
     --sjdbGTFfile annotation.gtf \
     --sjdbOverhang 90 --runThreadN 8

# 2) align + quantify, with Velocyto features on top of Gene
STAR --runMode alignReads \
     --genomeDir STAR_index \
     --readFilesIn R2.fastq.gz R1.fastq.gz \
     --readFilesCommand zcat \
     --soloType CB_UMI_Simple \
     --soloCBwhitelist 3M-february-2018.txt \
     --soloCBstart 1 --soloCBlen 16 --soloUMIstart 17 --soloUMIlen 12 \
     --soloFeatures Gene Velocyto \
     --outSAMtype BAM Unsorted \
     --runThreadN 8 \
     --outFileNamePrefix sample_

3M-february-2018.txt is the 10x v3 barcode whitelist (ships with Cell Ranger, cellranger-cs/*/lib/python/cellranger/barcodes/); for v2 chemistry use 737K-august-2016.txt and --soloCBlen 16 --soloUMIlen 10.

Output: sample_Solo.out/Velocyto/raw/{spliced,unspliced,ambiguous}.mtx plus matching barcodes.tsv / features.tsv in the same folder.

kb-python (0.30.2, --workflow nac)#

nac (“nascent and cDNA”) replaced the older lamanno workflow name in recent kb-python and writes ready-to-use layers directly into an .h5ad.

pip install kb-python==0.30.2

# 1) build the nac index (once per reference)
kb ref --workflow nac \
       -i index.idx -g t2g.txt \
       -f1 cdna.fa -f2 nascent.fa \
       -c1 cdna_t2c.txt -c2 nascent_t2c.txt \
       genome.fa annotation.gtf

# 2) pseudoalign + quantify one sample
kb count --workflow nac \
         -i index.idx -g t2g.txt \
         -c1 cdna_t2c.txt -c2 nascent_t2c.txt \
         -x 10xv3 -o out_dir \
         --h5ad \
         R1.fastq.gz R2.fastq.gz

Output: out_dir/counts_unfiltered/adata.h5ad, already containing adata.layers["mature"], ["nascent"], and ["ambiguous"] — this is the mature/nascent naming scATrans resolves automatically (see §2). Older kb-python (<0.27) used --workflow lamanno and wrote spliced.mtx / unspliced.mtx matrix files instead of an .h5ad.

alevin-fry (0.11.2) + pyroe (0.9.3), USA mode#

Fastest option for large cohorts; needs salmon and the alevin-fry binary in addition to the pyroe Python package.

pip install pyroe==0.9.3
# salmon and alevin-fry are separate binaries, e.g. via conda:
# conda install -c bioconda salmon=1.10.3 alevin-fry=0.11.2

# 1) build a spliced+intron ("splici") reference — read-length is your R2 length minus
#    flank-trim-length (91 - 5 = 86 here, for a 91bp 10x v3 cDNA read)
pyroe make-splici genome.fa annotation.gtf 91 splici_ref \
     --flank-trim-length 5 --filename-prefix splici

# 2) index and map
salmon index -t splici_ref/splici_fl86.fa -i splici_idx -p 8

salmon alevin -l ISR -i splici_idx \
     -1 R1.fastq.gz -2 R2.fastq.gz \
     --chromiumV3 --sketch -p 8 -o alevin_map

# 3) resolve unspliced/spliced/ambiguous (USA) counts
alevin-fry generate-permit-list -d fw -k -i alevin_map -o af_quant
alevin-fry collate -t 8 -i af_quant -r alevin_map
alevin-fry quant -t 8 -i af_quant -o af_quant_res \
     --tg-map splici_ref/splici_fl86_t2g_3col.tsv \
     --resolution cr-like --use-mtx

Use --chromium instead of --chromiumV3 for 10x v2 chemistry. Output: loaded in Python with pyroe.load_fry, not read directly (see §2).

2. Load the output into AnnData#

Each tool’s output lands in AnnData slightly differently.

velocyto (.loom) — layers are already named spliced / unspliced:

import anndata as ad

adata = ad.read_loom("SAMPLE_DIR/velocyto/sample.loom")
# adata.layers["spliced"], adata.layers["unspliced"], adata.layers["ambiguous"]

STARsolo (Velocyto/raw/*.mtx) — three separate Matrix Market files to stitch together:

import scanpy as sc
import pandas as pd

base = "sample_Solo.out/Velocyto/raw"
spliced = sc.read_mtx(f"{base}/spliced.mtx").T
unspliced = sc.read_mtx(f"{base}/unspliced.mtx").T

genes = pd.read_csv(f"{base}/features.tsv", header=None, sep="\t")[1].values
barcodes = pd.read_csv(f"{base}/barcodes.tsv", header=None)[0].values

adata = spliced
adata.var_names = genes
adata.obs_names = barcodes
adata.layers["spliced"] = spliced.X.copy()
adata.layers["unspliced"] = unspliced.X

kb-python (--workflow nac) — already a ready .h5ad, no renaming needed (scATrans resolves mature/nascent automatically):

adata = sc.read_h5ad("out_dir/counts_unfiltered/adata.h5ad")
# adata.layers["mature"], adata.layers["nascent"], adata.layers["ambiguous"]

alevin-fry (via pyroe):

from pyroe import load_fry

# output_format="velocity" -> layers "spliced" and "unspliced" directly
# (there is no "scVelo" format string; "velocity" is the correct one)
adata = load_fry("af_quant_res", output_format="velocity")

Merge into your QC’d AnnData#

Attach layers to the object you already filtered and clustered. Do not restart analysis from the raw velocyto / kb-python file.

The next cells build two tiny stand-ins — adata_main and adata_velocity — with a barcode-suffix mismatch, the usual reason this join returns almost nothing.

rng = np.random.default_rng(0)
n_cells, n_genes = 30, 15
gene_names = [f"Gene{i}" for i in range(n_genes)]

# adata_main: already QC'd, Cell Ranger multi-sample barcodes carry a "-1" suffix
barcodes_main = [f"{''.join(rng.choice(list('ACGT'), 16))}-1" for _ in range(n_cells)]
adata_main = ad.AnnData(
    X=rng.poisson(5, size=(n_cells, n_genes)).astype(float),
    obs=pd.DataFrame(
        {"condition": rng.choice(["Control", "Disease"], n_cells)}, index=barcodes_main
    ),
    var=pd.DataFrame(index=gene_names),
)

# adata_velocity: same cells/genes, but the quantifier's barcodes have no "-1" suffix
# (a very common real-world mismatch) and cover a slightly different gene set
barcodes_velocity = [b.replace("-1", "") for b in barcodes_main]
velocity_genes = gene_names[:-2] + ["DecoyGeneA", "DecoyGeneB"]
adata_velocity = ad.AnnData(
    X=rng.poisson(3, size=(n_cells, n_genes)).astype(float),
    obs=pd.DataFrame(index=barcodes_velocity),
    var=pd.DataFrame(index=velocity_genes),
)
adata_velocity.layers["unspliced"] = rng.poisson(1.2, size=(n_cells, n_genes)).astype(float)
adata_velocity.layers["spliced"] = adata_velocity.X.copy()

print("adata_main barcodes: ", adata_main.obs_names[:2].tolist())
print("adata_velocity barcodes:", adata_velocity.obs_names[:2].tolist())

adata_main barcodes:  ['TGGCCAAAATGTGGTG-1', 'GGGTCTGACTGATGTA-1']
adata_velocity barcodes: ['TGGCCAAAATGTGGTG', 'GGGTCTGACTGATGTA']

# Naive intersection BEFORE fixing the suffix mismatch
common_cells = adata_main.obs_names.intersection(adata_velocity.obs_names)
print(f"common cells without fixing barcodes: {len(common_cells)} / {n_cells}")

common cells without fixing barcodes: 0 / 30

Zero (or near-zero) overlap — exactly the silent failure mode to check for before merging. Fix the suffix, then merge properly.

# Fix the suffix, then merge
adata_velocity.obs_names = [b + "-1" for b in adata_velocity.obs_names]

common_cells = adata_main.obs_names.intersection(adata_velocity.obs_names)
common_genes = adata_main.var_names.intersection(adata_velocity.var_names)
print(f"common cells after fixing barcodes: {len(common_cells)} / {n_cells}")
print(
    f"common genes: {len(common_genes)} / {n_genes} "
    f"(velocity file also had {sorted(set(adata_velocity.var_names) - set(gene_names))})"
)

adata_merged = adata_main[common_cells, common_genes].copy()
adata_velocity_aligned = adata_velocity[common_cells, common_genes]
adata_merged.layers["spliced"] = adata_velocity_aligned.layers["spliced"]
adata_merged.layers["unspliced"] = adata_velocity_aligned.layers["unspliced"]

print(adata_merged)

common cells after fixing barcodes: 30 / 30
common genes: 13 / 15 (velocity file also had ['DecoyGeneA', 'DecoyGeneB'])
AnnData object with n_obs × n_vars = 30 × 13
    obs: 'condition'
    layers: 'spliced', 'unspliced'

scv.utils.merge(adata_main, adata_velocity) does the same intersection (plus scVelo bookkeeping) if you have pip install "scatrans[advanced]":

import scvelo as scv
scv.utils.merge(adata_main, adata_velocity)

Other empty-join causes:

  • Raw vs filtered barcodes. Quantifiers often emit the unfiltered whitelist. Intersection is fine if you only need cells already in adata_main.

  • ambiguous layer. scATrans does not use it. Leave it out. Do not add it to unspliced.

4. Sanity-check before running scATrans#

scat.qc.regime_diagnosis folds the global unspliced fraction into a reliability score that partition_de_by_mechanism uses automatically. Run it once here to make sure the merge above produced something sane before a full analysis.

frac = scat.qc.unspliced_global(adata_merged)
r = scat.qc.regime_diagnosis(adata_merged)
print(f"global unspliced fraction: {frac:.2f}")
print(f"regime: {r['regime']}, reliability: {r['reliability']:.2f}")
print(r["message"])

global unspliced fraction: 0.30
regime: ok, reliability: 1.00
unspliced fraction 30.1% is in the normal band; proxy not obviously corrupted.

A typical 10x 3′ library lands around 10–45% unspliced. Above ~50–70% usually means nuclear enrichment, gDNA, or swapped layers — see the FAQ. Next cell scales unspliced up to show that failure mode.

adata_bad = adata_merged.copy()
adata_bad.layers["unspliced"] = (
    adata_bad.layers["unspliced"] * 25
)  # simulate gDNA/nuclear contamination

frac_bad = scat.qc.unspliced_global(adata_bad)
r_bad = scat.qc.regime_diagnosis(adata_bad)
print(f"global unspliced fraction: {frac_bad:.2f}")
print(f"regime: {r_bad['regime']}, reliability: {r_bad['reliability']:.2f}")
print(r_bad["message"])

global unspliced fraction: 0.92
regime: high_unspliced, reliability: 0.00
unspliced fraction 91.5% is high (>= 45%): possible nuclear/gDNA contamination -> gamma fit and the nascent proxy may be unreliable; mechanism annotations down-weighted.

Reliability drops sharply — partition_de_by_mechanism would scale down mechanism_confidence accordingly rather than reporting mechanism labels it can’t back up.

Metabolic labeling (scNT-seq, sci-fate)#

If you already have new/old (or labeled/unlabeled) from the demultiplexing pipeline, rename them. Do not run velocyto on labeling data.

adata_labeling = adata_merged.copy()
adata_labeling.layers["old"] = adata_labeling.layers.pop("spliced")
adata_labeling.layers["new"] = adata_labeling.layers.pop("unspliced")

# --- the rename scATrans needs ---
adata_labeling.layers["spliced"] = adata_labeling.layers["old"]  # pre-existing / mature
adata_labeling.layers["unspliced"] = adata_labeling.layers["new"]  # newly synthesized / nascent

print(sorted(adata_labeling.layers.keys()))

['new', 'old', 'spliced', 'unspliced']

Labeling-based nascent fractions are a cleaner proxy than intron capture, but they still go through regime_diagnosis.

Next: Quickstart · Partition DE by mechanism — transcription vs stabilization.