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
)

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
)

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
)

6.2 Concatenate stringtie files

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

6.3 Cluster stringtie files

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

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
)

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
)

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
)

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/*stats

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