ATACseq 基本原理如下图所示。真核生物 DNA 与核小体结合形成染色质，结合的紧密程度是动态变化的。当 DNA 需要复制或转录时，结合变得松散让 DNA 暴露。暴露的 DNA 能被 Tn5 转座酶切割，与核小体紧密结合的 DNA 无法被切割。Tn5 转座酶切割 DNA 时插入测序接头，经过 PCR 扩增等步骤就完成了测序建库。
如下图所示，因为间隔的核小体数目不同，ATACseq 测序插入片段（fragments）分布会呈现多个峰的分布。在约 <100,200,400,600bp 处有峰，分别对应着 NER (nucleosome-free regions) 和 单、双、三核小体区域的 reads. 因此我认为常用的 2*150 的双端测序不适合 ATACseq. 因为最需要的 NER 区域 reads 是很短的，用 2*150 等于浪费许多的测序通量。
ATACseq 一般人种要求 50M 以上比对的 fragments (reads), 因为线粒体 DNA 没有核小体结合，测序出现大量线粒体数据是正常的，线粒体数据占比在 20-80% 都有可能。
ATACseq 分析流程如下图所示。质控、比对、peak calling 是必须的上游分析，Motif 等下游分析是个性化的。
准备工作
下载参考基因组、建立索引、移除不需要区域等。参考基因组在 GENCODE 下载。
一般用 bowtie2 比对，建立 bowtie2 索引。
bowtie2-build -f GRCh38.primary_assembly.genome.fa GRCh38
一般只分析常染色体，并移除 Blacklist 和线粒体。因此将基因组区域分为两个 bed 文件，一个是需要的区域一个是移除的区域。
# include bed
samtools faidx GRCh38.primary_assembly.genome.fa
awk -v FS="\t" -v OFS="\t" '{print $1 FS "0" FS ($2-1)}' GRCh38.primary_assembly.genome.fa.fai | \
grep -e "^chr[0-9]" > GRCh38.primary_assembly.genome.bed.temp
bedtools subtract -a GRCh38.primary_assembly.genome.bed.temp -b ${Blacklist} > GRCh38.primary_assembly.atac_included.bed
rm GRCh38.primary_assembly.genome.bed.temp
# exclude bed
# grep -v
awk -v FS="\t" -v OFS="\t" '{print $1 FS "0" FS ($2-1)}' GRCh38.primary_assembly.genome.fa.fai | \
grep -v -e "^chr[0-9]" > GRCh38.primary_assembly.genome.bed.temp
cat GRCh38.primary_assembly.genome.bed.temp ${Blacklist} | bedtools sort -i - | \
bedtools merge -i - > GRCh38.primary_assembly.atac_excluded.bed
rm GRCh38.primary_assembly.genome.bed.temp
fastq 质控
缠绕核小体 DNA 约 147bp 与相邻核小体连接的 DNA 约 20-90bp. 加上测序接头等约 135bp 长度会达到 200bp 左右，因此最后文库片段长度可能是 200-1000bp 左右，并且主要的部分在 600bp 一下，但 ATACseq 建库片段分布可能因为样本类型、细胞数量、处理过程等有关，也许文库片段分布有所差异。
原始 fastq 用 fastqc 生成质控报告，然后用 fastp 进行过滤。
fastqc -t 4 *fq.gz
fastp -i ${raw_dir}/${sample}_1.fq.gz -o ${clean_dir}/${sample}_R1.fq.gz \
-I ${raw_dir}/${sample}_2.fq.gz -O ${clean_dir}/${sample}_R2.fq.gz \
--compression 6 --json ${clean_dir}/${sample}_fastp.json --report_title ${sample}
用 multiqc 将质控报告汇总。
multiqc -o fastqc_report *fastqc.zip
multiqc -o fastp_report *fastp.json
bowtie2 比对
因为下一步 peak calling 用 genrich 支持多比对 reads 分析，所以 bowtie2 比对设置 -k 10 参数，-X 2000 允许插入片段最大 2000bp, 默认值 500. 注意 genrich 要求 -n 排序。
# bed 文件是移除了 blacklist 和线粒体的常染色体区域
bowtie2 --very-sensitive -k 10 -X 2000 --threads 6 -x ${GRCh38} \
-1 ${clean_dir}/${sample}_R1.fq.gz -2 ${clean_dir}/${sample}_R2.fq.gz \
-S ${bam_dir}/${sample}.sam
# -n
samtools view --threads 4 -F 524 -L ${atac_included.bed} -b ${bam_dir}/${sample}.sam \
| samtools sort -n --threads 4 -O BAM -o ${bam_dir}/${sample}_nSorted.bam
# -p
samtools sort --threads 4 -O BAM ${bam_dir}/${sample}_pSorted.bam ${bam_dir}/${sample}_nSorted.bam
gatk --java-options "-Xmx64G -Xms8G" MarkDuplicates --INPUT ${bam_dir}/${sample}_pSorted.bam \
--OUTPUT ${bam_dir}/${sample}_MD.bam --METRICS_FILE ${bam_dir}/${sample}_MD.txt
samtools index ${bam_dir}/${sample}_MD.bam
同时考虑有一些下游分析不适合用 -k 10 的比对，同步做一个没有 -k 参数设置的比对，并按照 -p 排序。
bowtie2 --very-sensitive -X 2000 --threads 6 -x ${GRCh38} \
-1 ${clean_dir}/${sample}_R1.fq.gz -2 ${clean_dir}/${sample}_R2.fq.gz \
-S ${bam_dir}/${sample}.nk.sam
samtools view --threads 4 -F 524 -q 30 -L ${atac_included.bed} -b ${bam_dir}/${sample}.nk.sam \
| samtools sort --threads 4 -O BAM -o ${bam_dir}/${sample}_Sorted.nk.bam
gatk --java-options "-Xmx64G -Xms8G" MarkDuplicates --INPUT ${bam_dir}/${sample}_Sorted.nk.bam \
--OUTPUT ${bam_dir}/${sample}_MD.nk.bam --METRICS_FILE ${bam_dir}/${sample}_md.metrics
samtools index ${bam_dir}/${sample}_MD.nk.bam
samtools view --threads 4 -F 1024 -b -o ${bam_dir}/${sample}_RD.nk.bam ${bam_dir}/${sample}_MD.nk.bam
samtools index ${bam_dir}/${sample}_RD.nk.bam
# 为了计算线粒体 Reads 占比，将原始 sam 进行排序和建立索引
samtools view --threads 4 -b ${bam_dir}/${sample}.nk.sam | samtools sort \
--threads 4 -O BAM -o ${bam_dir}/${sample}_Raw_Sorted.nk.bam
samtools index ${bam_dir}/${sample}_Raw_Sorted.nk.bam
samtools view 命令 -F 524 -q 30 移除不合格的比对。
$ samtools flags 524
0x20c 524 UNMAP,MUNMAP,QCFAIL
$ samtools flags 1024
0x400 1024 DUP
检查比对报告，可能 "aligned concordantly >1 times" 比例会比较高，原因一是 --very-sensitive 参数，让比对质量稍低于最佳比对的仍保留，第二是 ATACseq 数据线粒体含量高，线粒体有许多区域序列类似。
前面说了 ATACseq 插入片段应该在约 <100,200,400,600bp 大小有峰，用 gatk CollectInsertSizeMetrics 命令分析比对后插入片段分布。
gatk CollectInsertSizeMetrics -H ${bam_dir}/${sample}_InsertSize.pdf \
-I ${bam_dir}/${sample}_MD.bam -O ${bam_dir}/${sample}_InsertSize.txt
（可选步骤）移动 Reads
如下面示意图所示，Tn5 建库时会插入 9bp 的序列，所以比对到正链（+）的 reads 移动 +4 bp, 比对负链（-）的 reads 移动 -5 bp.
这个移动可以用 deeptools 的 alignmentSieve 命令完成。
alignmentSieve --numberOfProcessors 8 --ATACshift -b \
${bam_dir}/${sample}_nSorted.bam -o ${bam_dir}/${sample}_Shift.bam
一般的分析这点 reads 位移没影响，不需要进行这个步骤。
Genrich peak calling
Genrich peak calling 过程：
- 从比对结果计算实验样本和对照样本的每个位置覆盖度，并计算每一组 pileup 结果。
- 计算信号背景水平，如果有对照从对照样本计算，如果没有从实验样本计算。背景计算是总的 read/fragment/interval 长度除以基因组长度（从 SAM/BAM 表头计算），如果有
-e或-E等参数，会扣除相应长度。 - 每个位置计算 p 值和 q 值。
- 计算 p 值和 q 值显著的区域的 AUC.
- 合并相邻近区域并计算总 AUC 是否超过设置的阈值进行 call peak.
Genrich 的 -j 参数表示 ATAC 模式，默认调整 reads/fragments 5bp 位置。ATACseq 跟 ChIP-Seq 不同，一般没有 control input. 所以实验组和对照组分开用 Genrich 进行分析，得到各自的 peak 继续下游分析。
Genrich -t ${bam_dir}/${sample1}_nSorted.bam,${bam_dir}/${sample2}_nSorted.bam,${bam_dir}/${sample3}_nSorted.bam \
-o ${genrich_dir}/${group}.narrowPeak -f ${genrich_dir}/${group}_pq.bed \
-k ${genrich_dir}/${group}_pileup_p.bed -b ${genrich_dir}/${group}_reads.bed \
-E ${rm_bed} -m 30 -q 0.05 -a 500 -r -j
同一组的重复样品 -t 参数输入并且逗号分隔；-r 移除 PCR 重复；-m 过滤比对 MAPQ 值，Bowtie2 比对一般取 30 阈值；-b 将 reads/fragments 输出到 bed 文件；-f 输出每个区间每个样本及合并后的 p 值和显著性；-k 输出每个样本 pileups 和 p 值。
因为 -a 参数不好确定，可以先设置 -X 参数让 genrich 输出分析日志文件，但不进行 peak calling, 然后用不同的 -a 参数用 -P 参数模式进行 peak calling.
Genrich -t ${bam_dir}/${group}_1_nSorted.bam,${bam_dir}/${group}_2_nSorted.bam,${bam_dir}/${group}_3_nSorted.bam \
-f ${genrich_dir}/${group}_pq.bed -k ${genrich_dir}/${group}_pileup.bed \
-b ${genrich_dir}/${group}_reads.bed -E ${ex_bed} -m 30 -q 0.05 -r -j -X
# 用不同 -a 参数进行 peak calling
Genrich -f ${genrich_dir}/${group}_pq.bed -o ${genrich_dir}/${group}_a300.narrowPeak -a 300 -l 40 -q 0.05 -P
Genrich -f ${genrich_dir}/${group}_pq.bed -o ${genrich_dir}/${group}_a400.narrowPeak -a 400 -l 40 -q 0.05 -P
Genrich -f ${genrich_dir}/${group}_pq.bed -o ${genrich_dir}/${group}_a500.narrowPeak -a 500 -l 40 -q 0.05 -P
Genrich -f ${genrich_dir}/${group}_pq.bed -o ${genrich_dir}/${group}_a600.narrowPeak -a 600 -l 40 -q 0.05 -P
Genrich -f ${genrich_dir}/${group}_pq.bed -o ${genrich_dir}/${group}_a700.narrowPeak -a 700 -l 40 -q 0.05 -P
Genrich -f ${genrich_dir}/${group}_pq.bed -o ${genrich_dir}/${group}_a800.narrowPeak -a 800 -l 40 -q 0.05 -P
默认结果输出为 narrowPeak 格式。第五列 score 是平均 AUC (total AUC / bp) * 1000, 最大值 1000. 最后一列是 peak 信号最显著位置与 peak 起始距离，可视为 summit 位置。
chr1 9963 10500 peak_0 1000 . 2459.160400 11.532567 8.601989 93
chr1 11141 11518 peak_1 1000 . 1298.728394 10.535641 7.766460 207
chr1 13230 13959 peak_2 1000 . 1910.681885 10.124692 7.418398 290
FRiP 计算
Fraction of reads in peaks (FRiP) 指落入峰区域的 reads 占比，应该要高于 0.3 较好。这个指标非强制性，不同数据表现并不相同，计算也不用追求精确。
分别计算 BAM 文件总 fragments 数和处于 peak 区域的 fragments 数，再相除。
samtools view -c ${bam_dir}/${sample}_MD.nk.bam
samtools view -c -L ${genrich_dir}/${group}.narrowPeak ${bam_dir}/${sample}_MD.nk.bam
两命令得到的 fragments 数目相除便是 FRiP.
或者将峰文件转换为 SAF 文件，然后用 featureCounts 计算 counts/fragments. 其报告的 assigned reads 比例可视为 FRIP.
或者 DiffBind 分析时得到。
ChIPseeker 注释
ChIPseeker 对 peaks 进行基因组注释，这个注释是按照 peak 与基因位置关系给的，所以如果位置重合了，默认按照以下优先度分配：Promoter > 5UTR > 3UTR > Exon > Intron > Downstream > Intergenic.
ChIPseeker 读取数据后 annotatePeak 函数进行注释，注释结果用 as.data.frame 转换为数据框再保存到文件。参数 addFlankGeneInfo 和 flankDistance 报告所有在 peak 一定范围内的的基因。
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
library(org.Hs.eg.db)
library(ChIPseeker)
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
peaks <- readPeakFile(peak_path)
peak_anno <- annotatePeak(peaks, tssRegion = c(-3000, 3000), TxDb = txdb,
annoDb = "org.Hs.eg.db", level = "gene",
addFlankGeneInfo = TRUE, flankDistance = 1000)
DiffBind 差异分析
ATACseq 差异分析是比较困难的步骤，DiffBind 用得比较多，但它的结果也不一定那么可靠。DiffBind 原理是把数据当作 RNAseq 差异表达分析，只是把 RNAseq 分析的基因修改为 peak. 每一个 consensus peak 等于一个基因。明白了这个原理甚至可以自己手动去做，先从所有样本 peakset 取得 consensus peakset. 然后 featureCount 计算每个 peak reads 数，最后输入 DESeq2 进行差异分析。
DiffBind 输入是 csv 或 excel 格式的信息表，表里包含了分析的样品、数据路径、分组等等。
以下是表里可接受的列，只填自己需要填的列，缺少的列 DiffBind 会填充 NA.
- SampleID
- Tissue
- Factor
- Condition
- Treatment
- Replicate | Replicate number of sample
- bamReads | file path for bam file containing aligned reads for ChIP sample
- bamControl | ile path for bam file containing aligned reads for control sample
- Spikein
- ControlID
- Peaks | path for file containing peaks for sample
- PeakCaller
- raw | text file file; peak score is in fourth column
- bed
- narrow | default peak.format: narrowPeaks file
- macs | MACS .xls file
- swembl | SWEMBL .peaks file
- bayes | bayesPeak file
- peakset | peakset written out using pv.writepeakset
- fp4 | FindPeaks v4
- PeakFormat
- ScoreCol | column in peak files that contains peak scores
- LowerBetter | logical indicating that lower scores signify better peaks
- counts | file path for externally computed read counts
准备好信息表后 dba 函数读取信息创建 "dba" 对象。
library(DiffBind)
library(tidyverse)
info_path <- file.path(diffbind_dir, "kLEC.csv")
dba1 <- dba(sampleSheet = info_path)
第二步计算 counts 矩阵。DiffBind 首先从输入的 peaks 分析出 consensus peakset. 也可以自己通过 peaks 参数提供自己计算的 consensus peakset. 然后计算每个 peak 的 summit 位置，从 summit 向上下游延申相同长度（由 summits 参数决定）得到相同 peak 长度的 peakset 计算 counts 矩阵并下游分析。
默认 summits 是 200 最后每个 peak 长度是 401. ATAC-Seq peak 较小应设为 100 或更小。参数 bRemoveDuplicates 移除重复 reads. 但一些情况下软件计算的重复并不准确，比如双端测序，为了准确计算需要设置 bUseSummarizeOverlaps 参数为 TRUE.
dba2 <- dba.count(DBA = dba1, bRemoveDuplicates = TRUE, bUseSummarizeOverlaps = TRUE,
bParallel = TRUE, peaks = merged_peaks, summits = 100, filter = 1)
因为会进行过滤移除一些 peaks (如 counts 太少)，用 dba.peakset 函数获取最终 consensus peakset 和 binding affinity matrix. 后者由 dba.count 的 score 参数设定。这个函数也用于修改 DBA 的 peakset.
cons_peaks <- dba.peakset(DBA = dba2, bRetrieve = TRUE)
# 查看
> cons_peaks[1:2]
GRanges object with 2 ranges and 6 metadata columns:
seqnames ranges strand | CTR_1 CTR_2 CTR_3 TM_1 TM_2 TM_3
<Rle> <IRanges> <Rle> | <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
1 chr1 10019-10219 * | 305.641 320.698 261.559 692.019 855.200 699.159
2 chr1 11240-11440 * | 197.586 304.026 375.066 125.906 169.172 165.673
-------
seqinfo: 22 sequences from an unspecified genome; no seqlengths
DiffBind 有多种均一化方法，建议先查看帮助文档选择合适自己数据的方法。
均一化重要的一步是计算文库大小，DiffBind 支持用 full, RiP, background. full 指用完整文库大小，从 BAM 文件计算；RiP 指仅计算落入 consensus peaksets 的 reads; backgroud 指根据背景信号进行均一化，原理是 ATAC-Seq 得到的 peak 区域是比较小的，大约在 100-600bp. 如果将基因组分成大的 bin 比如 15000bp 大小，理论上这种大小 bin 的信号差异应该是技术误差而不是生物学差异。设置 background 参数后包含 consensus peakset 的染色体将被分 bin 计算背景信号差异。
dba3 <- dba.normalize(DBA = dba2, method = "DESeq2", library = "full")
根据实验设计的分组，分析差异的 peak.
dba_contrast <- dba.contrast(dba3, design = ~ Treatment,
contrast = c("Treatment", "CTR", "TM"))
dba_diff <- dba.analyze(dba_contrast, method = "DESeq2", bBlacklist = FALSE,
bGreylist = FALSE)
查看 contrast 信息和样品相关系数。
> dba.show(dba_diff, bContrasts = TRUE)
Factor Group Samples Group2 Samples2 DB.DESeq2
1 Treatment TM 3 CTR 3 34640
提取差异数据结果。从 dba.show 函数结果知道需要的 contrast 编号是 1.
diff_peaks1 <- dba.report(dba_diff, contrast = 1)
默认返回 Granges 结构数据，Conc 是所有样本平均 log reads 数，Fold 是 log2 差异倍数。
保存结果到 BED 和 CSV 格式。
# 转换到 tibble 包含了 P 值等全部信息
diff_peaks2 <- bind_cols(as_tibble(granges(diff_peaks1)), as_tibble(mcols(diff_peaks1)))
# 保存到文件
rtracklayer::export(diff_peaks1, con = diff_bed, format = "BED")
write_csv(diff_peaks2, diff_csv)
MEME-ChIP motif 分析
转录因子结合到开放的染色质区域调控基因表达，转录因子结合需要识别出特定的序列，这个特定序列称之为 motif, 结合的位点称为 TFBS (TF binding sites)。转录因子可以允许 TFBS 有一定的可塑性 （variations/flexible），所以 motif 序列不完全是固定死的。Motif 大部分不长，6-12 bp 左右。人类大约有 1600 个转录因子，其中超过 2/3 已鉴定了 motif.
motif 序列一般有两种模式，一是回文序列，比如 CACGTG 其反向互补序列也是 CACGTG. 二是两段保守序列被一段非保守序列分隔，这往往是因为结合的转录因子是二聚体，分别识别两段保守序列。
motif 的分析往往受假阳性困扰，这称为无效定理（Futility Theorem）。比如说往往在基因组序列中观察到大量的潜在转录因子结合位点，其中很少是真正起作用的，大部分预测的转录因子结合位点是无效的。所以 motif 分析后，要想办法进行人工筛选。
MEME-ChIP 主要执行以下步骤：
- 在输入序列的中间区域（默认 100bp）进行 motif 发现（MEME, STREME）
- CentriMo 分析哪些 motif 是富集在区域中心的
- Tomtom 分析他们与已知 motif 相似性
- 根据 motif 相似性对显著结果归类
- motif spacing analysis (SpaMo) (分析 motif 与结合在相邻位置 motif 的物理作用) (
-spamo-skip参数跳过这步分析) - 分析 motif 结合位置（FIMO），默认选取每一类别最显著的 motif 进行分析
MEME-ChIP 默认输入序列是约 500bp 长且中间 100bp 是 motif 区域。
ATAC-Seq 分析得到 peak 是长度不一的，因此取 summit 向两边各延申 250bp 得到 500bp 序列。
awk -v FS="\t" -v OFS="\t" '{midpos=$2+$10;print $1,midpos-250,midpos+250;}' \
${genrich_dir}/${group}.narrowPeak > ${meme_dir}/${group}.bed
bedtools getfasta -fo ${meme_dir}/${group}.fa -fi ${GRCh38} \
-bed ${meme_dir}/${group}.bed
考虑 motif 分析假阳性问题和转录因子一般结合到启动子区，只取注释到启动子区的 peaks 进行 motif 分析，是值得尝试的。上面命令取所有 peaks.
MEME-ChIP 将各步骤封装了，因此运行命令很简单：
meme-chip -meme-maxw 30 -meme-p 6 -oc ${meme_dir}/${group} \
-db ${meme_db} ${meme_dir}/${group}.fa
motif 数据库在 https://meme-suite.org/meme/db/motifs 下载。
MEME-ChIP 输出结果在 meme.html 查看；结果的汇总在 summary.tsv 文件；combined.meme 文件包含所有被鉴定出的 Motif.
motif E value 表示该 motif 多大概率是统计错误出现，而不是真的 motif. E 值越小说明发现的 motif 越大概率为真。
MEME-ChIP 的结果如果发现不方便批量提取结果，也可以自己单独运行它的子分析。比方说它 FIMO 只分析了部分 motif 而你需要分析全部的，那么可以单独进行 FIMO 分析。
fimo --parse-genomic-coord -oc ${fimo_dir} --bgfile ${meme_dir}/background --motif SP1_HUMAN.H11MO.1.A ${meme_db} ${meme_dir}/${group}.fa
可视化
可视化是个非常宽泛的概念，同时许多工具会自带一些可视化方法，这里介绍 deptools 和 EnrichedHeatmap 生成信号热图。
deeptools
deeptools 提供了 bam, bigwig 处理、QC、画图等许多工具。本教程介绍 plotHeatmap 画 ATACseq 信号热图。
第一步 bamCoverage 命令将 Bam 文件转换成 bigwig 文件。
bamCoverage --numberOfProcessors 8 --effectiveGenomeSize 2747877777 --normalizeUsing RPGC \
--outFileFormat bigwig -b ${bam_dir}/${sample}_MD.bam -o ${bw_dir}/${sample}.bw
参数 --effectiveGenomeSize 的设置值在 Effective Genome Size — deepTools 3.4.3 documentation 查看。
RPGC 计算方法：
Scale factor 计算方法：
(total number of mapped reads * fragment length) / effective genome size
第二步 computeMatrix 命令生成画图矩阵。但是先用 R 语言提取包含 peak 的启动子区域。注意输出的 BED 要包含 strand 列，computeMatrix 命令会正确处理 strand 信息。
library(rtracklayer)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
promoter <- promoters(genes(txdb), upstream = 3000, downstream = 3000)
peak <- import(peak_path, format = "narrowPeak")
overlap_index <- findOverlaps(promoter, peak)
keep_promoter <- unique(promoter[queryHits(overlap_index)])
export.bed(keep_promoter, bed_path)
从 bigwig 文件生成用于热图的矩阵，多个 bigwig 文件用空格分隔。
computeMatrix scale-regions --regionsFileName ${bw_dir}/${group}_promoter.bed --regionBodyLength 6000 \
--outFileName ${bw_dir}/${group}_matrix.gz --binSize 60 --numberOfProcessors 4 \
--scoreFileName ${bw_dir}/${group}_1.bw ${bw_dir}/${group}_2.bw ${bw_dir}/${group}_3.bw
最后 plotHeatmap 命令画图。
plotHeatmap --matrixFile ${bw_dir}/${group}_matrix.gz --outFileName ${bw_dir}/${group}_heatmap.pdf \
--zMin 0 --zMax 8 --xAxisLabel Promoter --startLabel 3Kb --endLabel 3Kb
EnrichedHeatmap
EnrichedHeatmap 基于 ComplexHeatmap 的信号富集热图 R 包。其原理是将 peakset 信号转换为矩阵并热图可视化。
library(rtracklayer)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
library(ChIPseeker)
library(EnrichedHeatmap)
library(circlize)
取得基因组启动子区域，注意用 genes(txdb) 而不是 txdb 限制为基因的启动子，否则条目将非常多。
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
promoter <- promoters(genes(txdb), upstream = 3000, downstream = 3000)
> promoter[1:3]
GRanges object with 3 ranges and 1 metadata column:
seqnames ranges strand | gene_id
<Rle> <IRanges> <Rle> | <character>
1 chr19 58359752-58365751 - | 1
10 chr8 18388282-18394281 + | 10
100 chr20 44649234-44655233 - | 100
-------
seqinfo: 595 sequences (1 circular) from hg38 genome
用 ChIPseeker 包的 readPeakFile 函数读取 narrowPeak 文件。或者如果是 rtracklayer 包的 import 函数，它支持读取更多的格式，比如 bigwig,bed 等。
peak <- readPeakFile(peakPath)
将启动子区 peakset 数据转换为画图矩阵。注意将没信号的启动子区数据移除。
mat <- normalizeToMatrix(peak, promoter, extend = 0, value_column = "V5", k = 100, mean_mode = "w0")
# 移除没信号启动子区
mat <- mat[rowSums(mat) > 0,]
读取后 narrowPeak 数据 "V5" 列为峰信号强度值，所以 value_column 参数设置 "V5".
参数 k 决定 bins/windows 数目（矩阵的列数）。参数 mean_mode 决定计算平均信号强度方法，可选以下方法：
Following illustrates different settings for mean_mode (note there is one signal region overlapping with other signals):
40 50 20 values in signal regions
++++++ +++ +++++ signal regions
30 values in signal region
++++++ signal region
================= a window (17bp), there are 4bp not overlapping to any signal regions.
4 6 3 3 overlap
absolute: (40 + 30 + 50 + 20)/4
weighted: (40*4 + 30*6 + 50*3 + 20*3)/(4 + 6 + 3 + 3)
w0: (40*4 + 30*6 + 50*3 + 20*3)/(4 + 6 + 3 + 3 + 4)
coverage: (40*4 + 30*6 + 50*3 + 20*3)/17
最后 EnrichedHeatmap 函数画图。
hm_color <- colorRamp2(breaks = c(0, 1000), colors = c("#fff5f0", "#cb1c1d"))
top_anno <- HeatmapAnnotation(enriched = anno_enriched(gp=gpar(col="black", lwd = 1.2)))
enrich_hm <- EnrichedHeatmap(mat = mat, col = hm_color, axis_name = c(-3000, 3000),
top_annotation = top_anno, show_heatmap_legend = FALSE, use_raster = FALSE)
参考资料
ATAC-Seq data analysis
ATAC-seq Guidelines - Harvard FAS Informatics
ATAC-seq 150bp reads
ChIP-seq Peak Annotation and Functional Analysis | Introduction to ChIP-Seq using high-performance computing
Library QC for ATAC-Seq and CUT&Tag AKA “Does My Library Look Okay?”
Make Enriched Heatmaps
Using EnrichedHeatmap for visualization of NGS experiments
An Introduction to the GenomicRanges Package
MEME-ChIP - MEME Suite
Ma, Wenxiu, William S. Noble, and Timothy L. Bailey. "Motif-based analysis of large nucleotide data sets using MEME-ChIP." Nature protocols 9.6 (2014): 1428-1450.
Yan, Feng, et al. "From reads to insight: a hitchhiker’s guide to ATAC-seq data analysis." Genome biology 21.1 (2020): 22.
Brouilette, Scott, et al. "A simple and novel method for RNA‐seq library preparation of single cell cDNA analysis by hyperactive Tn5 transposase." Developmental Dynamics 241.10 (2012): 1584-1590.
Gontarz, Paul, et al. "Comparison of differential accessibility analysis strategies for ATAC-seq data." Scientific reports 10.1 (2020): 1-13.
Yin, Senlin, et al. "Transcriptomic and open chromatin atlas of high-resolution anatomical regions in the rhesus macaque brain." Nature communications 11.1 (2020): 1-13.
Buenrostro, Jason D., et al. "ATAC‐seq: a method for assaying chromatin accessibility genome‐wide." Current protocols in molecular biology 109.1 (2015): 21-29.
Yan, Feng, et al. "From reads to insight: a hitchhiker’s guide to ATAC-seq data analysis." Genome biology 21.1 (2020): 22.
Buenrostro, Jason D., et al. "Transposition of native chromatin for multimodal regulatory analysis and personal epigenomics." Nature methods 10.12 (2013): 1213.
