Skip to content

Transcript assembly commands

Jump to bottom
Gemy George Kaithakottil edited this page Mar 27, 2026 · 26 revisions

1 Set the genome environment variable

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

2 Build Hisat2 index

(
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
)

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
)

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
)

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
)

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
)

4.2 Scallop sample SRR5956436

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

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
)

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
)

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