组蛋白修饰类型与功能¶

# === ChIP-seq数据分析流程 ===

# 1. 质控与比对
fastqc H3K27ac_chip.fastq.gz                    # 质控
trim_galore --paired H3K27ac_R1.fq.gz H3K27ac_R2.fq.gz  # 去接头

# 比对到参考基因组
bowtie2 -x hg38_index \
  -1 H3K27ac_R1_trimmed.fq.gz \
  -2 H3K27ac_R2_trimmed.fq.gz \
  -p 8 \                           # 8个线程
  --very-sensitive \                # 高灵敏度模式
  -S H3K27ac.sam                    # 输出SAM文件

# SAM转BAM，排序，去重复
samtools sort H3K27ac.sam -o H3K27ac.sorted.bam   # 排序
samtools index H3K27ac.sorted.bam                  # 建索引
picard MarkDuplicates \
  I=H3K27ac.sorted.bam \
  O=H3K27ac.dedup.bam \
  M=dup_metrics.txt \
  REMOVE_DUPLICATES=true            # 去PCR重复

# 2. 峰检测（MACS3）
macs3 callpeak \
  -t H3K27ac.dedup.bam \           # ChIP样本
  -c input.dedup.bam \              # Input对照
  -f BAMPE \                        # 双端测序
  -g hs \                           # 人类基因组大小
  -n H3K27ac \                      # 输出前缀
  --broad \                         # 宽峰模式（组蛋白修饰推荐）
  --broad-cutoff 0.1 \              # 宽峰q值阈值
  -q 0.05                           # 窄峰q值阈值

# 3. 生成BigWig信号文件（用于可视化）
bamCoverage \
  -b H3K27ac.dedup.bam \
  -o H3K27ac.bw \
  --normalizeUsing RPKM \           # RPKM归一化
  --binSize 10 \                    # 10bp分辨率
  -p 8                              # 8个线程

# === 组蛋白修饰下游分析 ===
import pandas as pd
import matplotlib.pyplot as plt

# 1. 读取MACS3的peak文件
peaks = pd.read_csv("H3K27ac_peaks.broadPeak", sep="\t", header=None,
                     names=["chr","start","end","name","score",
                            "strand","fc","pvalue","qvalue"])
print(f"H3K27ac峰数量: {len(peaks)}")
print(f"平均峰宽: {(peaks['end'] - peaks['start']).mean():.0f} bp")

# 2. 峰的基因组分布统计
# 需要用ChIPseeker或HOMER进行注释
# 这里展示结果解读框架
annotation_results = {
    "Promoter (<=1kb)": 35,    # 启动子区域占比
    "Promoter (1-2kb)": 8,
    "Intron": 30,              # 内含子（可能是增强子）
    "Intergenic": 20,          # 基因间区（远端增强子）
    "Exon": 5,
    "3'UTR": 2
}

plt.figure(figsize=(8, 6))
plt.pie(annotation_results.values(), labels=annotation_results.keys(),
        autopct='%1.1f%%')
plt.title("H3K27ac Peak Genomic Distribution")
plt.savefig("peak_annotation.png", dpi=150)

# 3. 不同组蛋白修饰的信号分布对比
# 用deepTools的computeMatrix和plotHeatmap
# computeMatrix reference-point -S H3K4me3.bw H3K27ac.bw H3K27me3.bw \
#   -R genes.bed --referencePoint TSS -a 3000 -b 3000 -o matrix.gz
# plotHeatmap -m matrix.gz -o heatmap.png --colorMap RdBu_r

# === 使用ChIPseeker进行峰注释 ===
library(ChIPseeker)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)

# 1. 读取peak文件
peak <- readPeakFile("H3K27ac_peaks.broadPeak")

# 2. 注释到基因组特征
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
peakAnno <- annotatePeak(peak, TxDb=txdb,
                          tssRegion=c(-3000, 3000))  # TSS上下游3kb

# 3. 可视化注释结果
plotAnnoPie(peakAnno)       # 基因组区域饼图
plotDistToTSS(peakAnno)     # 到TSS的距离分布