????大家好這里是有時(shí)候季更负乡，有時(shí)候月更的基因組學(xué)研究生虚吟。上一期記錄了ATAC數(shù)據(jù)從Fastq數(shù)據(jù)到bam文件的處理腳本寸认。
ATAC數(shù)據(jù)一鍵處理
基因組學(xué)研究生，公眾號：基因組學(xué)研究生一個(gè)不成熟的小腳本串慰，ATAC數(shù)據(jù)一鍵預(yù)處理從fastq到treated bam
????后來看到有同學(xué)后臺問我有沒有ATAC-seq可視化教程偏塞，所以簡單寫了寫，希望能幫到部分人~?
Ps：那位問我的同學(xué)邦鲫，由于我消息晚了幾天看到灸叼，所以回復(fù)不了了，希望這篇文章你能看到庆捺。
ATAC數(shù)據(jù)可視化
一古今、bam to peak????
假設(shè)我們已經(jīng)得到了bam文件，接下來首先使用MACS2計(jì)算的peak滔以，代碼非常簡單捉腥，如下：
macs2?callpeak?-t?bamfile1?bamfile2...?-n?output_file?-g?mm?-f?BAM?--nomodel??--shift?-100?--extsize?200二、Make Atlas?
????工具：bedtools
????為了使得數(shù)據(jù)之間具有可比性你画，我們需要制作一個(gè)矩陣抵碟，類似于基因表達(dá)矩陣(行為Gene,列為樣本)桃漾。ATAC矩陣則是行為Genome Region（即peak），列為樣本拟逮。但是由于每個(gè)樣本的開放Region不一樣撬统，所以我們需要制作一個(gè)包含所有Region的Region Atlas。
? ? MACS2將會輸出3個(gè)文件敦迄，其中包括一個(gè)后綴為“_summits.bed”的文件恋追，即peak的中心位置，我們使用這個(gè)文件制作總的Peak Atlas罚屋。
cat?A_summits.bed|awk?'BGEIN{OFS="\t"}{print?$1,$2-250,$3+250}' > A_summits_250bp.bedcat?A_summits_250bp.bed?B_summits_250bp.bed...?|?awk?'{if($2<0){print?$1?"\t"?0?"\t"?$3}?else?{print?$0}}'?|sortBed?-i?-?|bedtools?merge?-i?-?>?atlas三苦囱、計(jì)算各個(gè)樣本的Signal
????工具：deeptools
????獲得Atlas后，我們來計(jì)算每個(gè)樣本在每個(gè)Peak內(nèi)的信號值脾猛。
multiBamSummary BED-file --BED ${atlas_file} --bamfiles ${bamfile}*sort.bam -o ${out}all_count.npz --outRawCounts ${out}all_count.txt????這樣我們就得到了一個(gè)Genome Accessibility Matrix
四沿彭、可視化
????工具：R
????有了matrix，就可以做很多事情了尖滚『砹酰可以做差異分析，peak注釋漆弄，與RNA expr聯(lián)合分析等睦裳。
? ? ATAC熱圖就很好看，先簡單做個(gè)熱圖：
library(pheatmap) #?原始數(shù)據(jù)是有很多peak的撼唾，peak太多作圖很慢廉邑，我們直接去掉一些沒有變化的Region，#?假設(shè)你已經(jīng)read?table命名為countcount$std?<-?apply(count,1,sd)?# 計(jì)算標(biāo)準(zhǔn)差count.variable?<-?count[which(count$std>1),]?#直接截取標(biāo)準(zhǔn)差大于1的Regionpheatmap(count.variable[,1:6],cluster_cols?=?F,cluster_rows?=?T,show_rownames?=?F,scale?=?"row")#?可以根據(jù) pheatmap 參數(shù)調(diào)整圖片倒谷，這里選取前6列作圖????差異分析可以用DESeq2包蛛蒙，注釋可以用clusterProfiler, Chipseeker包。今天就記錄到這里吧渤愁，歡迎后臺回復(fù)和交流~
歡迎關(guān)注基因組學(xué)研究生
