補充RNA-seq流程
以前都是自己搭RNA-seq流程，雖然可以完成任務，但是數據量一多，批次多起來，就非常難管理。
既然別人提供了這么好的流程，那就要用起來，管理起來不是一般的輕松。
ENCODE-DCC/rna-seq-pipeline
安裝比較麻煩，沒有針對local的一鍵安裝，但我們可以借鑒Dockerfile文件里面的安裝方法。
也可以模仿chip-seq-pipeline2/scripts/install_conda_env.sh的安裝方法。
conda install -c anaconda ncurses
conda install -c r r-base
conda install -c conda-forge ghostscript
conda install -c anaconda zlib
conda install -c anaconda xz
conda install samtools==1.9
conda install -c bioconda star
conda install -c bioconda kallisto
conda install -c bioconda rsem
pip install pandas==0.24.2
pip install pysam==0.15.3
pip install qc-utils==19.8.1
mkdir manual_install_tools && cd manual_install_tools
#wget http://zlib.net/zlib-1.2.11.tar.gz && tar -xvf zlib-1.2.11.tar.gz
#cd zlib-1.2.11 && ./configure --prefix /home/lizhixin/softwares/rna-seq-pipeline/manual_install_tools/zlib-1.2.11 && make && make install && rm ../zlib-1.2.11.tar.gz
#wget https://tukaani.org/xz/xz-5.2.3.tar.gz && tar -xvf xz-5.2.3.tar.gz
#cd xz-5.2.3 && ./configure && make && make install && rm ../xz-5.2.3.tar.gz
#wget https://github.com/alexdobin/STAR/archive/2.5.1b.tar.gz && tar -xzf 2.5.1b.tar.gz
#cd STAR-2.5.1b && make STAR && rm ../2.5.1b.tar.gz
# wget https://github.com/pachterlab/kallisto/releases/download/v0.44.0/kallisto_linux-v0.44.0.tar.gz && tar -xzf kallisto_linux-v0.44.0.tar.gz
#git clone --branch 1.9 --single-branch https://github.com/samtools/samtools.git && \
# git clone -b 1.9 --single-branch git://github.com/samtools/htslib.git && \
# cd samtools && make && make install && cd ../ && rm -rf samtools* htslib*
#wget https://github.com/deweylab/RSEM/archive/v1.2.31.zip
# unzip v1.2.31.zip && rm v1.2.31.zip
# cd RSEM-1.2.31 && make
git clone https://github.com/ENCODE-DCC/kentutils_v385_bin_bulkrna.git
已經測試成功，等待大數據的考驗。
還需要下載REFERENCE。
kallisto index -i kallisto.hg19.self.build.idx Homo_sapiens.GRCh37.cdna.all.fa.gz
2020年05月22日
需要再把iPSC的數據跑完
發現每一個pipeline都需要准備一個獨立的conda環境，看到下面的安裝前的卸載命令了嗎！！！
安裝也需要一定的資源，需要登錄計算節點。
$ bash scripts/uninstall_conda_env.sh # uninstall it for clean-install $ bash scripts/install_conda_env.sh
在這之前需要：
https://github.com/ENCODE-DCC/atac-seq-pipeline/blob/master/docs/install_conda.md
准備一個單獨的文件夾，裝一下conda；
conda config --set auto_activate_base false
conda activate encode-atac-seq-pipeline
注意：GitHub上下載的pipeline需要做部分修改，細節參照下面。
現在管理程序需要消耗3個以上的core，所以需要投遞到稍微大一點的節點！
思來想去，還是覺得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