4A. Genome-based quantification (STAR + featureCounts)

Estimated time: 70 minutes

Overview

Questions

• How do we map RNA-seq reads to a reference genome?
• How do we determine library strandness before alignment?
• How do we quantify reads per gene using featureCounts?
• How do we submit mapping jobs on RCAC clusters with SLURM?

Objectives

• Check RNA-seq library strandness using Salmon.
• Index a reference genome for STAR.
• Map reads to the genome using STAR.
• Quantify gene level counts with featureCounts.
• Interpret basic alignment and counting statistics.

Introduction

---

In this episode we move from raw FASTQ files to genome based gene counts.
This workflow has three major steps:

1. Determine library strandness: STAR does not infer strandness automatically. featureCounts requires the correct setting.
2. Map reads to the genome using STAR: STAR is a splice aware aligner designed for RNA-seq data.
3. Count mapped reads per gene using featureCounts

This represents one of the two standard quantification strategies in RNA-seq.
The alternative method (B. transcript based quantification with Salmon or Kallisto) will be covered in Episode 4B.

Step 1: Checking RNA-seq library strandness

---

Callout

Why strandness matters

RNA-seq libraries can be:

• unstranded
  
• stranded forward
  
• stranded reverse

Stranded protocols are common in modern kits, but old GEO and SRA datasets are often unstranded.
featureCounts must be given the correct strand setting. If a stranded library is counted with the wrong direction, most reads land on the strand opposite their gene and are left unassigned (or counted to an antisense gene), so the assignment rate collapses. Counting a stranded library as unstranded loses less, but reads where genes on opposite strands overlap can no longer be assigned.

STAR does not infer strandness automatically, so we must determine it before alignment.

Animation (80 seconds): stranded vs unstranded RNA-seq. In a dUTP library, read 2 matches the RNA strand; Salmon reports this as ISR and featureCounts calls it -s 2. Our dataset turns out to be unstranded, so we count with -s 0, and reads where genes on opposite strands overlap end up among the Unassigned_Ambiguity reads in the featureCounts summary. Press Play to start, or open it in a new tab.

Detecting strandness with Salmon

Salmon can infer strandness directly from raw FASTQ files without any alignment.
We run Salmon with -l A, which tells it to explore all library types.

BASH

sinteractive -A rcac-rnaseq -q standby -p cpu -N 1 -n 4 --time=2:00:00

cd $SCRATCH/rnaseq-workshop
mkdir -p results/strand_check

module load biocontainers
module load salmon

salmon index --transcripts data/gencode.vM38.transcripts-clean.fa \
             --index data/salmon_index \
             --threads 4

salmon quant --index data/salmon_index \
             --libType A \
             --mates1 data/WT_Bcell_mock_rep1_R1.fastq.gz \
             --mates2 data/WT_Bcell_mock_rep1_R2.fastq.gz \
             --output results/strand_check \
             --threads 4

Building the Salmon index and running the check takes about 8 minutes on 4 cores. The important output is in results/strand_check/lib_format_counts.json:

What does this file look like?

JSON

{
    "read_files": "[ data/WT_Bcell_mock_rep1_R1.fastq.gz, data/WT_Bcell_mock_rep1_R2.fastq.gz]",
    "expected_format": "IU",
    "compatible_fragment_ratio": 1.0,
    "num_compatible_fragments": 11994681,
    "num_assigned_fragments": 11994681,
    "num_frags_with_concordant_consistent_mappings": 10472924,
    "num_frags_with_inconsistent_or_orphan_mappings": 2280227,
    "strand_mapping_bias": 0.5082214861866657,
    "MSF": 0,
    "OSF": 0,
    "ISF": 5322565,
    "MSR": 0,
    "OSR": 0,
    "ISR": 5150359,
    "SF": 1149290,
    "SR": 1130937,
    "MU": 0,
    "OU": 0,
    "IU": 0,
    "U": 0
}

Interpreting lib_format_counts.json

This JSON file reports Salmon’s automatic library type detection.

• "expected_format": "IU": This is Salmon’s conclusion. “I” means Inward (correct for paired-end reads) and “U” means Unstranded.
• "strand_mapping_bias" (0.508 here): This is the key evidence. A value near 0.5 (50%) indicates that reads mapped equally to both the sense and antisense strands, the definitive sign of an unstranded library.
• "ISF" and "ISR": These are the counts of inward-facing fragments whose first read maps to the forward strand of the transcript (ISF) or to the reverse strand (ISR): 5.3 and 5.2 million here. Because these values are almost equal, they confirm the ~50/50 split seen in the strand bias.

Common results:

Code	Meaning
IU	unstranded
ISR	stranded reverse (paired end)
ISF	stranded forward (paired end)
SR	stranded reverse (single end)
SF	stranded forward (single end)

Conclusion: The data is unstranded. We will record this and supply it to STAR/featureCounts.

Step 2: Preparing the reference genome for alignment

---

STAR requires a genome index. Indexing builds a searchable data structure that speeds up alignment.

Required input:

• genome FASTA file
• annotation file (GTF or GFF3)
• output directory for the index

Optionally, annotation improves splice junction detection and yields better alignments.

Callout

Why STAR?

• Splice-aware
  
• Extremely fast
  
• Produces high-quality alignments
  
• Well supported and widely used
  
• Compatible with downstream tools such as featureCounts, StringTie, and RSEM

Optional alternatives include HISAT2 and Bowtie2, but STAR remains the standard for high-depth Illumina RNA-seq.

Challenge

Exercise: Read the STAR manual

Check the “Generating genome indexes” section of the
STAR manual

Which options are required for indexing?

Show me the solution

Key options:

• --runMode genomeGenerate
  
• --genomeDir directory to store the index
  
• --genomeFastaFiles the FASTA file
  
• --sjdbGTFfile annotation (recommended)
  
• --runThreadN number of threads
  
• --sjdbOverhang read length minus one (100 is a safe default)

Indexing the Genome with STAR

Below is a minimal SLURM job script for building the STAR genome index. Create a file named index_genome.sh in $SCRATCH/rnaseq-workshop/scripts with this content:

BASH

#!/bin/bash
#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=48
#SBATCH --account=rcac-rnaseq
#SBATCH --qos=standby
#SBATCH --partition=cpu
#SBATCH --time=2:00:00
#SBATCH --job-name=star_index
#SBATCH --output=cluster-%x.%j.out
#SBATCH --error=cluster-%x.%j.err

module load biocontainers
module load star

data_dir="${SCRATCH}/rnaseq-workshop/data"
genome="${data_dir}/GRCm39.primary_assembly.genome.fa"
gtf="${data_dir}/gencode.vM38.primary_assembly.basic.annotation.gtf"
indexdir="${data_dir}/star_index"

STAR \
    --runMode genomeGenerate \
    --runThreadN ${SLURM_CPUS_ON_NODE} \
    --genomeDir ${indexdir} \
    --genomeFastaFiles ${genome} \
    --sjdbGTFfile ${gtf}

Submit the job from the scripts directory:

BASH

cd $SCRATCH/rnaseq-workshop/scripts
sbatch index_genome.sh

Indexing takes about 15 minutes and needs about 30 GB of memory. On Negishi, memory comes with the cores you request (there is no separate memory request), so the script asks for 48 cores to get enough memory for indexing a mammalian genome. Check progress with squeue -u $USER; the job log is written to cluster-star_index.<jobid>.out in the scripts directory.

Callout

Short on time during the workshop?

If your index job is still waiting in the queue when we reach the mapping step, use the prebuilt index from the completed results copy instead of your own:

BASH

ln -s /depot/workshop/data/rnaseq-workshop_results/data/star_index $SCRATCH/rnaseq-workshop/data/star_index

Run this only if $SCRATCH/rnaseq-workshop/data/star_index does not exist yet (cancel your index job first with scancel <jobid>). The mapping script below works the same with either index.

What will the output look like?

The index directory is now ready for use in mapping and should have contents like this: Location: $SCRATCH/rnaseq-workshop/data Contents of star_index/:

chrLength.txt
chrNameLength.txt
chrName.txt
chrStart.txt
exonGeTrInfo.tab
exonInfo.tab
geneInfo.tab
Genome
genomeParameters.txt
Log.out
SA
SAindex
sjdbInfo.txt
sjdbList.fromGTF.out.tab
sjdbList.out.tab
transcriptInfo.tab

Callout

Do I need the GTF for indexing?

No, but it is strongly recommended.
Including it improves splice junction detection and increases mapping accuracy.

Step 3: Mapping reads with STAR

---

Once the genome index is ready, we can align reads.

For aligning reads, important STAR options are:

• --runThreadN
• --genomeDir
• --readFilesIn
• --readFilesCommand zcat for gzipped FASTQ files
• --outSAMtype BAM SortedByCoordinate
• --outSAMunmapped Within

STAR’s defaults are tuned for mammalian genomes. If working with plants, fungi, or repetitive genomes, additional tuning may be required.

Challenge

Exercise: Write a STAR alignment command

Write a basic STAR alignment command using your genome index and paired FASTQ files.

Show me the solution

BASH

STAR \
    --runThreadN 12 \
    --genomeDir data/star_index \
    --readFilesIn data/SRR1234567_1.fastq.gz data/SRR1234567_2.fastq.gz \
    --readFilesCommand zcat \
    --outFileNamePrefix results/mapping/SRR1234567 \
    --outSAMtype BAM SortedByCoordinate \
    --outSAMunmapped Within

Mapping many FASTQ files with a SLURM array

When having multiple samples, array jobs simplify submission.

First, we will create samples.txt file with just the sample names

BASH

cd $SCRATCH/rnaseq-workshop/data
ls *_R1.fastq.gz | sed 's/_R1.fastq.gz//' > ${SCRATCH}/rnaseq-workshop/scripts/samples.txt

Save the array job below as $SCRATCH/rnaseq-workshop/scripts/map_reads.sh. Each array task reads its sample name from line ${SLURM_ARRAY_TASK_ID} of samples.txt, which the script looks for in the directory you submit from, so always submit it from scripts/.

BASH

#!/bin/bash
#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=20
#SBATCH --account=rcac-rnaseq
#SBATCH --qos=standby
#SBATCH --partition=cpu
#SBATCH --time=1:00:00
#SBATCH --job-name=read_mapping
#SBATCH --array=1-8
#SBATCH --output=cluster-%x.%j.out
#SBATCH --error=cluster-%x.%j.err

module load biocontainers
module load star

FASTQ_DIR="$SCRATCH/rnaseq-workshop/data"
GENOME_INDEX="$SCRATCH/rnaseq-workshop/data/star_index"
OUTDIR="$SCRATCH/rnaseq-workshop/results/mapping"

mkdir -p $OUTDIR

SAMPLE=$(sed -n "${SLURM_ARRAY_TASK_ID}p" samples.txt)

R1=${FASTQ_DIR}/${SAMPLE}_R1.fastq.gz
R2=${FASTQ_DIR}/${SAMPLE}_R2.fastq.gz

STAR \
    --runThreadN ${SLURM_CPUS_ON_NODE} \
    --genomeDir $GENOME_INDEX \
    --readFilesIn $R1 $R2 \
    --readFilesCommand zcat \
    --outFileNamePrefix $OUTDIR/$SAMPLE \
    --outSAMunmapped Within \
    --outSAMtype BAM SortedByCoordinate

Submit from the scripts directory:

BASH

cd $SCRATCH/rnaseq-workshop/scripts
sbatch map_reads.sh

Each sample takes about 2 to 3 minutes on 20 cores and up to 28 GB of memory. The eight tasks run in parallel when the queue allows.

Discussion

Why use an array job?

Think about this:

• We have eight samples, each with two FASTQ files
• Having eight separate slurm jobs means:
  • Copy-pasting the same command eight times
  • Higher chance of typos or mistakes
• if mapping with single slurm file, sequential execution = slower
• Samples can run independently

Why is a job array a good choice?

Arrays allow us to run identical commands across many inputs:

• All samples are processed the same way
  
• Runs in parallel
  
• Simplifies submission
  
• Reduces copy-paste errors
  
• Minimizes runtime cost on the cluster

Step 4: Assessing alignment quality

---

STAR produces several useful files, particularly:

• Log.final.out
  Contains mapping statistics, including % uniquely mapped.

• Aligned.sortedByCoord.out.bam
  Ready for downstream quantification.

To summarize alignment statistics across samples, we use MultiQC:

BASH

cd $SCRATCH/rnaseq-workshop
module load biocontainers
module load multiqc

multiqc results/mapping -o results/qc_alignment

Key metrics to review

• Uniquely mapped reads — typically 60–90% is good
  
• Multi-mapped reads — depends on genome, but >20% is concerning
  
• Unmapped reads — check reasons (mismatches? too short?)
  
• Mismatch rate — high values indicate quality or contamination issues

Challenge

Exercise: Examine your alignment metrics

Using your MultiQC alignment report:

1. Which sample has the highest uniquely mapped percentage?
   
2. Are any samples clear outliers?
   
3. Does mapped/unmapped reads appear consistent across samples?

Show me the solution

Interpretation for this dataset:

1. Uniquely mapped reads range from about 62 to 83 percent; mock_rep4 has the highest (83 percent).
2. IR_rep1 and IR_rep2 are lowest (62 and 64 percent), and IR_rep3, IR_rep4, and mock_rep4 highest (77 to 83 percent): the same split as the read quality groups in Episode 3. The two conditions map about equally well on average (70 percent for IR, 73 percent for mock). These rates are within the 60 to 90 percent range given above, so no sample fails, but the spread is worth recording. When one condition maps systematically worse, check whether the difference also shows up in the PCA in Episode 5A.
3. Most unmapped reads are classified as “too short” (1.4 to 4.6 million pairs per sample): STAR could align only part of the read pair. With 51 bp reads this usually reflects low quality read ends (see the FastQC report in Episode 3), rRNA or other sequence absent from the primary assembly, and reads spanning unannotated junctions.

Step 5: Quantifying gene counts with featureCounts

---

Once the reads have been aligned and sorted, we can quantify how many read pairs (fragments) map to each gene. This gives us gene-level counts that we will later use for differential expression. For alignment-based RNA-seq workflows, featureCounts is the recommended tool because it is fast, robust, and compatible with most gene annotation formats.

featureCounts uses two inputs:

1. A set of aligned BAM files (sorted by coordinate)
2. A GTF annotation file describing gene models

It assigns each read (or read pair) to a gene based on the exon boundaries.

Callout

Important: stranded or unstranded?

featureCounts needs to know whether your dataset is stranded.

• -s 0 unstranded
• -s 1 stranded forward
• -s 2 stranded reverse

From our Salmon strandness check, this dataset behaves as unstranded, so we will use -s 0.

Callout

Counting reads or read pairs?

With paired-end data, each fragment produces two reads. Counting fragments (read pairs) counts each cDNA molecule once; counting reads counts most fragments twice and makes the counts look larger and more precise than they are. -p tells featureCounts that the BAM files are paired-end. In Subread 2.0.1, the version on Negishi, -p alone counts fragments. From Subread 2.0.2 onward, -p alone counts reads, and you must add --countReadPairs to count fragments; check your version with featureCounts -v before reusing this script elsewhere.

Running featureCounts

We will create a new directory for count output:

BASH

cd $SCRATCH/rnaseq-workshop/scripts
mkdir -p $SCRATCH/rnaseq-workshop/results/counts

Below is an example SLURM script to count all BAM files at once. Save it as $SCRATCH/rnaseq-workshop/scripts/count_features.sh.

BASH

#!/bin/bash
#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=16
#SBATCH --account=rcac-rnaseq
#SBATCH --qos=standby
#SBATCH --partition=cpu
#SBATCH --time=0:30:00
#SBATCH --job-name=featurecounts
#SBATCH --output=cluster-%x.%j.out
#SBATCH --error=cluster-%x.%j.err

module load biocontainers
module load subread

GTF="$SCRATCH/rnaseq-workshop/data/gencode.vM38.primary_assembly.basic.annotation.gtf"
BAMS="$SCRATCH/rnaseq-workshop/results/mapping/*.bam"
output_dir="$SCRATCH/rnaseq-workshop/results/counts"

featureCounts \
  -T ${SLURM_CPUS_PER_TASK} \
  -p \
  -s 0 \
  -t exon \
  -g gene_id \
  -a ${GTF} \
  -o ${output_dir}/gene_counts.txt \
  ${BAMS}

Submit the job from the scripts directory (the previous block left you there):

BASH

sbatch count_features.sh

Counting all eight BAM files takes about 7 minutes.

gene_counts.txt will contain:

• one row per gene
• one column per sample

Callout

Generate a clean count matrix

BASH

cd $SCRATCH/rnaseq-workshop/results/counts
cut -f 1,7- gene_counts.txt |\
  sed 's+/scratch/negishi/'$USER'/rnaseq-workshop/results/mapping/++g' |\
  grep -v "^#" |\
  sed 's/Aligned.sortedByCoord.out.bam//g' \
  > gene_counts_clean.txt

Step 6: QC on counts

---

Before proceeding to differential expression, it is good practice to perform some QC on the count data. We can again use MultiQC to summarize featureCounts statistics:

BASH

cd $SCRATCH/rnaseq-workshop
mkdir -p results/qc_counts
module load biocontainers
module load multiqc
multiqc results/counts -o results/qc_counts

Challenge

Exercise: Examine your counting metrics

Using your MultiQC featureCounts report:

1. Which sample has the lowest percent assigned, and what are common reasons for lower assignment?
2. Which unassigned category is the largest, and what does it indicate?
3. Does the sample with the most assigned reads also have the highest percent assigned?
4. What does “Unassigned: No Features” mean, and why might it occur?
5. Are multi-mapped reads a concern for RNA-seq differential expression?

Show me the solution

Interpretation for this dataset. Assigned fragments are about 32 to 40 percent of the total. This looks low, but the total in the featureCounts summary includes every record in the BAM file: STAR kept unmapped pairs (--outSAMunmapped Within) and multimapping pairs, and featureCounts reports both as unassigned. Measured against uniquely mapped pairs only, the assignment rate is much higher.

1. The sample with the lowest assignment rate here is mock_rep2 (32 percent). Causes include low complexity, incomplete annotation, or more intronic reads.
2. It differs between samples: Unmapped is largest for IR_rep1 and IR_rep2 (about 25 percent), MultiMapping for IR_rep3, IR_rep4, and mock_rep3 (24 to 26 percent), and No Features for mock_rep1, mock_rep2, and mock_rep4 (23 to 35 percent). Unmapped pairs did not align at all; No Features pairs aligned outside annotated exons; MultiMapping pairs aligned to several loci and are not counted by default.
3. Here it does: IR_rep4 has both the most assigned pairs (9.6 million) and the highest percentage (40 percent). In general they need not agree: total assigned reads depend on sequencing depth, while percent assigned reflects library and mapping quality.
4. Reads mapped outside annotated exons, often from intronic regions, unannotated transcripts, or incomplete annotation.
5. Multi-mapping is expected for paralogs and repeats. It is usually fine as long as the rate is consistent across samples.

Summary

---

Key Points

• STAR is a fast splice-aware aligner widely used for RNA-seq.
• Genome indexing requires FASTA and optionally GTF for improved splice detection.
• Mapping is efficiently performed using SLURM array jobs.
• Alignment statistics (from Log.final.out + MultiQC) must be reviewed.
• Strandness should be determined before counting; Salmon infers it directly from a sample of FASTQ reads.
• Final BAM files are ready for downstream counting.
• featureCounts is used to obtain gene-level counts for differential expression.
• Exons are counted and summed per gene.
• The output is a simple matrix for downstream statistical analysis.