ChIP-seq | ATAC-seq | 数据分析流程

思来想去,还是觉得ENCODE的流程靠谱,所以又花了快一周来调试,终于排除万难,跑成功了。【2019年12月08日】

以下是ATAC生成的结果目录:

call-align         call-call_peak_pooled  call-filter           call-idr_ppr                    call-overlap      call-pool_ta_pr2              call-spr
call-align_mito    call-call_peak_ppr1    call-frac_mito        call-idr_pr                     call-overlap_ppr  call-qc_report                call-tss_enrich
call-annot_enrich  call-call_peak_ppr2    call-fraglen_stat_pe  call-jsd                        call-overlap_pr   call-read_genome_tsv          metadata.json
call-bam2ta        call-call_peak_pr1     call-gc_bias          call-macs2_signal_track         call-pool_ta      call-reproducibility_idr
call-call_peak     call-call_peak_pr2     call-idr              call-macs2_signal_track_pooled  call-pool_ta_pr1  call-reproducibility_overlap

  

马上测试ChIP-seq:

报错:

Error: package or namespace load failed for ‘caTools’:
 package ‘caTools’ was installed by an R version with different internals; it needs to be reinstalled for use with this R version
Error: package ‘caTools’ could not be loaded

用这个conda装了R3.6.1,装了几个包,可能把原来的包给覆盖了,只能重新装一个这个pipeline了,以后自己的包就不要装在默认的conda的目录里了。

必须为这两个流程准备独立的conda,不要用conda装任何R和python的包。

ChIP-seq也终于跑成功了:

call-align               call-call_peak_ppr2  call-idr_ppr                    call-pool_ta_pr1
call-align_ctl           call-call_peak_pr1   call-jsd                        call-pool_ta_pr2
call-align_R1            call-call_peak_pr2   call-macs2_signal_track         call-qc_report
call-bam2ta              call-choose_ctl      call-macs2_signal_track_pooled  call-read_genome_tsv
call-bam2ta_ctl          call-filter          call-overlap                    call-reproducibility_overlap
call-bam2ta_no_dedup_R1  call-filter_ctl      call-overlap_ppr                call-spr
call-call_peak           call-filter_R1       call-overlap_pr                 call-xcor
call-call_peak_pooled    call-fraglen_mean    call-pool_ta                    metadata.json
call-call_peak_ppr1      call-gc_bias         call-pool_ta_ctl

  

 

核心几点:

1. 用新的conda,以前所有的conda环境都要移除,~/.bash_profile,特别是R相关的

export R_LIBS_USER=/home/-/softwares/R_lib3
export R_LIBS=/home/-/softwares/R_lib3

2. 源码里面有两个文件的grep -P要改为grep -E

3. wdl文件里的默认资源设置要改,否则会溢出  

 

 

质量控制QC:

Quality control in ChIP-seq data Using the ChIPQC package 

ChIP-seq Data quality and Analysis Pipeline - 刘小乐大佬,Cistrome project

ChiLin: a comprehensive ChIP-seq and DNase-seq quality control and analysis pipeline


 

参考链接:

一篇文章学会ChIP-seq分析(上)

一篇文章学会ChIP-seq分析(下)

第10篇:ATAC-Seq、ChIP-Seq、RNA-Seq整合分析(本系列完结,内附目录)

 

基本步骤

  1. 用FastQC进行质控检测
  2. 用Trimmomatic进行质量过滤
  3. 用Bowtie2比对
  4. 用MACS2 call peaks

 

H3K27Ac + H3K4Me1 = Active enhancers
H3K27Ac + H3K4Me3 = Active promoters
H3K4Me1 - H3K27Ac = Inactive/Poised enhancers.

 


ENCODE已经有非常成熟的pipeline了,会用就行了。【已测试,均可运行,非常高效】

ENCODE-DCC

关于整个流程,也有非常详细的介绍。

ChIP-seq Data Standards and Processing Pipeline

ATAC-seq Data Standards and Prototype Processing Pipeline

还有比这更靠谱的pipeline吗,肯定是没了。但是开始学习时肯定不推荐用这些pipeline,封装得太好了,完全学不到什么精华知识。

安装完成后,建议自己构建小测试数据来测试流程,我测试了没问题。

 

这些流程的NB之处:

  • 分析随时能够续上,Resume pipeline,除非你有HPC的root权限;
  • 良好的组织形式,报告,结果文件,能看懂结果就好
  • 考虑得非常细,每一个细节,这就是ENCODE大型项目的态度(这就是为什么不推荐自己分析,这么多细节,搞透最少要一个月) 

缺点也非常明显:新手无法掌控这些流程

 

折腾了几天后的感想:

1. 这个流程设计得太复杂了,什么平台都想囊括在内,导致什么平台都没有真正做好;

2. 最核心的工具部署测试环节没了,导致软件的部署调试非常难,表面上看是部署成功了,但是可能跑了一天就报错终止了;

3. 流程恢复非常不智能,一个小报错就足以让你一天白跑;

4. 流程交互性不够,需要一个最基本的统计,整个流程跑了多少,多少成功,多少失败,最好有一个任务监控的报告monitor;

5. 这些流程设计之初就是为了ENCODE项目,只要能在他们机器上部署成功,能批量运行就行;后面又想服务大众,可兼容性没那么简单,每个机器本地环境不同,标准化部署很难;

 

看了GitHub的issue,才发现大家都有一样的问题,这个pipeline的根本还是顶层设计有问题,累死开发者,气死使用者

如果我是领导,我会这样拯救这个项目:

1. 独立出所有的云平台pipeline,因为云平台的配置部署简单稳定,不存在疑难杂症(问题是大部分实验室根本就没有普及云平台,用不到);

2. 单独做一份针对本地的PBS和SGE用户的原始脚本pipeline,这部分就是GitHub上面大家普遍抱怨的部分,这部分责任本来就不属于流程开发者,我就给你提供一个大框架,软件安装测试你们自己去做,不要把这部分责任揽到自己头上,吃力不讨好;

不要什么都想要,又什么都没做好;具体问题具体分析,适应形势,所有的一刀切最终都是悲剧。

 

另外提一下,流程部署最重要的就是稳定性,特别是同一批数据,绝对不能用不同的流程,甚至是软件版本、参数都不能变,否则这些数据之间就没有了可比性,一个实验室一旦搭建好了分析流程,一般都要用10年。

 

必须要手动修改这个流程了,否则只有放弃。

 

环境报错【错误的定位到了以前的conda里】:

from pysam.libchtslib import *
from matplotlib import pyplot as plt

ModuleNotFoundError: No module named 'matplotlib'
ImportError: libhts.so.1: cannot open shared object file: No such file or directory  

把conda重新安装一下,以前的conda全部清理掉,留下这个sandbox的conda

qsub: would exceed queue's generic per-user limit

队列的问题,It is possible that you tried submitting the job to a queue that has a max_queued limit—a specified maximum number of jobs that each user can submit to a queue (running or waiting)—and you have already reached that limit.

 

source activate encode-chip-seq-pipeline

json配置文件里的内存不要瞎设置,picard就保持默认就行,因为它默认只给8G。 

 

不要build genome了,直接下载,chip和ATAC都一样

bash scripts/download_genome_data.sh hg19 ~/softwares/chip-seq-pipeline2/db  

经查看,两个的基因组是一模一样的,不必重复下载,用同一个数据库,test数据可以顺利跑完。

diff scripts/download_genome_data.sh ../atac-seq-pipeline/scripts/download_genome_data.sh

 

用测试数据跑有个问题,数据太小没有结果,导致后面跑不下去,并不是流程有什么问题。

比如以下这些file在小数据集中是不会生成的。

bfilt.frip.qc

bfilt.narrowPeak.bb

bfilt.narrowPeak.hammock.gz 

就会报错:

ln: failed to access ‘/home/-/project2/analysis/ATAC-seq/encode-pipeline/encc/subsample/result/atac/a1db1b2f-64b6-4bce-b51b-2d5ee2d316bb/call-idr_ppr/execution/*.bfilt.narrowPeak.bb’: No such file or directory
ln: failed to access ‘/home/-/project2/analysis/ATAC-seq/encode-pipeline/encc/subsample/result/atac/a1db1b2f-64b6-4bce-b51b-2d5ee2d316bb/call-idr_ppr/execution/*.bfilt.narrowPeak.hammock.gz*’: No such file or directory
mkdir: cannot create directory ‘/home/-/project2/analysis/ATAC-seq/encode-pipeline/encc/subsample/result/atac/a1db1b2f-64b6-4bce-b51b-2d5ee2d316bb/call-idr_ppr/execution/glob-08ed81b9c4c9ccf6c3692d9ea29b11e0’: File exists
ln: failed to access ‘/home/-/project2/analysis/ATAC-seq/encode-pipeline/encc/subsample/result/atac/a1db1b2f-64b6-4bce-b51b-2d5ee2d316bb/call-idr_ppr/execution/*.bfilt.narrowPeak.hammock.gz*’: No such file or directory
ln: failed to access ‘/home/-/project2/analysis/ATAC-seq/encode-pipeline/encc/subsample/result/atac/a1db1b2f-64b6-4bce-b51b-2d5ee2d316bb/call-idr_ppr/execution/*.frip.qc’: No such file or directory

  

如果不是在云服务器或者有root权限,就无法使用数据库来存储任务进度,失败的任务只能重头再来,没法恢复。

所以对于那些没有root权限,在HPC上跑的人来说,提前的部署测试至关重要,不要直接上大数据开跑。

 

报错:=>> PBS: job killed: ncpus 2.35 exceeded limit 1 (sum)

不是所有的任务都可以在json里定义资源的实用,最本质的资源定义在流程的atac.wdl文件里。

因为我是统一队列,所以CPU和MEM可以统一,cpu=10,mem=40000

这个需要针对自己的集群设置一下。

 

好好分析一下,为什么--db file 这种形式的数据库无法恢复任务,为什么不稳定?

 

 


 

Linux安装问题:

grep: this version of PCRE is compiled without UTF support

grep -P 改为 grep -E就可以解决了,参考链接。【这种方法没有解决根本问题,流程里还有其他代码】

所以最根本的就是直接更新PCRE:【安装也解决不了问题,所以还是每一个grep改一下吧,改两处就行】

创建一个grep.test.sh测试grep是否正常工作:

conda env list | grep -E "\batac-seq-pipeline\s"  

 其次安装pcre

source activate encode-chip-seq-pipeline
conda install -c anaconda pcre
conda deactivate

source activate encode-atac-seq-pipeline
conda install -c anaconda pcre
conda deactivate

  

老版的conda用source activate,新版的用conda activate

source activate encode-chip-seq-pipeline
source activate encode-atac-seq-pipeline
source deactivate
conda deactivate

  

  

ChIP-seq的control问题:

Guide: Getting Started with ChIP-Seq

省钱了,以后chip-seq不用做control了!用机器学习替代ChIP-seq中的control 

suinleelab/AIControl.jl

ChIP-seq 分析原理

ChIP-seq阴阳-正负对照

 

ATAC-seq(都要用绝对路径,不然找不到文件;还需要指定输出目录,不然默认会输出到home目录)

echo "caper run ~/softwares/atac-seq-pipeline/atac.wdl -i ~/project2/ATAC-seq/ENCC/test.full.json --out-dir ~/project2/ATAC-seq/ENCC/result" | qsub -V -N ATAC -q large -l nodes=1:ppn=12,walltime=84:00:00,mem=120gb

 

关于平台的选择:

这些pipeline适配了超级多的平台,各种云平台,服务器,本地机。

 

对于普通用户肯定是在本地、PBS或SGE上跑。

  • 本地跑:需要非常多的核心,我12个核都跑不起来,所以看来最少要20个核。
  • PBS:还好我也有PBS集群,设置好queue即可并行,非常简单。

  

 


 

ChIP-seq练习数据

download.list

ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP017/SRP017311/SRR620204/SRR620204.sra
ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP017/SRP017311/SRR620205/SRR620205.sra
ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP017/SRP017311/SRR620206/SRR620206.sra
ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP017/SRP017311/SRR620207/SRR620207.sra
ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP017/SRP017311/SRR620208/SRR620208.sra
ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP017/SRP017311/SRR620209/SRR620209.sra

fastq-dump

ls *sra |while read id; do ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump –split-3 $id;done

 

Please find attached two sets of ENCC enhancers. Basically,

  • "encc-enhancer.bed" is enhancers defined with H3K27ac & H3K4me1 activity
  • "encc-enhancer-atac.bed" is enhancers defined with H3K27ac & H3K4me1 activity as well as open chromatin (ATAC-seq) signal summits.

The procedure to produce these files is as follow:

  1. Process ATAC-seq and ChIP-seq data using standard ENCODE pipeline (https://github.com/ENCODE-DCC/chip-seq-pipeline https://github.com/kundajelab/atac_dnase_pipelines) => Get narrowpeaks for ATAC-seq and ChIP-seq
  2. Following http://compbio.mit.edu/ChromHMM/ChromHMM_manual.pdf, binarized the peak files and run ChromHMM on histone marks and ATAC-seq
  3. Examine the resulting emission probability matrix (see emission.png), found that segment types E2 and E3 are more like enhancers (accessible or not)
  4. Extract E2+E3 as encc-enhancer.bed
  5. To get enhancer with unified length, further overlap encc-enhaner.bed with ATAC-seq summit (summit is the base with the highest level in a peak) and use +/-500 bp around the ATAC-seq summit as enhancers. => encc-enhancer-atac.bed

 

上一篇:flutter chip标签组件


下一篇:【原创】中断子系统-ARM GPIO中断处理流程