ATACseq 基本原理如下圖所示徙瓶。真核生物 DNA 與核小體結(jié)合形成染色質(zhì)，結(jié)合的緊密程度是動態(tài)變化的。當 DNA 需要復(fù)制或轉(zhuǎn)錄時笔咽，結(jié)合變得松散讓 DNA 暴露。暴露的 DNA 能被 Tn5 轉(zhuǎn)座酶切割霹期，與核小體緊密結(jié)合的 DNA 無法被切割叶组。Tn5 轉(zhuǎn)座酶切割 DNA 時插入測序接頭，經(jīng)過 PCR 擴增等步驟就完成了測序建庫历造。
如下圖所示甩十，因為間隔的核小體數(shù)目不同，ATACseq 測序插入片段（fragments）分布會呈現(xiàn)多個峰的分布吭产。在約 <100,200,400,600bp 處有峰侣监，分別對應(yīng)著 NER (nucleosome-free regions) 和 單、雙臣淤、三核小體區(qū)域的 reads. 因此我認為常用的 2*150 的雙端測序不適合 ATACseq. 因為最需要的 NER 區(qū)域 reads 是很短的橄霉，用 2*150 等于浪費許多的測序通量。
ATACseq 一般人種要求 50M 以上比對的 fragments (reads), 因為線粒體 DNA 沒有核小體結(jié)合荒典，測序出現(xiàn)大量線粒體數(shù)據(jù)是正常的酪劫，線粒體數(shù)據(jù)占比在 20-80% 都有可能吞鸭。
ATACseq 分析流程如下圖所示。質(zhì)控覆糟、比對刻剥、peak calling 是必須的上游分析，Motif 等下游分析是個性化的滩字。
準備工作
下載參考基因組造虏、建立索引、移除不需要區(qū)域等麦箍。參考基因組在 GENCODE 下載漓藕。
一般用 bowtie2 比對，建立 bowtie2 索引挟裂。
bowtie2-build -f GRCh38.primary_assembly.genome.fa GRCh38
一般只分析常染色體享钞，并移除 Blacklist 和線粒體。因此將基因組區(qū)域分為兩個 bed 文件诀蓉，一個是需要的區(qū)域一個是移除的區(qū)域栗竖。
# 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 質(zhì)控
纏繞核小體 DNA 約 147bp 與相鄰核小體連接的 DNA 約 20-90bp. 加上測序接頭等約 135bp 長度會達到 200bp 左右，因此最后文庫片段長度可能是 200-1000bp 左右渠啤，并且主要的部分在 600bp 一下狐肢，但 ATACseq 建庫片段分布可能因為樣本類型、細胞數(shù)量沥曹、處理過程等有關(guān)份名，也許文庫片段分布有所差異。
原始 fastq 用 fastqc 生成質(zhì)控報告妓美，然后用 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 將質(zhì)控報告匯總。
multiqc -o fastqc_report *fastqc.zip
multiqc -o fastp_report *fastp.json
bowtie2 比對
因為下一步 peak calling 用 genrich 支持多比對 reads 分析部脚，所以 bowtie2 比對設(shè)置 -k 10 參數(shù)想邦，-X 2000 允許插入片段最大 2000bp, 默認值 500. 注意 genrich 要求 -n 排序裤纹。
# bed 文件是移除了 blacklist 和線粒體的常染色體區(qū)域
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 參數(shù)設(shè)置的比對，并按照 -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 參數(shù)，讓比對質(zhì)量稍低于最佳比對的仍保留奸汇，第二是 ATACseq 數(shù)據(jù)線粒體含量高施符，線粒體有許多區(qū)域序列類似往声。
前面說了 ATACseq 插入片段應(yīng)該在約 <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 過程：
- 從比對結(jié)果計算實驗樣本和對照樣本的每個位置覆蓋度陆盘，并計算每一組 pileup 結(jié)果普筹。
- 計算信號背景水平，如果有對照從對照樣本計算隘马，如果沒有從實驗樣本計算太防。背景計算是總的 read/fragment/interval 長度除以基因組長度（從 SAM/BAM 表頭計算），如果有
-e或-E等參數(shù)酸员，會扣除相應(yīng)長度杏头。 - 每個位置計算 p 值和 q 值。
- 計算 p 值和 q 值顯著的區(qū)域的 AUC.
- 合并相鄰近區(qū)域并計算總 AUC 是否超過設(shè)置的閾值進行 call peak.
Genrich 的 -j 參數(shù)表示 ATAC 模式沸呐，默認調(diào)整 reads/fragments 5bp 位置醇王。ATACseq 跟 ChIP-Seq 不同，一般沒有 control input. 所以實驗組和對照組分開用 Genrich 進行分析崭添，得到各自的 peak 繼續(xù)下游分析寓娩。
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
同一組的重復(fù)樣品 -t 參數(shù)輸入并且逗號分隔；-r 移除 PCR 重復(fù)呼渣；-m 過濾比對 MAPQ 值棘伴，Bowtie2 比對一般取 30 閾值；-b 將 reads/fragments 輸出到 bed 文件屁置；-f 輸出每個區(qū)間每個樣本及合并后的 p 值和顯著性焊夸；-k 輸出每個樣本 pileups 和 p 值。
因為 -a 參數(shù)不好確定蓝角，可以先設(shè)置 -X 參數(shù)讓 genrich 輸出分析日志文件阱穗，但不進行 peak calling, 然后用不同的 -a 參數(shù)用 -P 參數(shù)模式進行 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 參數(shù)進行 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
默認結(jié)果輸出為 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) 指落入峰區(qū)域的 reads 占比，應(yīng)該要高于 0.3 較好患朱。這個指標非強制性鲁僚，不同數(shù)據(jù)表現(xiàn)并不相同，計算也不用追求精確。
分別計算 BAM 文件總 fragments 數(shù)和處于 peak 區(qū)域的 fragments 數(shù)冰沙，再相除侨艾。
samtools view -c ${bam_dir}/${sample}_MD.nk.bam
samtools view -c -L ${genrich_dir}/${group}.narrowPeak ${bam_dir}/${sample}_MD.nk.bam
兩命令得到的 fragments 數(shù)目相除便是 FRiP.
或者將峰文件轉(zhuǎn)換為 SAF 文件，然后用 featureCounts 計算 counts/fragments. 其報告的 assigned reads 比例可視為 FRIP.
或者 DiffBind 分析時得到拓挥。
ChIPseeker 注釋
ChIPseeker 對 peaks 進行基因組注釋蒋畜，這個注釋是按照 peak 與基因位置關(guān)系給的，所以如果位置重合了撞叽，默認按照以下優(yōu)先度分配：Promoter > 5UTR > 3UTR > Exon > Intron > Downstream > Intergenic.
ChIPseeker 讀取數(shù)據(jù)后 annotatePeak 函數(shù)進行注釋姻成，注釋結(jié)果用 as.data.frame 轉(zhuǎn)換為數(shù)據(jù)框再保存到文件。參數(shù) addFlankGeneInfo 和 flankDistance 報告所有在 peak 一定范圍內(nèi)的的基因愿棋。
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 用得比較多，但它的結(jié)果也不一定那么可靠糠雨。DiffBind 原理是把數(shù)據(jù)當作 RNAseq 差異表達分析才睹，只是把 RNAseq 分析的基因修改為 peak. 每一個 consensus peak 等于一個基因。明白了這個原理甚至可以自己手動去做甘邀，先從所有樣本 peakset 取得 consensus peakset. 然后 featureCount 計算每個 peak reads 數(shù)琅攘，最后輸入 DESeq2 進行差異分析。
DiffBind 輸入是 csv 或 excel 格式的信息表松邪，表里包含了分析的樣品坞琴、數(shù)據(jù)路徑、分組等等逗抑。
以下是表里可接受的列剧辐，只填自己需要填的列，缺少的列 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 函數(shù)讀取信息創(chuàng)建 "dba" 對象邮府。
library(DiffBind)
library(tidyverse)
info_path <- file.path(diffbind_dir, "kLEC.csv")
dba1 <- dba(sampleSheet = info_path)
第二步計算 counts 矩陣荧关。DiffBind 首先從輸入的 peaks 分析出 consensus peakset. 也可以自己通過 peaks 參數(shù)提供自己計算的 consensus peakset. 然后計算每個 peak 的 summit 位置，從 summit 向上下游延申相同長度（由 summits 參數(shù)決定）得到相同 peak 長度的 peakset 計算 counts 矩陣并下游分析褂傀。
默認 summits 是 200 最后每個 peak 長度是 401. ATAC-Seq peak 較小應(yīng)設(shè)為 100 或更小忍啤。參數(shù) bRemoveDuplicates 移除重復(fù) reads. 但一些情況下軟件計算的重復(fù)并不準確，比如雙端測序仙辟，為了準確計算需要設(shè)置 bUseSummarizeOverlaps 參數(shù)為 TRUE.
dba2 <- dba.count(DBA = dba1, bRemoveDuplicates = TRUE, bUseSummarizeOverlaps = TRUE,
bParallel = TRUE, peaks = merged_peaks, summits = 100, filter = 1)
因為會進行過濾移除一些 peaks (如 counts 太少)同波，用 dba.peakset 函數(shù)獲取最終 consensus peakset 和 binding affinity matrix. 后者由 dba.count 的 score 參數(shù)設(shè)定。這個函數(shù)也用于修改 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 有多種均一化方法欺嗤，建議先查看幫助文檔選擇合適自己數(shù)據(jù)的方法参萄。
均一化重要的一步是計算文庫大小，DiffBind 支持用 full, RiP, background. full 指用完整文庫大小煎饼，從 BAM 文件計算；RiP 指僅計算落入 consensus peaksets 的 reads; backgroud 指根據(jù)背景信號進行均一化校赤，原理是 ATAC-Seq 得到的 peak 區(qū)域是比較小的吆玖，大約在 100-600bp. 如果將基因組分成大的 bin 比如 15000bp 大小筒溃，理論上這種大小 bin 的信號差異應(yīng)該是技術(shù)誤差而不是生物學差異。設(shè)置 background 參數(shù)后包含 consensus peakset 的染色體將被分 bin 計算背景信號差異沾乘。
dba3 <- dba.normalize(DBA = dba2, method = "DESeq2", library = "full")
根據(jù)實驗設(shè)計的分組怜奖，分析差異的 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 信息和樣品相關(guān)系數(shù)。
> dba.show(dba_diff, bContrasts = TRUE)
Factor Group Samples Group2 Samples2 DB.DESeq2
1 Treatment TM 3 CTR 3 34640
提取差異數(shù)據(jù)結(jié)果翅阵。從 dba.show 函數(shù)結(jié)果知道需要的 contrast 編號是 1.
diff_peaks1 <- dba.report(dba_diff, contrast = 1)
默認返回 Granges 結(jié)構(gòu)數(shù)據(jù)歪玲，Conc 是所有樣本平均 log reads 數(shù)，F(xiàn)old 是 log2 差異倍數(shù)掷匠。
保存結(jié)果到 BED 和 CSV 格式滥崩。
# 轉(zhuǎn)換到 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 分析
轉(zhuǎn)錄因子結(jié)合到開放的染色質(zhì)區(qū)域調(diào)控基因表達，轉(zhuǎn)錄因子結(jié)合需要識別出特定的序列讹语，這個特定序列稱之為 motif, 結(jié)合的位點稱為 TFBS (TF binding sites)钙皮。轉(zhuǎn)錄因子可以允許 TFBS 有一定的可塑性 （variations/flexible），所以 motif 序列不完全是固定死的顽决。Motif 大部分不長短条，6-12 bp 左右。人類大約有 1600 個轉(zhuǎn)錄因子才菠，其中超過 2/3 已鑒定了 motif.
motif 序列一般有兩種模式茸时，一是回文序列，比如 CACGTG 其反向互補序列也是 CACGTG. 二是兩段保守序列被一段非保守序列分隔赋访，這往往是因為結(jié)合的轉(zhuǎn)錄因子是二聚體屹蚊，分別識別兩段保守序列。
motif 的分析往往受假陽性困擾进每，這稱為無效定理（Futility Theorem）汹粤。比如說往往在基因組序列中觀察到大量的潛在轉(zhuǎn)錄因子結(jié)合位點，其中很少是真正起作用的田晚，大部分預(yù)測的轉(zhuǎn)錄因子結(jié)合位點是無效的嘱兼。所以 motif 分析后，要想辦法進行人工篩選贤徒。
MEME-ChIP 主要執(zhí)行以下步驟：
- 在輸入序列的中間區(qū)域（默認 100bp）進行 motif 發(fā)現(xiàn)（MEME, STREME）
- CentriMo 分析哪些 motif 是富集在區(qū)域中心的
- Tomtom 分析他們與已知 motif 相似性
- 根據(jù) motif 相似性對顯著結(jié)果歸類
- motif spacing analysis (SpaMo) (分析 motif 與結(jié)合在相鄰位置 motif 的物理作用) (
-spamo-skip參數(shù)跳過這步分析) - 分析 motif 結(jié)合位置（FIMO）芹壕，默認選取每一類別最顯著的 motif 進行分析
MEME-ChIP 默認輸入序列是約 500bp 長且中間 100bp 是 motif 區(qū)域。
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 分析假陽性問題和轉(zhuǎn)錄因子一般結(jié)合到啟動子區(qū)，只取注釋到啟動子區(qū)的 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 數(shù)據(jù)庫在 https://meme-suite.org/meme/db/motifs 下載。
MEME-ChIP 輸出結(jié)果在 meme.html 查看潘明；結(jié)果的匯總在 summary.tsv 文件行剂；combined.meme 文件包含所有被鑒定出的 Motif.
motif E value 表示該 motif 多大概率是統(tǒng)計錯誤出現(xiàn)，而不是真的 motif. E 值越小說明發(fā)現(xiàn)的 motif 越大概率為真钳降。
MEME-ChIP 的結(jié)果如果發(fā)現(xiàn)不方便批量提取結(jié)果厚宰，也可以自己單獨運行它的子分析。比方說它 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 文件轉(zhuǎn)換成 bigwig 文件客税。
bamCoverage --numberOfProcessors 8 --effectiveGenomeSize 2747877777 --normalizeUsing RPGC \
--outFileFormat bigwig -b ${bam_dir}/${sample}_MD.bam -o ${bw_dir}/${sample}.bw
參數(shù) --effectiveGenomeSize 的設(shè)置值在 Effective Genome Size — deepTools 3.4.3 documentation 查看况褪。
RPGC 計算方法：
Scale factor 計算方法：
(total number of mapped reads * fragment length) / effective genome size
第二步 computeMatrix 命令生成畫圖矩陣。但是先用 R 語言提取包含 peak 的啟動子區(qū)域更耻。注意輸出的 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 信號轉(zhuǎn)換為矩陣并熱圖可視化锯七。
library(rtracklayer)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
library(ChIPseeker)
library(EnrichedHeatmap)
library(circlize)
取得基因組啟動子區(qū)域，注意用 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 函數(shù)讀取 narrowPeak 文件【匏或者如果是 rtracklayer 包的 import 函數(shù)噪猾，它支持讀取更多的格式，比如 bigwig,bed 等筑累。
peak <- readPeakFile(peakPath)
將啟動子區(qū) peakset 數(shù)據(jù)轉(zhuǎn)換為畫圖矩陣袱蜡。注意將沒信號的啟動子區(qū)數(shù)據(jù)移除。
mat <- normalizeToMatrix(peak, promoter, extend = 0, value_column = "V5", k = 100, mean_mode = "w0")
# 移除沒信號啟動子區(qū)
mat <- mat[rowSums(mat) > 0,]
讀取后 narrowPeak 數(shù)據(jù) "V5" 列為峰信號強度值慢宗，所以 value_column 參數(shù)設(shè)置 "V5".
參數(shù) k 決定 bins/windows 數(shù)目（矩陣的列數(shù)）坪蚁。參數(shù) 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 函數(shù)畫圖。
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.
