bio-chipseq-peak-annotation

## Version Compatibility Reference examples tested with: MACS3 3.0+, clusterProfiler 4.10+ Before using code patterns, verify installed versions match. If versions differ: - R: `packageVersion('<pkg>')` then `?function_name` to verify parameters If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying. # Peak Annotation with ChIPseeker **"Annotate my ChIP-seq peaks to genes"** → Assign peaks to genomic features (promoter, exon, intron, intergenic), find nearest genes, and calculate TSS distances. - R: `ChIPseeker::annotatePeak(peaks, TxDb=txdb)` ## Load Peaks and Annotations ```r 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 Peaks **Goal:** Assign each ChIP-seq peak to its nearest gene and genomic feature category. **Approach:** Use annotatePeak with a TxDb annotation database to classify peaks as promoter, exon, intron, or intergenic and retrieve the nearest gene symbol. ```r # Annotate with default settings peak_anno <- annotatePeak( peaks, TxDb = txdb, annoDb = 'org.Hs.eg.db' ) # View annotation summary peak_anno ``` ## Custom Promoter Definition ```r # 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' ) ``` ## Extract Annotated Data Frame ```r # 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) ``` ## Get Genes with Peaks in Promoter ```r # Filter for promoter peaks promoter_peaks <- anno_df[grep('Promoter', anno_df$annotation), ] # Get unique genes promoter_genes <- unique(promoter_peaks$SYMBOL) ``` ## Annotation Pie Chart ```r # Pie chart of genomic feature distribution plotAnnoPie(peak_anno) # Bar plot alternative plotAnnoBar(peak_anno) ``` ## Distance to TSS Plot ```r # Distribution of peaks relative to TSS plotDistToTSS(peak_anno, title = 'Distribution of peaks relative to TSS') ``` ## Compare Multiple Peak Sets **Goal:** Compare genomic feature distributions across multiple ChIP-seq experiments (e.g., different histone marks). **Approach:** Read and annotate each peak file separately, then use plotAnnoBar and plotDistToTSS on the annotation list for side-by-side comparison. ```r # 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) ``` ## Venn Diagram of Peak Overlap ```r # Find overlapping peaks genes_list <- lapply(anno_list, function(x) as.data.frame(x)$SYMBOL) vennplot(genes_list) ``` ## Coverage Plot ```r # Plot peak coverage around TSS covplot(peaks, weightCol = 'V5') # V5 is score column in narrowPeak ``` ## Profile Heatmap Around TSS **Goal:** Visualize the distribution of ChIP-seq signal around transcription start sites. **Approach:** Extract promoter regions from the TxDb, build a tag matrix of signal at those regions, and plot as a heatmap or average profile. ```r # 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') ``` ## Functional Enrichment of Peak Genes **Goal:** Determine which biological processes are enriched among genes with ChIP-seq peaks in their promoters. **Approach:** Extract Entrez IDs from annotated peaks and run GO enrichment analysis with clusterProfiler. ```r 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 ) ``` ## Seq2Gene - All Genes in Peak Regions ```r # Find all genes overlapping peak regions (not just nearest) genes_in_peaks <- seq2gene(peaks, tssRegion = c(-1000, 1000), flankDistance = 3000, TxDb = txdb) ``` ## Different Organisms ```r # 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) ``` ## Key Functions | Function | Purpose | |----------|---------| | readPeakFile | Read peak file (BED, narrowPeak) | | annotatePeak | Annotate peaks to genes | | plotAnnoPie | Pie chart of annotations | | plotAnnoBar | Bar plot of annotations | | plotDistToTSS | Distance to TSS distribution | | getPromoters | Get promoter regions | | getTagMatrix | Coverage matrix around regions | | tagHeatmap | Heatmap of signal | | plotAvgProf | Average profile plot | | seq2gene | Map peaks to all overlapping genes | ## Annotation Categories | Category | Description | |----------|-------------| | Promoter | Within tssRegion of TSS | | 5' UTR | 5' untranslated region | | 3' UTR | 3' untranslated region | | Exon | Coding exon | | Intron | Intronic region | | Downstream | Within 3kb downstream | | Distal Intergenic | Beyond gene regions | ## Related Skills - peak-calling - Generate peak files with MACS3 - differential-binding - Find differential peaks - pathway-analysis - Functional enrichment - chipseq-visualization - Additional visualizations