Hi-C技術通過交聯(lián)等實驗步驟獲得空間上相連的DNA片段，即物理位置上較遠的DNA片段之間的互作信息。根據(jù)染色體內部的互作概率顯著高于染色體之間的互作概率將不同的contig或者scaffold分成不同的染色體；根據(jù)在同一條染色體上，互作概率隨著互作距離的增加而減少將同條染色體的contig或者scaffold進行排序和定向。
HiC-Pro是2015年發(fā)表的一款處理Hi-C數(shù)據(jù)用來輔助組裝基因組的工具，也是現(xiàn)在主流應用的一款工具。記錄一下自己如何使用這個軟件的，之前的安裝記錄見另一篇文章Hic-Pro安裝。
一.文件準備
需要準備參考基因組fasta文件，bowtie2建立的索引文件，含有酶切片段信息的bed文件，基因組大小信息的sizes文件，被比對的Hi-C的reads序列文件。
1.下載Hi-C Pro的測試數(shù)據(jù)及hg19參考基因組
http://nservant.github.io/HiC-Pro/QUICKSTART.html#test-dataset
參考基因組GCF_000001405.25_GRCh37.p13_genomic.fna.gz
2.建立bowtie2索引
使用bowtie2-build建立索引
bowtie2-build -f GCF_000001405.25_GRCh37.p13_genomic.fna.gzGCF_000001405.25_GRCh37.p13_genomic
# -f 指定參考基因組,空格后面接索引文件前綴，索引文件前綴和基因組名稱最好一致；
結果生成6個bt2文件。
3. bed文件，使用digest_genome.py生成酶切片段文件。
digest_genome.py腳本在HiC-Pro-3.1.0/bin/utils目錄下，正常對自己的數(shù)據(jù)使用命令行為：
python3 /HiC-Pro-3.1.0/bin/utils/digest_genome.py? hifi_canu.contigs.fasta-r dpnii -o hifi_canu.contigs.bed；
#-r 指定自己的酶切種類名稱，大小寫都可以。
酶的名稱或序列，在代碼給了如下字典，所以此處-r寫^GATC或mboi(注意全小寫字母)都可以，^為酶切位點。
測試數(shù)據(jù)已經建好了bed文件，在HiC-Pro-3.1.0/annotation目錄下，HindIII_resfrag_hg19.bed。
第一列是序列名，這里不是原始輸入的序列名，而是HiC-Pro自行指定的序列名，第二列和第三列相關，第三列是程序檢測到的酶切位點對應序列的起始位置。例如，假設采用的酶是dnpii，那么參考基因組序列上GATC位點的起始位置。
4.sizes文件，記錄基因組中每條contig/scaffold/chromosome的長度信息的文件
正常要對自己的數(shù)據(jù)使用samtools建立索引，
samtools faidx hifi_canu.contigs.fasta，
awk '{print $1 "\t" $2}' hifi_canu.contigs.fasta.fai > hifi_canu.contigs.sizes，
測試數(shù)據(jù)已建好了sizes文件，同bed文件在同一目錄下，chrom_hg19.sizes。
這個文件記錄了genome.fasta中序列長度信息。
5.整理reads目錄結構
在HiCPro的源碼中只會讀入指定目錄的子目錄的文件，源碼如下：
所以，建立reads目錄結構如下：
注意：其中reads名稱默認必須是_R1.fastq.gz和_R2.fastq.gz結尾的。
6.運行主程序HiC-Pro前修改config-hicpro.txt文件
將HiCPro安裝目錄下的config-hicpro.txt文件拷貝到腳本運行目錄下，并進行修改。通常需修改的地方有以下幾個，其中用前半部分參數(shù)非常重要，通常每次運行都要修改的。所有的參數(shù)官網(wǎng)有給出解釋，地址鏈接為：http://nservant.github.io/HiC-Pro/MANUAL.html#manual
配置文件修改:
N_CPU，給定的CPU內存數(shù)，給的越多，運行的越快（根據(jù)服務器配置）；
BOWTIE2_IDX_PATH，填寫用bowtie2對reference建立索引所在目錄，是目錄，不是文件名，也不是前綴名；
REFERENCE_GENOME，填寫比對的參考基因組前綴；
注意：REFERENCE_GENOME一定要和bowtie2建立的索引對應上；
GENOME_SIZE，前面生成的sizes文件；
GENOME_FRAGMENT，前面用digest_genome.py程序生成的bed文件；
LIGATION_SITE，酶切位點末端補平再次連接后形成的嵌合序列，例如HindIII，則為AAGCTAGCTT；如果是MboI則序列為GATCGATC；如果是dnpii,則為
剩下一些參數(shù)也很重要，要根據(jù)自己的需求進行修改，
LOGFILE， log 文件名稱。
JOB_MEM, 內存大小,內存越大，計算越快；比如100gb或者1000M。
SORT_RAM = 1000M，我曾設置過3gb和200gb，結果都報錯。
MIN_MAPQ: 最低的質量分數(shù)，用于篩選，表示低于該MAPQ值會被過濾。
BOWTIE2_GLOBAL_OPTIONS:默認GLOBAL比對設置。
BOWTIE2_LOCAL_OPTIONS:默認LOCAL比對設置。
MIN_FRAG_SIZE: 最小的理論酶切片段大小,eg. 100。
MAX_FRAG_SIZE: 最大的理論酶切片段大小,eg 100000，這范圍越大，reads數(shù)越多。
MIN_INSERT_SIZE: 最小的文庫片段大小,eg.100。
MAX_INSERT_SIZE: 最大的文庫片段大小,eg.1000，HiC建庫的插入片段長度一般在300-500bp。
BIN_SIZE:需要生成的矩陣分辨率（bp)，也就是bin的長度，比如原文件中設置的20000 40000 150000500000 1000000，也就是20kb-1Mb之間，bin的值越小，分辨率越高，處理起來占用內存空間越大。
MATRIX_FORMAT：生成矩陣的形式，可選參數(shù)為：complete、asis、upper和lower，默認值為upper，表示保留上半部分。如果后續(xù)還要做A/B compartment分析的話，需修改成complete，否則會影響PCA分析的結果。
二.運行命令主程序
time /public/home/lq/software/HiC-Pro-3.1.0/bin/HiC-Pro \
-c/public/home/lq/script/denovo/HiC-Pro/test/config-hicpro.txt \
-i /public/home/lq/hlong/denovo/hicpro/test/02.reads \
-o /public/home/lq/software/HiC-Pro-3.1.0/test/hicpro_test
運行到Combine R1/R2 alignment files 時報錯了，后換成自己的數(shù)據(jù)就沒事了。
整個過程共11步，用了6h多，這還是我的小數(shù)據(jù)的測試用的時間。
三．結果解析
所有的輸出，都在自己運行軟件的時候指定的out文件夾下面，結果目錄如下:
bowtie_results:比對結果所在目錄；
hic_results:hic矩陣及分析結果所在目錄；
logs:存放分析日志；
rawdata：鏈接了原始數(shù)據(jù)；tmp：存放中間文件。
bowtie_results目錄下共有三個文件夾：
bwt2：存放合并后的bam文件和統(tǒng)計結果；
bwt2_global：存放全局比對結果；
bwt2_local：存放局部比對結果；
hic_results目錄下面有四個文件夾：
data：存放valid pair reads及其他無效數(shù)據(jù)文件；
matrix：存放不同分辨率矩陣文件；
pic：存放統(tǒng)計分析圖片；
stats：存放統(tǒng)計表；
hic_result/data目錄中的文件:
allVaildPairs:合并后的Valid pairs數(shù)據(jù)
DEPairs:Dangling end pairs數(shù)據(jù)
DumpPairs:實際片段長度和理論片段長度不同的數(shù)據(jù)
REPairs：酶切片段重新連接的pairs
FiltePairs:基于min/max insert/fragment size過濾的pairs，MAPQ過低的pairs；SCPairs：片段自連的pairs
hic_result/matrix目錄
matrix:存放不同分辨率矩陣文件, 分為raw和iced文件，raw: 初始的關聯(lián)矩陣iced:ice校正后的矩陣,供后續(xù)分析使用。
hic_result/pic目錄
plotHiCContactRanges_Example1.pdf有效互作中各類型比例圖；
plotHiCFragmentSize_Example1.pdf有效互作的片段大小分布圖；
plotMappingPairing_Example1.pd合并后雙端比對過濾結果圖；
hic_result/stats目錄:
1.LZ-3-15-1_allValidPairs.mergestat
這個文件主要記錄的是valid pairs中去除PCR duplication后，trans比對(比對到reference中不同序列)和cis比對(比對到reference中同一條序列)的情況。
其中valid_interaction與xx.mRSstat文件中一致；valid_interaction_rmdup表示去除PCR duplication后的valid interaction。
Valid interaction rmdup = Trans interaction + Cis interaction
2. LZ-3-15-1.mpairstat
這個文件主要記錄的是reads對的情況，包括
兩端均未比對上的reads pair（Unmapped_pairs）：3980707 ???? 16.746
只有一端比對上的reads pair(Pairs_with_singleton)：7666170????? 32.249
低質量的reads pair（Low_qual_pairs）：0????? 0.0
唯一比對reads pair(Unique_paired_alignments)：4945079 ???? 20.803
Unique paired alignments用于后續(xù)分析
3. LZ-3-15-1.mRSstat
這個文件主要記錄的是過濾掉的invalid Hi-C products，包括Dangling end pairs、Religation pairs、Self Cycle pairs、Dumped pairs等，如下圖所示
4. LZ-3-15-1_R1.mmapstat和LZ-3-15-1_R2.mmapstat
它們記錄了PE reads分開比對的結果。
以R1.mmapstat文件為例：
total_R1是總的R1 reads；
mapped_R1有由兩個部分組成，分別為：
第一步?(HiCPro稱為global alignment)比對上的reads pair(即global_R1)，
第二步比對(HiCPro稱為local alignment)比對上的reads對(即local_R1)。
具體關系如下：
Total reads = Unmapped pairs + Pairs with singleton + Low qual pairs + Unique paired alignments
四.報錯信息及修改
1.在進行到 Run quality checks for all samples時報錯了，進入Makefile文件中查看腳本181行內容，發(fā)現(xiàn)是在運行make_plots.sh腳本時出現(xiàn)錯誤，
參考Github上ISSUEs上鏈接：https://github.com/nservant/HiC-Pro/issues/163
發(fā)現(xiàn)我的dixon_2M.allValidPairs是空的，該文件在/hic_results/data/dixon_2M目錄下，原因有可能是因為自己下載的參考基因組的染色體編號和作者提供的bed文件不配套。換自己的數(shù)據(jù)后就沒事了
2.參數(shù)設置
后改成3gb
五.參考
Ref1: https://wap.sciencenet.cn/blog-2970729-1182259.html
主程序HiC-Pro可以單線程來跑（不使用--parallel選項），或者采用多線程來跑(使用--parallel選項)
/PATH/TO/HiC-Pro_2.11.1/bin/HiC-Pro --input?/PATH/TO/02.reads/?--output hicpro_output --conf?/PATH/TO/config-hicpro.txt?--parallel
?Ref2: https://blog.csdn.net/qq_50637636/article/details/120827491
Ref3: http://www.itdecent.cn/p/facec96ee6ac
Ref3: Nicolas Servant, Nelle Varoquaux, Bryan R. Lajoie. et al. HiC-Pro: an optimized and flexible pipeline for Hi-C data processing. Genome Biology. 2015.
本文使用 文章同步助手 同步
