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

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 annotation-guided:

(
mkdir -p IsoQuant/Mikado_Compare
for file in IsoQuant/PacBio_ref/pacbio_ref/*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 annotation-guided:

(
mkdir -p IsoQuant/Mikado_Compare
for file in IsoQuant/ONT_ref/ont_ref/*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

1.11 Other things to try

Re-run the annotation-guided commands with extra reporting options if you want known-vs-novel style output.
Add --threads 4 or another suitable value for faster runs.
Compare ont_ref vs pacbio_ref in a genome browser alongside gencode.v39.annotation.regionA.gtf.
Explore whether annotation-guided mode improves precision compared with de novo assembly.

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

These commands demonstrate how to assemble long-read transcript models with StringTie using the minimap2-derived BAM files produced during an IsoQuant-style workflow.

Aim of this exercise: compare three related strategies for long-read transcript reconstruction using PacBio alignments: assemble all reads together in one run, assemble each sample separately, and merge the per-sample assemblies into a non-redundant transcript set.

2.1 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.

2.2 Check StringTie help

stringtie --help

Notes

The -L option enables long-read processing mode.
The -G option can be used to guide assembly with a reference annotation.
The --merge mode combines multiple GTF files into a non-redundant set of transcripts.

2.3 Build a single BAM containing all PacBio reads

(
mkdir -p StringTie_LongRead/PacBio/all_reads
samtools merge -f \
    StringTie_LongRead/PacBio/all_reads/pacbio_all_reads.bam \
    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 && \
samtools index StringTie_LongRead/PacBio/all_reads/pacbio_all_reads.bam
)

Notes

This creates one combined PacBio BAM across all six samples.
Merging alignments lets you see what StringTie produces when all long-read evidence is assembled together in a single run.
The BAM is indexed so it can also be inspected later in a genome browser.

2.4 Assemble all PacBio reads together with StringTie long-read mode

(
mkdir -p StringTie_LongRead/PacBio/all_reads
stringtie \
    -L \
    -G $human_gff \
    -o StringTie_LongRead/PacBio/all_reads/pacbio_all_reads.stringtie.gtf \
    -l PACBIO_ALL \
    StringTie_LongRead/PacBio/all_reads/pacbio_all_reads.bam \
    > StringTie_LongRead/PacBio/all_reads/stringtie.log 2>&1
)

Notes

-L turns on long-read mode.
-G $human_gff guides the assembly with the human reference annotation.
-l PACBIO_ALL gives the output transcript IDs a distinctive prefix.

2.5 Assemble each PacBio sample separately with StringTie long-read mode

(
for bam in Inputs/Alignments/*pacbio*.bam; do
    sample=$(basename "$bam" .subsample_chr20regionA.aligned.bam)
    echo -e "\n\n## Running StringTie long-read assembly on $sample"
    mkdir -p StringTie_LongRead/PacBio/${sample}
    stringtie \
        -L \
        -G $human_gtf \
        -l ${sample}_STRG \
        -o StringTie_LongRead/PacBio/${sample}/${sample}.gtf \
        "$bam" \
        > StringTie_LongRead/PacBio/${sample}/stringtie.log 2>&1
done
)

Notes

This runs one assembly per PacBio BAM file.
Sample-specific labels are helpful because otherwise all assemblies would share the default STRG transcript prefix.
Per-sample assemblies make it easier to compare transcript content between individual samples before merging.

2.6 Collect the per-sample GTF files

(
mkdir -p StringTie_LongRead/PacBio/Assemblies
for gtf in StringTie_LongRead/PacBio/*/*.gtf; do
    sample=$(basename "$gtf")
    cp "$gtf" StringTie_LongRead/PacBio/Assemblies/${sample}
done
)

Notes

This copies the per-sample GTF files into a single directory.
Keeping them together makes the merge command easier to read and modify.

2.7 Merge the per-sample PacBio StringTie assemblies

(
mkdir -p StringTie_LongRead/PacBio/Merge
stringtie \
    --merge \
    -G $human_gff \
    -l PACBIO_MERGE \
    -o StringTie_LongRead/PacBio/Merge/pacbio_merged.stringtie.gtf \
    StringTie_LongRead/PacBio/Assemblies/*.gtf
)

Notes

--merge combines the input GTF files into a non-redundant transcript set.
Providing -G $human_gff includes the reference annotation during merge mode.
This merged assembly can be compared with the single combined-BAM assembly from section 9.4.

2.8 Quick transcript counts for the three PacBio strategies

grep -c $'\ttranscript\t' StringTie_LongRead/PacBio/all_reads/pacbio_all_reads.stringtie.gtf
grep -c $'\ttranscript\t' StringTie_LongRead/PacBio/*/sample_*.gtf
grep -c $'\ttranscript\t' StringTie_LongRead/PacBio/Merge/pacbio_merged.stringtie.gtf

Notes

This gives a quick first comparison of how many transcript models are produced by each approach.
Counts alone are not enough to judge quality, but they are a useful starting point.

2.9 Optional: compare the three PacBio StringTie outputs to the reference

(
mkdir -p StringTie_LongRead/PacBio/Mikado_Compare
for file in \
    StringTie_LongRead/PacBio/all_reads/pacbio_all_reads.stringtie.gtf \
    StringTie_LongRead/PacBio/Merge/pacbio_merged.stringtie.gtf \
    StringTie_LongRead/PacBio/*/sample_*.gtf; do
    echo -e "\n\n## $file"
    outfile=$(basename "$file" .gtf)
    mikado compare \
        -r $human_gff \
        -p "$file" \
        -o StringTie_LongRead/PacBio/Mikado_Compare/mikado_compare_ref_v_${outfile}
done
)

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

columntab < stringtie_longread_pacbio_compare_ref.tsv | less -S

Notes

This allows you to compare transcript-level sensitivity, precision and F1 across the different PacBio StringTie strategies.
It is particularly useful for asking whether a merged set improves recovery over the per-sample or pooled-BAM assemblies.

2.10 Try yourself: repeat the StringTie long-read workflow for ONT

Useful starting points:

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

stringtie -L -G $human_gff -o out.gtf some_ont.aligned.bam

Questions to explore:

Does ONT produce more transcript models than PacBio?
Does assembling all ONT reads together behave differently from merging per-sample ONT assemblies?
Are the ONT assemblies more likely to include novel or fragmented transcript models?
Does the merged ONT assembly improve transcript-level sensitivity, precision or F1 when compared with the reference?

2.11 Other things to try

Re-run the StringTie assemblies without -G and compare annotation-guided versus annotation-free long-read assembly.
Compare the PacBio StringTie output with the PacBio IsoQuant output for the same region.
Load the PacBio BAM files, the StringTie GTFs and gencode.v39.annotation.regionA.gtf into a genome browser and inspect a few loci by eye.
Try changing StringTie merge settings such as -f, -c, -F or -T and see whether the balance between sensitivity and precision changes.

Uh oh!

Long‐read transcript assembly with IsoQuant and StringTie

1.1 Set reference variables

1.2 Check IsoQuant help

1.3 PacBio full-length reads, annotation-guided

1.4 PacBio full-length reads, annotation-free

1.5 ONT reads, annotation-guided

1.6 ONT reads, annotation-free

1.7 Useful output files to inspect

1.8 Compare the final GTF files

1.9 Quick transcript counts

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

1.11 Other things to try

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

2.1 Check the available long-read BAM files

2.2 Check StringTie help

2.3 Build a single BAM containing all PacBio reads

2.4 Assemble all PacBio reads together with StringTie long-read mode

2.5 Assemble each PacBio sample separately with StringTie long-read mode

2.6 Collect the per-sample GTF files

2.7 Merge the per-sample PacBio StringTie assemblies

2.8 Quick transcript counts for the three PacBio strategies

2.9 Optional: compare the three PacBio StringTie outputs to the reference

2.10 Try yourself: repeat the StringTie long-read workflow for ONT

2.11 Other things to try

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Clone this wiki locally