Transcript assembly commands

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
)

The stringtie and scallop assemblies can be viewed in the Apollo browse

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
)

To quickly compare the transcript level sn, pr and F1 across samples use the following script

Inputs/Scripts/mikado_compare_stats_pivot.py 
    --rows t80 \
    --label_regex '(?<=mikado_compare_ref_v_).*(?=\.stats$)' \ 
    -o mikado_compare_ref_v_stringtie_scallop.tsv \
    Mikado_Compare/mikado_compare_ref_*stats

And view with

columntab <  mikado_compare_ref_v_stringtie_scallop.tsv |less -S

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
)

Notes

This extends the same comparison across all samples.
Useful for checking whether the effect of Portcullis filtering is consistent.

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
)

Notes

This creates pooled pass and fail junction sets across all samples.
Useful for building a consensus splice-junction resource.

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
)

Notes

This converts the pooled BED junction files into GFF-like output.
These files can be useful for loading into a genome browser or annotation workflow.

To clean up the directory to the original starting files

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

Notes

This removes generated files so the exercise can be rerun from a clean state.

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?

Extra suggestions for self-exploration

1. Check alignment quality in more detail

samtools flagstat Hisat2/SRR5956436/Hisat2.bam
samtools stats Hisat2/SRR5956436/Hisat2.bam | less

Ask yourself:

What fraction are properly paired?
How many secondary or supplementary alignments are present?

2. Compare transcript counts quickly

grep -c $'\ttranscript\t' Stringtie/SRR5956436/SRR5956436.gtf
grep -c $'\ttranscript\t' Scallop/SRR5956436/SRR5956436.gtf

3. Compare exon counts

grep -c $'\texon\t' Stringtie/SRR5956436/SRR5956436.gtf
grep -c $'\texon\t' Scallop/SRR5956436/SRR5956436.gtf

4. Look at a locus in a genome browser

Suggested tracks:

reference annotation GTF
HISAT2 BAM
StringTie GTF
Scallop GTF
Portcullis pass junction GFF
Portcullis fail junction GFF

Ask yourself:

Which assembler better matches the aligned read evidence?
Are fail junctions visibly weak or unsupported?

5. Try changing one parameter

For example:

change --max-intronlen
remove --dta
try a different sample
compare filtered vs unfiltered BAM assembly

6. Ask biological questions, not only technical ones

Are novel transcripts plausible isoforms, or likely artefacts?
Does combining samples help recover more complete gene models?
Does splice-junction filtering improve precision more than sensitivity?

Uh oh!

Transcript assembly commands

1 Set the genome environment variable

2 Build Hisat2 index

3 Perform HISAT2 alignment

3.1 Sample SRR5956436

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

3.3 Align all samples

4 Transcript assemblies

4.1 Stringtie sample SRR5956436

4.2 Scallop sample SRR5956436

4.3 Assemble all samples

Stringtie

Scallop

5 Mikado Compare

5.1 Compare stringtie vs scallop

5.2 Compare assemblies vs reference annotation

6 Combining assemblies

6.1 Stringtie --merge

6.2 Concatenate stringtie files

6.3 Cluster stringtie files

6.4 Compare vs reference annotation

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)

7.2 Portcullis full

7.3 Compare results unfiltered vs pass junctions (Sample SRR7947124)

7.4 Compare results unfiltered vs pass junctions (All samples)

7.5 Combine all pass or fail junction files

7.6 Convert to GFF

Other things to try

Extra suggestions for self-exploration

1. Check alignment quality in more detail

2. Compare transcript counts quickly

3. Compare exon counts

4. Look at a locus in a genome browser

5. Try changing one parameter

6. Ask biological questions, not only technical ones

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

1 Set the `genome` environment variable

6.1 Stringtie `--merge`