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

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 Portcullis full

for folder in Inputs/Reads/SRR*; do echo -ne "\n\n## Running portcullis on sample ${folder} ..." && cd 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.2 Compare results unfiltered vs pass junctions (Sample SRR5956436)

rm *fix.gff3; cp Hisat2/SRR5956436/portcullis_out/2-junc/portcullis_all.junctions.exon.gff3 SRR5956436_portcullis_all.junctions.exon.gff3 && cp Hisat2/SRR5956436/portcullis_out/3-filt/portcullis_filtered.pass.junctions.exon.gff3 SRR5956436_portcullis_filtered.pass.junctions.exon.gff3 && for i in *gff3; do cat $i | sed 's/match\t/transcript\t/g' | sed 's/match_part\t/exon\t/g' > $i.fix.gff3; done && mikado compare -r Inputs/Ref_Annotation/Athaliana_447_Araport11.gene_exons.regionA.gtf -p SRR5956436_portcullis_all.junctions.exon.gff3.fix.gff3 -o mikado_compare_portcullis_all && mikado compare -r Inputs/Ref_Annotation/Athaliana_447_Araport11.gene_exons.regionA.gtf -p SRR5956436_portcullis_filtered.pass.junctions.exon.gff3.fix.gff3 -o mikado_compare_portcullis_pass

7.3 Combine all pass or fail junction files

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

7.4 Convert to GFF

junctools convert -if bed -of egff -o Hisat2/portcullis_pass_junctions.gff Hisat2/portcullis_pass_junctions.bed && junctools convert -if bed -of egff -o Hisat2/portcullis_fail_junctools.gff 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