備注來源：
生信技能樹 [生信技能樹](javascript:void(0);) 2018-07-30
現(xiàn)在就可以跟我一起來一步步走通Hic數(shù)據(jù)分析流程啦本股，首先回憶一下準(zhǔn)備工作咯。
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安裝比較麻煩，我就不重復(fù)寫了走敌，大家仔細(xì)查看第3講：流程及軟件 。
很明顯，前面教程我們講到我們的Hic-pro安裝在：
source activate hic/home/zengjianming/biosoft/hicpro/bin/HiC-Pro_2.10.0/bin/HiC-Pro -h
但是使用的學(xué)員服務(wù)器沒有配置好，導(dǎo)致我的home目錄下面空間不足嚎莉，所以我挪動了data盤下面的空間，這樣就有28T的空間啦沛豌。
所以軟件挪了地方：
/data/training/zengjianming/biosoft/hicpro/bin/HiC-Pro_2.10.0/HiC-Pro/bin/HiC-Pro --help
然后我們這一講的重點是Hic-pro軟件的使用趋箩。
如果僅僅是看軟件的幫助文檔，會發(fā)現(xiàn)其實非常容易使用加派；
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個參數(shù)而已叫确，而最重要的當(dāng)然是配置文件的準(zhǔn)備咯。作者給出的示例文件不是一般的復(fù)雜芍锦；
## 記住 /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/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 bacteriacat >genome.sizeChromosome 4016942
測試數(shù)據(jù)
來自于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
復(fù)雜的是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
但是我們這里的測試數(shù)據(jù)是 新月柄桿菌 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 也是 用的內(nèi)切酶，其它還有 BamHI匈棘、EcoRI、HindIII析命、NdeI主卫、XhoI, 參考：https://www.biomart.cn/experiment/430/457/741/43956.htm
那么我需要針對特定酶和特定的參考基因組序列來設(shè)置好GENOME_FRAGMENT參數(shù)所需文件逃默。
代碼如下：
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這個參數(shù)，太尷尬：https://groups.google.com/forum/#!msg/hic-pro/PrX00iA0N1U/apzlw8ArBwAJ 也就是說每個樣本的fq文件都需要在自己獨立的文件夾下面瘩将。
因為軟件包含的步驟較多吟税，所以會比較耗時!
輸出結(jié)果解讀
所有的輸出，都在自己運行軟件的時候指定的out文件夾下面姿现。
共有bowtie_results和hic_results兩個結(jié)果文件夾肠仪。
其中bowtie_results文件夾下共有三個文件夾：bwt2、bwt2_global和bwt2_local备典，分別是序列比對結(jié)果异旧、染色體間關(guān)聯(lián)比對結(jié)果和染色體內(nèi)部關(guān)聯(lián)比對結(jié)果。其中bwt2文件夾下有一些數(shù)據(jù)統(tǒng)計結(jié)果的輸出文件提佣，如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文件夾下是各類結(jié)果數(shù)據(jù)統(tǒng)計倚喂；
matrix文件夾下分為 iced 和 raw 兩個文件夾每篷，分別標(biāo)準(zhǔn)化后的關(guān)聯(lián)矩陣和初始的關(guān)聯(lián)矩陣。
重點必須是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圖吧刑巧。
似乎是還停留在對比對結(jié)果的統(tǒng)計層面上喧兄！
