Lab 4a: Genome Alignment (STAR + Salmon)

Create the STAR index

cd ~/genomics

# Activate environment
conda activate genomics

# Create directory for STAR index
mkdir -p references/star_index

# Build STAR index for chromosome 11
# This takes 5-10 minutes and ~4 GB RAM
STAR --runMode genomeGenerate \
    --runThreadN 4 \
    --genomeDir references/star_index \
    --genomeFastaFiles references/Homo_sapiens.GRCh38.dna.chromosome.11.fa \
    --sjdbGTFfile references/Homo_sapiens.GRCh38.110.chr11.gtf \
    --sjdbOverhang 69 \
    --genomeSAindexNbases 11

# Check the generated index files
ls -lh references/star_index/

Align one sample

cd ~/genomics

# Align KO_1 sample
STAR --runThreadN 4 \
    --genomeDir references/star_index \
    --readFilesIn trimmed_data/KO_1_SRR10045016_1.trimmed.fastq.gz \
                  trimmed_data/KO_1_SRR10045016_2.trimmed.fastq.gz \
    --outFileNamePrefix alignment/KO_1_ \
    --outSAMattributes Standard \
    --readFilesCommand zcat \
    --outSAMtype BAM SortedByCoordinate \
    --quantMode GeneCounts

# Check output files
ls -lh alignment/KO_1_*

Align all samples

cd ~/genomics

# Define samples
SAMPLES="KO_1_SRR10045016 KO_2_SRR10045017 KO_3_SRR10045018 WT_1_SRR10045019 WT_2_SRR10045020 WT_3_SRR10045021"

# Align each sample
for SAMPLE in $SAMPLES; do
    echo "Aligning $SAMPLE..."

    STAR --runThreadN 4 \
        --genomeDir references/star_index \
        --readFilesIn trimmed_data/${SAMPLE}_1.trimmed.fastq.gz \
                      trimmed_data/${SAMPLE}_2.trimmed.fastq.gz \
        --outFileNamePrefix alignment/${SAMPLE}_ \
        --outSAMattributes Standard \
        --readFilesCommand zcat \
        --outSAMtype BAM SortedByCoordinate \
        --quantMode GeneCounts
        2>> logs/star.log
done

echo "Alignment complete!"
ls -lh alignment/*.bam

Review alignment log

cd ~/genomics

# View the alignment summary (Log.final.out)
cat alignment/KO_1_Log.final.out

# Key metrics to look for:
# - Uniquely mapped reads %: Should be >80% for good data
# - % of reads mapped to multiple loci: High values may indicate repetitive sequences
# - % of reads unmapped: Should be relatively low

Compare alignment rates across samples

Go through each file manually, see alignment/sampleName_Log.final.out and look for Sample,Input_Reads,Uniquely_Mapped,Unique%,Multi_Mapped,Unmapped rows

Index BAM files with samtools

cd ~/genomics/

# Index all BAM files (required for IGV and many tools)
for bam in alignment/*_Aligned.sortedByCoord.out.bam; do
    echo "Indexing $bam..."
    samtools index $bam
done

# Verify index files were created
ls -lh alignment/*.bai

Explore BAM file structure

# View BAM header
samtools view -H alignment/KO_1_Aligned.sortedByCoord.out.bam | head -20

# View first 5 aligned reads
samtools view alignment/KO_1_Aligned.sortedByCoord.out.bam | head -5

# Count total aligned reads
samtools view -c alignment/KO_1_Aligned.sortedByCoord.out.bam

# Count properly paired reads
samtools view -c -f 2 alignment/KO_1_Aligned.sortedByCoord.out.bam

Generate alignment statistics

# Detailed statistics with samtools flagstat
samtools flagstat alignment/KO_1_Aligned.sortedByCoord.out.bam

# Generate stats for all samples
for bam in alignment/*_Aligned.sortedByCoord.out.bam; do
    echo "--- $bam ---"
    samtools flagstat $bam | head -5
    echo ""
done

View in IGV (Integrative Genomics Viewer)

# Launch IGV (if installed via conda)
igv &

# In IGV:
# 1. Genomes > Load Genome from File > references/Homo_sapiens.GRCh38.dna.chromosome.11.fa
# 2. File > Load from File > references/references/Homo_sapiens.GRCh38.110.chr11.gtf
# 3. File > Load from File > alignment/KO_1_Aligned.sortedByCoord.out.bam
# 4. Navigate to a gene: e.g., search for "MDK" (Midkine gene on chr11)
# 5. Check reads mapping to each exon in "Midkine" gene
# 6. Load BAM file for a WT sample and see if there are any differences with the KO sample above