Mapping RNAseq reads

Mapping RNAseq reads should be performed with a dedicated, splice-aware aligner, such as STAR. This is an ultrafast RNAseq aligner that outperforms many other mappers. It uses a sequential maximum mappable seed search in uncompressed suffix arrays followed by seed clustering and a stitching procedure. In addition to unbiased de novo detection of canonical splice junctions (i.e., flanked by GU-AG that make up to 99% of splice sites in most organisms) STAR can discover non-canonical splice sites and chimeric (fusion) transcripts.

Normally, it is recommended to use a reference that has already annotated splice junctions or a 2-phase mapping approach with STAR. In the latter approach, in the 1^st pass STAR mapping of every sample, de novo information on the intron-exon (reference genome) and exon-exon (reference transcriptome) junctions is being collected and used to generate a new genome index for the 2^nd pass mapping. We are going to use the first approach where splice junctions are already available in an annotation file, in .gff format.

First we need to prepare the reference genome for alignment by generating genome index files from the .fasta file (and .gtf or .gff annotation files if available). We will use a subset of 548 contigs (a total of 75,170,146 bp) from the reference genome of Heliosperma pusillum v.1.0 (Szukala et al. 2023). Each contig may contain several genes.

For small genomes/transcriptomes --genomeSAindexNbases (a pre-indexing setting) must be scaled to min(14, log₂(Genome length)/2-1). The size of the bins for genome storage (--genomeChrBinNBits) needs to be also scaled to min(18, log₂(GenomeLength/NoReferences)). The setting –sjdbGTFfile gives the path to a GTF/GFF file with annotations, including splice junctions. Finally, --sjdbOverhang should be set to maximum read length -1.

In the next step we align the reads of each sample one by one to the reference genome. Here .fasta or .fastq files can be used, but they should be quality trimmed/filtered. To exemplify the step, we will map single end reads from just one individual - the reads have been already trimmed with trimmomatic.

STAR will write several output files, including alignments in .sam or .bam format, mapping summary statistics (a file ending in Log.final.out), splice junctions (a file ending with SJ.out.tab), unmapped reads, etc. There are many options that can be used, here we specify the most important ones.

The options explained:
--runThreadN: number of threads
--outFilterMismatchNoverLmax: report a mapping only if the ratio of mismatches to mapped length is less or equal to this value. Several other filtering values are available.