Long‐read transcript assembly with IsoQuant and StringTie

These commands demonstrate a simple long-read transcript assembly workflow using IsoQuant or StringTie on PacBio and Oxford Nanopore reads from a small human chromosome 20 region.

Aim of this exercise: assemble long-read transcript models from PacBio and ONT data, compare annotation-guided and annotation-free runs, and inspect how the two sequencing platforms differ in transcript recovery and novelty.

---

1.1 Set reference variables

cd /home/train/Annotation_workshop/Transcriptome_assembly

export human_genome=Inputs/Reference/Homo_sapiens.GRCh38.dna.primary_assembly_chr20regionA.fa
export human_gff=Inputs/Ref_Annotation/gencode.v39.annotation.regionA.gff

Notes

• human_genome stores the human region FASTA path.
• human_gff stores the matching GENCODE annotation for the same region.
• Using shell variables makes later commands easier to read and modify.

---

1.2 Check IsoQuant help

isoquant --help

Notes

• Useful for confirming that isoquant.py is on your PATH.
• Check out the various options

---

1.3 PacBio full-length reads, annotation-guided

(
mkdir -p IsoQuant/PacBio_ref
isoquant \
    --data_type pacbio \
    -o IsoQuant/PacBio_ref \
    -r $human_genome \
    -g $human_gff \
    --complete_genedb \
    --fastq \
        Inputs/Reads/sample_1_pacbio.subsample_chr20regionA.fastq \
        Inputs/Reads/sample_2_pacbio.subsample_chr20regionA.fastq \
        Inputs/Reads/sample_3_pacbio.subsample_chr20regionA.fastq \
        Inputs/Reads/sample_4_pacbio.subsample_chr20regionA.fastq \
        Inputs/Reads/sample_5_pacbio.subsample_chr20regionA.fastq \
        Inputs/Reads/sample_6_pacbio.subsample_chr20regionA.fastq \
    --fl_data \
    --stranded forward \
    --prefix pacbio_ref \
    --polya_trimmed all
)

Notes

• --fl_data tells IsoQuant to treat the reads as full-length transcript reads.
• -g $human_gff --complete_genedb provides the GENCODE annotation.
• --polya_trimmed all assumes all reads are already polyA-trimmed and properly oriented.

Try yourself

• Check out the IsoQuant outputs l -R IsoQuant/ see Output files for details
• IsoQuant aligns the reads with minimap check the bam alignments with samtools view -H to see the alignment options used

---

1.4 PacBio full-length reads, annotation-free

(
mkdir -p IsoQuant/PacBio_denovo
isoquant \
    --data_type pacbio \
    -o IsoQuant/PacBio_denovo \
    -r $human_genome \
    --fastq \
        Inputs/Reads/sample_1_pacbio.subsample_chr20regionA.fastq \
        Inputs/Reads/sample_2_pacbio.subsample_chr20regionA.fastq \
        Inputs/Reads/sample_3_pacbio.subsample_chr20regionA.fastq \
        Inputs/Reads/sample_4_pacbio.subsample_chr20regionA.fastq \
        Inputs/Reads/sample_5_pacbio.subsample_chr20regionA.fastq \
        Inputs/Reads/sample_6_pacbio.subsample_chr20regionA.fastq \
    --fl_data \
    --stranded forward \
    --prefix pacbio_denovo \
    --polya_trimmed all
)

Notes

• This run omits the annotation, so IsoQuant will build transcript models from the long-read evidence and genome alone.
• This is useful for comparing how much the annotation influences transcript discovery.

---

1.5 ONT reads, annotation-guided

(
mkdir -p IsoQuant/ONT_ref
isoquant \
    --data_type nanopore \
    -o IsoQuant/ONT_ref \
    -r $human_genome \
    -g $human_gff \
    --complete_genedb \
    --fastq \
        Inputs/Reads/HC1_ont.subsample_chr20regionA.fastq \
        Inputs/Reads/HC3_ont.subsample_chr20regionA.fastq \
        Inputs/Reads/HC4_ont.subsample_chr20regionA.fastq \
        Inputs/Reads/HT1_ont.subsample_chr20regionA.fastq \
        Inputs/Reads/HT5_ont.subsample_chr20regionA.fastq \
        Inputs/Reads/HT7_ont.subsample_chr20regionA.fastq \
    --stranded forward \
    --prefix ont_ref \
    --polya_trimmed none
)

Notes

• --polya_trimmed none means polyA trimming is not assumed to have been done already.
• This uses the reference annotation to guide transcript assignment and comparison.

---

1.6 ONT reads, annotation-free

(
mkdir -p IsoQuant/ONT_denovo
isoquant \
    --data_type nanopore \
    -o IsoQuant/ONT_denovo \
    -r $human_genome \
    --fastq \
        Inputs/Reads/HC1_ont.subsample_chr20regionA.fastq \
        Inputs/Reads/HC3_ont.subsample_chr20regionA.fastq \
        Inputs/Reads/HC4_ont.subsample_chr20regionA.fastq \
        Inputs/Reads/HT1_ont.subsample_chr20regionA.fastq \
        Inputs/Reads/HT5_ont.subsample_chr20regionA.fastq \
        Inputs/Reads/HT7_ont.subsample_chr20regionA.fastq \
    --stranded forward \
    --prefix ont_denovo \
    --polya_trimmed none
)

Notes

• This is the ONT equivalent of the de novo PacBio run.
• It allows comparison of annotation-guided versus annotation-free assembly.

---

1.7 Useful output files to inspect

find IsoQuant -maxdepth 3 -type f | less

Notes

• This gives a quick overview of the files written by each IsoQuant run.
• The output directory is controlled by -o and file prefixes by --prefix.

---

1.8 Compare the final GTF files

find IsoQuant -name "*.gtf" | sort

Notes

• This helps locate the transcript-model outputs from each run.
• A simple first comparison is annotation-guided versus de novo, within PacBio and ONT separately.
• Note the annotation-guided runs contain an additional SAMPLE_ID.extended_annotation.gtf GTF file with the entire reference annotation plus all discovered novel transcripts.

---

1.9 Quick transcript counts

(
find IsoQuant -name "*.gtf" | sort | while read -r file; do
    printf "%s\n" "$file"
    cut -f3 "$file" |sort -n |uniq -c |grep -v '#' |sort -n;
done
)

Try yourself

• Does annotation-guided mode report fewer or more transcript models than de novo mode?
• Does ONT produce more novel models than PacBio on this toy dataset?
• Are any differences driven mainly by mono-exonic predictions?

To list each IsoQuant GTF together with the total number of transcripts and the number of single-exon and multi-exon transcripts:

(
find IsoQuant -name "*.gtf" | sort | while read -r file; do
    printf "%s\n" "$file"
    printf "transcripts\t"
    grep -c $'\ttranscript\t' "$file"
    awk -F'\t' '
    $0 !~ /^#/ && $3=="exon" {
        n = split($9, x, ";")
        tid = ""
        for (i=1; i<=n; i++) {
            gsub(/^[ \t]+|[ \t]+$/, "", x[i])
            if (x[i] ~ /^transcript_id /) {
                sub(/^transcript_id "/, "", x[i])
                sub(/"$/, "", x[i])
                tid = x[i]
                break
            }
        }
        if (tid != "") tx[tid]++
    }
    END {
        single=0
        multi=0
        for (t in tx) {
            if (tx[t]==1) single++
            else if (tx[t]>1) multi++
        }
        print "single_exon_transcripts\t" single
        print "multi_exon_transcripts\t" multi
    }' "$file"
    echo
done
)

---

1.10 Optional: compare IsoQuant assemblies to the reference with Mikado Compare

First convert the IsoQuant output with gffread -T. This removes gene feature lines that do not contain a transcript_id, which can otherwise cause mikado compare to crash.

(
find IsoQuant -name "*.gtf" ! -name "*.gffread.gtf" | sort | while read -r file; do
    out="${file%.gtf}.gffread.gtf"
    gffread -T -o "$out" "$file"
    echo "$(basename "$file") -> $(basename "$out") DONE"
done
)

PacBio:

(
mkdir -p IsoQuant/Mikado_Compare
for file in IsoQuant/{PacBio_ref/pacbio_ref,PacBio_denovo/pacbio_denovo}/*gffread.gtf; do
    echo -e "\n\n## $file"
    outfile=$(basename ${file} .gtf)
    mikado compare \
        -r $human_gff \
        -p $file \
        -o IsoQuant/Mikado_Compare/mikado_compare_ref_v_${outfile}
done
)

ONT:

(
mkdir -p IsoQuant/Mikado_Compare
for file in IsoQuant/{ONT_ref/ont_ref,ONT_denovo/ont_denovo}/*gffread.gtf; do
    echo -e "\n\n## $file"
    outfile=$(basename ${file} .gtf)
    mikado compare \
        -r $human_gff \
        -p $file \
        -o IsoQuant/Mikado_Compare/mikado_compare_ref_v_${outfile}
done
)

To compare transcript-level sn, pr and F1 quickly:

(
Inputs/Scripts/mikado_compare_stats_pivot.py \
    --rows t80 \
    --label_regex '(?<=mikado_compare_ref_v_)[^/]+(?=\.stats$)' \
    -o isoquant_compare_ref.tsv \
    IsoQuant/Mikado_Compare/mikado_compare_ref_*.stats
)

columntab < isoquant_compare_ref.tsv | less -S

Notes

• ref.extended_annotation.gtf models include the entire reference annotation plus all discovered novel transcripts so sn is expectd to be 100%
• As expected providing the reference annotation to IsoQuant improves reconstruction i.e. compare denovo.transcript_models to ref.transcript_models

label	Sn	Pr	F1
ont_denovo.transcript_models	3.72	38.61	6.79
ont_ref.extended_annotation	100.00	96.32	98.13
ont_ref.transcript_models	11.55	75.00	20.01
pacbio_denovo.transcript_models	5.63	12.37	7.74
pacbio_ref.extended_annotation	100.00	72.18	83.84
pacbio_ref.transcript_models	8.02	17.21	10.94

---

2 Long-read transcript assembly with StringTie from minimap2 BAM files

These commands demonstrate a simple long-read transcript assembly workflow using StringTie on precomputed minimap2-aligned BAM files. The BAM files in Inputs/Alignments/ were generated elsewhere and can be used directly as input to StringTie, provided they are coordinate-sorted. StringTie uses alignment files as input, and -L enables long-read processing mode. (github.com)

Aim of this exercise: assemble transcripts from PacBio long-read alignments in three ways: using all PacBio BAMs together in one pooled run, assembling each sample separately, and merging the per-sample assemblies with stringtie --merge.

Check the available long-read BAM files

ls -lh Inputs/Alignments/*pacbio*.bam

Notes

• These BAM files are the long-read alignments that will be used as input to StringTie.
• StringTie expects coordinate-sorted RNA-seq alignments in BAM, SAM or CRAM format.
• For long-read transcript assembly, StringTie should be run with -L.

Check StringTie help

stringtie --help

Notes

• The -L option enables long-read processing mode.
• The -G option (not used in this example) can be used to guide assembly with a reference annotation.
• The --mix option (not used in this example) can be used with mixed reads i.e. both short and long read data alignments are expected
• The --merge mode combines multiple GTF files into a non-redundant set of transcripts.

---

2.1 Create output directories

mkdir -p StringTie_LongRead/{PacBio_all,PacBio_by_sample,PacBio_merge}

Notes

• This creates separate directories for the pooled assembly, per-sample assemblies, and merged assembly.

---

2.2 Assemble all PacBio reads together

(
stringtie -L \
    -l PBA \
    Inputs/Alignments/sample_1_pacbio.subsample_chr20regionA.aligned.bam \
    Inputs/Alignments/sample_2_pacbio.subsample_chr20regionA.aligned.bam \
    Inputs/Alignments/sample_3_pacbio.subsample_chr20regionA.aligned.bam \
    Inputs/Alignments/sample_4_pacbio.subsample_chr20regionA.aligned.bam \
    Inputs/Alignments/sample_5_pacbio.subsample_chr20regionA.aligned.bam \
    Inputs/Alignments/sample_6_pacbio.subsample_chr20regionA.aligned.bam \
    -o StringTie_LongRead/PacBio_all/pacbio_all.gtf
)

Notes

• -L enables long-read assembly mode in StringTie.
• -l PBA gives the assembled transcript IDs a distinctive pooled PacBio prefix.
• Multiple BAM files can be supplied directly for one pooled assembly, so there is no need to merge the BAMs first with samtools merge.
• This produces one combined transcript assembly from all PacBio alignments together.

---

2.3 Assemble each PacBio sample separately

(
for bam in Inputs/Alignments/sample_*_pacbio.subsample_chr20regionA.aligned.bam; do
    sample=$(basename "$bam" .subsample_chr20regionA.aligned.bam)
    mkdir -p StringTie_LongRead/PacBio_by_sample/${sample}
    stringtie -L \
        -l "${sample}_PBA" \
        "$bam" \
        -o StringTie_LongRead/PacBio_by_sample/${sample}/${sample}.gtf;
    echo "${sample}_PBA" DONE
done
)

Notes

• This runs one StringTie assembly per PacBio sample.
• -l "${sample}_PBA" gives each sample-specific assembly a distinctive transcript ID prefix.
• Keeping separate per-sample assemblies makes it possible to compare sample-specific transcript structures before merging.

---

2.4 Merge the per-sample PacBio StringTie assemblies

(
stringtie --merge \
    -l PBM \
    -o StringTie_LongRead/PacBio_merge/pacbio_merged.gtf \
    StringTie_LongRead/PacBio_by_sample/sample_*_pacbio/sample_*_pacbio.gtf
)

Notes

• stringtie --merge combines transcript models from multiple GTF files into one merged transcript set.
• -l PBM gives the merged transcript IDs a distinctive prefix.
• This is different from the pooled assembly above: here, each sample is assembled first, and the resulting GTF files are merged afterwards.
• Comparing the pooled assembly with the merged per-sample assembly is a useful way to explore how these approaches differ.

---

2.5 Locate the StringTie long-read GTF files

find StringTie_LongRead -name "*.gtf" | sort

Notes

• This lists the pooled assembly, all per-sample assemblies, and the merged assembly.

---

2.6 Count transcripts in each StringTie long-read assembly

(
find StringTie_LongRead -name "*.gtf" | sort | while read -r file; do
    printf "%s\t" "$file"
    grep -c $'\ttranscript\t' "$file"
done
)

Notes

• This provides a quick comparison of transcript counts across the different assembly strategies.

Try yourself

• How do the counts for pacbio_all.gtf and pacbio_merged.gtf compare?

---

2.7 Count single-exon and multi-exon transcripts in each StringTie long-read assembly

(
find StringTie_LongRead -name "*.gtf" | sort | while read -r file; do
    printf "%s\n" "$file"
    printf "transcripts\t"
    grep -c $'\ttranscript\t' "$file"
    awk -F'\t' '
    $0 !~ /^#/ && $3=="exon" {
        n = split($9, x, ";")
        tid = ""
        for (i=1; i<=n; i++) {
            gsub(/^[ \t]+|[ \t]+$/, "", x[i])
            if (x[i] ~ /^transcript_id /) {
                sub(/^transcript_id "/, "", x[i])
                sub(/"$/, "", x[i])
                tid = x[i]
                break
            }
        }
        if (tid != "") tx[tid]++
    }
    END {
        single=0
        multi=0
        for (t in tx) {
            if (tx[t]==1) single++
            else if (tx[t]>1) multi++
        }
        print "single_exon_transcripts\t" single
        print "multi_exon_transcripts\t" multi
    }' "$file"
    echo
done
)

Notes

• This helps distinguish whether differences in transcript count are mainly due to single-exon or multi-exon models.

Try Yourself

• Check out pacbio_all.gtf and pacbio_merged.gtf, try rerunning the merging step with different --merge options see stringtie -h

---

2.8 Optional: compare StringTie long-read assemblies to the reference annotation

(
mkdir -p StringTie_LongRead/Mikado_Compare
for file in $(find StringTie_LongRead -name "*.gtf" | sort); do
    echo -e "\n\n## $file"
    outfile=$(basename "${file}" .gtf)
    mikado compare \
        -r Inputs/Ref_Annotation/gencode.v39.annotation.regionA.gtf \
        -p "$file" \
        -o StringTie_LongRead/Mikado_Compare/mikado_compare_ref_v_${outfile}
done
)

Notes

• This compares the pooled, per-sample, and merged StringTie assemblies against the reference annotation.
• It can be useful for checking whether the pooled assembly or the merged per-sample assembly gives better agreement with the reference.

---

2.9 Try yourself

• Repeat the same workflow using the ONT BAM files in Inputs/Alignments/.
• Compare transcript counts between the pooled ONT assembly and the merged per-sample ONT assembly.
• Check whether ONT assemblies contain a higher or lower fraction of single-exon models than the PacBio assemblies.
• Compare the sn, pr and F1 for stringtie to the isoquant denovo models

Final Note

The following sn pr plots compare stringtie and IsoQuant assemblies (in denovo mode i.e. no reference annotation provided) for the pacbio samples, note this is from the full dataset not just the test region and the reference annotation was filtered to remove models with 0 nF1 i.e. not even partially reconstructable from the sample data. Stringtie | cat All is a simple concatenation of the indivdual sample files and is included to show the maximum achievable sn, this will contain redundant models and would need further filtering to be of practical use. Note on the full dataset removing single exon (monoexonic) models from the predictions and reference improves Isoquant performance see isoquant -h for options relating to this specifically

--report_novel_unspliced, -u REPORT_NOVEL_UNSPLICED"> report novel monoexonic transcripts (true/false), default: false for ONT, true for other data types

The --merge All run was with default parameters and modifying these will potentially improve sensitivity with likely a descrease in precision

A strategy of cherry picking from across all samples OR all samples + the combined runs potentially can improve sn and pr. See Mikado