Bio Chipseq Visualization

# Compute scale factor from spike-in reads (ChIP-Rx Drosophila or CUT&RUN E. coli)
SCALE=$(echo "scale=6; 1.0 / $SPIKE_IN_READS_M" | bc)  # 1 per million spike reads
bamCoverage -b chip.bam -o chip.bw --scaleFactor $SCALE --binSize 10
# DO NOT also pass --normalizeUsing; deepTools multiplies the two factors, reintroducing depth normalization

# Standard within-sample (CPM)
bamCoverage -b chip.bam -o chip.bw \
    --normalizeUsing CPM --binSize 10 \
    --extendReads 200 --numberOfProcessors 8

# Cross-sample at 1x genome coverage (ENCODE)
bamCoverage -b chip.bam -o chip.bw \
    --normalizeUsing RPGC --effectiveGenomeSize 2701495761 \
    --binSize 10 --extendReads 200

# ChIP vs Input log2 ratio (visualization of enrichment)
bamCompare -b1 chip.bam -b2 input.bam -o chip_vs_input.bw \
    --operation log2 --binSize 50 --extendReads 200 \
    --pseudocount 1 --skipZeroOverZero

# Compute matrix centered on TSS
computeMatrix reference-point \
    --referencePoint TSS \
    -b 3000 -a 3000 \
    -R genes.bed \
    -S chip.bw input.bw \
    -o matrix.gz \
    --outFileSortedRegions sorted_regions.bed \
    --numberOfProcessors 8 \
    --skipZeros

# Heatmap with k-means clustering (biology emerges from clusters)
plotHeatmap -m matrix.gz \
    -o heatmap.pdf \
    --kmeans 3 \
    --colorMap RdBu_r \
    --zMin -3 --zMax 3 \
    --refPointLabel TSS \
    --heatmapHeight 12 \
    --whatToShow 'heatmap and colorbar'

# Profile plot (average signal across regions)
plotProfile -m matrix.gz \
    -o profile.pdf \
    --perGroup \
    --plotTitle 'H3K4me3 around TSS'

computeMatrix scale-regions \
    -R genes.bed \
    -S chip.bw \
    -b 3000 -a 3000 \
    -m 5000 \
    -o matrix_genebody.gz \
    --numberOfProcessors 8

plotProfile -m matrix_genebody.gz -o genebody_profile.pdf --perGroup

multiBamSummary bins -b sample1.bam sample2.bam sample3.bam \
    --binSize 10000 -o results.npz \
    --numberOfProcessors 8

plotCorrelation -in results.npz \
    --corMethod spearman \
    --whatToPlot heatmap \
    --plotNumbers -o correlation.pdf \
    --outFileCorMatrix correlation.tab
# Replicates should correlate > 0.8 (narrow), > 0.6 (broad)

# tracks.ini
[x-axis]

[chip-h3k27ac]
file = h3k27ac.bw
color = darkblue
height = 3
title = H3K27ac

[chip-h3k4me3]
file = h3k4me3.bw
color = darkred
height = 3
title = H3K4me3

[peaks-narrowpeak]
file = peaks.narrowPeak
file_type = narrow_peak
color = black
height = 0.5
title = MACS peaks

[se-bed]
file = super_enhancers.bed
color = orange
height = 0.5
title = Super-enhancers

[genes]
file = genes.gtf
color = darkgreen
prefered_name = gene_name
height = 4

pyGenomeTracks --tracks tracks.ini --region chr1:1000000-1500000 -o region.pdf

library(Gviz)
library(GenomicRanges)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)

chr <- 'chr1'; start <- 1e6; end <- 1.1e6
itrack <- IdeogramTrack(genome = 'hg38', chromosome = chr)
gtrack <- GenomeAxisTrack()
dtrack <- DataTrack(range = 'sample.bw', genome = 'hg38',
                     type = 'histogram', name = 'ChIP', col.histogram = 'darkblue')
grtrack <- GeneRegionTrack(TxDb.Hsapiens.UCSC.hg38.knownGene,
                            genome = 'hg38', chromosome = chr, name = 'Genes')
plotTracks(list(itrack, gtrack, dtrack, grtrack), from = start, to = end, chromosome = chr)

library(EnrichedHeatmap)
library(rtracklayer)

# Normalize bigWig signal to a matrix around target sites
signal <- import('sample.bw')
tss <- promoters(txdb, upstream = 0, downstream = 1)
mat <- normalizeToMatrix(signal, tss, extend = 3000, mean_mode = 'w0', w = 50)

# Heatmap with customization
EnrichedHeatmap(mat, name = 'Signal', col = c('white', 'red'),
                top_annotation = HeatmapAnnotation(lines = anno_enriched()))

library(ChIPseeker)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)

peaks <- readPeakFile('peaks.narrowPeak')
promoter <- getPromoters(TxDb = TxDb.Hsapiens.UCSC.hg38.knownGene,
                          upstream = 3000, downstream = 3000)
tagMatrix <- getTagMatrix(peaks, windows = promoter)

# Tag heatmap and average profile
# tagHeatmap in ChIPseeker >= 1.36 takes palette (RColorBrewer name), not xlim/color;
# xlim is read from the tagMatrix window. plotAvgProf still uses xlim/conf.
tagHeatmap(tagMatrix, palette = 'Reds')
plotAvgProf(tagMatrix, xlim = c(-3000, 3000), conf = 0.95,
             xlab = 'Distance from TSS (bp)', ylab = 'Peak density')

# IGV batch script for reproducible screenshots
cat > igv.batch << 'EOF'
new
genome hg38
load chip.bw
load peaks.bed
load super_enhancers.bed
goto chr1:1000000-1100000
snapshot region1.png
goto chr2:50000000-51000000
snapshot region2.png
exit
EOF

igv.sh -b igv.batch