Lab 4: Read Alignment with STAR

Connect to Ibex and start an interactive session

# Connect to Ibex
ssh username@ilogin.ibex.kaust.edu.sa

# Start an interactive session with 8 CPUs and 8 GB RAM
# (chr11 STAR index is ~4 GB; never run analyses on the login node)
srun --cpus-per-task=8 --mem=8G --time=01:00:00 --pty /bin/bash

Load required modules

# Load alignment and BAM processing tools
module load star
module load samtools
module load picard

# Verify that all tools are loaded
module list

# Check versions
STAR --version
samtools --version | head -1
picard MarkDuplicates --version 2>&1 | head -1

Navigate to your workshop directory and create the output folder

# Navigate to your workshop directory on scratch
cd /ibex/user/$USER/workshop

# Create the directory for alignment outputs
mkdir -p aligned logs

# Confirm your directory structure
ls -la

/ibex/user/$USER/workshop/
├── raw_data/          # original full SRA FASTQs and SRR_Acc_List.txt
├── chr11_raw_data/    # chr11 FASTQs from Zenodo (KO/WT naming)
├── trimmed/           # fastp-trimmed chr11 FASTQs (*_trimmed_1/2.fastq.gz) from Lab 3
├── references/        # chr11 FASTA, GTF, and STAR index
├── aligned/           # will be created now
└── logs/              # SLURM log files

Copy the chr11 reference files from the shared course directory

cd /ibex/user/$USER/workshop

SHARED=/biocorelab/BIX/resources/genomes/workshops/human-chr11
mkdir -p references/chr11

# Chr11 genome FASTA and GTF
cp ${SHARED}/GRCh38.dna.chromosome.11.fa   references/chr11/
cp ${SHARED}/GRCh38.chr11.gtf  references/chr11/
cp ${SHARED}/GRCh38.cdna.chr11.fa  references/chr11/

# Verify
ls -lh references/chr11/

Build the STAR index for chr11 (interactive)

cd /ibex/user/$USER/workshop

mkdir -p references/star_index_chr11

STAR --runMode genomeGenerate \
    --genomeDir references/star_index_chr11 \
    --genomeFastaFiles references/chr11/GRCh38.chr11.fa \
    --sjdbGTFfile references/chr11/GRCh38.chr11.gtf \
    --runThreadN 8 \
    --genomeSAindexNbases 11

Index building takes approximately 3–5 minutes. Verify it completed:

ls -lh references/star_index_chr11/
# Expect: Genome, SA, SAindex, chrName.txt, genomeParameters.txt, etc.

Run STAR alignment for all samples

cd /ibex/user/$USER/workshop

for r1 in trimmed/*_trimmed_1.fastq.gz; do
    sample=$(basename "$r1" _trimmed_1.fastq.gz)
    echo "=== Aligning $sample ==="
    mkdir -p aligned/${sample}

    STAR --runThreadN 8 \
        --genomeDir references/star_index_chr11 \
        --readFilesIn trimmed/${sample}_trimmed_1.fastq.gz trimmed/${sample}_trimmed_2.fastq.gz \
        --readFilesCommand zcat \
        --outSAMtype BAM SortedByCoordinate \
        --outSAMattributes NH HI AS NM MD \
        --outFileNamePrefix aligned/${sample}/ \
        --outSAMstrandField intronMotif \
        --alignSJoverhangMin 8 \
        --alignSJDBoverhangMin 1 \
        --quantMode GeneCounts \
        --runMode alignReads

    echo "Done: $sample"
done

Understand the key STAR alignment parameters

Check expected output files for each sample

# Check the output directory for one sample
ls -lh aligned/KO_1_SRR10045016/

Inspect a BAM file from your alignment output

Use the BAM you produced in Part 3 to explore the format directly.

cd /ibex/user/$USER/workshop

BAM=aligned/KO_1_SRR10045016/Aligned.sortedByCoord.out.bam

# View the full BAM header (@HD, @SQ, @PG lines)
samtools view -H $BAM

# View the first 5 alignment records in SAM text format
samtools view $BAM | head -5

# Summary alignment statistics
samtools flagstat $BAM

# Count all reads in the BAM
samtools view -c $BAM

# Count only mapped reads (exclude FLAG bit 0x4 = unmapped)
samtools view -c -F 4 $BAM

Examine the CIGAR string and optional tags

Pick one alignment line and decode its CIGAR string to understand how the read spans a splice junction.

# Show one alignment record with all optional tags, formatted for readability
samtools view $BAM | head -1 | tr '\t' '\n' | cat -n

# Find reads that span a splice junction (N in CIGAR = intron skip)
samtools view $BAM | awk '$6 ~ /N/' | head -3

Index the BAM files with samtools

A BAM index (.bai file) allows tools like IGV and samtools to jump directly to any genomic region without reading the entire file.

cd /ibex/user/$USER/workshop

# Index the BAM file for every sample
for bam in aligned/*/Aligned.sortedByCoord.out.bam; do
    echo "Indexing $bam ..."
    samtools index "$bam"
done

# Verify — each BAM should now have a matching .bai file
ls -lh aligned/KO_1_SRR10045016/

Mark PCR duplicates with Picard

for bam in aligned/*/Aligned.sortedByCoord.out.bam; do
    sample=$(basename $(dirname "$bam"))
    echo "Marking duplicates for $sample ..."
    picard MarkDuplicates \
        I=aligned/${sample}/Aligned.sortedByCoord.out.bam \
        O=aligned/${sample}/marked_dup.bam \
        M=aligned/${sample}/dup_metrics.txt
    samtools index aligned/${sample}/marked_dup.bam
done

Inspect the Picard duplication metrics

# View the duplication report for one sample
cat aligned/KO_1_SRR10045016/dup_metrics.txt

Run samtools flagstat on all samples

cd /ibex/user/$USER/workshop

# Print flagstat output for every sample
for bam in aligned/*/Aligned.sortedByCoord.out.bam; do
    sample=$(basename $(dirname "$bam"))
    echo "=== $sample ==="
    samtools flagstat "$bam"
    echo ""
done

Inspect the STAR alignment summary log

# View the STAR summary log for the first sample
cat aligned/KO_1_SRR10045016/Log.final.out

# Loop to print key stats for all samples
for log in aligned/*/Log.final.out; do
    sample=$(basename $(dirname "$log"))
    echo "=== $sample ==="
    grep -E "Uniquely mapped|mapped to multiple|unmapped" "$log"
    echo ""
done

To run IGV locally, copy BAM files to your local machine or mount Ibex on your local machine

# Run these commands from your LOCAL terminal (not Ibex)
# Replace "username" and "SRRxxxxxxxx" with real values

scp username@ilogin.ibex.kaust.edu.sa:/ibex/user/username/workshop/aligned/KO_1_SRR10045016/Aligned.sortedByCoord.out.bam ./
scp username@ilogin.ibex.kaust.edu.sa:/ibex/user/username/workshop/aligned/KO_1_SRR10045016/Aligned.sortedByCoord.out.bam.bai ./

Download IGV from https://igv.org/doc/desktop/ if not already installed.

Load the genome and BAM file in IGV

1. In IGV, select Genomes → Load Genome from Server and choose Human (GRCh38/hg38)
2. Select File → Load from File and open the BAM file (Aligned.sortedByCoord.out.bam)
3. In the search box at the top, type GAPDH and press Enter to navigate to that gene

Navigate to the TARDBP gene locus

Since this dataset is a TDP-43 (TARDBP) study, navigate to the gene of interest:

# Type this gene name in the IGV search box to navigate:
TARDBP

# You can also navigate by coordinates:
# chr1:11,012,655-11,045,430  (GRCh38)