ATACSeq 开放染色质分析【高通量测序（NGS）原理】表观组学

生信 | ATAC-Seq基础分析+高级分析+多组学分析

2021-07-01  本文已影响0人  生信卷王

写在前面

ATAC-seq的首次提出：【文章原文】【scATAC原文】
ATAC-seq植物界的里程碑：【文章原文】【中文解读】
ATAC-seq针对不同植物物种进行研究：【文章原文】【中文解读】
ATAC-seq针对不同植物组织进行研究：【文章原文】【中文解读】
ATAC-seq针对不同处理条件进行研究：【文章原文】【中文解读】
ATAC-seq生信分析方法综述：【文章原文】【中文解读1】【中文解读2】
ATAC-seq基础原理一文了解：【一文了解ATAC-seq】

一、基础知识

1. 解密表观遗传学的三个方向与测序方法

[^xr2]:2. 名词解释

二、基础分析

1. 分析前的准备

1.1 下载ATAC测序数据（或者用你自己的数据也行）
# 提取100w条reads代码
ls *gz|while read id;do zcat $id|head -4000000|gzip > ./100/${id};done
image.png
1.2 下载tair10基因组与gff3注释文件
wget https://www.arabidopsis.org/download_files/Genes/TAIR10_genome_release/TAIR10_chromosome_files/TAIR10_chr_all.fas
wget https://www.arabidopsis.org/download_files/Genes/TAIR10_genome_release/TAIR10_gff3/TAIR10_GFF3_genes_transposons.gff
# 改名
mv TAIR10_chr_all.fas TAIR10_genome.fa
mv TAIR10_GFF3_genes_transposons.gff TAIR10.gff3
1.3 安装miniconda
1.4 所需软件安装（分析中所用到的软件都汇集于此，下文不再介绍如何下载）\color{red}{请先激活conda环境再安装！！！！} 
conda install -c bioconda -y fastqc
conda install -c bioconda -y trimmomatic
conda install -c bioconda -y bwa
conda install -c bioconda -y samtools
conda install -c bioconda -y picard
conda install -c bioconda -y macs2
conda install -c bioconda -y bedtools
conda install -c bioconda -y deeptools
conda install -c bioconda -y homer

\color{red}{后续分析是使用SRR7512044用做测试（如下），但思路和方法都可以借鉴}

测试数据

2. 测序质量检测

fastqc -t 8 -o ./ *.fastq.gz

-t: 线程
-o: 存放路径，不用指定前缀，默认为.fastq.gz前面的字段
*.fastq.gz: fastqc可以同时接受多个fastq.gz文件，因此采用正则表达式【*】表示全部

2. 测序质量检控制

# 定义一些变量，后面修改起来方便
filename=SRR7512044

trimmomatic PE -threads 10 \
./${filename}_1.fastq.gz ./${filename}_2.fastq.gz \
./${filename}_1.clean.fq.gz ./${filename}_1.drop.fq.gz \
./${filename}_2.clean.fq.gz ./${filename}_2.drop.fq.gz \
HEADCROP:15 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:24

PE: 双端模式。需要两个输入文件，正向测序序列和反向测序序列：sample_R1.fastq.gz和sample_R2.fastq.gz；
以及四个输出文件：sample_1.clean.fq.gz和sample_1.drop.fq.gz；sample_2.clean.fq.gz和sample_2.drop.fq.gz
-threads: 线程数
HEADCROP: 从 reads 的开头切掉指定数量的碱基，即从fastqc中看的
LEADING: 从 reads 的开头切除质量值低于阈值的碱基
TRAILING: 从 reads 的末尾开始切除质量值低于阈值的碱基
SLIDINGWINDOW: 从 reads 的 5' 端开始，进行滑窗质量过滤，切掉碱基质量平均值低于阈值的滑窗
MINLEN: 如果经过剪切后 reads 的长度低于阈值则丢弃这条 reads，即上面的计算方法：（76-15）×40%≈24

3. 建立索引与序列比对

bwa index -p /路径/bwaIndexTair /路径/TAIR10_genome.fa

-p: 索引前缀，后缀自动补充
TAIR10_genome.fa: 基因组文件

# 定义一些变量，后面用
filename=SRR7512044
index=/路径/bwaIndexTair

bwa mem -M -t 8 \
-R "@RG\tID:${filename}\tSM:${filename}\tLB:WXS\tPL:Illumina" \
$index \
${filename}_1.clean.fq.gz \
${filename}_2.clean.fq.gz \
| samtools sort -O bam -@ 10 -o ./${filename}.raw.bam
#无论干啥，生成bam文件后就要建立bam索引，并检查比对情况
samtools index -@ 10 -b ./${filename}.raw.bam
samtools flagstat ./${filename}.raw.bam > ./${filename}.raw.stat
cat ${filename}.raw.stat

-M: 将 shorter split hits 标记为次优，以兼容 Picard’s markDuplicates 软件
-R: STR 完整的read group的头部，可以用 '\t' 作为分隔符， 在输出的SAM文件中被解释为制表符TAB. read group 的ID，会被添加到输出文件的每一个read的头部，强推去看关于-R的解释【-R的使用】
-O: 表示输出的bam文件
-o: 输出的bam文件名
-t与@: 线程数

14608455 + 0 in total (QC-passed reads + QC-failed reads) reads总数
37967 + 0 secondary                                       出现比对到参考基因组多个位置的reads数
0 + 0 supplementary                                       可能存在嵌合的reads数
0 + 0 duplicates                                          重复的reads数
14590894 + 0 mapped (99.88% : N/A)                        比对到参考基因组上的reads数
14570488 + 0 paired in sequencing                         属于PE read的reads总数。
7285244 + 0 read1                                         PE read中Read 1 的reads 总数。
7285244 + 0 read2                                         PE read中Read 2 的reads 总数。
14507068 + 0 properly paired (99.56% : N/A)               完美比对的reads总数。PE两端reads比对到同一条序列，且根据比对结果推断的插入片段大小符合设置的阈值。
14551500 + 0 with itself and mate mapped                  PE两端reads都比对上参考序列的reads总数。
1427 + 0 singletons (0.01% : N/A)                         PE两端reads，其中一端比上，另一端没比上的reads总数。
26260 + 0 with mate mapped to a different chr             PE read中，两端分别比对到两条不同的序列的reads总数。
17346 + 0 with mate mapped to a different chr (mapQ>=5)   PE read中，两端分别比对到两条不同

4. 去除PCR重复并再次进行质控

picard MarkDuplicates \
REMOVE_DUPLICATES=true \
I=${filename}.raw.bam \
O=${filename}.rmdup.bam \
M=${filename}.rmdup.log
# 实时监测我们的数据发生了什么变化
samtools index ${filename}.rmdup.bam
samtools flagstat ${filename}.rmdup.bam > ./${filename}.rmdup.stat

REMOVE_DUPLICATES=true: 表示将检测出的PCR duplication直接去除；
I: 或者INPUT指定输入的bam文件；
O: 或者OUTPUT指定输出去除PCR duplication的bam文件；
M: 或者METRICS_FILE指定输出的matrix log文件。

samtools view -h -f 2 -q 30 ${filename}.rmdup.bam \
| grep -v -e "mitochondria" -e "*" -e "chloroplast" \
| samtools sort -O bam -@ 10 -o - > ${filename}.last.bam
# 实时监测我们的数据发生了什么变化
samtools index ${filename}.last.bam
samtools flagstat ${filename}.last.bam > ./${filename}.last.stat

-h: 输出的sam文件带header，默认不带
-f: 输出含有所有flag的reads
-q 比对的最低质量值，一般20/30都可以
其他参数介绍看 samtools 使用简述，需要说明的是由于比对到线粒体和叶绿体的reads不是我们关注的重点，因此使用grep剔除掉，当然了，做动物的话需要剔除什么视情况而定了。

5. 统计过滤后的结果中Insertsize长度分布

picard CollectInsertSizeMetrics \
I=./${filename}.last.bam \
O=./${filename}.last.insertsize \
H=./${filename}.last.hist.pdf

有关下图的介绍，请看【给你bam文件，你会画插入片段长度分布图吗？】

5. \color{red}{peaks\ calling}与FRiP_score/SPOT_score的计算

effective_genome_size=119481543

bedtools bamtobed -i ${filename}.last.bam > ${filename}.last.bed
macs2 callpeak \
-t ${filename}.last.bed \
-g ${effective_genome_size} \
--nomodel --shift -100 --extsize 200 \
-n ${filename}.last \
-q 0.01 --outdir ./peaks
shift 与 extsize 到底在干啥？ 移峰 
前几行是callpeak时的命令。后面是：
1.染色体号
2.peak起始位点
3.peak结束位点
4.peak区域长度
5.peak的峰值位点（summit position）
6.peak 峰值的高度（pileup height at peak summit, -log10(pvalue) for the peak summit）
7.peak的富集倍数（相对于random Poisson distribution with local lambda）

narrowPeak文件是BED6+4格式，可以上传到IGV查看
1.染色体号
2.peak起始位点
3.peak结束位点
4.peak name（样本名+_peak_1的后缀之类的）
5.-10*log10qvalue
6.正负链
7.fold change
8.-log10pvalue
9.-log10qvalue
10.相对于峰的峰顶位置（relative summit position to peak start）

BED格式的文件，包含peak的summits（峰顶）位置
1.染色体号
2.峰的峰顶起始位点（注意是峰顶，不是峰，他只是一个点）
3.峰的峰顶起始位点（因此二者就差一个碱基的位置）
4.峰顶（样本名+_peak_1的后缀之类的）
5.-log10pvalue

Reads=$(bedtools intersect -a ${filename}.last.bed -b ${filename}_peaks.narrowPeak |wc -l|awk '{print $1}')
totalReads=$(wc -l ${filename}.last.bed|awk '{print $1}') 
echo ${Reads} ${totalReads} ${filename} 'FRiP_score:' $(echo "scale=2;100*${Reads}/${totalReads}"|bc)'%'

二、高级分析

所谓的高级并不高大上，做了有可能也不能让你的文章提升档次，望周知！

1. IGV可视化

bamCoverage --bam ${filename}.last.bam -o ${filename}.last.bw \
    --binSize 10 \
    --normalizeUsing RPGC \
    --effectiveGenomeSize ${effective_genome_size} \
    --ignoreForNormalization chrX \
    --extendReads

注：如果发现没有这个包pysam.libchtslib，直接安装pip install numpydoc
-o: 生成bw文件
binSize: 自行设定覆盖度计算的窗口大小(bin)
normalizeUsing: 数据准化方法，有三种RPGC，RPKM，CPM，BPM
ignoreForNormalization: 设置你想忽略的染色体名称
extendReads: 让reads扩展至fragment的大小

IGV可视化bw、narrowpeak，bed，bam，gff文件

2. TSS/TES 可视化

awk -vFS='[\t=;]' -vOFS="\t" \
'{if ($3=="gene") print $1,$4,$5,$10,$6,$7}' \
TAIR10.gff|sed 's/Chr//g' > refseq.bed  
#由于tair官网上提供的基因组文件中染色人命名没有Chr，为保持一致，这里去掉

computeMatrix reference-point --referencePoint TSS  -p 15 \
    -b 2000 -a 2000 \
    -R ./refseq.bed \
    -S *.bw \
    -o TSS.gz \
    --skipZeros \
    --outFileSortedRegions Heatmap1sortedRegions.bed
plotHeatmap -m TSS.gz -out TSS.Heatmap.pdf --plotFileFormat pdf
plotProfile -m TSS.gz -out TSS.Profile.pdf --plotFileFormat pdf --perGroup

computeMatrix scale-regions  -p 15  \
    -b 2000 -a 2000 \
    -R ./refseq.bed \
    -S *.bw \
    --skipZeros -o body.gz
plotHeatmap -m body.gz -out body.Heatmap.pdf --plotFileFormat pdf
plotProfile -m body.gz -out body.Profile.pdf --plotFileFormat pdf --perGroup

3. 相关性可视化（依旧是使用deeptools）

multiBigwigSummary bins -b *.last.bw -o number.of.bins.npz

plotCorrelation -in number.of.bins.npz \
--corMethod pearson \
--skipZeros \
--plotTitle "Pearson Correlation of Average Scores Per Transcript" \
--plotFileFormat pdf \
--whatToPlot heatmap --colorMap RdYlBu --plotNumbers \
-o heatmap_PearsonCorr_bigwigScores.pdf \
--outFileCorMatrix PearsonCorr_bigwigScores.tab
heatmap 
plotCorrelation -in number.of.bins.npz \
--corMethod pearson \
--skipZeros \
--plotTitle "Pearson Correlation of Average Scores Per Transcript" \
--plotFileFormat pdf \
--whatToPlot scatterplot \
-o scatterplot_PearsonCorr_bigwigScores.pdf \
--outFileCorMatrix PearsonCorr_bigwigScores.tab
相关性点图

5. ChIPseeker对peaks进行注释和可视化

5.1 结构注释

5.1.1 准备工作：导入包，没有的安装吧。另外需要使用gff/gtf构建一个TxDb。因为官网上TxDb也不是每一个物种都有（目前共43个），因此我们自己手动构建。

#if (!requireNamespace("BiocManager", quietly = TRUE))
#    install.packages("BiocManager")
#BiocManager::install("clusterProfiler")

library(clusterProfiler)
library(ChIPseeker)
library(GenomicFeatures)
txdb <- makeTxDbFromGFF("public/bioinfoYu/genome/tair10/TAIR10.gff",
                        format="gff")    #也可以使用gtf
keytypes(txdb)    #感兴趣的话，可以用以下方法探索txdb都包含了什么内容
keys(txdb)

5.1.2 单个文件的结构注释（不推荐在这一步加OrgDB的内容，没用）

#读入单个summits文件
peaks <- readPeakFile("SRR7512044_summits.bed")
#结构注释
peakAnno <- annotatePeak(peaks,
                         TxDb=txdb,
                         tssRegion=c(-1000, 1000))

上面的tssRegion取值有说法的，因为判定"Promoter", "Exon", "Intron", "Intergenic"等结构的依据就是根据范围来的。比如定义1000以内的，则在1kbp之内规定为"Promoter (<=1kb)"，而1k-2k内定义为"Downstream (1-2kb)"。但你定义2000以内呢，则在1kbp之内规定为"Promoter (<=1kb)"，而1k-2k内定义为"Promoter (1-2kb)"。虽然，ChIPseeker会标记出来区域，但还是三思后再取范围。

#注释完，进行可视化，多种图可供选择
plotAnnoBar(peakAnno)
plotDistToTSS(peakAnno)
vennpie(peakAnno)
plotAnnoPie(peakAnno)
#install.packages("ggupset")
library(ggupset)
upsetplot(peakAnno)
#install.packages("ggimage")
library(ggimage)
upsetplot(peakAnno, vennpie=TRUE)

#最后将我们的注释结果转为数据框，便于查看
df <- as.data.frame(peakAnno)
#将注释到的基因提取出来（第14列），用于后续功能分析
gene <- df[,14]
提供可视化的方法

5.1.3 多个文件的结构注释

#一次也可以读入多个summits文件，使用list存储，然后使用lapply注释
files = list(peak1 = "peaks.bed", peak2 = ("peak.bed"))
peakAnnoList <- lapply(files, 
                       annotatePeak,
                       TxDb=txdb,
                       tssRegion=c(-2000, 2000))
plotAnnoBar(peakAnnoList)
plotDistToTSS(peakAnnoList)
读入多个summits文件的可视化

5.2 功能注释

5.2.1 可选支线
由于这里需要物种的OrgDb注释库，但并不是所有的物种都有啊！怎么办？戳看【构建自己的OrgDb库】。如果你足够幸运，研究的物种正好在官网仅提供的20个物种中存在，那就拿来用吧。除了以上这些都好说，还有是\color{red}{ID转换}，后面再说。
5.2.2 OrgDb下载
本次测试是用拟南芥的样本，在官网中可以直接将所需要的的OrgDb库下载到服务器上，这样就不用在R里面重复在网上获取浪费时间了。比如说拟南芥，就可以下载。当然也可以直接install了，看你喜欢哪种了。

wget http://www.bioconductor.org/packages/release/data/annotation/src/contrib/org.At.tair.db_3.13.0.tar.gz
tar -zvxf org.At.tair.db_3.13.0.tar.gz

5.2.3 导入OrgDb

#从本地下载OrgDb库（如果你直接下载到了R的lib里面，就直接library也可）
install.packages("/public/bioinfoYu/genome/tair10/org.At.tair.db", repos=NULL)
#另外一种方法
#if (!requireNamespace("BiocManager", quietly = TRUE))install.packages("BiocManager")
#BiocManager::install("org.At.tair.db")
#载入库
library(org.At.tair.db)
#感兴趣的话，看看org.At.tair.db包含的都是些啥
#查看包含的列名
keytypes(org.At.tair.db)
#查看包含的基因ID
keys(org.At.tair.db)
#随便选择一个基因和一些列名，看一看到第是什么东东，不过因为导入了其他包，下面命令报错，我至今还没有解决
select(org.At.tair.db,keys="AT1G01010",columns = c("ENTREZID","GO","SYMBOL"))

5.2.4 ID转换（非必须）
上面5.1.2已经得到peak附近的基因集【gene】了，如果需要ID转换，就是用Y叔的clusterProfiler中的bitr函数

#首先，去除版本号，如AT1G50040.1后面的.1，当然我们的测试集没有，因此不用做
#gene_rm <- gsub(".[0-9]+$","",gene)
#下面进行ID转换
keytypes(org.At.tair.db)
gene_last <- bitr(geneID = gene, 
                    fromType = "TAIR",
                    toType = c("ENTREZID", "SYMBOL", "GENENAME"),
                    OrgDb = org.At.tair.db)

# geneID    需要转换的ID
# fromType  当前ID类型
# toType    转换成什么ID，使用keytypes()查看有哪些类型
# OrgDb     注释数据库

5.2.5 GO功能注释

# 因为keyType不接受数据框类型的数据，因此我们转换为字符型
genelist <- gene.symbol[,2]  #注意我这里把ENTREZID单独拿出来了，所以下面自然是keyType = "ENTREZID",
go_result <- enrichGO(genelist,
                      keyType = "ENTREZID", 
                      OrgDb = org.At.tair.db, 
                      ont = "BP", #BP,MF,CC,ALL可选
                      pAdjustMethod = "BH", 
                      qvalueCutoff = 0.05, 
                      readable = FALSE)
#可视化，showCategory=20调整显示的条目多少
dotplot(go_result, showCategory=20)
#将结果转为数据框，可以选择导出
goresult <- as.data.frame(go_result)
GO可视化

5.2.6 KEGG通路注释

#geneClusters必须是数据框类型，ID必须是ENTREZID类型
class(gene.symbol)
compKEGG <- compareCluster(geneClusters=gene.symbol, 
                           fun = "enrichKEGG",
                           organism = "ath",    #在 http://www.genome.jp/kegg/catalog/org_list.html 上查看物种简写
                           pvalueCutoff  = 0.05, 
                           pAdjustMethod = "BH")
dotplot(compKEGG, showCategory = 20, title = "KEGG Pathway Enrichment Analysis")
KEGG可视化

6. homer寻找motif

# 提取peaks的位置信息文件
awk '{print $4"\t"$1"\t"$2"\t"$3"\t+"}' \
SRR7512044_peaks.narrowPeak > SRR7512044.homer_peaks.tmp
# 寻找motif
findMotifsGenome.pl SRR7512044.homer_peaks.tmp \
./TAIR10_genome.fa \  #指定参考基因组
./SRR7512044.motif \  #输出文件夹的名字
-size 200 -len 8,10,12

7. DiffBind找差异peak

首先根据我自己的理解将找差异peak的原理说一下：既然是找差异的peak，那首先就要保证不同样本间peak所在的区间基本一致才可以。因此第一步就是找到多个样本相同的peak区间。第二步，基于这些相同区域的peak开始找差异，这个地方的差异和传统的RNA-Seq找差异基因的原理一模一样，都是对落在峰区间的reads进行计数，数量不同的话差异就由此而来。但软件都帮我们做了，就不用考虑那么多了

SampleID：不能有重复
Factor：进行分组，后面要做差异分析的依据
Replicate：重复的编号
bamReads：最后生成的bam文件。全路径+bam文件名
Peaks：macs2生成的narrowPeak文件。全路径+narrowPeak文件名
PeakCaller：说明peak类型，使用narrowPeak文件则为narrow。使用xls文件，则用macs。使用bed文件，则用bed

# 安装包
if (!requireNamespace("BiocManager", quietly = TRUE))install.packages("BiocManager")
BiocManager::install("DiffBind")
library(DiffBind)
#导入数据，dba函数会将会将bam文件一同导入，因此csv中的路径一定要准确
tamoxifen <- dba(sampleSheet="atDiff.csv")
#查看样本间的相关性
plot(tamoxifen)

#计算每个样本中的reads数
tamoxifen <- dba.count(tamoxifen)
#简单查看计算结果
info <- dba.show(tamoxifen)
libsizes <- cbind(LibReads=info$Reads, FRiP=info$FRiP, PeakReads=round(info$Reads * info$FRiP))
rownames(libsizes) <- info$ID
libsizes

#统计结果
       LibReads FRiP PeakReads
cmt2_1  6720462 0.29   1948934
cmt2_2  7747491 0.33   2556672
wild_1  7078309 0.36   2548191
wild_2  8030752 0.34   2730456
ddcc_1  6370232 0.30   1911069
ddcc_2  6956429 0.35   2434750

#默认基于测序深度对数据进行标椎化
tamoxifen <- dba.normalize(tamoxifen)
norm <- dba.normalize(tamoxifen, bRetrieve=TRUE)
norm
#查看进行标准化的文库大小(library-size)，其实就是各样本文库大小总和的平均值，你可以计算验证一下
normlibs <- cbind(FullLibSize=norm$lib.sizes, NormFacs=norm$norm.factors, NormLibSize=round(norm$lib.sizes/norm$norm.factors))
rownames(normlibs) <- info$ID
normlibs

#标椎化结果
       FullLibSize  NormFacs NormLibSize
cmt2_1     6720462 0.9398443     7150612
cmt2_2     7747491 1.0834724     7150612
wild_1     7078309 0.9898885     7150612
wild_2     8030752 1.1230859     7150612
ddcc_1     6370232 0.8908652     7150612
ddcc_2     6956429 0.9728438     7150612

#分组，格式是表头在最前面，要分的组依次写在后面，只能两两比较，因此后面只能写两组，但可以多执行几次，都会追加到tamoxifen 中
#分组1，后面使用contrast=1单独查看
tamoxifen <- dba.contrast(tamoxifen, contrast=c("Factor","ddcc","wild"))
#分组2，后面使用contrast=2单独查看
tamoxifen <- dba.contrast(tamoxifen, contrast=c("Factor","cmt2","wild"))
#当然还可有分组3,4,5等，均可以使用contrast=number单独查看
tamoxifen
#按照分组分别进行差异分析，默认使用DESeq2进行计算，可以选择method = DBA_EDGER(edgR)，或者两个都要method = DBA_ALL_METHODS
tamoxifen <- dba.analyze(tamoxifen)
dba.show(tamoxifen, bContrasts=TRUE)
#查看差异分析的结果与导出为csv文件
#查看第1组差异分析的结果
tamoxifen.DB_1 <- dba.report(tamoxifen,contrast=1)
tamoxifen.DB_1
write.csv(tamoxifen.DB_1, file="diffBind_result.csv")

#查看第2组差异分析的结果
tamoxifen.DB_2 <- dba.report(tamoxifen,contrast=2)
tamoxifen.DB_2
write.csv(tamoxifen.DB_2, file="diffBind_result.csv")

##传统的上下调在找差异peak中称为“Gain”或“Loss”
#查看第1组的差异peak数量“Fold>0或<0”控制是“Gain”或“Loss”
sum(tamoxifen.DB$Fold>0,contrast=1)
sum(tamoxifen.DB$Fold<0,contrast=1)
#查看第2组的差异peak数量“Fold>0或<0”控制是“Gain”或“Loss”
sum(tamoxifen.DB$Fold>0,contrast=2)
sum(tamoxifen.DB$Fold<0,contrast=2)

#韦恩图可视化
dba.plotVenn(tamoxifen, contrast=1, bDB=TRUE, bGain=TRUE, bLoss=TRUE, bAll=FALSE)

#PCA
dba.plotPCA(tamoxifen) #这里可以使用所有样本进行PCA
dba.plotPCA(tamoxifen, contrast=1, label=DBA_FACTOR)#单独对分组1进行PCA

#MA plots
dba.plotMA(tamoxifen, contrast=1)

#Volcano plots
dba.plotVolcano(tamoxifen, contrast=1)

#Box plots
dba.plotBox(tamoxifen, contrast=1)

#Heatmap
hmap <- colorRampPalette(c("red", "black", "green"))(n = 13)
#对所有样本的所有的差异peak画热图
dba.plotHeatmap(tamoxifen,correlations=FALSE,scale="row",colScheme = hmap)
#单独对分组1中所有的差异peak画热图
dba.plotHeatmap(tamoxifen, contrast=1, correlations=FALSE, scale="row", colScheme = hmap)
dba.plotHeatmap(tamoxifen,correlations=FALSE,scale="row",colScheme = hmap)

#Profiling and Profile Heatmaps
dba.plotProfile(tamoxifen)

8. footprint分析

目前我还没有找到使用于所有物种的软件，HINT-ATAC用来做footprint分析的话，应该只支持不几个物种，比较鸡肋。ATACseqQC应该是一款不错的软件，最近还在探索。各位如有建议，欢迎提出。

三、ATAC-seq与多组学数据联合分析

转录因子ChIP-seq：由于大部分转录因子结合的是染色质开放区域，所以ATAC-seq的Peak可能和转录因子ChIP-seq的Peak存在部分重叠的情况，而且ATAC-seq得到的Peak长度往往更长，因此ATAC-seq数据和转录因子ChIP-seq数据可以相互验证。转录因子在ChIP-seq中独有的Peak暗示这个转录因子可能是结合在异染色质区域的驱动型转录因子（Pioneer TFs），驱动型转录因子随后招募染色质重塑复合体以及其它转录因子开始转录。另外，联合分析已经报道的ChIP-seq数据可以更准确地分析转录因子的足迹。
组蛋白修饰ChIP-seq：ATAC-seq数据同样可以和组蛋白修饰ChIP-seq数据进行联合分析，其中转录激活性修饰（H3K4me3，H3K4me1和H3K27ac等）与染色质开放程度呈正相关，转录抑制性修饰（H3K27me3）与染色质开放程度呈负相关。联合已知的增强子和启动子之间的相互作用数据也可以帮助构建调控网络。
RNA-seq：ATAC-seq数据可以通过联合分析RNA-seq数据来发现哪些差异表达的基因是受染色质可及性调控的，进一步可以推测这些差异表达的基因哪些是受开放染色质中具有motif和footprint的转录因子调控的，因此ATAC-seq与RNA-seq的联合分析有助于破译基因调控网络和细胞异质性。

上一篇下一篇

猜你喜欢

热点阅读