思来想去,还是觉得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整合分析(本系列完结,内附目录)
基本步骤
- 用FastQC进行质控检测
- 用Trimmomatic进行质量过滤
- 用Bowtie2比对
- 用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
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:
- 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
- Following http://compbio.mit.edu/ChromHMM/ChromHMM_manual.pdf, binarized the peak files and run ChromHMM on histone marks and ATAC-seq
- Examine the resulting emission probability matrix (see emission.png), found that segment types E2 and E3 are more like enhancers (accessible or not)
- Extract E2+E3 as encc-enhancer.bed
- 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