/bio-long-read-sequencing-isoseq-analysis

name: bio-long-read-sequencing-isoseq-analysis
description: Analyze PacBio Iso-Seq data for full-length isoform discovery and quantification. Use when characterizing transcript diversity or identifying novel splice variants.
tool_type: cli
primary_tool: IsoSeq3

Version Compatibility

Reference examples tested with: minimap2 2.26+, pandas 2.2+, pysam 0.22+, samtools 1.19+

Before using code patterns, verify installed versions match. If versions differ:

• Python: `pip show <package>` then `help(module.function)` to check signatures
• R: `packageVersion('<pkg>')` then `?function_name` to verify parameters
• CLI: `<tool> --version` then `<tool> --help` to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Iso-Seq Analysis

**"Analyze full-length isoforms from my Iso-Seq data"** → Process PacBio HiFi reads through CCS generation, primer removal, clustering, and isoform classification to discover novel transcript variants.

• CLI: `isoseq3 refine` → `isoseq3 cluster` → `pbmm2 align` → `sqanti3_qc.py`

IsoSeq3 Pipeline Overview

# Full pipeline: subreads -> HQ transcripts
# 1. CCS: Generate circular consensus sequences
# 2. Lima: Remove primers and demultiplex
# 3. Refine: Remove polyA and concatemers
# 4. Cluster: Group into isoforms
# 5. Polish: Generate high-quality consensus (optional with HiFi)

CCS Generation

# Generate CCS from subreads (skip if using HiFi reads)
ccs input.subreads.bam ccs.bam \
    --min-rq 0.9 \
    --min-passes 3 \
    --num-threads 32

# For HiFi reads, CCS is already done
# Start directly from HiFi reads

Primer Removal with Lima

# Iso-Seq specific primer removal
lima ccs.bam primers.fasta demux.bam \
    --isoseq \
    --peek-guess \
    --num-threads 16

# Output: demux.primer_5p--primer_3p.bam
# Lima reports also contain demux statistics

# Check lima report
cat demux.lima.summary

Primer File Format

>primer_5p
AAGCAGTGGTATCAACGCAGAGTACATGGG
>primer_3p
AAGCAGTGGTATCAACGCAGAGTAC

Refine Full-Length Reads

# Remove polyA tails and concatemers
isoseq3 refine demux.primer_5p--primer_3p.bam primers.fasta refined.bam \
    --require-polya \
    --min-polya-length 20

# Output: refined.bam (full-length non-chimeric reads)
# Also: refined.filter_summary.json

# Check refinement stats
cat refined.filter_summary.json | jq

Cluster Into Isoforms

# Cluster similar transcripts
isoseq3 cluster refined.bam clustered.bam \
    --verbose \
    --use-qvs \
    --num-threads 32

# Output files:
# - clustered.bam: Clustered transcripts
# - clustered.hq_transcripts.fasta: High-quality consensus
# - clustered.lq_transcripts.fasta: Low-quality consensus
# - clustered.cluster_report.csv: Cluster membership

Align to Reference

# Map HQ transcripts to reference genome
minimap2 -ax splice:hq \
    -uf \
    --secondary=no \
    reference.fa \
    clustered.hq_transcripts.fasta \
    | samtools sort -o aligned.bam

samtools index aligned.bam

# For downstream analysis
pbmm2 align reference.fa clustered.bam aligned.bam \
    --preset ISOSEQ \
    --sort

Collapse Redundant Isoforms

# Collapse mapped transcripts
isoseq3 collapse aligned.bam collapsed.gff

# Output:
# - collapsed.gff: Collapsed transcript models
# - collapsed.abundance.txt: Read counts per isoform
# - collapsed.group.txt: Isoform groupings

# Convert to GTF
gffread collapsed.gff -T -o collapsed.gtf

SQANTI3 Quality Control

# Classify isoforms against reference annotation
sqanti3_qc.py \
    clustered.hq_transcripts.fasta \
    reference.gtf \
    reference.fa \
    -o sqanti_output \
    --aligner_choice minimap2 \
    --cage_peak cage_peaks.bed \
    --polyA_motif_list polyA_motifs.txt \
    --cpus 16

# Key output files:
# - sqanti_output_classification.txt: Per-isoform metrics
# - sqanti_output_junctions.txt: Splice junction details
# - sqanti_output.params.txt: Run parameters

SQANTI3 Categories

| Category | Code | Description | |----------|------|-------------| | Full Splice Match | FSM | All junctions match reference | | Incomplete Splice Match | ISM | Subset of reference junctions | | Novel In Catalog | NIC | Novel combination of known junctions | | Novel Not in Catalog | NNC | Contains novel junction | | Antisense | AS | Overlaps gene on opposite strand | | Genic | G | Within gene but no junction match | | Intergenic | IR | Between genes | | Fusion | FU | Spans multiple genes |

SQANTI3 Filtering

# Filter artifacts using SQANTI3 rules
sqanti3_filter.py \
    sqanti_output_classification.txt \
    --isoforms clustered.hq_transcripts.fasta \
    --gtf collapsed.gtf \
    --faa predicted_proteins.faa \
    -o sqanti_filtered

# Custom filtering
python << 'EOF'
import pandas as pd

classification = pd.read_csv('sqanti_output_classification.txt', sep='\t')

# Keep FSM, ISM, NIC with evidence
keep = classification[
    (classification['structural_category'].isin(['full-splice_match', 'incomplete-splice_match', 'novel_in_catalog'])) &
    (classification['FL'] >= 2) &
    (classification['bite'] == 'FALSE')
]
keep['isoform'].to_csv('filtered_isoforms.txt', index=False, header=False)
EOF

Quantification with Pigeon

# PacBio's isoform quantification tool
pigeon classify \
    aligned.bam \
    reference.gtf \
    reference.fa \
    -o pigeon_output

# Produces count matrix and classification
pigeon report pigeon_output_classification.txt

TAMA for Annotation Merge

# Merge Iso-Seq with reference annotation
# First, convert to TAMA format
tama_format_convert.py \
    -i collapsed.gtf \
    -f gtf \
    -o isoseq.bed

# Create file list
echo -e "isoseq.bed\tcapped\t1\t1" > file_list.txt
echo -e "reference.bed\tcapped\t1\t2" >> file_list.txt

# Merge annotations
tama_merge.py \
    -f file_list.txt \
    -p me