ChIP-seq分析全流程¶

# 质控
fastqc -t 8 -o qc/ *.fastq.gz
multiqc qc/ -o multiqc_report/

# 接头去除
fastp --in1 chip_R1.fastq.gz --in2 chip_R2.fastq.gz \
  --out1 chip_clean_R1.fastq.gz --out2 chip_clean_R2.fastq.gz \
  --json chip_fastp.json --thread 8

# 同样处理Input对照
fastp --in1 input_R1.fastq.gz --in2 input_R2.fastq.gz \
  --out1 input_clean_R1.fastq.gz --out2 input_clean_R2.fastq.gz \
  --json input_fastp.json --thread 8

# 建索引（只需一次）
bowtie2-build /ref/GRCh38.fa /ref/bt2_hg38

# 比对（single-end）
bowtie2 -x /ref/bt2_hg38 \
  -U chip_clean.fastq.gz \
  -S chip.sam \
  -p 16 \
  --very-sensitive \
  --no-mixed \
  --no-discordant

# 比对（paired-end）
bowtie2 -x /ref/bt2_hg38 \
  -1 chip_clean_R1.fastq.gz \
  -2 chip_clean_R2.fastq.gz \
  -S chip.sam \
  -p 16 \
  --very-sensitive \
  -X 1000  # 最大fragment长度

# 转BAM + 排序
samtools view -bS -q 20 -F 1804 chip.sam | \
  samtools sort -@ 8 -o chip_sorted.bam -
samtools index chip_sorted.bam

# 去重
picard MarkDuplicates \
  I=chip_sorted.bam \
  O=chip_dedup.bam \
  M=chip_dup_metrics.txt \
  REMOVE_DUPLICATES=true
samtools index chip_dedup.bam

# 移除blacklist区域
bedtools intersect -v -abam chip_dedup.bam \
  -b /ref/hg38-blacklist.v2.bed > chip_final.bam
samtools index chip_final.bam

# 同样处理Input对照
bowtie2 -x /ref/bt2_hg38 -U input_clean.fastq.gz -S input.sam -p 16 --very-sensitive
samtools view -bS -q 20 -F 1804 input.sam | samtools sort -@ 8 -o input_sorted.bam -
picard MarkDuplicates I=input_sorted.bam O=input_dedup.bam M=input_dup.txt REMOVE_DUPLICATES=true
bedtools intersect -v -abam input_dedup.bam -b /ref/hg38-blacklist.v2.bed > input_final.bam
samtools index input_final.bam

# 比对统计
samtools flagstat chip_final.bam

# 注意：以下命令中macs2可替换为macs3（核心callpeak算法一致，结果几乎相同）
# MACS3新增功能：HMMRATAC（ATAC-seq专用）、cutoff-analysis、callvar（变异检测）
# 安装：pip install macs3

# === Narrow peak（转录因子）===
macs2 callpeak \
  -t chip_final.bam \        # treatment（ChIP样本）
  -c input_final.bam \       # control（Input对照）
  -f BAM \                    # 输入格式
  -g hs \                     # 基因组大小(hs=人,mm=鼠)
  -n TF_sample1 \             # 输出前缀
  --outdir peaks/ \
  -q 0.05 \                   # FDR阈值
  --keep-dup all \            # 已去重，保留所有reads
  --call-summits              # 输出summit坐标

# === Broad peak（组蛋白修饰 H3K27me3/H3K36me3）===
macs2 callpeak \
  -t chip_H3K27me3.bam \
  -c input_final.bam \
  -f BAM \
  -g hs \
  -n H3K27me3_sample1 \
  --outdir peaks/ \
  --broad \                   # broad peak模式
  --broad-cutoff 0.1 \       # broad peak FDR
  -q 0.05 \
  --keep-dup all

# === Paired-end数据 ===
macs2 callpeak \
  -t chip_final.bam \
  -c input_final.bam \
  -f BAMPE \                  # paired-end BAM
  -g hs \
  -n TF_PE_sample1 \
  --outdir peaks/ \
  -q 0.05

# === 生成BigWig信号轨迹（用于可视化）===
# 使用deeptools
bamCoverage -b chip_final.bam \
  -o chip_signal.bw \
  --binSize 10 \
  --normalizeUsing RPGC \
  --effectiveGenomeSize 2913022398 \
  --extendReads 200 \
  -p 16

# Input归一化的信号（log2ratio）
bamCompare -b1 chip_final.bam \
  -b2 input_final.bam \
  -o chip_vs_input.bw \
  --binSize 50 \
  --normalizeUsing RPGC \
  --effectiveGenomeSize 2913022398 \
  --scaleFactorsMethod readCount \
  -p 16

# 统计peak数量和FRiP
echo "Peak count: $(wc -l < peaks/TF_sample1_peaks.narrowPeak)"

# 计算FRiP（Fraction of Reads in Peaks）
total_reads=$(samtools view -c chip_final.bam)
reads_in_peaks=$(bedtools intersect -a chip_final.bam -b peaks/TF_sample1_peaks.narrowPeak -u | samtools view -c)
echo "FRiP: $(echo "scale=4; $reads_in_peaks/$total_reads" | bc)"

# === IDR分析 ===
# 需要每个重复独立call的peak（按p-value排序）
# 注意：IDR输入需要用-p (p-value cutoff较宽松，如0.01)以获得更多候选peak

macs2 callpeak -t rep1.bam -c input.bam -f BAM -g hs -n rep1 -p 0.01 --keep-dup all
macs2 callpeak -t rep2.bam -c input.bam -f BAM -g hs -n rep2 -p 0.01 --keep-dup all

# 按p-value排序
sort -k8,8nr rep1_peaks.narrowPeak > rep1_sorted.narrowPeak
sort -k8,8nr rep2_peaks.narrowPeak > rep2_sorted.narrowPeak

# 运行IDR
idr --samples rep1_sorted.narrowPeak rep2_sorted.narrowPeak \
    --input-file-type narrowPeak \
    --rank p.value \
    --output-file idr_output.txt \
    --plot \
    --log-output-file idr.log

# IDR < 0.05的peak是高置信度peak
awk '$5 >= 540' idr_output.txt | wc -l  # IDR 0.05对应-log10编码的540

# === deepTools质控可视化 ===
# 多样本相关性热图
multiBamSummary bins -b rep1.bam rep2.bam input.bam \
  -o results.npz -p 16
plotCorrelation -in results.npz --corMethod spearman \
  -o correlation_heatmap.pdf --plotFileFormat pdf

# 指纹图（Fingerprint plot）
plotFingerprint -b chip_rep1.bam chip_rep2.bam input.bam \
  --labels "ChIP_rep1" "ChIP_rep2" "Input" \
  -o fingerprint.pdf -p 16

# TSS处的信号富集
computeMatrix reference-point --referencePoint TSS \
  -b 3000 -a 3000 \
  -R /ref/genes.bed \
  -S chip_signal.bw \
  -o matrix.gz -p 16

plotHeatmap -m matrix.gz -o tss_heatmap.pdf
plotProfile -m matrix.gz -o tss_profile.pdf

library(DiffBind)

# 创建样本信息表
samples <- data.frame(
  SampleID = c("Ctrl_rep1", "Ctrl_rep2", "Treat_rep1", "Treat_rep2"),
  Condition = c("Control", "Control", "Treatment", "Treatment"),
  Replicate = c(1, 2, 1, 2),
  bamReads = c("ctrl_rep1.bam", "ctrl_rep2.bam", "treat_rep1.bam", "treat_rep2.bam"),
  bamControl = c("input.bam", "input.bam", "input.bam", "input.bam"),
  Peaks = c("ctrl_rep1_peaks.narrowPeak", "ctrl_rep2_peaks.narrowPeak",
            "treat_rep1_peaks.narrowPeak", "treat_rep2_peaks.narrowPeak"),
  PeakCaller = rep("narrow", 4)
)

# 创建DiffBind对象
dba_obj <- dba(sampleSheet = samples)

# 计算读数（在consensus peaks中计数）
# DiffBind 3.x默认使用summarizeOverlaps，无需显式设置bUseSummarizeOverlaps
dba_obj <- dba.count(dba_obj)

# 归一化
dba_obj <- dba.normalize(dba_obj)

# 建立对比（DiffBind 3.x API）
dba_obj <- dba.contrast(dba_obj, categories = DBA_CONDITION)

# 差异分析（使用DESeq2）
dba_obj <- dba.analyze(dba_obj, method = DBA_DESEQ2)

# 获取结果
db_results <- dba.report(dba_obj, th = 0.05)
print(db_results)

# 可视化
# MA plot
dba.plotMA(dba_obj)

# 火山图
dba.plotVolcano(dba_obj)

# PCA
dba.plotPCA(dba_obj, DBA_CONDITION)

# 热图
dba.plotHeatmap(dba_obj, contrast = 1, correlations = FALSE)

# 导出差异peak
export(db_results, "diff_binding_peaks.bed")

# === 安装HOMER ===
# conda install -c bioconda homer
# 或手动：
# perl /path/to/homer/configureHomer.pl -install homer
# perl /path/to/homer/configureHomer.pl -install hg38

# === Motif发现 ===
# 对peak summit区域（±100bp）做motif分析
findMotifsGenome.pl peaks/TF_sample1_summits.bed \
  hg38 \
  motif_results/ \
  -size 200 \          # summit周围±100bp
  -mask \              # 屏蔽重复序列
  -p 16 \             # 线程数
  -preparsedDir /tmp/homer_preparsed/

# 对narrowPeak全区域做motif
findMotifsGenome.pl peaks/TF_sample1_peaks.narrowPeak \
  hg38 \
  motif_full_results/ \
  -size given \        # 使用peak的实际大小
  -mask \
  -p 16

# 自定义背景（如：所有开放染色质区域）
findMotifsGenome.pl diff_peaks.bed \
  hg38 \
  motif_custom_bg/ \
  -size 200 \
  -bg all_peaks.bed \  # 自定义背景
  -mask -p 16

# === Peak注释（HOMER）===
annotatePeaks.pl peaks/TF_sample1_peaks.narrowPeak \
  hg38 \
  -annStats peak_annotation_stats.txt \
  > peak_annotations.txt

# 注释统计包含：promoter/exon/intron/intergenic比例

# === 基因组分布可视化（R: ChIPseeker）===

# R: ChIPseeker peak注释
library(ChIPseeker)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
library(org.Hs.eg.db)

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

# 读取peak文件
peaks <- readPeakFile("peaks/TF_sample1_peaks.narrowPeak")

# 注释
peakAnno <- annotatePeak(peaks, TxDb = txdb, annoDb = "org.Hs.eg.db",
                          tssRegion = c(-3000, 3000))

# 可视化
plotAnnoPie(peakAnno)           # 饼图
plotAnnoBar(peakAnno)           # 条形图
plotDistToTSS(peakAnno)         # 到TSS距离分布

# 导出注释结果
anno_df <- as.data.frame(peakAnno)
write.csv(anno_df, "peak_annotation_table.csv")

# GO/KEGG富集分析（peak关联的基因）
library(clusterProfiler)
genes <- anno_df$geneId[!is.na(anno_df$geneId)]
ego <- enrichGO(gene = genes, OrgDb = org.Hs.eg.db, ont = "BP", readable = TRUE)
dotplot(ego, showCategory = 20)

# ===== 完整ChIP-seq流水线 =====

SAMPLE="TF_rep1"
INPUT="Input_rep1"
REF="/ref/GRCh38.fa"
BT2_IDX="/ref/bt2_hg38"
BLACKLIST="/ref/hg38-blacklist.v2.bed"
THREADS=16

# 1. 质控+trim
fastp --in1 ${SAMPLE}.fastq.gz --out1 ${SAMPLE}_clean.fastq.gz --thread 8 --json ${SAMPLE}.json
fastp --in1 ${INPUT}.fastq.gz --out1 ${INPUT}_clean.fastq.gz --thread 8 --json ${INPUT}.json

# 2. 比对
bowtie2 -x ${BT2_IDX} -U ${SAMPLE}_clean.fastq.gz -p ${THREADS} --very-sensitive | \
  samtools view -bS -q 20 -F 1804 | samtools sort -@ 8 -o ${SAMPLE}_sorted.bam -
bowtie2 -x ${BT2_IDX} -U ${INPUT}_clean.fastq.gz -p ${THREADS} --very-sensitive | \
  samtools view -bS -q 20 -F 1804 | samtools sort -@ 8 -o ${INPUT}_sorted.bam -

# 3. 去重
picard MarkDuplicates I=${SAMPLE}_sorted.bam O=${SAMPLE}_dedup.bam M=${SAMPLE}_dup.txt REMOVE_DUPLICATES=true
picard MarkDuplicates I=${INPUT}_sorted.bam O=${INPUT}_dedup.bam M=${INPUT}_dup.txt REMOVE_DUPLICATES=true

# 4. 去blacklist
bedtools intersect -v -abam ${SAMPLE}_dedup.bam -b ${BLACKLIST} > ${SAMPLE}_final.bam
bedtools intersect -v -abam ${INPUT}_dedup.bam -b ${BLACKLIST} > ${INPUT}_final.bam
samtools index ${SAMPLE}_final.bam; samtools index ${INPUT}_final.bam

# 5. Peak calling
macs2 callpeak -t ${SAMPLE}_final.bam -c ${INPUT}_final.bam \
  -f BAM -g hs -n ${SAMPLE} --outdir peaks/ -q 0.05 --call-summits --keep-dup all

# 6. BigWig
bamCoverage -b ${SAMPLE}_final.bam -o ${SAMPLE}.bw --binSize 10 \
  --normalizeUsing RPGC --effectiveGenomeSize 2913022398 --extendReads 200 -p ${THREADS}

# 7. Motif
findMotifsGenome.pl peaks/${SAMPLE}_summits.bed hg38 motifs/${SAMPLE}/ -size 200 -mask -p ${THREADS}

# 8. 质控统计
echo "Total reads: $(samtools view -c ${SAMPLE}_final.bam)"
echo "Peaks: $(wc -l < peaks/${SAMPLE}_peaks.narrowPeak)"
echo "FRiP: $(echo "scale=4; $(bedtools intersect -a ${SAMPLE}_final.bam -b peaks/${SAMPLE}_peaks.narrowPeak -u | samtools view -c) / $(samtools view -c ${SAMPLE}_final.bam)" | bc)"