Skip to content

Transcript assembly commands

Jump to bottom
swarbred edited this page Apr 14, 2026 · 26 revisions

These commands demonstrate a simple RNA-seq transcript assembly workflow using HISAT2, StringTie, Scallop, Mikado Compare, and Portcullis.

Aim of this exercise: map RNA-seq reads to the genome, assemble transcripts with two different tools, compare those assemblies to each other and to a reference annotation, and explore how splice-junction filtering affects transcript evidence.

1 Set the genome environment variable

cd /home/train/Annotation_workshop/Transcriptome_assembly
export genome=Inputs/Reference/Athaliana_447_TAIR10_Chr3_clean.fa

Notes

  • cd moves into the expected working directory.
  • export genome=... stores the genome path in a shell variable so it can be reused later.
  • This makes later commands easier to read and edit.

2 Build Hisat2 index

hisat2-build -h
(
rm -rf Hisat2 Stringtie Scallop Assemblies Mikado_Compare Portcullis mikado_compare* *portcullis*
mkdir -p Hisat2/db
hisat2-build \
    $genome \
    Hisat2/db/Athaliana_447_TAIR10_Chr3_clean \
    > Hisat2/index.log 2>&1
)

Notes

  • rm -rf ... removes previous outputs so the exercise starts cleanly.
  • hisat2-build creates an index from the genome FASTA.
  • The index basename here is Hisat2/db/Athaliana_447_TAIR10_Chr3_clean.
  • The log is written to Hisat2/index.log.

3 Perform HISAT2 alignment

3.1 Sample SRR5956436

(
mkdir -p Hisat2/SRR5956436
hisat2 \
    --max-intronlen 50000 \
    --rna-strandness RF \
    --dta -p 2 \
    -x Hisat2/db/Athaliana_447_TAIR10_Chr3_clean \
    -1 Inputs/Reads/SRR5956436/*R1.*fastq.gz \
    -2 Inputs/Reads/SRR5956436/*R2.*fastq.gz \
    -S Hisat2/SRR5956436/Hisat2.sam
)

Notes

  • This aligns one paired-end RNA-seq sample to the genome.
  • -x points to the HISAT2 index basename.
  • --max-intronlen 50000 sets the maximum intron size HISAT2 will consider.
  • --rna-strandness RF specifies the stranded library type.
  • --dta means the alignments are tailored for downstream transcript assembly.
  • -S writes SAM output to a file.

3.2 Sample SRR5956436 - redirect alignment output, sort, convert to bam format and index the bam

(
mkdir -p Hisat2/SRR5956436
hisat2 \
    --max-intronlen 50000 \
    --rna-strandness RF \
    --dta -p 2 \
    -x Hisat2/db/Athaliana_447_TAIR10_Chr3_clean \
    -1 Inputs/Reads/SRR5956436/*R1.*fastq.gz \
    -2 Inputs/Reads/SRR5956436/*R2.*fastq.gz |
samtools sort -@ 2 -o Hisat2/SRR5956436/Hisat2.bam && \
samtools index Hisat2/SRR5956436/Hisat2.bam
)

Notes

  • This pipes HISAT2 output directly into samtools sort.
  • The result is a coordinate-sorted BAM file.
  • samtools index creates a .bai index so the BAM can be accessed by genomic position.

3.3 Align all samples

(
for folder in Inputs/Reads/SR*; do
    echo -ne "\n\n## Running Hisat2 on sample ${folder} ... " && \
    folder=$(basename $folder) && mkdir -p Hisat2/${folder} && \
    hisat2 \
        --max-intronlen 50000  --rna-strandness RF \
        --dta -p 2 \
        -x Hisat2/db/Athaliana_447_TAIR10_Chr3_clean \
        -1 Inputs/Reads/${folder}/*R1.*fastq.gz \
        -2 Inputs/Reads/${folder}/*R2.*fastq.gz |
    samtools sort -@ 2 -o Hisat2/${folder}/Hisat2.bam
done
)

Notes

  • This applies the same alignment procedure to all samples.
  • Each sample gets its own output folder in Hisat2/.

4 Transcript assemblies

mkdir -p {Stringtie,Scallop}/SRR5956436

4.1 Stringtie sample SRR5956436

(
stringtie Hisat2/SRR5956436/Hisat2.bam \
    -o Stringtie/SRR5956436/SRR5956436.gtf \
    > Stringtie/SRR5956436/stringtie.log 2>&1
)

Notes

  • StringTie assembles transcripts from the aligned reads in the BAM file.
  • The output is a GTF file containing transcript models.

4.2 Scallop sample SRR5956436

(
scallop -i Hisat2/SRR5956436/Hisat2.bam \
    -o Scallop/SRR5956436/SRR5956436.gtf \
    > Scallop/SRR5956436/scallop.log 2>&1
)

Notes

  • Scallop is another transcript assembler.
  • It uses the same BAM input, but may produce different transcript models from StringTie.

Try yourself

  • Does Scallop predict more transcripts than StringTie?
  • Are the extra predictions mostly multi-exon or single-exon?

4.3 Assemble all samples

Stringtie

(
for folder in Inputs/Reads/SRR*; do
    echo -ne "\n\n## Running stringtie on sample ${folder} ... " && \
    folder=$(basename ${folder}) && mkdir -p Assemblies Stringtie/${folder} && \
    stringtie \
        Hisat2/${folder}/Hisat2.bam \
        -l ${folder}_STRG \
        -o Stringtie/${folder}/${folder}.gtf \
        > Stringtie/${folder}/stringtie.log 2>&1 && \
    cp Stringtie/${folder}/${folder}.gtf Assemblies/stringtie-${folder}.gtf && \
    echo done
done
)

Notes

  • -l ${folder}_STRG adds a sample-specific label prefix to transcript IDs (without this assemblies would share IDs).
  • GTF files are copied into Assemblies/ for easier downstream comparison.

Scallop

(
for folder in Inputs/Reads/SRR*; do
    echo -ne "\n\n## Running scallop on sample ${folder} ... " && \
    folder=$(basename ${folder}) && mkdir -p Assemblies Scallop/${folder} && \
    scallop \
        -i Hisat2/${folder}/Hisat2.bam \
        -o Scallop/${folder}/${folder}.gtf \
        > Scallop/${folder}/scallop.log 2>&1 && \
    cp Scallop/${folder}/${folder}.gtf Assemblies/scallop-${folder}.gtf && \
    echo done
done
)

Try yourself

  • Compare transcript counts across samples.
  • Which assembler seems more conservative?

5 Mikado Compare

5.1 Compare stringtie vs scallop

(
mkdir -p Mikado_Compare
mikado compare \
    -r Assemblies/stringtie-SRR5956436.gtf \
    -p Assemblies/scallop-SRR5956436.gtf \
    -o Mikado_Compare/mikado_compare_stringtie-SRR5956436_v_scallop-SRR5956436
)

Notes

  • This compares transcript models from the two assemblers for the same sample.
  • Useful for showing how different assemblers may make different choices.

5.2 Compare assemblies vs reference annotation

Stringtie assemblies vs reference annotation

(
mkdir -p Mikado_Compare
for file in Assemblies/str*gtf; do
    echo -e "\n\n## $file" && \
    outfile=$(basename ${file} .gtf) && \
    mikado compare \
        -r Inputs/Ref_Annotation/Athaliana_447_Araport11.gene_exons.regionA.gtf \
        -p $file \
        -o Mikado_Compare/mikado_compare_ref_v_${outfile}
done
)

Scallop assemblies vs reference annotation

(
mkdir -p Mikado_Compare
for file in Assemblies/sca*gtf; do
    echo -e "\n\n## $file" && \
    outfile=$(basename ${file} .gtf) && \
    mikado compare \
        -r Inputs/Ref_Annotation/Athaliana_447_Araport11.gene_exons.regionA.gtf \
        -p $file \
        -o Mikado_Compare/mikado_compare_ref_v_${outfile}
done
)

Notes

  • Each assembly is compared against the reference annotation.
  • This helps assess how well the assemblies recover known transcript structures.

Try yourself

  • Which assembler matches the reference better overall?
  • Do some samples recover more of the annotation than others?

6 Combining assemblies

6.1 Stringtie --merge 

(
assemblies=(Assemblies/str*gtf)

echo -ne "\n\n## Running stringtie --merge on ${assemblies[@]}\n"
mkdir -p Assemblies/Stringtie_Merge
stringtie \
    --merge "${assemblies[@]}" \
    -o Assemblies/Stringtie_Merge/merged_stringtie.gtf
)

Notes

  • stringtie --merge combines transcript models from multiple samples into one merged set.

Try yourself

  • Does the merged set contain more transcripts than any single sample?
  • Check out the merge usage mode in the stringtie help (stringtie -h)

6.2 Concatenate stringtie files

cat Assemblies/str*gtf | gffread -T -o Assemblies/Stringtie_Merge/all_stringtie.gtf

Notes

  • This concatenates all StringTie GTFs and converts them into a single GTF.
  • It does not intelligently merge transcript structures.

6.3 Cluster stringtie files

cat Assemblies/str*gtf | gffread -T -M -K -Q -o Assemblies/Stringtie_Merge/all_stringtie_clustered.gtf

Notes

  • This is a more aggressive attempt to cluster or collapse related transcript models.
  • It is useful for comparing different ways of combining assemblies.

6.4 Compare vs reference annotation

(
mkdir -p Assemblies/Stringtie_Merge/Mikado_Compare
for file in Assemblies/Stringtie_Merge/*gtf; do
    echo -e "\n\n## $file"
    outfile=$(basename ${file} .gtf)
    mikado compare \
        -r Inputs/Ref_Annotation/Athaliana_447_Araport11.gene_exons.regionA.gtf \
        -p $file \
        -o Assemblies/Stringtie_Merge/Mikado_Compare/mikado_compare_ref_v_${outfile}
done
)

Notes

  • This compares the merged, concatenated, and clustered transcript sets to the reference annotation.

Try yourself

  • Does merging improve reference recovery?
  • Does merging with the default parameters lose some sensitivity compared to the concatenated transcripts ?
  • Try rerunning the merging with different merging options e.g. lowering -F or -T and rerunning the comparison to the reference, do we increase sn, pr or F1?

7 Portcullis

7.1 Realign reads omitting the --dta option (it is generally sensible for annotation purposes to run with the --dta option we are only removing this option for demonstration purposes as this is a small toy dataset)

(
for folder in Inputs/Reads/SR*; do
    echo -e "\n\n## Running Hisat2 on sample ${folder} ... " && \
    folder=$(basename $folder) && \
    mkdir -p Portcullis/Hisat2/${folder} && \
    hisat2 \
        --max-intronlen 50000 \
        --rna-strandness RF -p 2 \
        -x Hisat2/db/Athaliana_447_TAIR10_Chr3_clean \
        -1 Inputs/Reads/${folder}/*R1.*fastq.gz \
        -2 Inputs/Reads/${folder}/*R2.*fastq.gz |
    samtools sort -@ 2 -o Portcullis/Hisat2/${folder}/Hisat2.bam
done
)

Notes

  • --dta is omitted here only for demonstration purposes.
  • In a real annotation workflow, using --dta is usually sensible for transcript assembly.
  • Here the main aim is to generate alignments to explore Portcullis junction filtering.

7.2 Portcullis full

(
for folder in Inputs/Reads/SRR*; do
    echo -ne "\n\n## Running portcullis on sample ${folder} ... " && \
    cd Portcullis/Hisat2/$(basename ${folder}) && \
    rm -rf portcullis_out && \
    portcullis full \
        --keep_temp --exon_gff --save_bad \
        ../../../Inputs/Reference/Athaliana_447_TAIR10_Chr3_clean.fa \
        Hisat2.bam > portcullis.log 2>&1 && \
        cd ../../.. && \
        echo done
done
)

Notes

  • Portcullis evaluates splice junctions and separates likely good junctions from likely false ones.
  • --save_bad keeps failed junctions for later comparison.

7.3 Compare results unfiltered vs pass junctions (Sample SRR7947124)

(
for i in Portcullis/Hisat2/SRR7947124; do
    sample=$(basename $i)
    echo -e "\n\n## Processing .... " $sample
    rm -rf $i/Mikado_compare
    mkdir -p $i/Mikado_compare
    cp $i/portcullis_out/2-junc/portcullis_all.junctions.exon.gff3 \
        $i/Mikado_compare/${sample}_portcullis_all.junctions.exon.gff3
    cp $i/portcullis_out/3-filt/portcullis_filtered.pass.junctions.exon.gff3 \
        $i/Mikado_compare/${sample}_portcullis_filtered.pass.junctions.exon.gff3
    for f in $i/Mikado_compare/*gff3; do
        cat $f | sed -e 's/match\t/transcript\t/g' -e 's/match_part\t/exon\t/g' > $i/Mikado_compare/$(basename $f .gff3).fix.gff3
        mikado compare \
            -r Inputs/Ref_Annotation/Athaliana_447_Araport11.gene_exons.regionA.gtf \
            -p $i/Mikado_compare/$(basename $f .gff3).fix.gff3 \
            -o $i/Mikado_compare/mikado_compare_$(basename $f)
    done
done
)
paste Portcullis/Hisat2/SRR7947124/Mikado_compare/mikado_compare_SRR7947124_portcullis_all.junctions.exon.gff3.stats \
Portcullis/Hisat2/SRR7947124/Mikado_compare/mikado_compare_SRR7947124_portcullis_filtered.pass.junctions.exon.gff3.stats

Notes

  • This compares all detected junctions versus only pass junctions for one sample.
  • The sed step changes feature names so Mikado Compare can treat them more like transcript/exon features.

Try yourself

  • Does using only pass junctions improve agreement with the reference?
  • Is the improvement mainly in precision, sensitivity, or both?

7.4 Compare results unfiltered vs pass junctions (All samples)

(
for i in Portcullis/Hisat2/*; do
    sample=$(basename $i)
    echo -e "\n\n## Processing .... " $sample
    rm -rf $i/Mikado_compare
    mkdir -p $i/Mikado_compare
    cp $i/portcullis_out/2-junc/portcullis_all.junctions.exon.gff3 \
        $i/Mikado_compare/${sample}_portcullis_all.junctions.exon.gff3
    cp $i/portcullis_out/3-filt/portcullis_filtered.pass.junctions.exon.gff3 \
        $i/Mikado_compare/${sample}_portcullis_filtered.pass.junctions.exon.gff3
    for f in $i/Mikado_compare/*gff3; do
        cat $f | sed -e 's/match\t/transcript\t/g' -e 's/match_part\t/exon\t/g' > $i/Mikado_compare/$(basename $f .gff3).fix.gff3
        mikado compare \
            -r Inputs/Ref_Annotation/Athaliana_447_Araport11.gene_exons.regionA.gtf \
            -p $i/Mikado_compare/$(basename $f .gff3).fix.gff3 \
            -o $i/Mikado_compare/mikado_compare_$(basename $f)
    done
done
)

7.5 Combine all pass or fail junction files

(
junctools set --operator max \
    -o Portcullis/Hisat2/portcullis_pass_junctions.bed \
    union Portcullis/Hisat2/*/portcullis_out/3-filt/portcullis_filtered.pass.junctions.bed && \
junctools set --operator max \
    -o Portcullis/Hisat2/portcullis_fail_junctions.bed \
    union Portcullis/Hisat2/*/portcullis_out/3-filt/portcullis_filtered.fail.junctions.bed
)

7.6 Convert to GFF

(
junctools convert -if bed -of egff \
    -o Portcullis/Hisat2/portcullis_pass_junctions.gff \
    Portcullis/Hisat2/portcullis_pass_junctions.bed && \
junctools convert -if bed -of egff \
    -o Portcullis/Hisat2/portcullis_fail_junctions.gff \
    Portcullis/Hisat2/portcullis_fail_junctions.bed
)

To clean up the directory to the original starting files

rm -rf Stringtie Scallop Mikado_Compare Assemblies SRR5956436* mikado_compare_* Hisat2 Portcullis

Other things to try, have a go at filtering the Portcullis/Hisat2/SRR7947124/Hisat2.bam file using the pass junctions (.tab file) with portcullis bamfilt. Try assembling the filtered bam with stringtie and compare it to the results of assembling the unfiltered bam. Does this show any improvement?

Clone this wiki locally