1.1 Assembling transcripts based on RNA-seq data

RNA-Seq data is in general very useful in annotation projects as the data usually comes from the actual organism you are studying and thus avoids the danger of introducing errors caused by differences in gene structure between your study organism and other species.

Important points to remember before starting working with RNA-Seq:

• Check if sequence data is paired-end or single-end. Last generation of sequenced short reads (since 2013) are almost all paired. Anyway, it is important to check that information, which will be useful for the tools used in the next steps.

• Check if the sequence data are stranded. This information will be useful for the tools used in the next steps. (It is recommended to use stranded data to avoid transcript fusion during the transcript assembly process. That gives more reliable results.)

• Left / L / forward / 1 have identical meaning. It is the same for Right / R /Reverse / 2

module load bioinfo-tools
module load FastQC/0.11.5

mkdir fastqc
fastqc ./raw/ERR305399_1.fastq.gz -o fastqc/

scp user@rackham.uppmax.uu.se:[path-to-report.html] .

export PERL5LIB=$PERL5LIB:./GAAS/annotation/
./GAAS/annotation/Tools/bin/fastq_guessMyFormat.pl -i ./raw/ERR305399_1.fastq.gz

2.1 Trimmomatic

Trimmomatic performs a variety of useful trimming tasks for illumina paired-end and single ended data.The selection of trimming steps and their associated parameters are supplied on the command line.

mkdir trimmomatic
module load bioinfo-tools
module load trimmomatic/0.36

The following command line will perform the following: * Remove adapters (ILLUMINACLIP:TruSeq3-PE.fa:2:30:10) * Remove leading low quality or N bases (below quality 3) (LEADING:3) * Remove trailing low quality or N bases (below quality 3) (TRAILING:3) * Scan the read with a 4-base wide sliding window, cutting when the average quality per base drops below 15 (SLIDINGWINDOW:4:15) * Drop reads below the 36 bases long (MINLEN:36)

java -jar /sw/apps/bioinfo/trimmomatic/0.36/milou/trimmomatic-0.36.jar PE -threads 5 -phred33 ./raw/ERR305399_1.fastq.gz ./raw/ERR305399_2.fastq.gz trimmomatic/ERR305399.left_paired.fastq.gz trimmomatic/ERR305399.left_unpaired.fastq.gz trimmomatic/ERR305399.right_paired.fastq.gz trimmomatic/ERR305399.right_unpaired.fastq.gz ILLUMINACLIP:/sw/apps/bioinfo/trimmomatic/0.36/milou/adapters/TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

2.2 HISAT2

Once the reads have been trimmed, we use HISAT2 to align the RNA-seq reads to a genome in order to identify exon-exon splice junctions. HISAT2 is a fast and sensitive alignment program for mapping next-generation sequencing reads (whole-genome, transcriptome, and exome sequencing data) against a reference genome.

First you need to build an index of your chromosome 4

mkdir index

module load HISAT2
module load samtools/1.8
hisat2-build /proj/uppstore2017171/courses/RNAseqWorkshop/downloads/assembly_annotation/chromosome/chr4.fa ./index/chr4_index

Then you can run HISAT2 : * --phred33 Input qualities are ASCII chars equal to the Phred quality plus 33. This is also called the “Phred+33” encoding, which is used by the very latest Illumina pipelines (you checked it with the script fastq_guessMyFormat.pl)
* --rna-strandness <string> For single-end reads, use F or R. F means a read corresponds to a transcript. R means a read corresponds to the reverse complemented counterpart of a transcript. For paired-end reads, use either FR or RF. (RF means fr-firststrand see here for more explanation).
* --novel-splicesite-outfile <path> In this mode, HISAT2 reports a list of splice sites in the file : chromosome name genomic position of the flanking base on the left side of an intron genomic position of the flanking base on the right strand (+, -, and .) ‘.’ indicates an unknown strand for non-canonical splice sites.

hisat2 --phred33 --rna-strandness RF --novel-splicesite-outfile hisat2/splicesite.txt -S hisat2/accepted_hits.sam -p 5 -x index/chr4_index -1 trimmomatic/ERR305399.left_paired.fastq.gz -2 trimmomatic/ERR305399.right_paired.fastq.gz

Finally you need to change the SAM into BAM file and to sort it in order for StringTie to use the bam file to assemble the read into transcripts.

samtools view -bS -o hisat2/accepted_hits.bam hisat2/accepted_hits.sam
samtools sort -o hisat2/accepted_hits.sorted.bam hisat2/accepted_hits.bam

2.3 StringTie

StringTie is a fast and highly efficient assembler of RNA-Seq alignments into potential transcripts. It uses a novel network flow algorithm as well as an optional de novo assembly step to assemble and quantitate full-length transcripts representing multiple splice variants for each gene locus. You can add as input an annotation from GTF/GFF3 file to calculate TPM and FPKM values.

module load bioinfo-tools
module load StringTie

stringtie hisat2/accepted_hits.sorted.bam -o stringtie/transcripts.gtf

When done you can find your results in the directory outdir. The file transcripts.gtf includes your assembled transcripts.

You could now also visualise all this information using a genome browser, such as IGV. IGV requires a genome fasta file and any number of annotation files in GTF or GFF3 format (note that GFF3 formatted file tend to look a bit weird in IGV sometimes).

Transfer the GTF files to your computer using SCP:

scp user@rackham.uppmax.uu.se:[path-to-transcripts.gtf] .

Looking at your results, are you happy with the default values of StringTie (which we used in this exercise) or is there something you would like to change?

3.1 Assembly quality

There are different ways of assessing the quality of your assembly, you will find some of them here.

We will run BUSCO to check the the quality of the assembly. BUSCO provides measures for quantitative assessment of genome assembly, gene set, and transcriptome completeness (what we are going to do here). Genes that make up the BUSCO sets for each major lineage are selected from orthologous groups with genes present as single-copy orthologs in at least 90% of the species in the chosen branch of tree of life.

For the Trinity results :

module load BUSCO/3.0.2b

/sw/apps/bioinfo/BUSCO/3.0.2b/rackham/bin/run_BUSCO.py -i /proj/uppstore2017171/courses/RNAseqWorkshop/downloads/assembly_annotation/RNAseq/trinity/Trinity.fasta -o busco_trinity -l $BUSCO_LINEAGE_SETS/arthropoda_odb9 -m tran -c 5

Busco will take 30 to run so you can check the results in /proj/uppstore2017171/courses/RNAseqWorkshop/downloads/assembly_annotation/RNAseq/busco_trinity

For the guided assembly results, you need first to extract the transcript sequences from the gtf transcript file :

module load BioPerl

export PERL5LIB=$PERL5LIB:./GAAS/annotation/
./GAAS/annotation/Tools/bin/gff3_sp_extract_sequences.pl --cdna -g transcripts.gtf -f /proj/uppstore2017171/courses/RNAseqWorkshop/downloads/assembly_annotation/chromosome/chr4.fa -o transcripts_stringtie.fa

Then you can run BUSCO again :

/sw/apps/bioinfo/BUSCO/3.0.2b/rackham/bin/run_BUSCO.py -i stringtie/transcripts_stringtie.fa -o busco_stringtie -l $BUSCO_LINEAGE_SETS/arthropoda_odb9 -m tran -c 5

Compare the two BUSCO results, what do you think happened for StringTie?

3.2 What’s next?

Now you are ready either to annotate your RNA-Seq or you can use then to do the genome annotation.

For the de-novo assembly you can use the Trinity.fasta file obtained. For the genome-guided assembly you can either use the StringTie results transcripts.gtf, but you will often need to reformat it into a gff file.

End of document