本文可在 http://xuzhougeng.top/archives/hic-pro-an-optimized-and-flexible-pipeline-for-hi-c-data-processing 免費閱讀
HiC-Pro是一個高效率的Hi-C數(shù)據(jù)預處理工具,能夠應(yīng)用于dilution Hi-C, in situ Hi-C, DNase Hi-C, Micro-C, capture-C, capture Hi-C 和 HiChip 這些數(shù)據(jù)。
HiC-Pro的工作流程如下, 簡單的說就是先雙端測序各自比對鸭叙,然后進行合并,根據(jù)合并的結(jié)果篩選有效配對棠众。之后有效配對用于構(gòu)建contact maps.
安裝方法
HiC-Pro依賴于如下的軟件
- Bowtie2(>2.2.2), 用于序列比對
- Python2.7, 并安裝 Pysam, bx-python, numpy, scipy
- R, RColorBrewer + ggplot2
- samtools > 1.1
- GNU sort, 支持 -V, 按照version進行排序
為了保證環(huán)境的干凈檀训,我用conda進行了一個環(huán)境進行安裝
conda create -y -n hic-pro python=2.7 pysam bx-python numpy scipy samtools bowtie2
conda activate hic-pro
以2.11.1版本為例進行介紹砰嘁,最新的版本在https://github.com/nservant/HiC-Pro/releases檢查
wget https://github.com/nservant/HiC-Pro/archive/v2.11.1.tar.gz
tar -zxvf v2.11.1.tar.gz
cd HiC-Pro-2.11.1
因為我希望把HiC-Pro安裝到~/opt/bisofot
下,所以我需要修改當前目錄下的config-install.txt
中的PREFIX部分市殷,
PREFIX = /home/xzg/opt/bisofot
如果服務(wù)器支持任務(wù)投遞愕撰,可以修改CLUSTER_SYS部分, 設(shè)置為TORQUE, SGE, SLURM 或 LSF,
我的miniconda的安裝目錄是~/miniconda3
, 所以hic-pro環(huán)境的實際路徑是~/miniconda3/envs/hic-pro
make configure
make install
安裝結(jié)束之后被丧,/home/xzg/opt/bisofot
文件夾下就出現(xiàn)了HiC-Pro_2.11.1
盟戏,之后的軟件調(diào)用方式為
~/opt/biosfot/HiC-Pro_2.11.1/bin/HiC-Pro -h
如果沒有出現(xiàn)Error 就說明安裝成功了。
處理流程
讓我們新建一個項目文件夾甥桂,以一個測試數(shù)據(jù)集為例進行介紹柿究。
下載測試數(shù)據(jù)并解壓縮,該數(shù)據(jù)來自于Dixon et al. 2012 黄选, 使用HindIII 進行酶切
mkdir -p hic-pro && cd hic-pro
wget https://zerkalo.curie.fr/partage/HiC-Pro/HiCPro_testdata.tar.gz && tar -zxvf HiCPro_testdata.tar.gz
之后將測序結(jié)果移動或者軟連接到fastq文件夾下
mkdir -p fastq
mv test_data/* fastq
ls fastq
# dixon_2M dixon_2M_2
創(chuàng)建注釋文件
為了處理原始數(shù)據(jù)蝇摸,HiC-Pro需要三個注釋文件
- BED文件,記錄可能的酶切位點
- table文件办陷,記錄每條contig/scaffold/chromosome的長度
- bowtie2索引
其中BED文件和table文件必須要放在HiC-Pro_2.11.1/annotations
目錄下貌夕,該文件夾下已經(jīng)有了人類hg19和小鼠mm10。 我們以GRCh38為例, 介紹如何創(chuàng)建這三個注釋信息
# 切換目錄
cd ~/opt/biosfot/HiC-Pro_2.11.1/annotation
# 下載GRCh38的序列
wget ftp://ftp.ensembl.org/pub/release-97/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
gunzip Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
# BED
~/opt/biosoft/HiC-Pro_2.11.1/bin/utils/digest_genome.py -r HindIII -o GRCh38_HindIII.bed Homo_sapiens.GRCh38.dna.primary_assembly.fa
# chromosome size
seqkit fx2tab -nl Homo_sapiens.GRCh38.dna.primary_assembly.fa | awk '{print $1"\t"$2}' > GRCh38.chrom.size
# bowtie2 index
bowtie2-build --threads 20 Homo_sapiens.GRCh38.dna.primary_assembly.fa GRCh38
配置HiC-Pro
拷貝HiC-Pro的配置文件到項目文件夾下
cp ~/opt/biosfot/HiC-Pro_2.11.1/config-hicpro.txt .
修改配置文件config-hicpro.txt
# 線程數(shù)
N_CPU = 80
# bowtie2索引, 絕對路徑
BOWTIE2_IDX_PATH = /home/xzg/opt/biosfot/HiC-Pro_2.11.1/annotation
# bowtie2索引時的前綴
REFERENCE_GENOME = GRCh38
# 參考基因組各染色體長度
GENOME_SIZE = GRCh38.chrom.size
# 酶切位點
GENOME_FRAGMENT = GRCh38_HindIII.bed
LIGATION_SITE = AAGCTAGCTT
對于LIGATION_SITE民镜,不同酶切位點對應(yīng)的序列為HindIII(AAGCTAGCTT), MboI(GATCGATC) , DpnII(GATCGATC), NcoI(CCATGCATGG)啡专。
我們要修改的參數(shù)其實就是上面幾個。當然該配置文件還有許多參數(shù)可以修改制圈,具體見https://github.com/nservant/HiC-Pro/blob/master/doc/MANUAL.md
運行如下代碼们童,啟動分析項目
~/opt/biosoft/HiC-Pro_2.11.1/bin/HiC-Pro -i fastq -o results -c config-hicpro.txt
HiC-Pro會新建一個工作目錄,results, 之后會遍歷fastq目錄鲸鹦,尋找其中的fastq文件慧库,將其軟連接到results下的rawdata, 之后就開始用bowtie2比對以及后續(xù)的分析奸远。
測試數(shù)據(jù)代碼運行到
Run ICE Normalization
就中斷了劲适,可能是用的參考基因組和原來的教程(hg19)不一樣, 不過我自己的數(shù)據(jù)集是沒有問題的惧笛。
最后的結(jié)果如下:
$ tree -L 2 results
results
├── bowtie_results # 比對之后的輸出
│ ├── bwt2
│ ├── bwt2_global
│ └── bwt2_local
├── config-hicpro.txt
├── hic_results
│ ├── data # 有效配對
│ ├── matrix # contact maps
│ ├── pic # 可視化質(zhì)控信息
│ └── stats #文字版質(zhì)控信息
├── logs # 各種日志
│ ├── dixon_2M
│ └── dixon_2M_2
├── rawdata -> /home/xzg/project/Tutorial/hic-pro/fastq
└── tmp
一個關(guān)鍵的結(jié)果就是data文件里的以.validPairs結(jié)尾的文件葛菇,有7+1列甘磨,