Transcript assembly commands

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 -e "\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 && mkdir -p 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 && mkdir -p 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 vs reference annotation

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

mkdir -p Mikado_Compare && for file in Assemblies/sca*gtf; do outfile=$(basename ${file} .gtf) && echo $file && 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="" && for i in Assemblies/str*gtf; do assemblies="$assemblies $i"; done; 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 outfile=$(basename ${file} .gtf) && echo $file && 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 "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 's/match\t/transcript\t/g' | sed '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 "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 's/match\t/transcript\t/g' | sed '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_junctools.gff Portcullis/Hisat2/portcullis_fail_junctions.bed

To clean up 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?