最近有粉絲在我們《生信技能樹(shù)》公眾號(hào)后臺(tái)吐槽說(shuō)某公司給他們測(cè)了ATAC-seq，只拿到差異peak，想要不差異的peak居然被告之是額外分析項(xiàng)目，要加錢(qián)。各種巧立名目的費(fèi)用讓他害怕，還不如直接找公司要來(lái)fastq測(cè)序數(shù)據(jù)，找我們從頭開(kāi)始分析。
一條龍服務(wù)，一個(gè)ATAC-seq項(xiàng)目的標(biāo)準(zhǔn)分析僅收費(fèi)1600。同樣的我把這個(gè)《ATAC-seq》任務(wù)安排給了學(xué)徒，感謝學(xué)徒在這個(gè)春節(jié)假期還兢兢業(yè)業(yè)完成任務(wù)！
下面是學(xué)徒的探索
環(huán)境搭建如果是全新服務(wù)器或者全新用戶，首先需要安裝conda（最適合初學(xué)者的軟件管理解決方案）： #一路yes下去
wget https://repo./miniconda/Miniconda3-latest-Linux-x86_64.sh
bash Miniconda3-4.6.14-Linux-x86_64.sh
source ~/.bashrc
然后使用conda安裝一些軟件或者軟件環(huán)境，比如下載測(cè)序數(shù)據(jù)文件的aspera軟件環(huán)境： conda create -n download -y
conda activate download
conda install -y -c hcc aspera-cli
which ascp
## 一定要搞清楚你的軟件被conda安裝在哪
ls -lh ~/miniconda3/envs/download/etc/asperaweb_id_dsa.openssh
還有ATAC-SEQ數(shù)據(jù)分析流程的相關(guān)軟件： ## 安裝好conda后需要設(shè)置鏡像。
conda config --add channels https://mirrors.tuna./anaconda/pkgs/free
conda config --add channels https://mirrors.tuna./anaconda/cloud/conda-forge
conda config --add channels https://mirrors.tuna./anaconda/cloud/bioconda
conda config --set show_channel_urls yes
conda create -n atac -y python=2 bwa
conda info --envs
# 可以用search先進(jìn)行檢索
conda search trim_galore
## 保證所有的軟件都是安裝在 atac 這個(gè)環(huán)境下面
conda activate atac
conda install -y trim-galore bedtools deeptools homer meme macs2 bowtie bowtie2 sambamba
conda search samtools
conda install -y samtools=1.11
然后構(gòu)建工作目錄架構(gòu)： # 注意組織好自己的項(xiàng)目
mkdir -p ~/project/atac/
cd ~/project/atac/
mkdir {sra,raw,clean,align,peaks,motif,qc}
cd raw
取決于個(gè)人習(xí)慣。 實(shí)戰(zhàn)數(shù)據(jù)準(zhǔn)備參考：使用ebi數(shù)據(jù)庫(kù)直接下載fastq測(cè)序數(shù)據(jù) , 需要自行配置好，然后去EBI里面搜索到的 fq.txt 路徑文件： #一次性下載所有的 fastq.gz樣本
dsa=$HOME/miniconda3/envs/download/etc/asperaweb_id_dsa.openssh
ls -lh $dsa
# conda activate download
# 自己搭建好 download 這個(gè) conda 的小環(huán)境哦。
x=_1
y=_2
for id in {73..80}
do
ascp -QT -l 300m -P33001 -i $dsa \
era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/SRR123/0$id/SRR123031$id/SRR123031$id$x.fastq.gz .
ascp -QT -l 300m -P33001 -i $dsa \
era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/SRR123/0$id/SRR123031$id/SRR123031$id$y.fastq.gz .
done
把上面的代碼存為代碼文件：download.sh ，然后使用下面的命令放在后臺(tái)下載即可： conda activate download
nohup bash download.sh &
得到的文件如下： 1.2G 10月 9 22:07 SRR12303173_1.fastq.gz
1.2G 10月 9 22:12 SRR12303173_2.fastq.gz
1.1G 10月 9 22:17 SRR12303174_1.fastq.gz
1.2G 10月 9 22:24 SRR12303175_1.fastq.gz
1.3G 10月 9 22:30 SRR12303175_2.fastq.gz
1.1G 10月 10 09:54 SRR12303176_1.fastq.gz
1.1G 10月 10 09:56 SRR12303176_2.fastq.gz
1.4G 10月 10 10:02 SRR12303177_1.fastq.gz
1.4G 10月 10 10:05 SRR12303177_2.fastq.gz
1.2G 10月 10 10:09 SRR12303178_1.fastq.gz
1.2G 10月 10 10:13 SRR12303178_2.fastq.gz
1.2G 10月 10 10:18 SRR12303179_1.fastq.gz
1.2G 10月 10 10:21 SRR12303179_2.fastq.gz
1.3G 10月 10 10:26 SRR12303180_1.fastq.gz
1.3G 10月 10 10:35 SRR12303180_2.fastq.gz
可以看到，aspera下載的時(shí)候，中間11個(gè)小時(shí)任務(wù)終止了，是我自己重新跑了aspera下載，續(xù)起來(lái)了的。而且如果你仔細(xì)看會(huì)發(fā)現(xiàn) SRR12303174這個(gè)樣品只有R1的fq文件缺失了R2，也是需要重新單獨(dú)下載。 conda activate download
dsa=$HOME/miniconda3/envs/download/etc/asperaweb_id_dsa.openssh
id=74;y=_2
ascp -QT -l 300m -P33001 -i $dsa \
era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/SRR123/0$id/SRR123031$id/SRR123031$id$y.fastq.gz .
## SRR12303174_2.fastq.gz 100% 1094MB 87.3Mb/s 01:27
## Completed: 1120670K bytes transferred in 87 seconds
## (104410K bits/sec), in 1 file.
## 全部文件下載完畢后，使用下面的命令檢查一下fq.gz文件是否完整
gzip -t *gz
只有項(xiàng)目的fq數(shù)據(jù)全部準(zhǔn)備而且確認(rèn)無(wú)誤后才能進(jìn)行下一步！ 測(cè)序數(shù)據(jù)的質(zhì)量控制這里選擇trim_galore軟件，自動(dòng)批量運(yùn)行： mkdir -p ~/project/atac/
cd ~/project/atac/
conda activate atac
trim_galore --help
for id in {73..80}
do
nohup trim_galore -q 25 --phred33 --length 35 -e 0.1 --stringency 4 --paired \
-o clean raw/SRR123031$id*.fastq.gz &
done
得到的文件如下： 828M 10月 10 17:00 clean/SRR12303173_1_val_1.fq.gz
836M 10月 10 17:00 clean/SRR12303173_2_val_2.fq.gz
797M 10月 10 18:09 clean/SRR12303174_1_val_1.fq.gz
808M 10月 10 18:09 clean/SRR12303174_2_val_2.fq.gz
900M 10月 10 17:03 clean/SRR12303175_1_val_1.fq.gz
917M 10月 10 17:03 clean/SRR12303175_2_val_2.fq.gz
787M 10月 10 16:58 clean/SRR12303176_1_val_1.fq.gz
794M 10月 10 16:58 clean/SRR12303176_2_val_2.fq.gz
992M 10月 10 17:13 clean/SRR12303177_1_val_1.fq.gz
1006M 10月 10 17:13 clean/SRR12303177_2_val_2.fq.gz
845M 10月 10 17:01 clean/SRR12303178_1_val_1.fq.gz
855M 10月 10 17:01 clean/SRR12303178_2_val_2.fq.gz
840M 10月 10 17:01 clean/SRR12303179_1_val_1.fq.gz
847M 10月 10 17:01 clean/SRR12303179_2_val_2.fq.gz
945M 10月 10 17:06 clean/SRR12303180_1_val_1.fq.gz
963M 10月 10 17:06 clean/SRR12303180_2_val_2.fq.gz
這個(gè)過(guò)濾還是有點(diǎn)狠的，之前1.3G現(xiàn)在都小于1G了。實(shí)際上可以走fastqc+multiqc的質(zhì)控看過(guò)濾前后的具體情況。 數(shù)據(jù)比對(duì)到參考基因組1、mm10小鼠參考基因組的下載#下載
mkdir -p ~/project/atac/ref
cd ~/project/atac/ref
nohup wget http://hgdownload.soe./goldenPath/mm10/bigZips/mm10.fa.gz &
#解壓
gunzip mm10.fa.gz
2、bowtie2-build構(gòu)建參考基因組索引文件這一步會(huì)生成6個(gè)索引文件,這一步耗時(shí)比較常。也可以自行下載對(duì)應(yīng)的參考基因組索引。 conda activate atac
nohup bowtie2-build mm10.fa mm10 1>log 2>&1 &
得到的文件如下： 848M 10月 10 17:20 mm10.1.bt2
633M 10月 10 17:20 mm10.2.bt2
6.0K 10月 10 16:36 mm10.3.bt2
633M 10月 10 16:36 mm10.4.bt2
2.6G 1月 23 2020 mm10.fa
848M 10月 10 18:05 mm10.rev.1.bt2
633M 10月 10 18:05 mm10.rev.2.bt2
3、bowtie2進(jìn)行批量比對(duì)首先制作配置文件： cd ~/project/atac/align
ls $HOME/project/atac/clean/*_1.fq.gz > 1
ls $HOME/project/atac/clean/*_2.fq.gz > 2
ls $HOME/project/atac/clean/*_2.fq.gz |cut -d"/" -f 7|cut -d"_" -f 1 > 0
paste 0 1 2 > config.clean ## 供mapping使用的配置文件
然后創(chuàng)建含有如下內(nèi)容的腳本（aligh.sh）：
## 相對(duì)目錄需要理解
bowtie2_index=$HOME/project/atac/ref/mm10
## 一定要搞清楚自己的bowtie2軟件安裝在哪里，以及自己的索引文件在什么地方?。?！
#source activate atac
cat config.clean |while read id;
do echo $id
arr=($id)
fq2=${arr[2]}
fq1=${arr[1]}
sample=${arr[0]}
## 比對(duì)過(guò)程15分鐘一個(gè)樣本
bowtie2 -p 5 --very-sensitive -X 2000 -x $bowtie2_index -1 $fq1 -2 $fq2 |samtools sort -O bam -@ 5 -o - > ${sample}.raw.bam
samtools index ${sample}.raw.bam
bedtools bamtobed -i ${sample}.raw.bam > ${sample}.raw.bed
samtools flagstat ${sample}.raw.bam > ${sample}.raw.stat
# https://github.com/biod/sambamba/issues/177
sambamba markdup --overflow-list-size 600000 --tmpdir='./' -r ${sample}.raw.bam ${sample}.rmdup.bam
samtools index ${sample}.rmdup.bam
## ref:https://www./p/170294/
## Calculate %mtDNA:
mtReads=$(samtools idxstats ${sample}.rmdup.bam | grep 'chrM' | cut -f 3)
totalReads=$(samtools idxstats ${sample}.rmdup.bam | awk '{SUM += $3} END {print SUM}')
echo '==> mtDNA Content:' $(bc <<< "scale=2;100*$mtReads/$totalReads")'%'
samtools flagstat ${sample}.rmdup.bam > ${sample}.rmdup.stat
samtools view -h -f 2 -q 30 ${sample}.rmdup.bam |grep -v chrM |samtools sort -O bam -@ 5 -o - > ${sample}.last.bam
samtools index ${sample}.last.bam
samtools flagstat ${sample}.last.bam > ${sample}.last.stat
bedtools bamtobed -i ${sample}.last.bam > ${sample}.bed
done
提交腳本的代碼是： conda activate atac
nohup bash aligh.sh 1>log 2>&1 &
全部運(yùn)行完畢后輸出非常多文件。 首先看bam文件，如下： 1.1G 10月 11 15:49 SRR12303173.last.bam
1.8G 10月 10 23:01 SRR12303173.raw.bam
1.3G 10月 11 15:48 SRR12303173.rmdup.bam
823M 10月 11 16:00 SRR12303174.last.bam
1.7G 10月 11 00:24 SRR12303174.raw.bam
976M 10月 11 15:59 SRR12303174.rmdup.bam
1.4G 10月 11 16:11 SRR12303175.last.bam
2.2G 10月 11 02:26 SRR12303175.raw.bam
1.6G 10月 11 16:09 SRR12303175.rmdup.bam
1.2G 10月 11 16:23 SRR12303176.last.bam
1.8G 10月 11 03:52 SRR12303176.raw.bam
1.4G 10月 11 16:21 SRR12303176.rmdup.bam
1.8G 10月 11 16:37 SRR12303177.last.bam
2.5G 10月 11 06:16 SRR12303177.raw.bam
2.1G 10月 11 16:35 SRR12303177.rmdup.bam
1.2G 10月 11 16:50 SRR12303178.last.bam
1.9G 10月 11 07:53 SRR12303178.raw.bam
1.4G 10月 11 16:48 SRR12303178.rmdup.bam
1.2G 10月 11 17:02 SRR12303179.last.bam
1.9G 10月 11 09:35 SRR12303179.raw.bam
1.4G 10月 11 17:00 SRR12303179.rmdup.bam
1.7G 10月 11 17:16 SRR12303180.last.bam
2.4G 10月 11 11:51 SRR12303180.raw.bam
1.9G 10月 11 17:14 SRR12303180.rmdup.bam
每個(gè)樣品分別會(huì)輸出3個(gè)bam文件，測(cè)序數(shù)據(jù)比對(duì)的bam，以及去除PCR重復(fù)后的bam，以及去除線粒體reads后的bam文件。 查看log日志，發(fā)現(xiàn)這些樣本的線粒體含量是： ==> mtDNA Content: 1.81%
==> mtDNA Content: 3.72%
==> mtDNA Content: 1.88%
==> mtDNA Content: 1.98%
==> mtDNA Content: 2.16%
==> mtDNA Content: 3.78%
==> mtDNA Content: 2.11%
==> mtDNA Content: 2.17%
因?yàn)槲覀兪鞘紫热コ齈CR重復(fù)然后計(jì)算線粒體含量，其實(shí)是不準(zhǔn)確的。 比對(duì)后的bam文件的統(tǒng)計(jì)測(cè)序文庫(kù)復(fù)雜度的檢驗(yàn)一個(gè)簡(jiǎn)單的含有awk腳本的shell腳本，代碼如下： ls *.last.bam|while read id;
do
bedtools bamtobed -bedpe -i $id | \
awk 'BEGIN{OFS="\t"}{print $1,$2,$4,$6,$9,$10}' | sort | uniq -c | \
awk 'BEGIN{mt=0;m0=0;m1=0;m2=0} ($1==1){m1=m1+1} ($1==2){m2=m2+1} {m0=m0+1} {mt=mt+$1} END{m1_m2=-1.0; if(m2>0) m1_m2=m1/m2;printf "%d\t%d\t%d\t%d\t%f\t%f\t%f\n",mt,m0,m1,m2,m0/mt,m1/m0,m1_m2}' > ${id%%.*}.nodup.pbc.qc;
done
腳本制作好了后命名為： conda activate atac
nohup bash stat_qc.sh &
Library complexity measures計(jì)算結(jié)果如下，...nodup.pbc.qc文件格式為： TotalReadPairs
DistinctReadPairs
OneReadPair
TwoReadPairs
NRF=Distinct/Total
PBC1=OnePair/Distinct
PBC2=OnePair/TwoPair
針對(duì)NRF、PBC1、PBC2這幾個(gè)指標(biāo)，ENCODE官網(wǎng)提供了標(biāo)準(zhǔn). 計(jì)算結(jié)果顯示NRF、PBC1、PBC2的值都非常完美，說(shuō)明我們進(jìn)行過(guò)濾和PCR去重的bam文件質(zhì)量上沒(méi)有問(wèn)題，可以用于后續(xù)的分析。
前面的步驟是為了輸出 last.bam 的文件，需要首先轉(zhuǎn)為tagAlign，然后作為macs的輸入文件去找peaks，拿到peaks后進(jìn)行注釋。 另外，后面的步驟我們換了一個(gè)課題，但是分析內(nèi)容是一致的，我把a(bǔ)spera下載的代碼同樣的共享在這里哈： x=_1
y=_2
for id in {93,98,99}
do
ascp -QT -l 300m -P33001 -i $dir era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/SRR126/0$id/SRR126920$id/SRR126920$id$x.fastq.gz .
ascp -QT -l 300m -P33001 -i $dir era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/SRR126/0$id/SRR126920$id/SRR126920$id$y.fastq.gz .
done
五、生成tagAlign格式文件1. Convert PE BAM to tagAlign- 對(duì)于單端序列。直接用bed格式就可以；對(duì)于雙端序列，推薦用bedpe格式。這兩種格式都可以稱(chēng)之為tagAlign，可以作為macs的輸入文件。
- tagAligen格式相比bam，文件大小會(huì)小很多，更加方便文件的讀取。在轉(zhuǎn)換得到tagAlign格式之后，我們就可以很容易的將坐標(biāo)進(jìn)行偏移
nohup ls *nodup.srt.name.bam|while read id; do bedtools bamtobed -bedpe -mate1 -i $id | gzip -nc > ${id%%.*}.nodup.srt.name.bedpe.gz;done &
#含有chrM的染色體的TagAlign文件
nohup ls *.nodup.srt.name.bedpe.gz | while read id; do zcat $id | awk 'BEGIN{OFS="\t"}{printf "%s\t%s\t%s\tN\t1000\t%s\n%s\t%s\t%s\tN\t1000\t%s\n",$1,$2,$3,$9,$4,$5,$6,$10}' | gzip -nc > ${id%%.*}.nodup.srt.name.tagAlign.gz; done &
#去除chrM的染色體的TagAlign文件
nohup ls *nodup.srt.name.bedpe.gz|while read id; do zcat $id | grep -P -v "^chrM" | awk 'BEGIN{OFS="\t"}{printf "%s\t%s\t%s\tN\t1000\t%s\n%s\t%s\t%s\tN\t1000\t%s\n",$1,$2,$3,$9,$4,$5,$6,$10}' | gzip -nc > ${id%%.*}.nodup.nomit.srt.name.tagAlign.gz; done
2. Stand Cross Correlation analysis| 用于評(píng)估ATAC/Chip實(shí)驗(yàn)質(zhì)量好壞的一個(gè)重要指標(biāo) |
NREADS=25000000
nohup ls *.nodup.srt.name.bedpe.gz | while read id; do zcat $id | grep -v “chrM” | shuf -n ${NREADS} --random-source=<(openssl enc -aes-256-ctr -pass pass:$(zcat -f ${id%%.*}.nodup.srt.name.tagAlign.gz | wc -c) -nosalt </dev/zero 2>/dev/null) | awk 'BEGIN{OFS="\t"}{print $1,$2,$3,"N","1000",$9}' | gzip -nc > ${id%%.*}.nodup.nomit.srt.name.$((NREADS / 1000000)).tagAlign.gz; done &
#命令最終會(huì)生成交叉相關(guān)質(zhì)量評(píng)估文件，*.cc.qc文件中會(huì)輸出包含11列的信息，重點(diǎn)關(guān)注9-11列的信息，cc.plot.pdf文件相當(dāng)于*.cc.qc文件的可視化
nohup ls *$((NREADS / 1000000)).tagAlign.gz | while read id; do Rscript $(which run_spp.R) -c=$id -p=10 -filtchr=chrM -savp=${id%%.*}.cc.plot.pdf -out=${id%%.*}.cc.qc; done &
#質(zhì)控結(jié)果查看，主要看NSC,RSC,Quality tag三個(gè)值即輸出文件的第9列，第10列，第11列。
ls *.cc.qc|while read id; do cat $id | awk '{print $9, "\t", $10, "\t", $11}';done
 
Normalized strand cross-correlation coefficent (NSC)：NSC是最大交叉相關(guān)值除以背景交叉相關(guān)的比率(所有可能的鏈轉(zhuǎn)移的最小交叉相關(guān)值)。NSC值越大表明富集效果越好，NSC值低于1.1表明較弱的富集，小于1表示無(wú)富集。NSC值稍微低于1.05，有較低的信噪比或很少的峰，這肯能是生物學(xué)真實(shí)現(xiàn)象，比如有的因子在特定組織類(lèi)型中只有很少的結(jié)合位點(diǎn)；也可能確實(shí)是數(shù)據(jù)質(zhì)量差。 Relative strand cross-correlation coefficient (RSC)：RSC是片段長(zhǎng)度相關(guān)值減去背景相關(guān)值除以phantom-peak相關(guān)值減去背景相關(guān)值。RSC的最小值可能是0，表示無(wú)信號(hào)；富集好的實(shí)驗(yàn)RSC值大于1；低于1表示質(zhì)量低。 QualityTag: Quality tag based on thresholded RSC (codes: -2:veryLow,-1:Low,0:Medium,1:High,2:veryHigh) 查看交叉相關(guān)性質(zhì)量評(píng)估結(jié)果，主要看NSC,RSC,Quality tag三個(gè)值，這三個(gè)值分別對(duì)應(yīng)輸出文件的第9列，第10列，第11列。
六、Call Peaks1、去除線粒體基因的TagAlign格式文件進(jìn)行shift操作，輸入macs2軟件去callpeaksmooth_window=150 # default
shiftsize=$(( -$smooth_window/2 ))
pval_thresh=0.01
nohup ls *nodup.nomit.srt.name.tagAlign.gz | while read id; \
do macs2 callpeak \
-t $id -f BED -n "${id%%.*}" -g mm -p $pval_thresh \
--shift $shiftsize --extsize $smooth_window --nomodel -B --SPMR --keep-dup all --call-summits; \
done &
2、去除ENCODE列入黑名單的區(qū)域- 去除黑名單的bed文件，用于后續(xù)的peaks注釋
BLACKLIST=/home/gongyuqi/project/ATAC/mm10.blacklist.bed.gz
#*_summits.bed為macs2軟件callpeak的結(jié)果文件之一
nohup ls *_summits.bed | while read id; do bedtools intersect -a $id -b $BLACKLIST -v > ${id%%.*}_filt_blacklist.bed; done &
#查看過(guò)濾黑名單的區(qū)域前后的bed文件的peaks數(shù)
ls *summits.bed|while read id; do cat $id |wc -l >>summits.bed.txt;done
ls *summits_filt_blacklist.bed|while read id; do cat $id |wc -l >>summits_filt_blacklist.bed.txt;done
past summits.bed.txt summits_filt_blacklist.bed.txt
 - 去除黑名單的narrowPeaks文件，用于后續(xù)的IDR評(píng)估
#使用IDR需要先對(duì)MACS2的結(jié)果文件narrowPeak根據(jù)-log10(p-value)進(jìn)行排序,-log10(p-value)在第八列。
# Sort by Col8 in descending order and replace long peak names in Column 4 with Peak_<peakRank>
#*_peaks.narrowPeak為macs2軟件callpeak的結(jié)果文件之一
NPEAKS=300000
ls *_peaks.narrowPeak | while read id; do sort -k 8gr,8gr $id | awk 'BEGIN{OFS="\t"}{$4="Peak_"NR ; print $0}' | head -n ${NPEAKS} | gzip -nc > ${id%%_*}.narrowPeak.gz; done
BLACKLIST=../BLACKLIST/mm10.blacklist.bed.gz
#生成不壓縮文件
ls *.narrowPeak.gz | while read id; do bedtools intersect -v -a $id -b ${BLACKLIST} | awk 'BEGIN{OFS="\t"} {if ($5>1000) $5=1000; print $0}' | grep -P 'chr[\dXY]+[ \t]' > ${id%%.*}.narrowPeak.filt_blacklist; done
#生成壓縮文件
#ls *.narrowPeak.gz | while read id; do bedtools intersect -v -a $id -b ${BLACKLIST} | awk 'BEGIN{OFS="\t"} {if ($5>1000) $5=1000; print $0}' | grep -P 'chr[\dXY]+[ \t]' | gzip -nc > ${id%%.*}.narrowPeak.filt_blacklist.gz; done
3、Irreproducibility Discovery Rate (IDR)評(píng)估| 用于評(píng)估重復(fù)樣本間peaks一致性的重要指標(biāo) |
首先生成narrowPeak_sample.txt文件，方便后續(xù)循環(huán)處理，生成文件內(nèi)容如下： nohup cat narrowPeak_sample.txt | while read id
do
arr=(${id})
Rep1=${arr[0]}
Rep2=${arr[1]}
sample=${Rep1%%.*}_${Rep2%%.*}_idr
idr --samples $Rep1 $Rep2 --input-file-type narrowPeak -o $sample --plot
done &
- SRR12692092.filt_blacklist.narrowPeak SRR12692093.filt_blacklist.narrowPeak
- 沒(méi)有通過(guò)IDR閾值的顯示為紅色
  - BRM014-10uM_24h_wt （樣本處理情況）
- SRR12692098.filt_blacklist.narrowPeak SRR12692099.filt_blacklist.narrowPeak
- 沒(méi)有通過(guò)IDR閾值的顯示為紅色
  IDR評(píng)估會(huì)同時(shí)考慮peaks間的overlap和富集倍數(shù)的一致性。通過(guò)IDR閾值（0.05）的占比越大，說(shuō)明重復(fù)樣本間peaks一致性越好。從idr的分析結(jié)果看，我們的測(cè)試數(shù)據(jù)還可以的呢。
IDR評(píng)估相關(guān)參考資料： 重復(fù)樣本的處理——IDR 4、Fraction of reads in peaks (FRiP)評(píng)估| 反映樣本富集效果好壞的評(píng)價(jià)指標(biāo) |
#生成bed文件
nohup ls *.nodup.bam|while read id;do (bedtools bamtobed -i $id >${id%%.*}.nodup.bed) ;done &
#批量計(jì)算FRiP
ls *_summits_filt_blacklist.bed|while read id;
do
echo $id
bed=${id%%_*}.nodup.bed
Reads=$(bedtools intersect -a $bed -b $id |wc -l|awk '{print $1}')
totalReads=$(wc -l $bed|awk '{print $1}')
echo $Reads $totalReads
echo '==> FRiP value:' $(bc <<< "scale=2;100*$Reads/$totalReads")'%'
done
FRiP值在5%以上算比較好的。但也不絕對(duì)，這是個(gè)軟閾值，可以作為一個(gè)參考。  FRiP評(píng)估相關(guān)參考資料： https://www.jianshu.com/p/09e05bcd6981?utm_campaign=maleskine&utm_content=note&utm_medium=seo_notes&utm_source=recommendation 七、Peak annotation1、Feature Distributionsetwd("path to bed file")
library(ChIPpeakAnno)
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
library(org.Mm.eg.db)
library(BiocInstaller)
library(ChIPseeker)
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
files = list(DMSO_24h_wt_rep1 = "SRR12692092_summits_filt_blacklist.bed",
DMSO_24h_wt_rep2 = "SRR12692093_summits_filt_blacklist.bed",
BRM014_10uM_24h_wt_rep1 = "SRR12692098_summits_filt_blacklist.bed",
BRM014_10uM_24h_wt_rep2 = "SRR12692099_summits_filt_blacklist.bed")
#匯總所有樣本
#plotAnnoBar和plotDistToTSS這兩個(gè)柱狀圖都支持多個(gè)數(shù)據(jù)同時(shí)展示
peakAnnoList <- lapply(files, annotatePeak,
TxDb=txdb,
tssRegion=c(-3000, 3000))
plotAnnoBar(peakAnnoList,title = " Feature Distribution")
plotDistToTSS(peakAnnoList,title = " Feature Distribution relative to TSS")
#例舉單個(gè)樣本
peakAnno <- annotatePeak(files[[1]],# 分別改成2或者3或者4即可，分別對(duì)應(yīng)四個(gè)文件
tssRegion=c(-3000, 3000),
TxDb=txdb,
annoDb="org.Mm.eg.db")
plotAnnoPie(peakAnnoLipeakAnnost)
upsetplot(peakAnno, vennpie=TRUE)
2、查看peaks在全基因組上的分布#輸入文件的準(zhǔn)備
DMSO_24h_wt_rep1<-read.csv("SRR12692092_summits_filt_blacklist.csv")
DMSO_24h_wt_rep1<-DMSO_24h_wt_rep1[,-4]
DMSO_24h_wt_rep2<-read.csv("SRR12692093_summits_filt_blacklist.csv")
DMSO_24h_wt_rep2<-DMSO_24h_wt_rep2[,-4]
BRM014_10uM_24h_wt_rep1<-read.csv("SRR12692098_summits_filt_blacklist.csv")
BRM014_10uM_24h_wt_rep1<-BRM014_10uM_24h_wt_rep1[,-4]
BRM014_10uM_24h_wt_rep2<-read.csv("SRR12692099_summits_filt_blacklist.csv")
BRM014_10uM_24h_wt_rep2<-BRM014_10uM_24h_wt_rep2[,-4]
#以DMSO_24h_wt_rep1為例
set.seed(123)
circos.initializeWithIdeogram(plotType = c("axis", "labels"))
circos.track(ylim = c(0, 1), panel.fun = function(x, y) {
chr = CELL_META$sector.index
xlim = CELL_META$xlim
ylim = CELL_META$ylim
circos.rect(xlim[1], 0, xlim[2], 1)
}, track.height = 0.15, bg.border = NA, bg.col=rainbow(24))
text(0, 0, "DMSO_24h_wt_rep1", cex = 1.5)
circos.genomicDensity(DMSO_24h_wt_rep1, col = c("#000080"), track.height = 0.2)
circos.clear()
看到這樣的結(jié)果，第一反應(yīng)就是————為什么兩種處理情況下染色體開(kāi)放程度那么像??？難道我代碼有問(wèn)題??？經(jīng)過(guò)反復(fù)檢查驗(yàn)證（將一個(gè)樣本chr1上的peaks都刪掉，再次運(yùn)行上述代碼，就會(huì)發(fā)現(xiàn)顯著的改變），可以確定分析上是沒(méi)有問(wèn)題的。這兩種處理導(dǎo)致的差異可能不是很顯著。再加上20萬(wàn)+的peaks放在這個(gè)小小的circos圖上展示，有些差異會(huì)被掩蓋掉。就如在做TSS富集分析的時(shí)候，單獨(dú)看TSS前后3Kb區(qū)域，可以看到有兩個(gè)峰，但在看TSS-genebody-TSE區(qū)域是，TSS處相對(duì)微弱的那個(gè)峰就被掩蓋掉了。 
3、拿到每個(gè)樣本中peaks對(duì)應(yīng)得基因名
這一步非常重要，拿到基因名就可以根據(jù)課題需要進(jìn)行差異分析等 #以DMSO_24h_wt樣本為例
#replicate 1
peakAnno_DMSO_24h_wt_rep1 <- annotatePeak(files[[1]],
tssRegion=c(-3000, 3000),
TxDb=txdb,
annoDb="org.Mm.eg.db")
genelist_DMSO_24h_wt_rep1_uniqe<-as.data.frame(unique(peakAnno_DMSO_24h_wt_rep1@anno@elementMetadata@listData[["SYMBOL"]]))
colnames(genelist_DMSO_24h_wt_rep1_uniqe)<-"symbol"
#replicate 2
peakAnno_DMSO_24h_wt_rep2 <- annotatePeak(files[[2]],
tssRegion=c(-3000, 3000),
TxDb=txdb,
annoDb="org.Mm.eg.db")
genelist_DMSO_24h_wt_rep2_uniqe<-as.data.frame(unique(peakAnno_DMSO_24h_wt_rep2@anno@elementMetadata@listData[["SYMBOL"]]))
colnames(genelist_DMSO_24h_wt_rep2_uniqe)<-"symbol"
#重復(fù)樣本間共同的開(kāi)放基因
venn.diagram(
x=list(
DMSO_24h_wt_rep1=genelist_DMSO_24h_wt_rep1_uniqe$symbol,
DMSO_24h_wt_rep2=genelist_DMSO_24h_wt_rep2_uniqe$symbol
),
filename = "DMSO_24h_wt.png",
lty="dotted",
lwd=3,
col="transparent",
fill=c("darkorchid1","cornflowerblue"),
alpha=0.5,
label.col=c("darkorchid1","white","darkblue") ,
cex=1,
fontfamily="serif",
fontface="bold",
cat.default.pos="text",
cat.col=c("darkorchid1","darkblue"),
cat.cex=0.6,
cat.fontfamily="serif",
cat.dist=c(0.3,0.3),
cat.pos=0
)
#查看各組內(nèi)樣本間的overlapping reads：DMSO_24h_wt， BRM014_10uM_24h_wt；
#以及組間peaks的異同情況：DMSO_24h_wt vs. BRM014_10uM_24h_wt
#代碼類(lèi)似上面的，就不一一展示了
從下圖可以看出，不管是組間還是組內(nèi)，差異的peaks數(shù)目都不是很多了，這一點(diǎn)也驗(yàn)證了我們上面做的再全基因組范圍內(nèi)查看peaks的分布結(jié)果。 網(wǎng)頁(yè)工具絕對(duì)是完成不了這樣的命令行數(shù)據(jù)分析哦這個(gè)是基于Linux的ngs數(shù)據(jù)的上游處理，目前沒(méi)有成熟的網(wǎng)頁(yè)工具支持這樣的分析。其實(shí)呢，如果你有時(shí)間請(qǐng)務(wù)必學(xué)習(xí)編程基礎(chǔ)，自由自在的探索海量的公共數(shù)據(jù)，輔助你的科研，那么： 如果你沒(méi)有時(shí)間從頭開(kāi)始學(xué)編程，也可以委托專(zhuān)業(yè)的團(tuán)隊(duì)付費(fèi)拿到同樣的數(shù)據(jù)分析， 比如我們。一條龍服務(wù)，一個(gè)簡(jiǎn)單的ATAC-seq項(xiàng)目的標(biāo)準(zhǔn)分析（從fq文件到peaks的注釋?zhuān)﹥H收費(fèi)1600，而且是可以拿到全部的數(shù)據(jù)和代碼哦！ 如果TAC-seq項(xiàng)目實(shí)驗(yàn)設(shè)計(jì)比較復(fù)雜，比如多個(gè)實(shí)驗(yàn)條件多個(gè)時(shí)間點(diǎn)，需要做差異分析或者時(shí)序分析，費(fèi)用會(huì)比較高昂，請(qǐng)謹(jǐn)慎聯(lián)系我們哈！ - 需要自己讀文獻(xiàn)篩選合適的數(shù)據(jù)集
- 提供1個(gè)小時(shí)左右的一對(duì)一講解轉(zhuǎn)錄組數(shù)據(jù)處理背景知識(shí)。
如果需要委托，直接在我們《生信技能樹(shù)》公眾號(hào)留言即可，我們會(huì)安排合適的生信工程師對(duì)接具體的項(xiàng)目。
|
