最近一周开始学习ATAC数据分析，目前只是尝试了简单分析，新手小白，如果有小失误欢迎留言~
1.准备环境
创建conda环境（如果已有成熟环境不需再次创建）
conda create -n atac -y python=2 bwa
conda info --envs
激活环境
conda activate atac
2.安装所需要的包
conda install -c bioconda -y fastqc
conda install -c bioconda -y trimmomatic
conda install -y hisat2
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
conda install -y sambamba
conda install -y deeptools
3.对下机数据rowdata进行简单质控--可选
执行命令
fastqc -t 8 -o /public/home/sss/ss/ATAC/raw ATAC1-2_FKDL210215276-1a_1.fq.gz
fastqc -t 8 -o /public/home/sss/ss/ATAC/rawATAC1-2_FKDL210215276-1a_2.fq.gz
说明
-t: 线程
-o: 存放路径，不用指定前缀，默认为.fastq.gz前面的字段
*.fastq.gz:fastqc可以同时接受多个fastq.gz文件，因此采用正则表达式*表示全部
结果会产生两个文件：
image
image
将HTLM文件下载到本地打开
image
注： 接头含量太多需要过滤以及去接头
4.使用trim_galore去除接头，trim_galore可以自动识别接头
下载trim_galore
curl -fsSL[https://github.com/FelixKrueger/TrimGalore/archive/0.6.6.tar.gz](https://github.com/FelixKrueger/TrimGalore/archive/0.6.6.tar.gz)-o trim_galore.tar.gz
tar xvzf trim_galore.tar.gz
######################################
去除NexteraPE
trim_galore -q 20 --phred33 --stringency 3 --length 20 -e 0.1 --paired./ATAC1-2_FKDL210215276-1a_1.fq.gz ./ATAC1-2_FKDL210215276-1a_2.fq.gz -o ./
产生新的.gz文件
image
image
5.对去除接头后的数据质检
执行命令
fastqc -t 8 -o /public/home/sss/ss/ATAC/raw ATAC1-2_FKDL210215276-1a_1_val_1.fq.gz \fastqc -t 8 -o/public/home/sss/ss/ATAC/raw ATAC1-2_FKDL210215276-1a_2_val_2.fq.gz
产生两个文件：
image
image
将htlm文件下载并打开:
image
目前已没有接头序列，重复序列比较多，后面比对后需要进一步去重
6.使用hisat2比对
最重要建立参考基因组的索引
mkdir reference && cd reference
参考基因组
wget http://rice.plantbiology.msu.edu/pub/data/Eukaryotic_Projects/o_sativa/annotation_dbs/pseudomolecules/version_7.0/all.dir/all.con
注释文件
wget http://rice.plantbiology.msu.edu/pub/data/Eukaryotic_Projects/o_sativa/annotation_dbs/pseudomolecules/version_7.0/all.dir/all.gff3
Hisat2 建立索引
hisat2-build -p 10 reference/all.con reference/Osativa
建立后的索引文件
image
建立align.sh文件
clean data
/public/home/sss/software/anaconda3/envs/chipseq/bin/hisat2 --dta -t-x /public/home/sss/ss/genome/NIP -1/public/home/sss/ss/ATAC/cleandata/ATAC1-2_FKDL210215276-1a_2_val_2.fq.gz -2 /public/home/sss/sss/ATAC/cleandata/NIP-1_FRAS210228372-2r_2.clean.fq.gz | /public/home/sss/anaconda2/envs/chipseq/bin/samtools sort -@ 8 -O bam -o-/public/home/sss/ss/ATAC/cleandata/ATAC_1-2.bam
注意：由于sam文件较大，因此我们直接跳过sam，使用samtools转换为排序后的bam文件，注意【|】不要漏掉，且排序是依据reads比对到基因组的位置排的
-O: 表示输出的bam文件
-o: 输出的bam文件名
-t与@: 线程数
比对后产生：
image
7.比对后的序列去重和排序
选用sambamba来去重复
sambamba markdup -r ATAC_1-2.bam ATAC_1-2.sambamba.rmdup.bam
smarttool提取可靠的比对结果
samtools view -f 2 -q 30 -o ATAC_1-2.rmdup.q30.bamATAC_1-2.sambamba.rmdup.bam
smarttool对比对结果排序
samtools sort -O bam -@ 4 -o ./ATAC_1-2.fraw.bam ATAC_1-2.rmdup.q30.bam
产生的bam 文件:
image
8.计算插入片段长度，FRiP值，IDR计算重复情况
提取片段的长度
从sam文件中求得insert size：在“read mapped in proper pair”的前提下，取第九列
samtools view ATAC_1-2.fraw.bam |awk -F'\t' 'function abs(x){return ((x < 0.0) ? -x : x)} {print$1"\t"abs($9)}' | sort | uniq | cut -f2 > fragment.length.txt
提取出片段长度，R代码如下：
image
得到结果图:
image
4
FRiP_score的计算
Reads=$(bedtools intersect -a ATAC_1-2_summits.bed -bATAC_1-2_peaks.narrowPeak |wc -l|awk '{print $1}')
totalReads=$(wc -l ATAC_1-2_summits.bed|awk '{print $1}')
echo ${Reads} ${totalReads} ATAC_1-2 'FRiP_score:' $(echo"scale=2;100*${Reads}/${totalReads}"|bc)'%'
9.使用过滤好的数据call peak
samtools构建index
samtools index ATAC_1-2.fraw.bam
bedtools转换为bed 文件
bedtools bamtobed -i ATAC_1-2.fraw.bam > ATAC_1-2.bed
使用macs2 callpeak
macs2 callpeak -t ATAC_1-2.bed -g380699722 --nomodel --shift -100 --extsize 200 -n ATAC_1-2 --outdir ./peak
得到的文件：
image
10.deeptools 进行TSS可视化
将 bam 文件转化为 bw 文件
bamCoverage -p 8 --bam ATAC_1-2.sambamba.rmdup.bam --binSize 10--centerReads --smoothLength 14 --normalizeUsing RPKM -ocoverageATAC_1-2.bigwig
bamCoverage -p 8 --bam ATAC_1-2.fraw.bam --binSize 10 --centerReads--smoothLength 14 --normalizeUsing RPKM -o coverageATAC_1-2.fraw.bigwig
建立.deeptools_TSS.sh文件
计算reference-point---需要基因组的Refseq文件--下载或者自己创制
mkdir -p /public/home/sss/ss/ATAC/tss
cd /public/home/sss/ss/ATAC/tss
source activate atac
computeMatrix reference-point --referencePoint TSS -p 15 \
-b 2000 -a 2000 \
-S /public/home/sss/ss/ATAC/cleandata/coverageATAC_1-2.bw \
-R /public/home/sss/ss/ATAC/NIP.deptools.bed \
--skipZeros \
-o matrix1_test_TSS.gz \
--outFileSortedRegions regions1_test_genes.bed
## both plotHeatmap andplotProfile will use the output from computeMatrix
plotHeatmap -m matrix1_test_TSS.gz -out test_Heatmap.png
plotHeatmap -m matrix1_test_TSS.gz -out test_Heatmap.pdf --plotFileFormat pdf --dpi 720
plotProfile -m matrix1_test_TSS.gz -out test_Profile.png
plotProfile -m matrix1_test_TSS.gz -out test_Profile.pdf --plotFileFormat pdf --perGroup --dpi 720
##########################################################################################
结果文件：
image
image
11.peaks注释
结构注释----会将peak所落在基因组上的区域结构注释出来，比如说启动子区域，UTR区域，内含子区域等等。同时，也会将与peak最临近的基因注释出来，非常好用
使用gff/gtf构建一个TxDb
R代码：
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install()
BiocManager::install(c("GenomicFeatures","AnnotationDbi"))
#install.packages("AnnotationDbi")
library(clusterProfiler)
library(ChIPseeker)
library(GenomicFeatures)
library(AnnotationDbi)
txdb <- makeTxDbFromGFF("D:/生信学习/ATAC/ricegenome/NIP.gtf",
format="gtf") #也可以使用gtf
keytypes(txdb) #感兴趣的话，可以用以下方法探索txdb都包含了什么内容
keys(txdb)
读入单个summits文件
peaks <- readPeakFile("D:/生信学习/ATAC/ATAC_1-2_summits.bed")
结构注释
peakAnno <- annotatePeak(peaks,
TxDb=txdb,
tssRegion=c(-1000, 1000))
#注释完，进行可视化，多种图可供选择
plotAnnoBar(peakAnno)
plotDistToTSS(peakAnno)
vennpie(peakAnno)
plotAnnoPie(peakAnno)
install.packages("ggupset")
library(ggupset)
upsetplot(peakAnno)
install.packages("ggimage")
library(ggimage)
upsetplot(peakAnno, vennpie=TRUE)
library(ggplot2)
最后将我们的注释结果转为数据框，便于查看
df <- as.data.frame(peakAnno)
将注释到的基因提取出来（第14列），用于后续功能分析
gene <- df[,14]
head(gene)
获得结果文件;
image
image
功能注释-----需要物种的OrgDb注释库
R代码：
#GO 富集分析
library(AnnotationHub)
hub <- AnnotationHub()
query(hub, "Oryza sativa")
nip_org <- hub[['AH96213']] ##下载目标物种到org数据
keytypes(nip_org)#查看可使用的key
length(keys(nip_org))
##我们可以使用riceidconverter来转换ID，甚至于直接检索GO注释：
require(clusterProfiler)
head(gene)
gene_SYMBOL <- riceidconverter::RiceIDConvert(gene,fromType='MSU',toType = 'SYMBOL')
head(gene_SYMBOL)
gene_SYMBOL1 <- as.character(gene_SYMBOL$SYMBOL)
head(gene_SYMBOL1)
gene_SYMBOL2 <- gene_SYMBOL1[grep("LOC*" ,gene_SYMBOL1)]
C4_3DAP_go = enrichGO(gene_SYMBOL2,keyType = "SYMBOL",OrgDb=nip_org, pvalueCutoff=1, qvalueCutoff=1)
dotplot(C4_3DAP_go,showCategory = 20)
cnetplot(C4_3DAP_go,circular=T, ###画为圈图
colorEdge=T,showCategory =15)
cnetplot(C4_3DAP_go, showCategory = 5)
write.table(C4_3DAP_go, 'C4_3DAP_go.txt', sep = '\t', quote = FALSE,row.names = FALSE)
kegg
install.packages("R.utils")
R.utils::setOption("clusterProfiler.download.method",'auto')
head(gene)
gene_list <- bitr(geneID = gene_SYMBOL1,
fromType ="SYMBOL",
toType =c("ENTREZID"),
OrgDb = nip_org)
head(gene_list)
genelist2<-(gene_list$ENTREZID)
head(genelist2)
compKEGG <- enrichKEGG(genelist2,organism="osa",keyType="kegg",pvalueCutoff=0.05,pAdjustMethod="BH",qvalueCutoff=0.2)
dotplot(compKEGG, showCategory = 20, title = "KEGG Pathway EnrichmentAnalysis")
################################################################
获得结果文件：
image
image
image
12.homer寻找motif
使用homer 进行 motif分析
conda 环境下安装 homor环境
conda install -c bioconda -y homer
awk '{print $4"\t"$1"\t"$2"\t"$3"\t+"}'\
ATAC_1-2_peaks.narrowPeak > ATAC_1-2_peaks.tmp
findMotifsGenome.pl ATAC_1-2_peaks.tmp \
/public/home/sss/ss/genome/NIP.fa \
./ATAC_1-2.motif \
-size 200 -len 8,10,12
#########################################
image
产生的文件：
image
将该目录下所有文件下载到本地，打开homerResults.html 文件，即可得到如下结果展示：
image
输出结果按照p-value排序，最后一列是一个链接到motif文件的超链接，可以从这个文件中找到包含此motif的其他序列。
在Best Match/Details 列中，HOMER将会展示与denovo motif最匹配的已知motif（该列还包含有一个名为More information 的超链接，
点击后会出现以下页面：该页面包含motif的一些基本信息，如链接到motfi文件的超链接，且可以生成motif logo的pdf格式。
参考文章链接：
原文链接：https://blog.csdn.net/u012110870/article/details/102804623
原文链接： https://www.jianshu.com/p/9aa719faa4b5
原文链接：https://cloud.tencent.com/developer/article/1360799
原文链接：[https://www.jianshu.com/p/9a31f5f01e7b
