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

⚠️ Windows users:

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
)

⚠️ Windows users:

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
)

⚠️ Windows users:

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
)

⚠️ Windows users:

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
)

⚠️ Windows users:

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
)

⚠️ Windows users:

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?
Check out isoquant option --report_novel_unspliced and note the default difference in behaviour for ONT and PacBio data

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
)

⚠️ Windows users:

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
)

⚠️ Windows users:

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
)

⚠️ Windows users:

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
)

⚠️ Windows users:

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
)

⚠️ Windows users:

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
)

⚠️ Windows users:

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
)

⚠️ Windows users:

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
)

⚠️ Windows users:

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
)

⚠️ Windows users:

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
)

⚠️ Windows users:

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
)

⚠️ Windows users:

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

To clean up to just the original files

cd /home/train/Annotation_workshop/Transcriptome_assembly
rm -rf -v !("Inputs"|"Example_output")

⚠️ Windows users:

cd /home/train/Annotation_workshop/Transcriptome_assembly && rm -rf -v !("Inputs"|"Example_output")

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

Mono_multi_exon_Exon_level_lenient precision_sensitivity

Mono_multi_exon_Gene_level_80_base_F1 precision_sensitivity

Only_multi_exon_Gene_level_80_base_F1 precision_sensitivity

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

To clean up to just the original files

Final Note

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