Aligning Reads to a Reference

Trimmed reads are just sequences until you place them on the genome. RNA-seq alignment has to be splice-aware, because a read can span an exon-exon junction and map to two exons separated by an intron. HISAT2 and STAR are the two standard spliced aligners; this lesson runs HISAT2 end to end to produce a sorted, indexed BAM that is ready for quantification.

Build the genome index

# Build a HISAT2 index from the reference genome FASTA
hisat2-build -p 4 GRCh38.primary_assembly.genome.fa grch38

This writes eight files named grch38.1.ht2 through grch38.8.ht2. The value grch38 is the index basename you pass to the aligner, not a single file on disk. STAR is a faster alternative but needs far more memory (roughly 30 GB for the human genome); it builds an index with STAR --runMode genomeGenerate and can emit a sorted BAM directly using --outSAMtype BAM SortedByCoordinate.

Align to the reference

hisat2 -p 4 -x grch38 \
  -1 sample_R1.trimmed.fastq.gz \
  -2 sample_R2.trimmed.fastq.gz \
  --summary-file sample.hisat2.log \
  | samtools sort -@ 4 -o sample.sorted.bam -

20000000 reads; of these:
  20000000 (100.00%) were paired; of these:
    1650000 (8.25%) aligned concordantly 0 times
    16200000 (81.00%) aligned concordantly exactly 1 time
    2150000 (10.75%) aligned concordantly >1 times
    ----
    1650000 pairs aligned concordantly 0 times; of these:
      210000 (12.73%) aligned discordantly 1 time
    ----
    1440000 pairs aligned 0 times concordantly or discordantly; of these:
      2880000 mates make up the pairs; of these:
        1900000 (65.97%) aligned 0 times
        720000 (25.00%) aligned exactly 1 time
        260000 (9.03%) aligned >1 times
95.25% overall alignment rate

HISAT2 streams SAM to stdout, so piping into samtools sort avoids ever writing an uncompressed SAM file, and the summary above also lands in sample.hisat2.log. The last line, the overall alignment rate, is the most useful number: a clean mammalian library is usually 90 percent or higher, while under about 70 percent points to leftover adapters, the wrong reference or species, or contamination from an organism that is not in your reference. Genomic-DNA and most rRNA reads still map to the genome, so they pad intronic and intergenic counts rather than dragging this rate down; you catch those with downstream QC, not this number. In the flagstat report that follows, primary mapped should equal this rate, properly paired counts mates that mapped in the expected orientation and distance, and singletons are reads whose mate failed to map.

Index and inspect

samtools index sample.sorted.bam
samtools flagstat sample.sorted.bam

40260000 + 0 in total (QC-passed reads + QC-failed reads)
40000000 + 0 primary
260000 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
38360000 + 0 mapped (95.28% : N/A)
38100000 + 0 primary mapped (95.25% : N/A)
40000000 + 0 paired in sequencing
20000000 + 0 read1
20000000 + 0 read2
36700000 + 0 properly paired (91.75% : N/A)
37120000 + 0 with itself and mate mapped
980000 + 0 singletons (2.45% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)