HiC数据分析实战之Hic-pro
通过第一讲：三维基因组学习笔记，我们了解了3D基因组研究范围，然后根据我在生信技能树发布的生信工程师标准提炼出基础技能，也就是第二讲：生信基础技能 。最后提炼出了数据分析流程，并且安装好了对应的软件，也就是第3讲：流程及软件 。不过中间我还插播了一个文献解读 。
现在才正式开始数据处理实战，其中实战的测试数据，参考基因组以及对应的软件安装都是在第3讲：流程及软件 。看懂了这些准备工作，现在就可以跟我一起来一步步走通Hic数据分析流程啦，首先回忆一下准备工作咯。
conda create -n hic python=2 bowtie2
source activate hic
conda search hiclab
conda install -y sra-tools samtools
conda install -y bedtools trim-galore
source activate hic
conda install -y numpy scipy matplotlib h5py cython numexpr statsmodels scikit-learn pandas
cd ~
ln -s /data/training/zengjianming/ data
mkdir -p /home/zengjianming/data/biosoft/hiclib
cd /home/zengjianming/data/biosoft/hiclib
pip install https://bitbucket.org/mirnylab/mirnylib/get/tip.tar.gz
pip install https://bitbucket.org/mirnylab/hiclib/get/tip.tar.gz ## 17.7MB , 44kB/s
mkdir -p /home/zengjianming/data/biosoft/hicpro
cd /home/zengjianming/data/biosoft/hicpro
git clone https://github.com/nservant/HiC-Pro.git ## 46.95 MiB | 134.00 KiB/s, done
cd HiC-Pro/
mkdir /home/zengjianming/data/biosoft/hicpro/bin
mkdir -p /data/training/zengjianming/biosoft/hicpro/bin/HiC-Pro_2.10.0
mkdir -p /home/zengjianming/bin
再次了解 Hic-pro 软件
Hic-pro安装比较麻烦，我就不重复写了，大家仔细查看第3讲：流程及软件 。
很明显，前面教程我们讲到我们的Hic-pro安装在：
source activate hic
/home/zengjianming/biosoft/hicpro/bin/HiC-Pro_2.10.0/bin/HiC-Pro -h
但是使用的学员服务器没有配置好，导致我的home目录下面空间不足，所以我挪动了data盘下面的空间，这样就有28T的空间啦。
所以软件挪了地方：
/data/training/zengjianming/biosoft/hicpro/bin/HiC-Pro_2.10.0/HiC-Pro/bin/HiC-Pro --help
然后我们这一讲的重点是Hic-pro软件的使用。
如果仅仅是看软件的帮助文档，会发现其实非常容易使用；
usage : HiC-Pro -i INPUT -o OUTPUT -c CONFIG [-s ANALYSIS_STEP] [-p] [-h] [-v]
-i|--input INPUT : input data folder; Must contains a folder per sample with input files
-o|--output OUTPUT : output folder
-c|--conf CONFIG : configuration file for Hi-C processing
重点就是3个参数而已，而最重要的当然是配置文件的准备咯。作者给出的示例文件不是一般的复杂；
# 记住 /data/training/zengjianming/ 目录就是 /home/zengjianming/data 目录。
cat /data/training/zengjianming/biosoft/hicpro/bin/HiC-Pro_2.10.0/HiC-Pro/config-hicpro.txt
mkdir -p ~/data/project/hic/
cd ~/data/project/hic/
参考基因组
mkdir -p ~/data/project/hic/ref
cd ~/data/project/hic/ref
wget ftp://ftp.ensemblgenomes.org/pub/bacteria/release-40/fasta/bacteria_20_collection/caulobacter_crescentus_na1000/dna/Caulobacter_crescentus_na1000.ASM2200v1.dna.toplevel.fa.gz
gunzip Caulobacter_crescentus_na1000.ASM2200v1.dna.toplevel.fa.gz
source activate hic
bowtie2-build Caulobacter_crescentus_na1000.ASM2200v1.dna.toplevel.fa bacteria
cat >genome.size
Chromosome 4016942
测试数据
来自于Tung B. K. Le et al. Science 2013 ：https://www.ncbi.nlm.nih.gov/sra/?term=srr824846
mkdir -p ~/data/project/hic/fq/s1/
cd ~/data/project/hic/fq/s1/
858M Jul 3 16:21 SRR824846_Q20L10_1.fastq.gz
857M Jul 3 16:22 SRR824846_Q20L10_2.fastq.gz
运行Hic-pro软件
上面说过使用这个软件需要修改一个配置文件， 简单的包括：
BOWTIE2_IDX_PATH = # bowtie2建立的索引所在的路径，记住绝对路径
REFERENCE_GENOME = # bowtie2建立的索引
GENOME_SIZE = # 一个文件记录着参考基因组中每条序列的大小，格式为chr01 10000000
复杂的是GENOME_FRAGMENT =，也就是 HiC消化片段位点文件，这个文件在每个HIC流程中都需要生成，一般HiC建库的酶为hindiii(A^AGCTT ) 或者dnpiii 生成命令为 /PATH/HiC-Pro-master/bin/utils/digest_genome.py -r hindiii -o Refgenome.fasta
在软件例子里面是：
GENOME_FRAGMENT = HindIII_resfrag_hg19.bed
LIGATION_SITE = AAGCTAGCTT
但是我们这里的测试数据是 新月柄杆菌 Caulobacter crescentus , 文章里面 Le TB et al., "High-resolution mapping of the spatial organization of a bacterial chromosome.", Science, 2013 Oct 24;342(6159):731-4 提到的是限制性酶切位点NcoI 也是 用的内切酶，其它还有 BamHI、EcoRI、HindIII、NdeI、XhoI, 参考：https://www.biomart.cn/experiment/430/457/741/43956.htm
那么我需要针对特定酶和特定的参考基因组序列来设置好GENOME_FRAGMENT参数所需文件。
代码如下：
mkdir -p ~/data/project/hic/digest
cd ~/data/project/hic/digest
bin=/data/training/zengjianming/biosoft/hicpro/bin/HiC-Pro_2.10.0/HiC-Pro/bin/utils/digest_genome.py
$bin -r C^CATGG -o bacteria.bed ../ref/Caulobacter_crescentus_na1000.ASM2200v1.dna.toplevel.fa
得到的 bacteria.bed 文件如下；
Chromosome 0 3188 HIC_Chromosome_1 0 +
Chromosome 3188 3551 HIC_Chromosome_2 0 +
Chromosome 3551 3737 HIC_Chromosome_3 0 +
Chromosome 3737 4500 HIC_Chromosome_4 0 +
Chromosome 4500 5578 HIC_Chromosome_5 0 +
Chromosome 5578 10834 HIC_Chromosome_6 0 +
Chromosome 10834 10840 HIC_Chromosome_7 0 +
Chromosome 10840 19183 HIC_Chromosome_8 0 +
Chromosome 19183 19381 HIC_Chromosome_9 0 +
Chromosome 19381 21967 HIC_Chromosome_10 0 +
然后修改配置文件后即可成功运行hicpro软件啦。
cd ~/data/project/hic/
source activate hic
cp /data/training/zengjianming/biosoft/hicpro/bin/HiC-Pro_2.10.0/HiC-Pro/config-hicpro.txt ./
bin=/data/training/zengjianming/biosoft/hicpro/bin/HiC-Pro_2.10.0/HiC-Pro/bin/HiC-Pro
nohup $bin -i fq -o out -c config-hicpro.txt 1> run.log 2>&1 &
这个配置文件问题简单，重要的地方如下：
BOWTIE2_IDX_PATH = /home/zengjianming/data/project/hic/ref/
REFERENCE_GENOME = bacteria
GENOME_SIZE = /home/zengjianming/data/project/hic/ref/genome.size
PAIR1_EXT = SRR824846_Q20L10_1
PAIR2_EXT = SRR824846_Q20L10_2
GENOME_FRAGMENT = /home/zengjianming/data/project/hic/digest/bacteria.bed
LIGATION_SITE = CCATGG
但是 我遇到的报错居然是在 input的fastq文件的folder这个参数，太尴尬：https://groups.google.com/forum/#!msg/hic-pro/PrX00iA0N1U/apzlw8ArBwAJ 也就是说每个样本的fq文件都需要在自己独立的文件夹下面。
因为软件包含的步骤较多，所以会比较耗时!
输出结果解读
所有的输出，都在自己运行软件的时候指定的out文件夹下面。
共有bowtie_results和hic_results两个结果文件夹。
其中bowtie_results文件夹下共有三个文件夹：bwt2、bwt2_global和bwt2_local，分别是序列比对结果、染色体间关联比对结果和染色体内部关联比对结果。其中bwt2文件夹下有一些数据统计结果的输出文件，如mpairstat文件（如下）。
Total_pairs_processed 23546961 100.0
Unmapped_pairs 22197 0.094
Low_qual_pairs 0 0.0
Unique_paired_alignments 17744193 75.357
Multiple_pairs_alignments 462162 1.963
Pairs_with_singleton 5318409 22.586
Low_qual_singleton 0 0.0
Unique_singleton_alignments 0 0.0
Multiple_singleton_alignments 0 0.0
Reported_pairs 17744193 75.357
hic_results文件夹下共有data、matrix以及pic三个文件夹。
data文件夹下是比对上的有效序列对，文件末尾为_allValidPairs；
pic文件夹下是各类结果数据统计；
matrix文件夹下分为 iced 和 raw 两个文件夹，分别标准化后的关联矩阵和初始的关联矩阵。
重点必须是hic_results文件，如下：
hic_results/
|-- [ 24] data
| `-- [4.0K] s1
| |-- [504M] _bacteria.bwt2pairs.DEPairs
| |-- [ 93K] _bacteria.bwt2pairs.DumpPairs
| |-- [ 83M] _bacteria.bwt2pairs.REPairs
| |-- [ 303] _bacteria.bwt2pairs.RSstat
| |-- [ 73M] _bacteria.bwt2pairs.SCPairs
| |-- [ 0] _bacteria.bwt2pairs.SinglePairs
| |-- [1.2G] _bacteria.bwt2pairs.validPairs
| |-- [1.2G] s1_allValidPairs
| |-- [ 150] s1_allValidPairs.mergestat
| `-- [ 486] s1.mRSstat
|-- [ 24] matrix
| `-- [ 41] s1
| |-- [ 99] iced
| | |-- [ 85] 1000000
| | | |-- [ 275] s1_1000000_iced.matrix
| | | `-- [ 125] s1_1000000_iced.matrix.biases
| | |-- [ 83] 150000
| | | |-- [6.6K] s1_150000_iced.matrix
| | | `-- [ 675] s1_150000_iced.matrix.biases
| | |-- [ 81] 20000
| | | |-- [345K] s1_20000_iced.matrix
| | | `-- [4.9K] s1_20000_iced.matrix.biases
| | |-- [ 81] 40000
| | | |-- [ 86K] s1_40000_iced.matrix
| | | `-- [2.4K] s1_40000_iced.matrix.biases
| | `-- [ 83] 500000
| | |-- [ 792] s1_500000_iced.matrix
| | `-- [ 225] s1_500000_iced.matrix.biases
| `-- [ 99] raw
| |-- [ 99] 1000000
| | |-- [ 139] s1_1000000_abs.bed
| | |-- [ 165] s1_1000000.matrix
| | `-- [ 18] s1_1000000_ord.bed -> s1_1000000_abs.bed
| |-- [ 96] 150000
| | |-- [ 783] s1_150000_abs.bed
| | |-- [4.1K] s1_150000.matrix
| | `-- [ 17] s1_150000_ord.bed -> s1_150000_abs.bed
| |-- [ 93] 20000
| | |-- [5.9K] s1_20000_abs.bed
| | |-- [216K] s1_20000.matrix
| | `-- [ 16] s1_20000_ord.bed -> s1_20000_abs.bed
| |-- [ 93] 40000
| | |-- [2.9K] s1_40000_abs.bed
| | |-- [ 53K] s1_40000.matrix
| | `-- [ 16] s1_40000_ord.bed -> s1_40000_abs.bed
| `-- [ 96] 500000
| |-- [ 253] s1_500000_abs.bed
| |-- [ 469] s1_500000.matrix
| `-- [ 17] s1_500000_ord.bed -> s1_500000_abs.bed
`-- [ 24] pic
`-- [ 188] s1
|-- [4.9K] plotHiCContactRanges_s1.pdf
|-- [5.2K] plotHiCFragment_s1.pdf
|-- [6.3K] plotHiCFragmentSize_s1.pdf
|-- [5.0K] plotMappingPairing_s1.pdf
`-- [5.1K] plotMapping_s1.pdf
18 directories, 40 files
因为这个菌参考基因组很小，才4M左右，所以取500K，就只有8个bin了，其实意义不大。
随便秀一个pdf图吧。
似乎是还停留在对比对结果的统计层面上！
后续再更新结果的生物学意义解读，以及后续的compartment以及TAD和loop的分析，敬请期待哦！
里面提到的各种文献，综述，测试数据，软件，待这个三维基因组更新完毕会统一打包发放哦，请大家不要着急。