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的分析，敬请期待哦！

里面提到的各种文献，综述，测试数据，软件，待这个三维基因组更新完毕会统一打包发放哦，请大家不要着急。