染色质可及性ATAC-seq全流程¶

# FastQC质控
fastqc -t 8 sample_R1.fastq.gz sample_R2.fastq.gz -o qc_raw/

# Trim Galore（自动检测Nextera接头）
trim_galore --paired --nextera --quality 20 --length 20 \
  --cores 4 --fastqc \
  sample_R1.fastq.gz sample_R2.fastq.gz \
  -o trimmed/

# 或手动指定Nextera接头
trim_galore --paired \
  -a CTGTCTCTTATACACATCT -a2 CTGTCTCTTATACACATCT \
  --quality 20 --length 20 \
  --cores 4 \
  sample_R1.fastq.gz sample_R2.fastq.gz \
  -o trimmed/

# Bowtie2索引（如未建立）
bowtie2-build --threads 16 genome.fa genome_bt2_idx

# Bowtie2比对
bowtie2 -x genome_bt2_idx \
  -1 trimmed/sample_R1_val_1.fq.gz \
  -2 trimmed/sample_R2_val_2.fq.gz \
  --very-sensitive \
  -X 2000 \
  --no-mixed --no-discordant \
  --threads 16 \
  2> align_stats.txt \
  | samtools sort -@ 8 -o aligned/sample.bam

samtools index aligned/sample.bam

# 查看比对统计
cat align_stats.txt

# 标记并去除PCR重复
picard MarkDuplicates \
  I=aligned/sample.bam \
  O=aligned/sample_dedup.bam \
  M=aligned/sample_dup_metrics.txt \
  REMOVE_DUPLICATES=true

# 多步过滤
samtools view -@ 8 -b -f 2 -F 1804 -q 30 \
  aligned/sample_dedup.bam \
  > aligned/sample_filtered.bam

# Flag解释:
# -f 2: 保留properly paired
# -F 1804: 去除unmapped(4), mate unmapped(8), not primary(256), 
#           fails QC(512), duplicate(1024)
# -q 30: MAPQ >= 30

samtools index aligned/sample_filtered.bam

# 统计过滤前后reads数
echo "原始比对reads: $(samtools view -c aligned/sample.bam)"
echo "过滤后reads: $(samtools view -c aligned/sample_filtered.bam)"

# 去除线粒体reads
samtools view -@ 8 -b aligned/sample_filtered.bam \
  $(samtools view -H aligned/sample_filtered.bam | grep "^@SQ" | \
    awk '{print $2}' | sed 's/SN://' | grep -v "chrM" | tr '\n' ' ') \
  > aligned/sample_noMT.bam

# 或者用samtools idxstats计算线粒体比例后直接过滤
echo "线粒体reads比例:"
samtools idxstats aligned/sample_filtered.bam | \
  awk '{if($1=="chrM") mt=$3; total+=$3} END{print mt/total*100 "%"}'

# 去除黑名单区域
bedtools intersect -v -abam aligned/sample_noMT.bam \
  -b encode_blacklist_hg38.bed \
  > aligned/sample_final.bam

samtools index aligned/sample_final.bam

echo "最终分析reads: $(samtools view -c aligned/sample_final.bam)"

# 使用alignmentSieve（deepTools）进行Tn5偏移校正
alignmentSieve \
  --bam aligned/sample_final.bam \
  --outFile aligned/sample_shifted.bam \
  --ATACshift \
  --numberOfProcessors 16

samtools sort -@ 8 -o aligned/sample_shifted_sorted.bam aligned/sample_shifted.bam
samtools index aligned/sample_shifted_sorted.bam

# 或使用自定义脚本（基于pysam）
python3 << 'EOF'
import pysam

inbam = pysam.AlignmentFile("aligned/sample_final.bam", "rb")
outbam = pysam.AlignmentFile("aligned/sample_shifted.bam", "wb", template=inbam)

for read in inbam:
    if read.is_unmapped or read.is_secondary or read.is_supplementary:
        continue
    if read.is_reverse:
        read.reference_start = read.reference_start - 5
    else:
        read.reference_start = read.reference_start + 4
    outbam.write(read)

inbam.close()
outbam.close()
EOF

# 使用Picard计算插入片段大小
picard CollectInsertSizeMetrics \
  I=aligned/sample_final.bam \
  O=qc/insert_size_metrics.txt \
  H=qc/insert_size_histogram.pdf

# 使用ATACseqQC（R包，更全面的QC）

library(ATACseqQC)
library(ggplot2)

# 片段大小分布
fragSize <- fragSizeDist("aligned/sample_final.bam", "Sample1")

# 或直接从BAM提取
library(Rsamtools)
param <- ScanBamParam(what="isize")
bam <- scanBam("aligned/sample_final.bam", param=param)
frag_sizes <- abs(bam[[1]]$isize)
frag_sizes <- frag_sizes[frag_sizes > 0 & frag_sizes < 1000]

ggplot(data.frame(size=frag_sizes), aes(x=size)) +
  geom_histogram(bins=500, fill="steelblue", color="white", linewidth=0.1) +
  geom_vline(xintercept=c(150, 200, 400), linetype="dashed", color="red") +
  labs(x="Fragment size (bp)", y="Count",
       title="ATAC-seq Fragment Size Distribution") +
  xlim(0, 800) +
  theme_minimal(base_size=14)
ggsave("qc/fragment_size_distribution.pdf", width=8, height=5)

# TSS富集图
library(ChIPseeker)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene

# 读取peaks用于TSS注释（先做peak calling再回来看）
peaks <- readPeakFile("peaks/sample_peaks.narrowPeak")
peakAnno <- annotatePeak(peaks, TxDb=txdb, tssRegion=c(-3000, 3000))
plotAnnoPie(peakAnno)

# MACS2 peak calling（推荐：BAMPE模式）
macs2 callpeak \
  -t aligned/sample_final.bam \
  -f BAMPE \
  -g hs \
  -n sample \
  --outdir peaks/ \
  --nomodel \
  --keep-dup all \
  -B --SPMR \
  -q 0.01

# 输出文件：
# peaks/sample_peaks.narrowPeak  - peak坐标和信号
# peaks/sample_summits.bed       - peak顶点
# peaks/sample_treat_pileup.bdg  - 信号track

# HMMRATAC（替代方案）
HMMRATAC \
  -b aligned/sample_final.bam \
  -i aligned/sample_final.bam.bai \
  -g genome.genome \
  -o peaks/sample_hmmratac

# 统计
echo "Peak数: $(wc -l < peaks/sample_peaks.narrowPeak)"
echo "Peak中位数宽度: $(awk '{print $3-$2}' peaks/sample_peaks.narrowPeak | sort -n | awk 'NR==int(NR/2)')"

# 生成bigWig用于可视化
bamCoverage \
  -b aligned/sample_final.bam \
  -o tracks/sample.bw \
  --normalizeUsing RPGC \
  --effectiveGenomeSize 2913022398 \
  --binSize 10 \
  -p 16

library(DiffBind)

# 创建DiffBind样本表
samples <- data.frame(
  SampleID = c("Control_1","Control_2","Control_3","Treatment_1","Treatment_2","Treatment_3"),
  Condition = c(rep("Control",3), rep("Treatment",3)),
  bamReads = paste0("aligned/", c("ctrl1","ctrl2","ctrl3","treat1","treat2","treat3"), "_final.bam"),
  Peaks = paste0("peaks/", c("ctrl1","ctrl2","ctrl3","treat1","treat2","treat3"), "_peaks.narrowPeak"),
  PeakCaller = rep("macs", 6)
)

# DiffBind分析
dba_obj <- dba(sampleSheet=samples)
dba_obj <- dba.count(dba_obj, bUseSummarizeOverlaps=TRUE)
dba_obj <- dba.normalize(dba_obj)
dba_obj <- dba.contrast(dba_obj, contrast=c("Condition","Treatment","Control"))
dba_obj <- dba.analyze(dba_obj, method=DBA_DESEQ2)

# 获取差异peaks
diff_peaks <- dba.report(dba_obj, th=0.05, fold=1)
cat("差异peaks数:", length(diff_peaks), "\n")
cat("  开放增加:", sum(diff_peaks$Fold > 0), "\n")
cat("  开放减少:", sum(diff_peaks$Fold < 0), "\n")

# MA图
dba.plotMA(dba_obj)
# 火山图
dba.plotVolcano(dba_obj)

# 保存差异peaks
export(diff_peaks, "results/differential_peaks.bed")

# TOBIAS流程
# 1. Tn5偏好校正
TOBIAS ATACorrect \
  --bam aligned/sample_final.bam \
  --genome genome.fa \
  --peaks peaks/sample_peaks.narrowPeak \
  --outdir tobias_output/ \
  --cores 16

# 2. 足迹评分
TOBIAS ScoreBigwig \
  --signal tobias_output/sample_final_corrected.bw \
  --regions peaks/sample_peaks.narrowPeak \
  --output tobias_output/sample_footprints.bw \
  --cores 16

# 3. TF结合检测（差异比较）
TOBIAS BINDetect \
  --motifs JASPAR2024_CORE_vertebrates.meme \
  --signals treatment_footprints.bw control_footprints.bw \
  --genome genome.fa \
  --peaks consensus_peaks.bed \
  --outdir tobias_output/bindetect \
  --cores 16 \
  --cond_names Treatment Control

# HOMER motif分析
# 安装HOMER后需要下载motif数据库
# perl configureHomer.pl -install hg38

findMotifsGenome.pl \
  peaks/sample_summits.bed \
  hg38 \
  motif_results/ \
  -size 200 \
  -mask \
  -p 16

# 对差异peaks分别分析
# 开放增加的peaks
findMotifsGenome.pl \
  results/peaks_gained.bed \
  hg38 \
  motif_gained/ \
  -size 200 -mask -p 16

# 开放减少的peaks
findMotifsGenome.pl \
  results/peaks_lost.bed \
  hg38 \
  motif_lost/ \
  -size 200 -mask -p 16

# MEME-ChIP（替代方案）
# 先提取peak序列
bedtools getfasta -fi genome.fa -bed peaks/sample_summits_200bp.bed > peak_seqs.fa
meme-chip -meme-maxw 20 -meme-nmotifs 10 peak_seqs.fa -o meme_output/

#!/bin/bash
# ====================================================
# ATAC-seq完整分析流程
# ====================================================

set -euo pipefail

SAMPLE="sample1"
THREADS=16
GENOME="ref/hg38.fa"
BT2_INDEX="ref/hg38_bt2"
BLACKLIST="ref/hg38_blacklist_v2.bed"
GENOME_SIZE="2913022398"  # hg38 effective genome size

R1="raw/${SAMPLE}_R1.fastq.gz"
R2="raw/${SAMPLE}_R2.fastq.gz"
OUTDIR="results/${SAMPLE}"

mkdir -p ${OUTDIR}/{qc,trimmed,aligned,peaks,tracks,motifs}

# ==== Step 1: 质控 ====
echo "[Step 1] 质控与接头去除..."
trim_galore --paired --nextera --quality 20 --length 20 \
  --cores 4 --fastqc \
  ${R1} ${R2} \
  -o ${OUTDIR}/trimmed/

# ==== Step 2: 比对 ====
echo "[Step 2] Bowtie2比对..."
bowtie2 -x ${BT2_INDEX} \
  -1 ${OUTDIR}/trimmed/${SAMPLE}_R1_val_1.fq.gz \
  -2 ${OUTDIR}/trimmed/${SAMPLE}_R2_val_2.fq.gz \
  --very-sensitive -X 2000 \
  --no-mixed --no-discordant \
  --threads ${THREADS} \
  2> ${OUTDIR}/qc/bowtie2_stats.txt \
  | samtools sort -@ 8 -o ${OUTDIR}/aligned/raw.bam

samtools index ${OUTDIR}/aligned/raw.bam

# ==== Step 3: 去重和过滤 ====
echo "[Step 3] 去重和过滤..."
picard MarkDuplicates \
  I=${OUTDIR}/aligned/raw.bam \
  O=${OUTDIR}/aligned/dedup.bam \
  M=${OUTDIR}/qc/dup_metrics.txt \
  REMOVE_DUPLICATES=true

samtools view -@ 8 -b -f 2 -F 1804 -q 30 \
  ${OUTDIR}/aligned/dedup.bam \
  > ${OUTDIR}/aligned/filtered.bam

# ==== Step 4: 去MT和黑名单 ====
echo "[Step 4] 去除线粒体和黑名单..."
# 计算MT比例
MT_PCT=$(samtools idxstats ${OUTDIR}/aligned/filtered.bam | \
  awk '{if($1=="chrM") mt=$3; total+=$3} END{printf "%.1f", mt/total*100}')
echo "线粒体reads比例: ${MT_PCT}%"

# 去除chrM
samtools view -@ 8 -b ${OUTDIR}/aligned/filtered.bam \
  $(seq 1 22 | sed 's/^/chr/' | tr '\n' ' ') chrX chrY \
  > ${OUTDIR}/aligned/noMT.bam

# 去除黑名单
bedtools intersect -v -abam ${OUTDIR}/aligned/noMT.bam \
  -b ${BLACKLIST} \
  > ${OUTDIR}/aligned/final.bam

samtools index ${OUTDIR}/aligned/final.bam
FINAL_READS=$(samtools view -c ${OUTDIR}/aligned/final.bam)
echo "最终有效reads: ${FINAL_READS}"

# ==== Step 5: Tn5偏移校正 ====
echo "[Step 5] Tn5偏移校正..."
alignmentSieve \
  --bam ${OUTDIR}/aligned/final.bam \
  --outFile ${OUTDIR}/aligned/shifted.bam \
  --ATACshift \
  --numberOfProcessors ${THREADS}

samtools sort -@ 8 -o ${OUTDIR}/aligned/shifted_sorted.bam ${OUTDIR}/aligned/shifted.bam
samtools index ${OUTDIR}/aligned/shifted_sorted.bam

# ==== Step 6: 片段大小分布 ====
echo "[Step 6] 片段大小分布..."
picard CollectInsertSizeMetrics \
  I=${OUTDIR}/aligned/final.bam \
  O=${OUTDIR}/qc/insert_metrics.txt \
  H=${OUTDIR}/qc/insert_histogram.pdf

# ==== Step 7: Peak calling ====
echo "[Step 7] MACS2 Peak calling..."
macs2 callpeak \
  -t ${OUTDIR}/aligned/final.bam \
  -f BAMPE -g hs \
  -n ${SAMPLE} \
  --outdir ${OUTDIR}/peaks/ \
  --nomodel --keep-dup all \
  -B --SPMR -q 0.01

echo "Peak数: $(wc -l < ${OUTDIR}/peaks/${SAMPLE}_peaks.narrowPeak)"

# ==== Step 8: BigWig生成 ====
echo "[Step 8] 生成信号track..."
bamCoverage \
  -b ${OUTDIR}/aligned/final.bam \
  -o ${OUTDIR}/tracks/${SAMPLE}.bw \
  --normalizeUsing RPGC \
  --effectiveGenomeSize ${GENOME_SIZE} \
  --binSize 10 -p ${THREADS}

# ==== Step 9: Motif分析 ====
echo "[Step 9] Motif富集分析..."
findMotifsGenome.pl \
  ${OUTDIR}/peaks/${SAMPLE}_summits.bed \
  hg38 ${OUTDIR}/motifs/ \
  -size 200 -mask -p ${THREADS}

# ==== QC汇总 ====
echo ""
echo "========== QC Summary =========="
echo "样本: ${SAMPLE}"
echo "比对率: $(grep 'overall' ${OUTDIR}/qc/bowtie2_stats.txt)"
echo "重复率: $(grep 'PERCENT' ${OUTDIR}/qc/dup_metrics.txt | cut -f9)"
echo "线粒体: ${MT_PCT}%"
echo "有效reads: ${FINAL_READS}"
echo "Peak数: $(wc -l < ${OUTDIR}/peaks/${SAMPLE}_peaks.narrowPeak)"
echo "================================"