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

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.

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 fraction of single-exon models than the PacBio assemblies.
Compare the PacBio and ONT results in a genome browser alongside the reference annotation.

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

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

2.1 Create output directories

2.2 Assemble all PacBio reads together

2.3 Assemble each PacBio sample separately

2.4 Merge the per-sample PacBio StringTie assemblies

2.5 Locate the StringTie long-read GTF files

2.6 Count transcripts in each StringTie long-read assembly

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

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

2.9 Try yourself

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