Hi-C技術(shù)通過(guò)交聯(lián)等實(shí)驗(yàn)步驟獲得空間上相連的DNA片段企垦，即物理位置上較遠(yuǎn)的DNA片段之間的互作信息勘畔。根據(jù)染色體內(nèi)部的互作概率顯著高于染色體之間的互作概率將不同的contig或者scaffold分成不同的染色體琐馆；根據(jù)在同一條染色體上慰枕，互作概率隨著互作距離的增加而減少將同條染色體的contig或者scaffold進(jìn)行排序和定向台汇。
HiC-Pro是2015年發(fā)表的一款處理Hi-C數(shù)據(jù)用來(lái)輔助組裝基因組的工具，也是現(xiàn)在主流應(yīng)用的一款工具稼钩。記錄一下自己如何使用這個(gè)軟件的顾稀，之前的安裝記錄見(jiàn)另一篇文章Hic-Pro安裝。
一.文件準(zhǔn)備
需要準(zhǔn)備參考基因組fasta文件坝撑，bowtie2建立的索引文件静秆，含有酶切片段信息的bed文件，基因組大小信息的sizes文件巡李，被比對(duì)的Hi-C的reads序列文件抚笔。
1.下載Hi-C Pro的測(cè)試數(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 指定參考基因組,空格后面接索引文件前綴，索引文件前綴和基因組名稱(chēng)最好一致侨拦；
結(jié)果生成6個(gè)bt2文件殊橙。
3. bed文件，使用digest_genome.py生成酶切片段文件阳谍。
digest_genome.py腳本在HiC-Pro-3.1.0/bin/utils目錄下蛀柴，正常對(duì)自己的數(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 指定自己的酶切種類(lèi)名稱(chēng)矫夯，大小寫(xiě)都可以。
酶的名稱(chēng)或序列吊洼，在代碼給了如下字典训貌，所以此處-r寫(xiě)^GATC或mboi(注意全小寫(xiě)字母)都可以，^為酶切位點(diǎn)。
測(cè)試數(shù)據(jù)已經(jīng)建好了bed文件递沪，在HiC-Pro-3.1.0/annotation目錄下豺鼻，HindIII_resfrag_hg19.bed。
第一列是序列名款慨，這里不是原始輸入的序列名儒飒，而是HiC-Pro自行指定的序列名，第二列和第三列相關(guān)檩奠，第三列是程序檢測(cè)到的酶切位點(diǎn)對(duì)應(yīng)序列的起始位置桩了。例如，假設(shè)采用的酶是dnpii埠戳，那么參考基因組序列上GATC位點(diǎn)的起始位置井誉。
4.sizes文件，記錄基因組中每條contig/scaffold/chromosome的長(zhǎng)度信息的文件
正常要對(duì)自己的數(shù)據(jù)使用samtools建立索引整胃，
samtools faidx hifi_canu.contigs.fasta颗圣，
awk '{print $1 "\t" $2}' hifi_canu.contigs.fasta.fai > hifi_canu.contigs.sizes，
測(cè)試數(shù)據(jù)已建好了sizes文件屁使，同bed文件在同一目錄下在岂，chrom_hg19.sizes。
這個(gè)文件記錄了genome.fasta中序列長(zhǎng)度信息蛮寂。
5.整理reads目錄結(jié)構(gòu)
在HiCPro的源碼中只會(huì)讀入指定目錄的子目錄的文件洁段，源碼如下：
所以，建立reads目錄結(jié)構(gòu)如下：
注意：其中reads名稱(chēng)默認(rèn)必須是_R1.fastq.gz和_R2.fastq.gz結(jié)尾的共郭。
6.運(yùn)行主程序HiC-Pro前修改config-hicpro.txt文件
將HiCPro安裝目錄下的config-hicpro.txt文件拷貝到腳本運(yùn)行目錄下祠丝，并進(jìn)行修改。通常需修改的地方有以下幾個(gè)除嘹，其中用前半部分參數(shù)非常重要写半，通常每次運(yùn)行都要修改的。所有的參數(shù)官網(wǎng)有給出解釋?zhuān)刂锋溄訛椋篽ttp://nservant.github.io/HiC-Pro/MANUAL.html#manual
配置文件修改:
N_CPU尉咕，給定的CPU內(nèi)存數(shù)叠蝇，給的越多，運(yùn)行的越快（根據(jù)服務(wù)器配置）年缎；
BOWTIE2_IDX_PATH悔捶，填寫(xiě)用bowtie2對(duì)reference建立索引所在目錄，是目錄单芜，不是文件名蜕该，也不是前綴名；
REFERENCE_GENOME洲鸠，填寫(xiě)比對(duì)的參考基因組前綴堂淡；
注意：REFERENCE_GENOME一定要和bowtie2建立的索引對(duì)應(yīng)上馋缅；
GENOME_SIZE，前面生成的sizes文件绢淀；
GENOME_FRAGMENT萤悴，前面用digest_genome.py程序生成的bed文件；
LIGATION_SITE皆的，酶切位點(diǎn)末端補(bǔ)平再次連接后形成的嵌合序列覆履，例如HindIII，則為AAGCTAGCTT费薄；如果是MboI則序列為GATCGATC内狗；如果是dnpii,則為
剩下一些參數(shù)也很重要，要根據(jù)自己的需求進(jìn)行修改义锥，
LOGFILE柳沙， log 文件名稱(chēng)。
JOB_MEM, 內(nèi)存大小,內(nèi)存越大拌倍，計(jì)算越快赂鲤；比如100gb或者1000M。
SORT_RAM = 1000M柱恤，我曾設(shè)置過(guò)3gb和200gb数初，結(jié)果都報(bào)錯(cuò)。
MIN_MAPQ: 最低的質(zhì)量分?jǐn)?shù)梗顺，用于篩選泡孩，表示低于該MAPQ值會(huì)被過(guò)濾。
BOWTIE2_GLOBAL_OPTIONS:默認(rèn)GLOBAL比對(duì)設(shè)置寺谤。
BOWTIE2_LOCAL_OPTIONS:默認(rèn)LOCAL比對(duì)設(shè)置仑鸥。
MIN_FRAG_SIZE: 最小的理論酶切片段大小,eg. 100。
MAX_FRAG_SIZE: 最大的理論酶切片段大小,eg 100000变屁，這范圍越大眼俊，reads數(shù)越多。
MIN_INSERT_SIZE: 最小的文庫(kù)片段大小,eg.100粟关。
MAX_INSERT_SIZE: 最大的文庫(kù)片段大小,eg.1000疮胖，HiC建庫(kù)的插入片段長(zhǎng)度一般在300-500bp。
BIN_SIZE:需要生成的矩陣分辨率（bp)闷板，也就是bin的長(zhǎng)度澎灸，比如原文件中設(shè)置的20000 40000 150000500000 1000000，也就是20kb-1Mb之間遮晚，bin的值越小性昭，分辨率越高，處理起來(lái)占用內(nèi)存空間越大鹏漆。
MATRIX_FORMAT：生成矩陣的形式巩梢，可選參數(shù)為：complete创泄、asis艺玲、upper和lower括蝠，默認(rèn)值為upper，表示保留上半部分饭聚。如果后續(xù)還要做A/B compartment分析的話(huà)忌警，需修改成complete，否則會(huì)影響PCA分析的結(jié)果秒梳。
二.運(yùn)行命令主程序
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
運(yùn)行到Combine R1/R2 alignment files 時(shí)報(bào)錯(cuò)了法绵，后換成自己的數(shù)據(jù)就沒(méi)事了。
整個(gè)過(guò)程共11步酪碘，用了6h多朋譬，這還是我的小數(shù)據(jù)的測(cè)試用的時(shí)間。
三．結(jié)果解析
所有的輸出兴垦，都在自己運(yùn)行軟件的時(shí)候指定的out文件夾下面徙赢，結(jié)果目錄如下:
bowtie_results:比對(duì)結(jié)果所在目錄；
hic_results:hic矩陣及分析結(jié)果所在目錄探越；
logs:存放分析日志狡赐；
rawdata：鏈接了原始數(shù)據(jù)；tmp：存放中間文件钦幔。
bowtie_results目錄下共有三個(gè)文件夾：
bwt2：存放合并后的bam文件和統(tǒng)計(jì)結(jié)果枕屉；
bwt2_global：存放全局比對(duì)結(jié)果；
bwt2_local：存放局部比對(duì)結(jié)果鲤氢；
hic_results目錄下面有四個(gè)文件夾：
data：存放valid pair reads及其他無(wú)效數(shù)據(jù)文件搀擂；
matrix：存放不同分辨率矩陣文件；
pic：存放統(tǒng)計(jì)分析圖片卷玉；
stats：存放統(tǒng)計(jì)表哥倔；
hic_result/data目錄中的文件:
allVaildPairs:合并后的Valid pairs數(shù)據(jù)
DEPairs:Dangling end pairs數(shù)據(jù)
DumpPairs:實(shí)際片段長(zhǎng)度和理論片段長(zhǎng)度不同的數(shù)據(jù)
REPairs：酶切片段重新連接的pairs
FiltePairs:基于min/max insert/fragment size過(guò)濾的pairs，MAPQ過(guò)低的pairs揍庄；SCPairs：片段自連的pairs
hic_result/matrix目錄
matrix:存放不同分辨率矩陣文件, 分為raw和iced文件咆蒿，raw: 初始的關(guān)聯(lián)矩陣iced:ice校正后的矩陣,供后續(xù)分析使用。
hic_result/pic目錄
plotHiCContactRanges_Example1.pdf有效互作中各類(lèi)型比例圖蚂子；
plotHiCFragmentSize_Example1.pdf有效互作的片段大小分布圖沃测；
plotMappingPairing_Example1.pd合并后雙端比對(duì)過(guò)濾結(jié)果圖；
hic_result/stats目錄:
1.LZ-3-15-1_allValidPairs.mergestat
這個(gè)文件主要記錄的是valid pairs中去除PCR duplication后食茎，trans比對(duì)(比對(duì)到reference中不同序列)和cis比對(duì)(比對(duì)到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
這個(gè)文件主要記錄的是reads對(duì)的情況附迷，包括
兩端均未比對(duì)上的reads pair（Unmapped_pairs）：3980707 ???? 16.746
只有一端比對(duì)上的reads pair(Pairs_with_singleton)：7666170????? 32.249
低質(zhì)量的reads pair（Low_qual_pairs）：0????? 0.0
唯一比對(duì)reads pair(Unique_paired_alignments)：4945079 ???? 20.803
Unique paired alignments用于后續(xù)分析
3. LZ-3-15-1.mRSstat
這個(gè)文件主要記錄的是過(guò)濾掉的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分開(kāi)比對(duì)的結(jié)果稻据。
以R1.mmapstat文件為例：
total_R1是總的R1 reads艾猜；
mapped_R1有由兩個(gè)部分組成，分別為：
第一步?(HiCPro稱(chēng)為global alignment)比對(duì)上的reads pair(即global_R1)捻悯，
第二步比對(duì)(HiCPro稱(chēng)為local alignment)比對(duì)上的reads對(duì)(即local_R1)匆赃。
具體關(guān)系如下：
Total reads = Unmapped pairs + Pairs with singleton + Low qual pairs + Unique paired alignments
四.報(bào)錯(cuò)信息及修改
1.在進(jìn)行到 Run quality checks for all samples時(shí)報(bào)錯(cuò)了，進(jìn)入Makefile文件中查看腳本181行內(nèi)容今缚，發(fā)現(xiàn)是在運(yùn)行make_plots.sh腳本時(shí)出現(xiàn)錯(cuò)誤算柳，
參考Github上ISSUEs上鏈接：https://github.com/nservant/HiC-Pro/issues/163
發(fā)現(xiàn)我的dixon_2M.allValidPairs是空的，該文件在/hic_results/data/dixon_2M目錄下姓言，原因有可能是因?yàn)樽约合螺d的參考基因組的染色體編號(hào)和作者提供的bed文件不配套瞬项。換自己的數(shù)據(jù)后就沒(méi)事了
2.參數(shù)設(shè)置
后改成3gb
五.參考
Ref1: https://wap.sciencenet.cn/blog-2970729-1182259.html
主程序HiC-Pro可以單線(xiàn)程來(lái)跑（不使用--parallel選項(xiàng)），或者采用多線(xiàn)程來(lái)跑(使用--parallel選項(xiàng))
/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.reibang.com/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.
本文使用 文章同步助手 同步
