ATAC-seq分析实操生信技能树健明教程

放在最前面的参考链接：给学徒的ATAC-seq数据实战（附上收费视频)

数据来源的文章：

  • The landscape of accessible chromatin in mammalian preimplantation embryos. Nature 2016 Jun 30;534(7609):652-7. PMID: 27309802


    image.png
数据的GEO号：GSE66581
配置环境之软件的安装，这里首推通过conda来创建一个project专属的环境
  • 可以无脑复制下面这段代码
# https://mirrors.tuna.tsinghua.edu.cn/help/anaconda/
# https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/ 
wget https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda3-latest-Linux-x86_64.sh
bash Miniconda3-latest-Linux-x86_64.sh
source ~/.bashrc 
## 安装好conda后需要设置镜像。
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda
conda config --set show_channel_urls yes

conda  create -n atac -y   python=2 bwa
conda info --envs
source activate atac
# 可以用search先进行检索
conda search trim_galore
## 保证所有的软件都是安装在 wes 这个环境下面
conda install -y sra-tools  
conda install -y trim-galore  samtools bedtools
conda install -y deeptools homer  meme
conda install -y macs2 bowtie bowtie2 
conda install -y  multiqc 
conda install -y  sambamba
数据的下载，由于原文数据太多，这里选取了四组数据来进行练习
  • 创建文件config.sra内容如下(fastq-dump时候使用，以便将sra文件转换成fastq文件时候加上我们所需要的样品名称):
2-cell-1 SRR2927015
2-cell-2 SRR2927016
2-cell-5 SRR3545580
2-cell-4 SRR2927018
  • 创建srr.list文件（里面为我们所需要下载的数据的SRR号）
SRR2927015
SRR2927016
SRR3545580
SRR2927018
数据下载
source activate atac 
mkdir -p  ~/project/atac/
cd ~/project/atac/
mkdir {sra,raw,clean,align,peaks,motif,qc}
cd sra 
cat srr.list |while read id;do ( nohup  prefetch $id & );done
## 默认下载目录：~/ncbi/public/sra/
  • 下载完成后数据大小如下：
-rw-r--r-- 1 4.2G Nov 20  2015 SRR2927015.sra
-rw-r--r-- 1 5.5G Nov 20  2015 SRR2927016.sra
-rw-r--r-- 1 2.0G Nov 20  2015 SRR2927018.sra
-rw-r--r-- 1 7.0G May 20  2016 SRR3545580.sra

第一步：将sra文件转换成fastq文件

mv ~/ncbi/public/sra/SRR* sra/
  • 创建sh文件01_fastq-dump.sh，内容如下：
## 下面需要用循环
## cd ~/project/atac/
source activate atac
dump=fastq-dump
analysis_dir=raw
# mkdir -p $analysis_dir # 由于之前已经创建过了，所以这里就无需创建了
## 下面用到的 config.sra 文件，就是上面自行制作的。

# $fastq-dump sra/SRR2927015.sra  --gzip --split-3  -A 2-cell-1 -O clean/
cat config.sra |while read id;
do echo $id
arr=($id) # 这里可以类似看成获得矩阵
srr=${arr[1]} # 这里表示提取矩阵的第二列，即SRR号
sample=${arr[0]} # 这里表示提取矩阵的第一列，即样本名称
#  测序数据的sra转fasq
nohup $dump -A  $sample -O $analysis_dir  --gzip --split-3  sra/$srr.sra &
done
  • 然后运行, 有集群的切勿在登陆节点运行，要切换到计算节点（不知道为什么我自己的不能通过提交任务qsub来运行）
sh 01_fastq-dump.sh

第二步，测序数据的过滤

  • 手动创建一个包含fastq文件路径的config.raw文本文件，第一列随意填，充数，第二列为fastq1的路径，第二列为fastq2的路径。
1 ~/project/atac/raw/2-cell-1_1.fastq.gz ~/project/atac/raw/2-cell-1_2.fastq.gz
2 ~/project/atac/raw/2-cell-2_1.fastq.gz ~/project/atac/raw/2-cell-2_2.fastq.gz
3 ~/project/atac/raw/2-cell-4_1.fastq.gz ~/project/atac/raw/2-cell-4_2.fastq.gz
4 ~/project/atac/raw/2-cell-5_1.fastq.gz ~/project/atac/raw/2-cell-5_2.fastq.gz
  • 创建02_trim_galore.sh文件，内容如下
cd ~/project/atac/
# mkdir -p clean 
source activate atac
# trim_galore -q 25 --phred33 --length 35 -e 0.1 --stringency 4 --paired -o clean/ raw/2-cell-1_1.fastq.gz raw/2-cell-1_2.fastq.gz
cat config.raw  |while read id;
do echo $id
arr=($id)
fq2=${arr[2]}
fq1=${arr[1]}
sample=${arr[0]}
nohup  trim_galore -q 25 --phred33 --length 35 -e 0.1 --stringency 4 --paired -o  clean  $fq1   $fq2  &
done
ps -ef |grep trim

第三步，数据质量的检测

  • 创建03_fastqc_multiqc.sh文本文件，内容如下
mkdir -p qc
cd ~/JMJ705_ChIP-seq/qc
mkdir -p clean
fastqc -t 5  ../clean/*gz -o clean/
mkdir -p raw
fastqc -t 5  ../raw/*gz -o raw/
  • 使用multiqc进行合并质检结果
cd raw/
multiqc ./*zip

cd clean/
multiqc ./*zip

第四步，比对

  • 健明说：比对需要的index，看清楚物种，根据对应的软件来构建，这里直接用bowtie2进行比对和统计比对率, 需要提前下载参考基因组然后使用命令构建索引，或者直接就下载索引文件：下载小鼠参考基因组的索引和注释文件, 这里用常用的mm10
  • 下载索引文件
# 索引大小为3.2GB， 不建议自己下载基因组构建，可以直接下载索引文件，代码如下：
mkdir referece && cd reference
wget -4 -q ftp://ftp.ccb.jhu.edu/pub/data/bowtie2_indexes/mm10.zip
unzip mm10.zip
  • 解压后文件大小
-rw-r--r-- 1  848M May  3  2012 mm10.rev.1.bt2
-rw-r--r-- 1  633M May  3  2012 mm10.rev.2.bt2
-rw-r--r-- 1  848M May  2  2012 mm10.1.bt2
-rw-r--r-- 1  633M May  2  2012 mm10.2.bt2
-rw-r--r-- 1  6.0K May  2  2012 mm10.3.bt2
-rw-r--r-- 1  633M May  2  2012 mm10.4.bt2
  • 取少量数据进行比对，测试流程
zcat clean/2-cell-1_1_val_1.fq.gz |head -10000 > test_file/test1.fq
zcat clean/2-cell-1_2_val_2.fq.gz |head -10000 > test_file/test2.fq

bowtie2 -x ~/project/atac/referece/mm10 -1 test1.fq  -2 test2.fq | samtools sort -@ 5 -O bam -o test.bam
# 三种不同的去重复软件
# 这里选用sambamba来去重复
sambamba markdup -r test.bam  test.sambamba.rmdup.bam
samtools flagstat test.sambamba.rmdup.bam

samtools flagstat test.sambamba.rmdup.bam
samtools flagstat test.bam

## 接下来只保留两条reads要比对到同一条染色体(Proper paired) ，还有高质量的比对结果(Mapping quality>=30)
## 顺便过滤 线粒体reads
samtools view -h -f 2 -q 30  test.sambamba.rmdup.bam |grep -v chrM| samtools sort  -O bam  -@ 5 -o - > test.last.bam
bedtools bamtobed -i test.last.bam  > test.bed
  • 测试无错误，可进行下一步分析, 创建04_bowtie2_align_sambamba_markdup.sh文件，基本copy健明的就行，改改文件路径，内容如下：
cd ~/project/atac/align

ls ~/project/atac/clean/*_1.fq.gz > 1
ls ~/project/atac/clean/*_2.fq.gz > 2
ls ~/project/atac/clean/*_2.fq.gz |cut -d"/" -f 8|cut -d"_" -f 1  > 0 ## 这里最好自己逐步分开运行一下检测0里面的结果是否是你的sample名称
paste 0 1 2  > config.clean ## 供mapping使用的配置文件

## 相对目录需要理解
bowtie2_index=~/project/atac/referece/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]}
## 比对过程15分钟一个样本
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.biostars.org/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

  • 其中bowtie2比对加入了-X 2000 参数，是最大插入片段，宽泛的插入片段范围(10-1000bp)

  • 得到的bam文件如下：

-rw-r--r-- 1  523M Oct 12 07:15 ./2-cell-5.last.bam
-rw-r--r-- 1  899M Oct 12 07:14 ./2-cell-5.rmdup.bam
-rw-r--r-- 1  5.5G Oct 12 06:51 ./2-cell-5.raw.bam
-rw-r--r-- 1  427M Oct 12 03:23 ./2-cell-4.last.bam
-rw-r--r-- 1  586M Oct 12 03:23 ./2-cell-4.rmdup.bam
-rw-r--r-- 1  1.8G Oct 12 03:17 ./2-cell-4.raw.bam
-rw-r--r-- 1  678M Oct 12 02:22 ./2-cell-2.last.bam
-rw-r--r-- 1  1.1G Oct 12 02:20 ./2-cell-2.rmdup.bam
-rw-r--r-- 1  4.6G Oct 12 02:00 ./2-cell-2.raw.bam
-rw-r--r-- 1  490M Oct 11 23:02 ./2-cell-1.last.bam
-rw-r--r-- 1  776M Oct 11 23:01 ./2-cell-1.rmdup.bam
-rw-r--r-- 1  3.7G Oct 11 22:48 ./2-cell-1.raw.bam
  • 上述脚本的步骤都可以拆分运行，比如bam文件构建index或者转为bed的:
ls *.last.bam|xargs -i samtools index {} 
ls *.last.bam|while read id;do (bedtools bamtobed -i $id >${id%%.*}.bed) ;done
ls *.raw.bam|while read id;do (nohup bedtools bamtobed -i $id >${id%%.*}.raw.bed & ) ;done
  • 最后得到的bed文件是
-rw-r--r-- 1  254M Oct 12 07:16 ./2-cell-5.bed
-rw-r--r-- 1  203M Oct 12 03:24 ./2-cell-4.bed
-rw-r--r-- 1  338M Oct 12 02:23 ./2-cell-2.bed
-rw-r--r-- 1  237M Oct 11 23:03 ./2-cell-1.bed

第五步，使用macs2找peaks

# macs2 callpeak -t 2-cell-1.bed  -g mm --nomodel --shift -100 --extsize 200  -n 2-cell-1 --outdir ../peaks/
cd ~/project/atac/peaks/
ls *.bed | while read id ;do (macs2 callpeak -t $id  -g mm --nomodel --shift  -100 --extsize 200  -n ${id%%.*} --outdir ./) ;done

-rw-r--r-- 1  1.1M Oct 12 14:35 2-cell-5_peaks.narrowPeak
-rw-r--r-- 1  690K Oct 12 14:35 2-cell-5_summits.bed
-rw-r--r-- 1  1.2M Oct 12 14:35 2-cell-5_peaks.xls
-rw-r--r-- 1  368K Oct 12 14:35 2-cell-4_peaks.narrowPeak
-rw-r--r-- 1  418K Oct 12 14:35 2-cell-4_peaks.xls
-rw-r--r-- 1  247K Oct 12 14:35 2-cell-4_summits.bed
-rw-r--r-- 1  1.2M Oct 12 14:35 2-cell-2_peaks.narrowPeak
-rw-r--r-- 1  1.4M Oct 12 14:35 2-cell-2_peaks.xls
-rw-r--r-- 1  805K Oct 12 14:35 2-cell-2_summits.bed
-rw-r--r-- 1  634K Oct 12 14:34 2-cell-1_peaks.narrowPeak
-rw-r--r-- 1  720K Oct 12 14:34 2-cell-1_peaks.xls
-rw-r--r-- 1  425K Oct 12 14:34 2-cell-1_summits.bed

第六步，计算插入片段长度，FRiP值，IDR计算重复情况

统计indel插入长度的分布

  • 看 bam文件第9列，在R里面统计绘图 bam文件详解
    image.png
  • 提取bam文件的第九列, 创建一个config.last_bam文件，里面内容包含bam文件的名称
2-cell-1.last.bam 2-cell-1.last
2-cell-2.last.bam 2-cell-2.last
2-cell-4.last.bam 2-cell-4.last
2-cell-5.last.bam 2-cell-5.last
  • 然后创建了提取bam文件的第九列indel插入长度信息的sh文件 indel_length.sh，内容如下：
cat config.last_bam |while read id;
do
arr=($id)
sample=${arr[0]}
sample_name=${arr[1]}
samtools view $sample |awk '{print $9}'  > ${sample_name}_length.txt
done
  • 然后我们得到四个bam文件的indel插入长度信息
-rw-r--r-- 1   24M Oct 12 15:27 2-cell-5.last_length.txt
-rw-r--r-- 1   19M Oct 12 15:27 2-cell-4.last_length.txt
-rw-r--r-- 1   32M Oct 12 15:26 2-cell-2.last_length.txt
-rw-r--r-- 1   22M Oct 12 15:26 2-cell-1.last_length.txt

cmd=commandArgs(trailingOnly=TRUE); 
input=cmd[1]; output=cmd[2]; 
a=abs(as.numeric(read.table(input)[,1])); 
png(file=output);
hist(a,
main="Insertion Size distribution",
ylab="Read Count",xlab="Insert Size",
xaxt="n",
breaks=seq(0,max(a),by=10)
); 

axis(side=1,
at=seq(0,max(a),by=100),
labels=seq(0,max(a),by=100)
);

dev.off()
  • 准备一个用于R语言批量绘制indel分布的文本输入文件config.indel_length_distribution
2-cell-1.last_length.txt 2-cell-1.last_length
2-cell-2.last_length.txt 2-cell-2.last_length
2-cell-4.last_length.txt 2-cell-4.last_length
2-cell-5.last_length.txt 2-cell-5.last_length
  • 有了上面的文件就可以批量检验bam文件进行出图。创建批量运行的shell脚本indel_length_distribution.sh
cat config.indel_length_distribution  |while read id;
do
arr=($id)
input=${arr[0]}
output=${arr[1]}
Rscript indel_length_distribution.R $input $output
done
image.png

FRiP值的计算：fraction of reads in called peak regions

  • 单个举例
bedtools intersect -a 2-cell-1.bed -b 2-cell-1_peaks.narrowPeak |wc -l
148210
wc  -l 2-cell-1.bed
5105850 
# 故2-cell-1的FRiP为
148210/5105850 = 0.0292
  • 批量计算 FRiP, 创建sh文件01_FRiP.sh
cd ~/project/atac/peaks
ls *narrowPeak|while  read id;
do
echo $id
bed=$(basename $id "_peaks.narrowPeak").bed
#ls  -lh $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
  • 运行结果为
$ sh 01_FRiP.sh 
2-cell-1_peaks.narrowPeak
148210 5105850
==> FRiP value: 2.90%
2-cell-2_peaks.narrowPeak
320407 7292154
==> FRiP value: 4.39%
2-cell-4_peaks.narrowPeak
90850 4399720
==> FRiP value: 2.06%
2-cell-5_peaks.narrowPeak
258988 5466482
==> FRiP value: 4.73%
  • 健明记录：

    • Fraction of reads in peaks (FRiP) - Fraction of all mapped reads that fall into the called peak regions, i.e. usable reads in significantly enriched peaks divided by all usable reads. In general, FRiP scores correlate positively with the number of regions. (Landt et al, Genome Research Sept. 2012, 22(9): 1813–1831)

    • 文章其它指标：https://www.nature.com/articles/sdata2016109/tables/4

可以使用R包看不同peaks文件的overlap情况

  • 注意这里由于R版本是3.5.1，所以需要GCC版本要大于等于4.8，由于本服务器系统版本为4.7，所以需更改GCC版本，通过module命令调用更高的GCC版本
  • 装包不成功，暂时不管
source /public/home/software/.bashrc
module load GCC/5.4.0-2.26
source activate atac
  • 将narrowPeak文件传入到本地，使用本地R进行可视化
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/") 
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
source("http://bioconductor.org/biocLite.R") 
library('BiocInstaller')
# biocLite("ChIPpeakAnno")
# biocLite("ChIPseeker")
library(ChIPseeker)
library(ChIPpeakAnno)
setwd("E://desktop/sept/ATAC-seq_practice/find_peaks_overlaping/")
list.files('./',"*.narrowPeak")
tmp = lapply(list.files('./',"*.narrowPeak"),function(x){
  return(readPeakFile(file.path('./', x)))
  })
tmp
ol <- findOverlapsOfPeaks(tmp[[1]],tmp[[2]]) # 这里选取的是第一个文件和第二个文件，即cell.1_peak_1和cell.2_peak
png('overlapVenn.png')
makeVennDiagram(ol)
dev.off()
image.png
  • 也可以使用专业软件，IDR 来进行计算出来，同时考虑peaks间的overlap，和富集倍数的一致性
conda  create -n py3 -y   python=3 idr
conda activate py3
idr -h 
idr --samples  2-cell-1_peaks.narrowPeak 2-cell-2_peaks.narrowPeak  --plot

第七步，deeptools的可视化

cd  ~/project/atac/deeptools_result
#source activate atac # 由于原本电脑存在deeptools所以就没必要激活了
#ls  *.bam  |xargs -i samtools index {} 
ls *last.bam |while read id;do
nohup bamCoverage -p 5 --normalizeUsingRPKM -b $id -o ${id%%.*}.last.bw &
done

# cd dup 
# ls  *.bam  |xargs -i samtools index {} 
# ls *.bam |while read id;do
# nohup bamCoverage --normalizeUsing CPM -b $id -o ${id%%.*}.rm.bw & 
# done

curl 'http://genome.ucsc.edu/cgi-bin/hgTables?hgsid=646311755_P0RcOBvAQnWZSzQz2fQfBiPPSBen&boolshad.hgta_printCustomTrackHeaders=0&hgta_ctName=tb_ncbiRefSeq&hgta_ctDesc=table+browser+query+on+ncbiRefSeq&hgta_ctVis=pack&hgta_ctUrl=&fbQual=whole&fbUpBases=200&fbExonBases=0&fbIntronBases=0&fbDownBases=200&hgta_doGetBed=get+BED' >mm10.refseq.bed

查看TSS附件信号强度：创建07_deeptools_TSS.sh

## both -R and -S can accept multiple files 
mkdir -p  ~/project/atac/tss
cd   ~/project/atac/tss 
# source activate atac # 由于我这里自己系统有就没调用了
computeMatrix reference-point  --referencePoint TSS  -p 15  \
-b 10000 -a 10000    \
-R ~/project/atac/mm10_Refgene/Refseq.bed  \
-S ~/project/atac/deeptools_result/*.bw  \
--skipZeros  -o matrix1_test_TSS.gz  \
--outFileSortedRegions regions1_test_genes.bed

##     both plotHeatmap and plotProfile 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.png

    image.png

    image.png
  • 查看基因body的信号强度,创建07_deeptools_Body.sh
#source activate atac
mkdir Body
cd ~/project/atac/Body
computeMatrix scale-regions  -p 15  \
-R ~/project/atac/mm10_Refgene/Refseq.bed  \
-S ~/project/atac/deeptools_result/*.bw  \
-b 10000 -a 10000  \
--skipZeros -o matrix1_test_body.gz
# plotHeatmap -m matrix1_test_body.gz  -out ExampleHeatmap1.png
plotHeatmap -m matrix1_test_body.gz  -out test_body_Heatmap.png
plotProfile -m matrix1_test_body.gz  -out test_body_Profile.png
plotProfile -m matrix1_test_body.gz -out test_Body_Profile.pdf --plotFileFormat pdf --perGroup --dpi 720

第八步：peaks注释

options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/") 
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
source("http://bioconductor.org/biocLite.R") 
library('BiocInstaller')
biocLite("ChIPpeakAnno")
library(ChIPpeakAnno)
setwd("E://desktop/sept/ATAC-seq_practice/peaks_annotaion/")
biocLite("TxDb.Mmusculus.UCSC.mm10.knownGene")
biocLite("org.Mm.eg.db")
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
promoter <- getPromoters(TxDb=txdb, 
                         upstream=3000, downstream=3000)
files = list(cell_1_summits = "2-cell-1_summits.bed", cell_2_summits = "2-cell-2_summits.bed",
                   cell_4_summits = "2-cell-4_summits.bed", cell_5_summits = "2-cell-5_summits.bed")
peakAnno <- annotatePeak(files[[1]], # 分别改成2或者3或者4即可，分别对应四个文件
                         tssRegion=c(-3000, 3000),
                         TxDb=txdb, annoDb="org.Hs.eg.db")
plotAnnoPie(peakAnno)
image.png
  • plotAnnoBar来展示
plotAnnoBar(peakAnno)
image.png
  • vennpie来展示
vennpie(peakAnno)
image.png
  • upsetplot来展示
upsetplot(peakAnno)
upsetplot(peakAnno, vennpie=TRUE)
image.png

vennpie=TRUE
  • ChIPseeker还注释了最近的基因，peak离最近基因的距离分布是什么样子的？ChIPseeker提供了plotDistToTSS函数来画这个分布：
plotDistToTSS(peakAnno,
              title="Distribution of transcription factor-binding loci\nrelative to TSS")
image.png
  • plotAnnoBar和plotDistToTSS这两个柱状图都支持多个数据同时展示，方便比较，比如：
peakAnnoList <- lapply(files, annotatePeak, 
                       TxDb=txdb,tssRegion=c(-3000, 3000))
plotAnnoBar(peakAnnoList)
image.png
  • 批量绘制距离TSS的百分比图
plotDistToTSS(peakAnnoList)
image.png
  • ChIPseeker还提供了一个vennplot函数，比如我想看注释的最近基因在不同样本中的overlap：
genes <- lapply(peakAnnoList, function(i) 
    as.data.frame(i)$geneId)
vennplot(genes[2:4], by='Vennerable')
  • ChIPseeker还提供了一个vennplot函数，比如我想看注释的最近基因在不同样本中的overlap：
genes <- lapply(peakAnnoList, function(i) 
    as.data.frame(i)$geneId)
vennplot(genes[2:4], by='Vennerable')

使用hommer进行注释

perl ~/miniconda3/envs/atac/share/homer-4.9.1-6/configureHomer.pl  -install mm10  # 此不能再计算节点运行，需先在登陆节点运行下载
## 保证数据库下载是OK
# ls -lh  ~/miniconda3/envs/atac/share/homer-4.9.1-5/data/genomes  
# source activate atac
mkdir hommer_anno
cp peaks/*narrowPeak hommer_anno/
cd   ~/project/atac/hommer_anno  
ls *.narrowPeak |while read id;
do 
echo $id
awk '{print $4"\t"$1"\t"$2"\t"$3"\t+"}' $id >${id%%.*}.homer_peaks.tmp
annotatePeaks.pl  {id%%.*}.homer_peaks.tmp mm10  1>${id%%.*}.peakAnn.xls
  2>${id%%.*}.annLog.txt
done
  • 然后将peakAnn.xls文件导入到本地，通过excel透视功能进行可视化
  • 用到的函数COUNTIF(), SUM()
  • 这里还用到一个用来上传gif文件，生成链接的在线网页 http://thyrsi.com/
  • 操作gif动图
image.png

第九步，motif寻找及注释

  • 创建08_hommer_motif.sh 文件
# mkdir -p  ~/project/atac/motif
cd   ~/project/atac/motif
# source activate atac
ls ../peaks/*.narrowPeak |while read id;
do
file=$(basename $id )
sample=${file%%.*}
echo $sample 
awk '{print $4"\t"$1"\t"$2"\t"$3"\t+"}' $id > ${sample}.homer_peaks.tmp
nohup findMotifsGenome.pl ${sample}.homer_peaks.tmp  mm10 ${sample}_motifDir -len 8,10,12  &
done
homerResults
knownResults

第十步，差异peaks分析

  • diffbind
  • DESeq2等后续进行分析

写在最后的话，通过此流程的实践，再次学到很多新的知识点以及脚本操作。

  • 批量运行命令以及可视化
  • 简单R脚本程序的编写
最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 203,547评论 6赞 477
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 85,399评论 2赞 381
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 150,428评论 0赞 337
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 54,599评论 1赞 274
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 63,612评论 5赞 365
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 48,577评论 1赞 281
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 37,941评论 3赞 395
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 36,603评论 0赞 258
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 40,852评论 1赞 297
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,605评论 2赞 321
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,693评论 1赞 329
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 33,375评论 4赞 318
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 38,955评论 3赞 307
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,936评论 0赞 19
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,172评论 1赞 259
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 43,970评论 2赞 349
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 42,414评论 2赞 342

推荐阅读更多精彩内容

  • 诗：天地又增了一分别离
    最是那一关车门的冰冷 暗淡了漫天的繁星 道一声珍重 沙扬娜拉 扬诗守礼，终不过差 再见了，我的青春年华 压抑呼吸，...
    天马行空云飞扬阅读 181评论 0赞 2
  • 一件医疗美容案件代理词
    代 理 词 尊敬的审判长、审判员 辽宁海星律师事务所根据相关法律规定接受本案原告的委托，指派我担任原告与被告人身损...
    口无遮拦的人阅读 1,182评论 0赞 0
  • 十笑
    第一笑 秋分时节的邂逅 莞尔一笑 思绪便如同席慕蓉的诗句， 寂静，欢喜 第二笑 漫天落叶的偶遇 相视一笑 心脏跟随...
    计无寒阅读 596评论 0赞 0
  • 放假了
    4月5日，阴有小雨。今天是小长假的第一天，今天二五班去竹泉村春游去了，我们没有参加，原因是16年去过，再...
    改变脾气阅读 247评论 0赞 0