最近一周開始學(xué)習(xí)ATAC數(shù)據(jù)分析，目前只是嘗試了簡(jiǎn)單分析，新手小白，如果有小失誤歡迎留言~
1.準(zhǔn)備環(huán)境
創(chuàng)建conda環(huán)境（如果已有成熟環(huán)境不需再次創(chuàng)建）
conda create -n atac -y python=2 bwa
conda info --envs
激活環(huán)境
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.對(duì)下機(jī)數(shù)據(jù)rowdata進(jìn)行簡(jiǎn)單質(zhì)控--可選
執(zhí)行命令
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
說(shuō)明
-t: 線程
-o: 存放路徑，不用指定前綴，默認(rèn)為.fastq.gz前面的字段
*.fastq.gz:fastqc可以同時(shí)接受多個(gè)fastq.gz文件，因此采用正則表達(dá)式*表示全部
結(jié)果會(huì)產(chǎn)生兩個(gè)文件：
image
image
將HTLM文件下載到本地打開
image
注： 接頭含量太多需要過(guò)濾以及去接頭
4.使用trim_galore去除接頭，trim_galore可以自動(dòng)識(shí)別接頭
下載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 ./
產(chǎn)生新的.gz文件
image
image
5.對(duì)去除接頭后的數(shù)據(jù)質(zhì)檢
執(zhí)行命令
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
產(chǎn)生兩個(gè)文件：
image
image
將htlm文件下載并打開:
image
目前已沒有接頭序列，重復(fù)序列比較多，后面比對(duì)后需要進(jìn)一步去重
6.使用hisat2比對(duì)
最重要建立參考基因組的索引
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文件較大，因此我們直接跳過(guò)sam，使用samtools轉(zhuǎn)換為排序后的bam文件，注意【|】不要漏掉，且排序是依據(jù)reads比對(duì)到基因組的位置排的
-O: 表示輸出的bam文件
-o: 輸出的bam文件名
-t與@: 線程數(shù)
比對(duì)后產(chǎn)生：
image
7.比對(duì)后的序列去重和排序
選用sambamba來(lái)去重復(fù)
sambamba markdup -r ATAC_1-2.bam ATAC_1-2.sambamba.rmdup.bam
smarttool提取可靠的比對(duì)結(jié)果
samtools view -f 2 -q 30 -o ATAC_1-2.rmdup.q30.bamATAC_1-2.sambamba.rmdup.bam
smarttool對(duì)比對(duì)結(jié)果排序
samtools sort -O bam -@ 4 -o ./ATAC_1-2.fraw.bam ATAC_1-2.rmdup.q30.bam
產(chǎn)生的bam 文件:
image
8.計(jì)算插入片段長(zhǎng)度，F(xiàn)RiP值，IDR計(jì)算重復(fù)情況
提取片段的長(zhǎng)度
從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
提取出片段長(zhǎng)度，R代碼如下：
image
得到結(jié)果圖:
image
4
FRiP_score的計(jì)算
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.使用過(guò)濾好的數(shù)據(jù)call peak
samtools構(gòu)建index
samtools index ATAC_1-2.fraw.bam
bedtools轉(zhuǎn)換為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 進(jìn)行TSS可視化
將 bam 文件轉(zhuǎn)化為 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文件
計(jì)算reference-point---需要基因組的Refseq文件--下載或者自己創(chuàng)制
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
##########################################################################################
結(jié)果文件：
image
image
11.peaks注釋
結(jié)構(gòu)注釋----會(huì)將peak所落在基因組上的區(qū)域結(jié)構(gòu)注釋出來(lái)，比如說(shuō)啟動(dòng)子區(qū)域，UTR區(qū)域，內(nèi)含子區(qū)域等等。同時(shí)，也會(huì)將與peak最臨近的基因注釋出來(lái)，非常好用
使用gff/gtf構(gòu)建一個(gè)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:/生信學(xué)習(xí)/ATAC/ricegenome/NIP.gtf",
format="gtf") #也可以使用gtf
keytypes(txdb) #感興趣的話，可以用以下方法探索txdb都包含了什么內(nèi)容
keys(txdb)
讀入單個(gè)summits文件
peaks <- readPeakFile("D:/生信學(xué)習(xí)/ATAC/ATAC_1-2_summits.bed")
結(jié)構(gòu)注釋
peakAnno <- annotatePeak(peaks,
TxDb=txdb,
tssRegion=c(-1000, 1000))
#注釋完，進(jìn)行可視化，多種圖可供選擇
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)
最后將我們的注釋結(jié)果轉(zhuǎn)為數(shù)據(jù)框，便于查看
df <- as.data.frame(peakAnno)
將注釋到的基因提取出來(lái)（第14列），用于后續(xù)功能分析
gene <- df[,14]
head(gene)
獲得結(jié)果文件;
image
image
功能注釋-----需要物種的OrgDb注釋庫(kù)
R代碼：
#GO 富集分析
library(AnnotationHub)
hub <- AnnotationHub()
query(hub, "Oryza sativa")
nip_org <- hub[['AH96213']] ##下載目標(biāo)物種到org數(shù)據(jù)
keytypes(nip_org)#查看可使用的key
length(keys(nip_org))
##我們可以使用riceidconverter來(lái)轉(zhuǎn)換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")
################################################################
獲得結(jié)果文件：
image
image
image
12.homer尋找motif
使用homer 進(jìn)行 motif分析
conda 環(huán)境下安裝 homor環(huán)境
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
產(chǎn)生的文件：
image
將該目錄下所有文件下載到本地，打開homerResults.html 文件，即可得到如下結(jié)果展示：
image
輸出結(jié)果按照p-value排序，最后一列是一個(gè)鏈接到motif文件的超鏈接，可以從這個(gè)文件中找到包含此motif的其他序列。
在Best Match/Details 列中，HOMER將會(huì)展示與denovo motif最匹配的已知motif（該列還包含有一個(gè)名為More information 的超鏈接，
點(diǎn)擊后會(huì)出現(xiàn)以下頁(yè)面：該頁(yè)面包含motif的一些基本信息，如鏈接到motfi文件的超鏈接，且可以生成motif logo的pdf格式。
參考文章鏈接：
原文鏈接：https://blog.csdn.net/u012110870/article/details/102804623
原文鏈接： http://www.itdecent.cn/p/9aa719faa4b5
原文鏈接：https://cloud.tencent.com/developer/article/1360799
原文鏈接：[http://www.itdecent.cn/p/9a31f5f01e7b
