• Added split-parallel-merge workflow for barcode & UMI calling on large BAM files (scCirRL_extract_reads, scCirRL_collect_ref_bc, scCirRL_assign_bc_from_tsv, scCirRL_merge_bc_umi)

scCirRL-seq is a computational pipeline designed for single-cell RNA-seq long-read data.

It mainly consists of three parts:

1. cell barcode/UMI calling: a standalone module for barcode/UMI calling and gene/transcript quantification
2. alternative splicing analysis: modules for identification of cell type-specific and allele-specific splicing
3. visualization of gene and transcript expression in single cells

• Haplotype-resolved full-length transcriptome analysis in single cells by scCirRL-seq
  • Updates (v0.0.1)
  • What is scCirRL-seq
  • Table of Contents
  • Installation
    • Operating system and python version
    • Via conda (install locally for now)
    • From source files
  • 0. Preprocessing
    • 0.1 Mapping
    • 0.2 (optional) Mark potential chimeric long reads
    • 0.3 Transcript identification
  • 1. Barcode & UMI calling
    • 1.1 Input
    • 1.2 Command
    • 1.3 Output
    • 1.4 Parallel barcode & UMI calling (for large BAM files)
  • 2. Cell clustering and annotation
  • 3. Cell type-specific splicing analysis
    • 3.1 Input
    • 3.2 Command
    • 3.3 Output
  • 4. Allele-specific splicing analysis
    • 4.1 Phase long reads with whatshap
    • 4.2 Identify allele-specific splicing
    • 4.3 Identify disease-associated GWAS SNPs from allele-specific spliced genes
  • 5. Visualization
    • 5.1. UMAP plot of gene VIM and EPCAM
    • 5.2. UMAP plot of CD44 transcripts
    • 5.3. UMAP plot of both CD44 gene and transcripts

scCirRL-seq was written in python3 and tested on Linux/Unix systems, so it may not work well with python2 and/or other systems(Windows/Mac).

git clone git@github.com:Xinglab/scCirRL-seq.git
cd scCirRL-seq
conda create --prefix ./conda_env
conda activate ./conda_env
conda install -c conda-forge -c bioconda python=3.8 --file conda_requirements.txt
python setup.py install

git clone git@github.com:Xinglab/scCirRL-seq.git
cd scCirRL-seq && pip install .

For mapping, we recommend using minimap2 in RNA splice mode. Any other long-read RNA-seq alignment tools can also be used here.

• Input

  • long_read.fq/fa: 1D long reads or RCA consensus sequences in fastq/fasta format
  • ref.fa: reference genome
  • (optional) anno.gtf: gene annotation file in GTF format
  • (optional) n_threads: number of threads to use

• Command

# 1. convert GTF to BED12 using paftools.js from minimap2
paftools.js gff2bed anno.gtf > anno_junc.bed12

# 2. mapping with minimap2
minimap2 ref.fa long_read.fq/fa         \
         --junc-bed anno_junc.bed12     \
         -Yax splice -ub -k14 -w4       \
         --sam-hit-only --secondary=no  \
         -t n_threads -o long_read.sam

# 3. convert sam to sorted bam
samtools view long_read.sam -b long_read.bam
samtools sort long_read.bam -@ n_threads -o long_read.sorted.bam

Optionally, you can use scCirRL_split_chimeric_read to mark potential chimeric long reads, which may rescue some reads from secondary alignments.

scCirRL_split_chimeric_read long_read.sorted.bam long_read.split.sorted.bam

For transcript identification, we recommend using ESPRESSO(≥1.3.1), other tools like Bambu or IsoQuant can also be used.

• Input
  • long_read.sorted.bam: sorted long-read alignment file in BAM format
  • ref.fa: reference genome file in FASTA format
  • anno.gtf: gene annotation file in GTF format
  • esp_output_dir: output directory
• Output
  • esp_output_dir/esp_cmpt.tsv: compatible isoforms for all reads
  • esp_output_dir/for_esp_input_N2_R0_updated.gtf: updated GTF annotation
• Command

# 1. create tab-separated input file for ESPRESSO
mkdir esp_output_dir 2> /dev/null
in_tsv=esp_output_dir/for_esp_input.tsv
abs_path=$(realpath long_read.sorted.bam)
base_name=$(basename long_read.sorted.bam)
echo -e "$abs_path\t$base_name" > $in_tsv

# 2. ESPRESSO S step
perl ESPRESSO_S.pl -L esp_output_dir/for_esp_input.tsv  \
                   -F ref.fa -A anno.gtf            \
                   -M notFilterOutchrM -T n_threads \
                   -O esp_output_dir

# 3. ESPRESSO C step
perl ESPRESSO_C.pl -I esp_output_dir -F ref.fa \
                   -X 0 -T n_threads

# 4. ESPRESSO Q step
perl ESPRESSO_Q.pl -L esp_output_dir/for_esp_input.updated \
                   -A anno.gtf -T n_threads            \
                   -V esp_output_dir/esp_cmpt.tsv

Note that -V esp_cmpt.tsv is optional in ESPRESSO Q step, but it is required if you want scCirRL-seq to output gene/transcript quantification information.

scCirRL-seq identifies barcode and UMI from sorted alignment BAM file of single-cell long reads without requiring reference barcode from short-read data.

• Required:
  • long_read.sorted.bam: sorted long-read BAM (recommend using minimap2)
• Optional:
  • read_isoform_compatible.tsv: tabular file of compatible isoforms for all reads, generated by ESPRESSO(≥1.3.1), Bambu or IsoQuant
  • updated.gtf: updated GTF annotation, generated by ESPRESSO(≥1.3.1), Bambu or IsoQuant
  • annotation.gtf: reference annotation GTF file, to retrieve gene names if gene names are not provided in updated.gtf
  • cell_barcode.tsv: reference cell barcode list. If provided, scCirRL-seq will directly use it to guide the barcode calling, only long reads with cell barcodes in the provided list will be kept
  • barcode sequence length (default: 16)
  • max. allowed edit distance between barcode and reference barcode (default: 2)
  • UMI sequence length (default: 12)
  • max. allowed edit distance between PCR-duplicated UMIs (default: 1)

scCirRL long_read.sorted.bam \
        output_dir           \
        -t updated.gtf       \
        -m read_isoform_compatible.tsv \
        -g annotation.gtf

• bc_umi.bam: BAM file with barcode/UMI information for each long read, BAM tags: CB and UB, for barcode and UMI. Note that only long read with barcode/UMI called are kept
• bc_umi.tsv: tabular file with barcode/UMI/gene/transcript information for all barcode-called long reads
• expression_matrix/: folder containing 10X format sparse expression matrix at both gene and transcript level, can be directly parsed using standard single-cell analysis tools, like Seurat/Azimuth
• ref_bc.tsv: reference cell barcode list identified from long reads. Note that if cell barcode list was provided to scCirRL-seq, ref_bc.tsv file will not be generated
• high_qual_bc_umi.rank.tsv: cell barcode list ranked by unique UMI count based on all high-quality long reads. Note that if cell barcode list was provided to scCirRL-seq, high_qual_bc_umi.rank.tsv file will not be generated

Example of bc_umi.tsv:

For large BAM files, the barcode and UMI calling step can be split into parallel jobs and then merged. This six-step workflow produces identical outputs to the main scCirRL command.

Overview

The intermediate reads TSV produced in step 1 stores the pre-extracted barcode/UMI candidate information per read (qname, mapping position, bc, umi, bu, is_perfect, ts_tag), so steps 2–6 never need to re-open the original BAM (except for the optional BAM output in step 4).

Step 1 — Extract reads (run in parallel, one job per chromosome/region)

# split by chromosome
# Note: chrM (MT) is required to be included for single-cell analysis, as it is used to estimate the ambient RNA contamination level
for CHR in chr1 chr2 chr3 chr4 chr5 chr6 chr7 chr8 chr9 chr10 \
           chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19 \
           chr20 chr21 chr22 chrX chrY chrM; do
    scCirRL_extract_reads long_read.sorted.bam \
                          ${CHR}.reads.tsv.gz  \
                          -r ${CHR} &
done
wait

Step 2 — Collect reference barcodes (single pass over all TSVs, runs once)

scCirRL_collect_ref_bc chr*.reads.tsv.gz \
                        ref_bc.tsv       \
                        --high-qual-bc-rank high_qual_bc.rank.tsv

If a reference barcode list is already available (e.g. from 10X Cell Ranger):

scCirRL_collect_ref_bc chr*.reads.tsv.gz ref_bc.tsv -l barcodes.tsv.gz

Step 3 — Split read–isoform compatible TSV by chromosome (runs once)

When a read–isoform compatible TSV is available (from ESPRESSO or Bambu), split it into per-chromosome files so that step 4 can filter each chunk on-the-fly without scanning the whole file.

scCirRL_split_cmpt_tsv read_isoform_compatible.tsv \
                        updated.gtf                 \
                        cmpt_by_chrom/              \
                        --prefix cmpt
# produces: cmpt_by_chrom/cmpt.chr1.tsv, cmpt.chr2.tsv, ..., cmpt.unassigned.tsv

Step 3 can be skipped if no read–isoform compatible TSV is available; in that case omit -m in step 4.

Step 4 — Assign barcodes & UMIs, tag BAM (run in parallel)

for CHR in chr1 chr2 ... chrX chrY chrM; do
    scCirRL_assign_bc_from_tsv                            \
        ${CHR}.reads.tsv.gz ref_bc.tsv                    \
        ${CHR}.bc_umi.tsv                                 \
        --region  ${CHR}                                  \
        --in-bam  long_read.sorted.bam                    \
        --out-bam ${CHR}.bc_umi.bam                       \
        -g annotation.gtf                                 \
        -t updated.gtf                                    \
        -m cmpt_by_chrom/cmpt.${CHR}.tsv &
done
wait

Each chunk BAM is indexed automatically after writing.

Note: Omit --in-bam and --out-bam if you only need the bc_umi.tsv output.

Step 5 — Merge

scCirRL_merge_bc_umi bc_umi.tsv                  \
    --tsv-chunks chr*.bc_umi.tsv                 \
    --bam-chunks chr*.bc_umi.bam                 \
    --out-bam    bc_umi.sorted.bam

Omit --bam-chunks and --out-bam if no BAM output was requested in step 4.

Step 6 — Generate expression matrix

scCirRL_make_expression_matrix bc_umi.tsv \
    expression_matrix/

Output files are written under expression_matrix/gene/ and expression_matrix/transcript/, each containing matrix.mtx.gz, features.tsv.gz, and barcodes.tsv.gz in 10X sparse format.

The merged bc_umi.tsv and bc_umi.sorted.bam from step 5 are drop-in replacements for the outputs of the main scCirRL command and are fully compatible with all downstream scCirRL-seq tools.

All the downstream single-cell long-read analyses rely on the cell type clustering result, which could be accomplished by running Seurat on the gene expression matrix.

Annotating the cell clusters is the process of assigning cell types to each cell cluster. This could be performed manually or based on known reference annotation like Azimuth.

For example, for human peripheral blood mononuclear cells (PBMC), the clustering and annotation result can be obtained by mapping to the Azimuth human PBMC reference dataset.

We provide an example R script for human PBMC data. For data from other species/tissues, this needs to be done manually by users.

Note that not having cell clusters annotated with cell types does not affect the downstream splicing analyses by scCirRL-seq. You can simply provide scCirRL-seq with the cluster IDs for cell barcodes.

Here is an example of the output file for this step, bc_to_cell_type.tsv:

Or, if no cell type annotation was available:

scCirRL-seq performs pairwise comparison to identify differentially spliced genes/transcripts between each two cell types/clusters.

• Required
  • expression_matrix/transcript: transcript count matrix directory, generated by scCirRL-seq
  • bc_to_cell_type.tsv: list of barcode and corresponding cell type or cluster ID, generated by Seurat/Azimuth
  • out_prefix: prefix of output files
• Optional
  • min. number of reads to keep a transcript for each cell type (default: 10)
  • min. number of transcripts to keep a gene for each cell type (default: 2)
  • list of genes of interest. If provided, all genes will still be processed, but only differentially spliced genes that are in this list will be output

scCirRL_cell_type_specific_splicing expression_matrix/transcript \
                                    bc_to_cell_type.tsv          \
                                    out_prefix

• *_cell_specific_spliced_genes.tsv: differentially spliced genes between 2 cell types. Gene FDR ≤0.05, max. delta ratio ≥0.05
• *_cell_specific_spliced_transcripts.tsv: differentially spliced transcripts between 2 cell types. Gene FDR ≤0.05, transcript p value ≤0.05, delta ratio ≥0.05

Example of *_cell_specific_spliced_genes.tsv:

Example of *_cell_specific_spliced_transcripts.tsv:

With additional whole-genome sequencing data, scCirRL-seq can identify allele-specific splicing within each cell type. The allele-specific splicing analysis mainly consists of three steps:

1. Phase long reads with whatshap
2. Identify allele-specific spliced genes/transcripts within each cell type
3. Search for disease-associated GWAS SNPs from identified allele-specific spliced genes

• Input
  • wgs_phased.vcf: WGS phased VCF file, generated from WGS short-read data using GATK and SHAPEIT
  • ref.fa: reference genome fasta file
  • bc_umi.bam: BAM file with barcode/UMI information, generated by scCirRL-seq
• Command

# sort and index `bc_umi.bam`
samtools sort bc_umi.bam -o bc_umi_sorted.bam
samtools index bc_umi_sorted.bam

# identify haplotype
whatshap haplotag wgs_phased.vcf bc_umi_sorted.bam \
                  --reference ref.fa               \
                  --ignore-read-groups             \
                  --skip-missing-contigs           \
                  -o bc_umi_hap.bam                \
                  --output-haplotag-list hap_list.tsv

• Output
  • bc_umi_hap.bam: BAM file with barcode/UMI and additional haplotype information
  • hap_list.tsv: tabular file containing haplotype information for barcode-called long reads

• Input

  • bc_umi.tsv: tabular file with barcode/UMI/gene/transcript information, generated by scCirRL-seq
  • hap_list.tsv: tabular file containing haplotype information, generated by whatshap, see above
  • bc_to_cell_type.tsv: list of barcode and corresponding cell type or cluster ID, generated by Seurat/Azimuth

• Command

scCirRL_allele_specific_splicing bc_umi.tsv \
                                 bc_to_cell_type.tsv \
                                 hap_list.tsv \
                                 output_prefix

• Output
  • *_allele_specific_spliced_genes.tsv: allele-specific spliced genes. Gene FDR ≤0.05, max. delta ratio ≥0.05
  • *_allele_specific_spliced_transcripts.tsv: allele-specific spliced transcripts. Gene FDR ≤0.05, transcript p value ≤0.05, delta ratio ≥0.05

Example of *_allele_specific_spliced_genes.tsv:

Example of *_allele_specific_spliced_transcripts.tsv:

• Input

  • *_allele_specific_spliced_genes.tsv: allele-specific spliced genes, generated by scCirRL_allele_specific_splicing
  • wgs_phased.vcf: WGS phased VCF file, generated from WGS short-read data using GATK and SHAPEIT
  • annotation.gtf: annotation GTF file
  • gwas_catalog_efo.tsv: GWAS catalog database with EFO term
  • ld.duckdb: linkage disequilibrium (LD) score database in duckdb format

• Command

scCirRL_gene_with_gwas_disease -g allele_spliced_genes.tsv \
                               annotation.gtf       \
                               wgs_phased.vcf       \
                               gwas_catalog_efo.tsv \
                               ld.duckdb            \
                               output_prefix

• Output
  • *_gwas_efo_term_stat.tsv: GWAS disease-associated genes, with traits described in EFO term
  • *_gwas_snps_ld: directory of LD scores and GWAS traits for all SNPs, each gene has one separate .ld file

Example of *_gwas_efo_term_stat.tsv:

Example of *_gwas_snps_ld/gene.ld:

Given genes/transcripts of interest (cell-type-specific splicing/allele-specific spciling), scCirRL-seq also offers visualization of your results in the single-cell UMAP space, based on the Seurat R package.

# gene/transcript expression folder generated by scCirRL-seq
gene_mtx <- "output_dir/expression_matrix/gene"
transcript_mtx <- "output_dir/expression_matrix/transcript"
# script for UMAP plot
source("visualization.R)

scCirRL_umap_plot(gene_mtx_dir = gene_mtx,
                  # `feature_mtx_dir` needs to be only one folder/object or the same size as `feature_list`
                  feature_mtx_dir = gene_mtx,
                  feature_list = c("VIM", "EPCAM"))

scCirRL_umap_plot(gene_mtx_dir = gene_mtx,
                  # `feature_mtx_dir` needs to be only one folder/object or the same size as `feature_list`
                  feature_mtx_dir = transcript_mtx,
                  feature_list = c("ENST00000433892", "ENST00000263398"))

scCirRL_umap_plot(gene_mtx_dir = gene_mtx,
                  # `feature_mtx_dir` needs to be only one folder/object or the same size as `feature_list`
                  feature_mtx_dir = c(gene_mtx, transcript_mtx, transcript_mtx),
                  feature_list = c("CD44", "ENST00000433892", "ENST00000263398"))

Note that, for gene_mtx_dir and feature_mtx_dir, you can also provide with Seurat object instead of a matrix folder.

gene_obj = make_seurat_obj(gene_mtx)
transcript_obj = make_seurat_obj(transcript_mtx)

scCirRL_umap_plot(gene_mtx_dir = gene_obj,
                  # `feature_mtx_dir` needs to be only one folder/object or same size as `feature_list`
                  feature_mtx_dir = c(gene_obj, transcript_obj, transcript_obj),
                  feature_list = c("CD44", "ENST00000433892", "ENST00000263398"))