备注来源：
生信技能树 [生信技能树](javascript:void(0);) 2018-07-30
现在就可以跟我一起来一步步走通Hic数据分析流程啦，首先回忆一下准备工作咯。
conda create -n hic python=2 bowtie2source activate hicconda search hiclabconda install -y sra-tools samtools conda install -y bedtools trim-galoresource activate hicconda install -y numpy scipy matplotlib h5py cython numexpr statsmodels scikit-learn pandas cd ~ ln -s /data/training/zengjianming/ datamkdir -p /home/zengjianming/data/biosoft/hiclibcd /home/zengjianming/data/biosoft/hiclibpip install https://bitbucket.org/mirnylab/mirnylib/get/tip.tar.gzpip install https://bitbucket.org/mirnylab/hiclib/get/tip.tar.gz ## 17.7MB , 44kB/smkdir -p /home/zengjianming/data/biosoft/hicprocd /home/zengjianming/data/biosoft/hicpro git clone https://github.com/nservant/HiC-Pro.git ## 46.95 MiB | 134.00 KiB/s, donecd HiC-Pro/mkdir /home/zengjianming/data/biosoft/hicpro/binmkdir -p /data/training/zengjianming/biosoft/hicpro/bin/HiC-Pro_2.10.0mkdir -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.txtmkdir -p ~/data/project/hic/cd ~/data/project/hic/
参考基因组
mkdir -p ~/data/project/hic/refcd ~/data/project/hic/refwget ftp://ftp.ensemblgenomes.org/pub/bacteria/release-40/fasta/bacteria_20_collection/caulobacter_crescentus_na1000/dna/Caulobacter_crescentus_na1000.ASM2200v1.dna.toplevel.fa.gzgunzip Caulobacter_crescentus_na1000.ASM2200v1.dna.toplevel.fa.gzsource activate hicbowtie2-build Caulobacter_crescentus_na1000.ASM2200v1.dna.toplevel.fa bacteriacat >genome.sizeChromosome 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.gz857M 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.bedLIGATION_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/digestcd ~/data/project/hic/digestbin=/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 hiccp /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-Pronohup $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 = bacteriaGENOME_SIZE = /home/zengjianming/data/project/hic/ref/genome.sizePAIR1_EXT = SRR824846_Q20L10_1PAIR2_EXT = SRR824846_Q20L10_2GENOME_FRAGMENT = /home/zengjianming/data/project/hic/digest/bacteria.bedLIGATION_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.0Unmapped_pairs 22197 0.094Low_qual_pairs 0 0.0Unique_paired_alignments 17744193 75.357Multiple_pairs_alignments 462162 1.963Pairs_with_singleton 5318409 22.586Low_qual_singleton 0 0.0Unique_singleton_alignments 0 0.0Multiple_singleton_alignments 0 0.0Reported_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.pdf18 directories, 40 files
因为这个菌参考基因组很小，才4M左右，所以取500K，就只有8个bin了，其实意义不大。
随便秀一个pdf图吧。
image
似乎是还停留在对比对结果的统计层面上！
网友评论