long-read-sequencing/isoseq-analysis/SKILL.md
Discovers, classifies, filters, and quantifies full-length transcript isoforms from PacBio Iso-Seq/Kinnex (HiFi) and Oxford Nanopore (cDNA/direct-RNA) long reads, using the isoseq+pigeon pipeline, SQANTI3, and ONT tools (IsoQuant, FLAIR, Bambu, StringTie2). Covers why a novel isoform is an artifact until proven otherwise (RT template-switching, intra-priming, and 5' degradation manufacture junctions and truncations), the SQANTI3 structural categories and their trust order, the Kinnex skera-split step, orthogonal CAGE/poly-A/short-read-junction validation, and why long-read isoform quantification needs EM. Use when building a full-length isoform catalog, classifying/filtering long-read transcripts, running Iso-Seq or ONT cDNA/dRNA analysis, or judging novel-isoform reliability.
npx skillsauth add GPTomics/bioSkills bio-long-read-sequencing-isoseq-analysisInstall this skill globally with one command. Works with Claude Code, Cursor, and Windsurf.
3 of 9 scanners reported clean
Some scanners were skipped, did not run, or reported a non-clean status. Review each row below.
Reference examples tested with: isoseq 4.3+, pigeon 1.2+, SQANTI3 5.2+, pbmm2 1.13+, minimap2 2.28+, IsoQuant 3.4+.
Before using code patterns, verify installed versions match. If versions differ:
<tool> --version then <tool> --help to confirm flagspip show <pkg> / packageVersion('<pkg>') for SQANTI3/IsoQuant/BambuResults depend on inputs that outlive the binary version - record them:
isoseq3 -> isoseq in v4; the classifier pigeon is a separate binary.If code throws an error, introspect the installed tool (isoseq --help, pigeon --help, sqanti3_qc.py --help) and adapt the example to the actual API rather than retrying.
"Find the isoforms in my long-read RNA data" -> Build a full-length isoform catalog, then classify and filter it against the reference with orthogonal end/junction support - because discovery without curation is a catalog of artifacts.
isoseq refine ... && isoseq collapse ... && pigeon classify ... && pigeon filter ... (PacBio), IsoQuant/FLAIR/Bambu (ONT)RT template-switching, intra-priming on genomic poly-A, and 5' RNA degradation actively MANUFACTURE novel junctions and truncated isoforms. So the classification + filter + orthogonal validation IS the analysis, not a QC postscript. Invert the posture from "I discovered N novel isoforms" to "I curated N novel isoforms that survived artifact filtering." Three consequences:
Reference comparison is junction-chain based. NIC > NNC in trust, always; ISM is a diagnostic, not a discovery.
| Category (field value) | Meaning | Trust |
|------------------------|---------|-------|
| FSM (full-splice_match) | every internal junction matches a reference transcript; ends may differ | highest (known); ends still need CAGE/polyA |
| ISM (incomplete-splice_match) | junction subset of a reference (fewer 5' exons) | low - the 5'-degradation/RT-dropoff signature; trust only with CAGE |
| NIC (novel_in_catalog) | novel combination of KNOWN splice sites | high among novels - RT-switching cannot fake a NIC |
| NNC (novel_not_in_catalog) | >=1 genuinely novel splice site | lower - where junction artifacts concentrate; needs canonical/short-read support |
| genic / genic_intron | overlaps introns/exons; within an intron | low - pre-mRNA / gDNA carryover |
| fusion | spans >=2 genes | RT-chimera until proven by short-read split reads |
| intergenic / antisense | no gene overlap / antisense | novel-gene candidate or artifact; needs ORF/CAGE/conservation |
Mono-exon transcripts have no junctions to validate and are the false-discovery sink (intra-priming + gDNA run unchecked) - require ORF + CAGE + polyA + conservation before belief.
| Data / goal | Tool | Why |
|-------------|------|-----|
| PacBio Iso-Seq/Kinnex, turnkey | isoseq + pigeon | native PacBio collapse + SQANTI-style classify/filter, SMRT Link integrated |
| Any long-read transcriptome, full curation | SQANTI3 | structural classification + ~50 QC descriptors + rules/ML filter + rescue; PacBio and ONT |
| ONT bulk discovery + quantification | IsoQuant | intron-graph; lowest novel FP rate among ONT tools |
| ONT, want built-in differential splicing | FLAIR | align -> correct junctions -> collapse -> diffSplice |
| Quantification with a precision knob | Bambu | NDR (novel discovery rate) calibrates precision; R/Bioconductor |
| Genome-guided assembly / hybrid short+long | StringTie2 -L (--mix) | fast long-read transcript assembly |
| ONT single-cell long-read isoforms | FLAMES | single-cell/spatial full-length isoforms |
| Differential isoform usage (DTU/DTE) | -> alternative-splicing | this skill yields the filtered set + counts and hands off |
PacBio Iso-Seq and ONT cDNA sequence reverse-transcribed cDNA (modifications erased; strand from primers); ONT direct-RNA sequences native RNA (true strand, poly-A length, modifications preserved, lower accuracy). Match the minimap2 preset to the chemistry:
minimap2 -ax splice ref.fa ont_cdna.fq # ONT cDNA (orient first with pychopper)
minimap2 -ax splice -uf -k14 ref.fa drna.fq # ONT direct RNA (stranded -> -uf, small k)
minimap2 -ax splice:hq -uf ref.fa hifi.fa # PacBio HiFi (or pbmm2 --preset ISOSEQ)
-uf forces the forward transcript strand - correct for stranded dRNA/Iso-Seq, wrong for unoriented ONT PCR-cDNA (orient with pychopper first).
# 0. Kinnex (MAS-seq) ONLY: deconcatenate the array into segmented reads FIRST
skera split movie.hifi_reads.bam mas_adapters.fasta movie.segmented.bam # skip for classic Iso-Seq
# 1. Remove cDNA primers; 2. produce FLNC (full-length non-chimeric)
lima movie.segmented.bam primers.fasta movie.fl.bam --isoseq --peek-guess
isoseq refine movie.fl.5p--3p.bam primers.fasta movie.flnc.bam --require-polya
# 3. cluster (reference-free) or skip and align FLNC directly; 4. map; 5. collapse to isoforms
isoseq cluster2 movie.flnc.bam clustered.bam # cluster2 scales to large sets
pbmm2 align --preset ISOSEQ --sort ref.fa clustered.bam mapped.bam
isoseq collapse --do-not-collapse-extra-5exons mapped.bam movie.flnc.bam collapsed.gff
# collapsed.flnc_count.txt = FLNC molecules per isoform = the real DEPTH metric
# 6. classify + filter with pigeon (needs the collapsed.sorted.gff after prepare, NOT a BAM)
pigeon prepare collapsed.gff # sorts the transcript GFF
pigeon prepare annotation.gtf ref.fa # sorts the annotation -> annotation.sorted.gtf, indexes genome
pigeon classify collapsed.sorted.gff annotation.sorted.gtf ref.fa \
--fl collapsed.flnc_count.txt --cage-peak cage.refTSS.bed --poly-a polyA.motif.list
pigeon filter collapsed_classification.txt --isoforms collapsed.sorted.gff
pigeon report --exclude-singletons collapsed_classification.filtered_lite_classification.txt saturation.txt
pigeon is PacBio's productized SQANTI3 (classify/filter, NOT a quantifier). Substitute SQANTI3 itself for the full descriptor set, ML filter, rescue module, and ONT support:
sqanti3_qc.py collapsed.gff annotation.gtf ref.fa --CAGE_peak cage.bed --polyA_motif_list polyA.txt \
--short_reads short_reads_fofn.txt # isoforms positional defaults to GTF/GFF; add --fasta for FASTA input
sqanti3_filter.py rules collapsed_classification.txt # or: sqanti3_filter.py ml ...
Trigger: reporting incomplete-splice-match transcripts as discoveries. Mechanism: 5' RNA degradation truncates FSMs into ISMs. Symptom: inflated novel/ISM fraction tracking RNA quality, not biology. Fix: treat ISM fraction as an integrity QC; keep ISMs only with CAGE 5' support.
Trigger: trusting 3' ends without poly-A validation. Mechanism: oligo-dT primes on a genomic internal A-stretch. Symptom: spurious short/mono-exon transcripts; perc_A_downstream_TTS >59%. Fix: SQANTI3/pigeon filter on downstream genomic A-content and poly-A motif; --require-polya alone does NOT catch this.
Trigger: treating NNC like NIC. Mechanism: novel splice sites are where RT template-switching and mapping artifacts land. Symptom: novel junctions absent from short-read data. Fix: require canonical junctions or short-read SJ coverage; prefer NIC.
Trigger: pigeon classify mapped.bam .... Mechanism: pigeon classifies the collapsed.sorted.gff after pigeon prepare, not an alignment. Symptom: wrong-input error. Fix: isoseq collapse -> pigeon prepare -> pigeon classify.
Trigger: raw isoform counts as abundance. Mechanism: discovery is depth-unsaturated; truncated reads are multi-isoform-compatible. Symptom: deeper libraries "have more isoforms"; double-counted abundance. Fix: rarefaction curve (--exclude-singletons); EM quantification (Bambu/IsoQuant/NanoCount), not raw FLNC counts.
| Threshold | Source | Rationale |
|-----------|--------|-----------|
| perc_A_downstream_TTS > 59-60% = intra-priming | SQANTI (Tardaguila 2018) | genomic A-rich window means the poly-A was internal, not the real tail |
| novel junction trusted if canonical OR short-read cov >= 3 | SQANTI3 rules filter | a single criterion for RT-switch/NNC artifacts |
| ML filter needs >= 250 Reference-Match FSM | SQANTI3 | enough true-positive labels to train; else falls back to rules |
| exclude singletons (1-FLNC) for saturation | pigeon report | singletons are the dominant unreliable novel bucket |
| FLNC count = depth metric | isoseq collapse | independently sequenced full-length molecules, before clustering/dedup |
| Error / symptom | Cause | Solution |
|-----------------|-------|----------|
| isoseq3: command not found | renamed in v4 | use isoseq (subcommands unchanged) |
| pigeon classify wrong input | fed a BAM | give the collapsed.sorted.gff after pigeon prepare |
| Huge novel-isoform count | filter skipped/underpowered | run pigeon/SQANTI3 filter with CAGE/polyA/short-read support |
| Many mono-exon novels | intra-priming / gDNA carryover | filter on poly-A; require ORF/CAGE for mono-exon |
| Wrong-strand spliced alignment | -uf on unoriented cDNA | orient with pychopper, or drop -uf for cDNA |
| Isoform counts not comparable across samples | depth-unsaturated discovery | EM quantification + rarefaction curve |
-uf)development
Installs 425 bioinformatics skills covering sequence analysis, RNA-seq, single-cell, variant calling, metagenomics, structural biology, and 56 more categories. Use when setting up bioinformatics capabilities or when a bioinformatics task requires specialized skills not yet installed.
testing
Chains a somatic (tumor-normal) SNV/indel and structural-variant pipeline end to end with GATK Mutect2 (or Strelka2), wiring the somatic-specific machinery - panel-of-normals and gnomAD germline-resource priors, GetPileupSummaries/CalculateContamination, and LearnReadOrientationModel FFPE/oxoG orientation-bias filtering fed into FilterMutectCalls. Use when calling somatic mutations from a tumor-normal pair (or tumor-only with PoN caveats), deciding which artifact filter removes which class of false positive, reasoning about VAF/purity/ploidy and clonal-vs-subclonal detection, adding somatic SV/CNV or TMB/MSI/signatures, or routing variants to AMP/ASCO/CAP tier and oncogenicity interpretation (never germline ACMG).
development
End-to-end pooled and single-cell CRISPR screen analysis from FASTQ to hit genes. Orchestrates library design QC, guide counting, six-stage screen QC (plasmid Gini, replicate Pearson, CEGv2 PR-AUC, copy-number artifact), method-appropriate hit calling across MAGeCK RRA/MLE, BAGEL2, drugZ, JACKS, and Chronos, cancer-cell-line copy-number correction (CRISPRcleanR / Chronos), batch correction for multi-batch screens, and the specialized branches for combinatorial paralog screens, single-cell Perturb-seq, base-editor variant-function screens, prime-editor screens, and in vivo bottleneck-aware screens. Use when analyzing any pooled CRISPR screen end-to-end, matching the hit-calling method to the experimental design, integrating copy-number correction into the pipeline, or branching the workflow for single-cell, combinatorial, base-editor, prime-editor, or in vivo variants.
development
Transcribe DNA to RNA and translate to protein using Biopython, with NCBI codon-table selection, CDS validation, and six-frame ORF finding. Use when converting a CDS or ORF to its amino-acid sequence, selecting a non-standard (mitochondrial, bacterial, ciliate) genetic code, validating a coding sequence, or scanning all reading frames.