library(ChIPseeker)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
library(org.Hs.eg.db)

txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene

# Read peaks from MACS3
peaks <- readPeakFile('sample_peaks.narrowPeak')

# Annotate with default settings
peak_anno <- annotatePeak(
    peaks,
    TxDb = txdb,
    annoDb = 'org.Hs.eg.db'
)

# View annotation summary
peak_anno

# Define promoter region (-3kb to +3kb from TSS)
peak_anno <- annotatePeak(
    peaks,
    TxDb = txdb,
    tssRegion = c(-3000, 3000),  # Promoter definition
    annoDb = 'org.Hs.eg.db'
)

# Convert to data frame
anno_df <- as.data.frame(peak_anno)

# Key columns: seqnames, start, end, annotation, distanceToTSS, SYMBOL, GENENAME
head(anno_df)

# Export to CSV
write.csv(anno_df, 'annotated_peaks.csv', row.names = FALSE)

# Filter for promoter peaks
promoter_peaks <- anno_df[grep('Promoter', anno_df$annotation), ]

# Get unique genes
promoter_genes <- unique(promoter_peaks$SYMBOL)

# Pie chart of genomic feature distribution
plotAnnoPie(peak_anno)

# Bar plot alternative
plotAnnoBar(peak_anno)

# Distribution of peaks relative to TSS
plotDistToTSS(peak_anno, title = 'Distribution of peaks relative to TSS')

# Read multiple peak files
peak_files <- list(
    H3K4me3 = 'H3K4me3_peaks.narrowPeak',
    H3K27ac = 'H3K27ac_peaks.narrowPeak',
    H3K27me3 = 'H3K27me3_peaks.broadPeak'
)

peak_list <- lapply(peak_files, readPeakFile)

# Annotate all
anno_list <- lapply(peak_list, annotatePeak, TxDb = txdb, annoDb = 'org.Hs.eg.db')

# Compare annotations
plotAnnoBar(anno_list)
plotDistToTSS(anno_list)

# Find overlapping peaks
genes_list <- lapply(anno_list, function(x) as.data.frame(x)$SYMBOL)
vennplot(genes_list)

# Plot peak coverage around TSS
covplot(peaks, weightCol = 'V5')  # V5 is score column in narrowPeak

# Get promoter coordinates
promoter <- getPromoters(TxDb = txdb, upstream = 3000, downstream = 3000)

# Get tag matrix
tagMatrix <- getTagMatrix(peaks, windows = promoter)

# Plot heatmap
tagHeatmap(tagMatrix, xlim = c(-3000, 3000), color = 'red')

# Average profile
plotAvgProf(tagMatrix, xlim = c(-3000, 3000), xlab = 'Distance from TSS')

library(clusterProfiler)

# Get genes from peaks
genes <- unique(anno_df$ENTREZID)

# GO enrichment
ego <- enrichGO(
    gene = genes,
    OrgDb = org.Hs.eg.db,
    ont = 'BP',
    pAdjustMethod = 'BH',
    pvalueCutoff = 0.05
)

# Find all genes overlapping peak regions (not just nearest)
genes_in_peaks <- seq2gene(peaks, tssRegion = c(-1000, 1000), flankDistance = 3000, TxDb = txdb)

# Mouse
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
library(org.Mm.eg.db)
peak_anno_mm <- annotatePeak(peaks, TxDb = TxDb.Mmusculus.UCSC.mm10.knownGene, annoDb = 'org.Mm.eg.db')

# Zebrafish
library(TxDb.Drerio.UCSC.danRer11.refGene)
library(org.Dr.eg.db)