Transcript assembly commands

1 Set the genome environment variable

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 -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 -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 -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 && 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 . && cp Hisat2/SRR5956436/portcullis_out/3-filt/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 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 portcullis_filtered.pass.junctions.exon.gff3.fix.gff3 -o mikado_compare_portcullis_pass