ChIP-seq & ATAC-seq

一、ChIP-seq vs ATAC-seq

二、ChIPseeker：peak 註解標準工具

# BiocManager::install(c('ChIPseeker','TxDb.Hsapiens.UCSC.hg38.knownGene','org.Hs.eg.db'))
# library(ChIPseeker); library(TxDb.Hsapiens.UCSC.hg38.knownGene)
# library(clusterProfiler); library(org.Hs.eg.db)

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

# 1) 讀 narrowPeak / BED 檔
# peaks <- readPeakFile('peaks/sample.narrowPeak')
# peaks      # GRanges

# 2) 註解 peak：每個 peak 對應到最近的基因 + genomic feature
# anno <- annotatePeak(peaks,
#                      TxDb = txdb,
#                      tssRegion = c(-3000, 3000),    # promoter 定義
#                      annoDb = 'org.Hs.eg.db',
#                      level = 'gene')
# anno_df <- as.data.frame(anno)
# head(anno_df[, c('seqnames','start','end','annotation','geneId','SYMBOL','distanceToTSS')])

# 3) 視覺化 / Visualize
# plotAnnoPie(anno)        # promoter / intron / intergenic / ... 圓餅
# plotAnnoBar(anno)
# vennpie(anno)
# upsetplot(anno)
# plotDistToTSS(anno)      # peak 距 TSS 距離分布

# 4) 多樣本比較 (e.g. 多個 ChIP)
# files <- list(s1 = 'peaks/s1.narrowPeak', s2 = 'peaks/s2.narrowPeak')
# annoList <- lapply(files, annotatePeak, TxDb = txdb, tssRegion = c(-3000, 3000))
# plotAnnoBar(annoList)

# 5) Peak heatmap：所有 peak 在 TSS 附近 ±3kb 的訊號分布
# promoter <- getPromoters(TxDb = txdb, upstream = 3000, downstream = 3000)
# tagMatrix <- getTagMatrix(peaks, windows = promoter)
# tagHeatmap(tagMatrix, xlim = c(-3000, 3000), color = 'firebrick')
# plotAvgProf(tagMatrix, xlim = c(-3000, 3000))

# 6) GO/KEGG 富集 (拿 anno 中的 geneId)
# entrez_ids <- unique(as.data.frame(anno)$geneId)
# ego <- enrichGO(entrez_ids, OrgDb = org.Hs.eg.db, ont = 'BP', readable = TRUE)
# dotplot(ego)

cat('ChIPseeker demo — install Bioc packages above to run on real peaks.\n')

三、DiffBind：差異結合分析

# BiocManager::install('DiffBind')
# library(DiffBind)

# 1) 準備 sample sheet (CSV，必填欄：SampleID, Condition, bamReads, Peaks, PeakCaller)
# samples <- read.csv('samples.csv')
# head(samples)
# #   SampleID  Tissue   Condition  Replicate  bamReads          Peaks                PeakCaller
# #   ctrl_1    HeLa     control    1          bam/ctrl_1.bam    peaks/ctrl_1.narrow  narrow

# 2) 建立 dba 物件
# dba <- dba(sampleSheet = samples)
# dba       # 摘要

# 3) 計算每個 peak 的 read count
# dba <- dba.count(dba, summits = 250)         # 以 peak summit ±250bp 為窗

# 4) 設定比較
# dba <- dba.contrast(dba, categories = DBA_CONDITION, minMembers = 2)

# 5) 執行差異分析（內部用 DESeq2）
# dba <- dba.analyze(dba)
# dba.show(dba, bContrasts = TRUE)

# 6) 取結果 / Extract results
# res <- dba.report(dba, contrast = 1, th = 0.05, fold = 1)
# res    # GRanges of differentially bound peaks
# as.data.frame(res)

# 7) 視覺化
# dba.plotMA(dba)            # MA plot
# dba.plotVolcano(dba)       # Volcano
# dba.plotHeatmap(dba, contrast = 1, correlations = FALSE)
# dba.plotPCA(dba, attributes = DBA_CONDITION)

cat('DiffBind demo — needs BAM + peak files + sample sheet to actually run.\n')

四、Motif 分析

# BiocManager::install(c('motifmatchr','JASPAR2024','TFBSTools',
#                        'BSgenome.Hsapiens.UCSC.hg38'))

# library(motifmatchr); library(JASPAR2024); library(TFBSTools)
# library(BSgenome.Hsapiens.UCSC.hg38)

# 1) 從 JASPAR 拿 motif PWMs
# opts <- list(species = 9606, collection = 'CORE')   # 9606 = human
# pwms <- getMatrixSet(JASPAR2024, opts)
# length(pwms)

# 2) 對 peaks 做 motif 匹配
# bs <- BSgenome.Hsapiens.UCSC.hg38
# motif_matches <- matchMotifs(pwms, peaks, genome = bs)
# motif_matches  # SummarizedExperiment of TRUE/FALSE matches

# 3) 計算每個 peak 中匹配到的 motif 數量 + 對 TF activity (chromVAR)
# library(chromVAR)
# dev <- computeDeviations(object = se_atac,         # ATAC counts SE
#                          annotations = motif_matches)
# variability <- computeVariability(dev)
# head(variability[order(-variability$variability), ])

cat('Motif demo — install motifmatchr + JASPAR2024 + chromVAR for full workflow.\n')

五、在 R 中畫 IGV-like 基因組軌道

# BiocManager::install(c('Gviz','rtracklayer'))
# library(Gviz); library(rtracklayer)
# library(TxDb.Hsapiens.UCSC.hg38.knownGene)

# # 區域定義 / Define region
# chr <- 'chr17'; start <- 7670000; end <- 7700000

# # 基因軌道 / Gene track
# txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
# gene_track <- GeneRegionTrack(txdb, chromosome = chr, start = start, end = end,
#                                name = 'Genes', transcriptAnnotation = 'symbol')

# # bigWig 訊號軌道 (e.g. ChIP-seq, ATAC-seq)
# bw1 <- DataTrack(range = 'data/sample1.bw', chromosome = chr,
#                  start = start, end = end, name = 'ChIP', type = 'h')
# bw2 <- DataTrack(range = 'data/sample2.bw', chromosome = chr,
#                  start = start, end = end, name = 'ATAC', type = 'h')

# # Peak 軌道
# peaks <- import('peaks/sample.narrowPeak')
# peak_track <- AnnotationTrack(peaks, chromosome = chr,
#                                start = start, end = end, name = 'Peaks')

# # 座標軸 + 染色體圖示
# axis  <- GenomeAxisTrack()
# ideo  <- IdeogramTrack(genome = 'hg38', chromosome = chr)

# # 一次畫所有軌道
# plotTracks(list(ideo, axis, gene_track, bw1, bw2, peak_track),
#            from = start, to = end)

# # 進階替代：karyoploteR、ggbio
cat('Gviz demo — needs bigWig / peak files. Replace paths to actually run.\n')

六、scATAC-seq 工具：Signac / ArchR

# 詳見 scRNA-seq tutorial 對應章節（參考 sister 教程）
# 標準流程概念：
#   1) Cell Ranger ATAC → fragments.tsv.gz + filtered_peak_bc_matrix
#   2) Signac::CreateChromatinAssay → 整合進 Seurat
#   3) QC (TSS enrichment, nucleosome signal)
#   4) LSI dimensionality reduction → UMAP → clustering
#   5) GeneActivity matrix → 與 scRNA 整合
#   6) Motif analysis (chromVAR) → cell-type TF activity

cat('scATAC demo — see Signac vignette for full single-cell pipeline.\n')