1. Generating Normalized BigWig Files with deepTools bamCoverage

To compare signal tracks across samples with varying sequencing depths, raw BAM files must be normalized. We use deepTools bamCoverage to generate BigWig tracks normalized to 1x genome coverage ($1\times \text{depth} = \text{BPM}$ / Reads Per Genomic Content):

for sample in HD_Rep1 HD_Rep2 Control_Rep1 Control_Rep2; do
    echo "Generating normalized BigWig for ${sample}..."
    bamCoverage \
        -b 03_filtered/${sample}.clean.bam \
        -o 07_visualization/${sample}.bw \
        --binSize 20 \
        --normalizeUsing BPM \
        --effectiveGenomeSize 2913022398 \
        --ignoreForNormalization chrX chrY chrM \
        --numberOfProcessors 8
done

2. Computing the TSS Signal Matrix with computeMatrix

We calculate H3K4me3 signal density within $\pm 3\text{kb}$ windows surrounding all human Transcription Start Sites (TSS):

# Compute signal matrix around TSS reference points
computeMatrix reference-point \
    --referencePoint TSS \
    -b 3000 -a 3000 \
    -R gencode.v38.annotation.gtf \
    -S 07_visualization/Control_Rep1.bw 07_visualization/Control_Rep2.bw \
       07_visualization/HD_Rep1.bw 07_visualization/HD_Rep2.bw \
    --skipZeros \
    -o 07_visualization/matrix_tss.gz \
    -p 8

3. Plotting Publication-Grade Heatmaps with plotHeatmap

plotHeatmap \
    -m 07_visualization/matrix_tss.gz \
    -out 07_visualization/h3k4me3_tss_heatmap.pdf \
    --colorMap RdYlBu_r \
    --samplesLabel "Control Rep1" "Control Rep2" "HD Rep1" "HD Rep2" \
    --regionsLabel "All Genes" \
    --kmeans 3

plotProfile \
    -m 07_visualization/matrix_tss.gz \
    -out 07_visualization/h3k4me3_tss_composite_profile.pdf \
    --samplesLabel "Control 1" "Control 2" "HD 1" "HD 2" \
    --plotTitle "H3K4me3 Enrichment at Promoters (HD vs Control)"