个人hic分析流程搭建1—使用hic-pro

个人hic分析流程搭建1

PDX-HiC_processingScripts+PDXHiC_supplemental+TNBC-project参考

1，数据处理前置任务

（1）安装hic-pro软件：

见我个人博客https://blog.csdn.net/weixin_62528784/article/details/142055404，上游分析也可以使用Juicer

（2）获取上游hic fastq数据：

从公共数据库GEO获取数据，见脚本说明https://github.com/MaybeBio/HiC_pipeline/blob/main/upstream/prefetch_fasterq.sh

（3）依据hic-pro软件需求获取3个注释文件：

首先获取参考基因组fa文件，hg38或者hg19都可以，按照自己研究需要；
此处使用hg19为例，建议到https://hgdownload.soe.ucsc.edu/goldenPath/hg19/下载最新的版本；
因为官方要求全程处理的chr信息和bowtie处理的index索引文件一致，所以我强烈建议只提取chr1-22+XY的fa数据作为参考基因组，其他非常规染色体数据全部舍弃，

具体实现方法有很多，seqkit，samtools faidx等，

①含有酶切片段信息的bed文件:

重点在于获取限制酶的信息，如果实验有提供信息就好，如果没有可以google到biostars上查询如何仅仅凭借原始fq文件判断限制酶信息，
有使用回帖ref看切割位点，有使用qc报告看k-mers，也有脚本估计的；
例如参考：https://mp.weixin.qq.com/s/82x7pNrd0t28-G0ofcIvuA(只是举个例子，尽量使用权威lab的脚本)，使用对应脚本见https://github.com/MaybeBio/TNBC-project/blob/main/1%2CHiC-pro upstream processing/predicted_Hic_restriction_enzyme.py。

确定了酶切信息之后，参考https://github.com/nservant/HiC-Pro/blob/master/doc/FAQ.md
FAQ页面有：
这里有如何产生bed文件方法，大意是版本更新之后有一个专门的函数，即utility digest_genome.py用来产生bed文件，

具体使用方法跳转到utility页面中查看：
https://github.com/nservant/HiC-Pro/blob/master/doc/UTILS.md

可以看到官方已经内置了一些单酶切的情况处理：
HindIII, DpnII, BglII and MboI，可以使用酶切位点文件（^指示），或者指明对应的酶；
当然还有多酶切，就比较复杂了，我本人曾经遇到过使用arima的酶切：

参考https://github.com/nservant/HiC-Pro/blob/master/doc/FAQ.md

https://github.com/nservant/HiC-Pro/issues/202

https://github.com/nservant/HiC-Pro/issues/328

总之按照官网指示以及issue中的回复，照猫画虎按图索骥处理自己的限制酶即可，
然后按照官网指示以及issue中使用者的经验，我执行了以下命令：

②染色体size的table文件：

可以从UCSC基因组浏览器中获取，但是个人建议还是依据自己处理之后的参考基因组fa文件重新处理（以hg19.fa为例），使用samtools+awk处理，


可以参考https://mp.weixin.qq.com/s/Z5kdd8M62TUGpQaz1NM6mA，

③bowtie2软件建立的索引

参考https://bowtie-bio.sourceforge.net/bowtie2/manual.shtml


实际操作时-f参数实际上就是指定参考基因组fa文件路径，也就是前置处理中去除了异常chr的hg19.fa的路径，生成文件的前缀<bt2_base>参数也可以提供一个文件夹路径的前缀，相当于将生成文件移动到某个文件夹中。总之生成6个文件（4个bt2,2个rev bt2），

2，配置并运行hic-pro

（1）明确命令运行所需配置参数以及文件：

基本上只要关注 -i、-o、-c 参数即可，其他参数都是并行、分步用的，本人一般直接一步运行。

-i参数：

所提供的文件夹嵌套3层，最外层提供给-i参数，第1子层是你自己定义的样本名称（注意这个要命名好，后续输出结果每个样本文件也是以这个为前缀命名），第2子层是具体样本的fq文件放置。
例如：

-o参数：
输出文件夹

-i参数：
运行hic-pro所需要的细节参数配置文件，
参考https://github.com/MaybeBio/TNBC-project/blob/main/1%2CHiC-pro upstream processing/config-hicpro.txt以及https://github.com/MaybeBio/TNBC-project/blob/main/1%2CHiC-pro upstream processing/config-hicpro_mb231.txt，

以上部分无需修改；

这里主要就是依据任务需求以及服务器的配置修改参数N_CPU、SORT_RAM、JOB_MEM参数，其他的暂时没有修改的必要，可以参考其他博客公众号推文、或者issue经验进行参考修改。

data这里就是-i参数中第3层输入的fq文件的后缀要对应，即fq文件_R1/2.fastq（或_R1/2.fastq.gz，解压与否无所谓）对应后缀_R1/2。

主要修改参数：

BOWTIE2_IDX_PATH=参考基因组建立索引index的路径位置，即6个bt2文件所在的位置（文件夹前缀，不带\）
REFERENCE_GENOME=注意上面这个参数对应，当时候bowtie2-build命名的前缀是什么，这里就用什么，也就是bt2文件的前缀，所以建议bowtie2-build建立index时前缀就命名为hg38、hg19之类，然后这里也就能够直接填hg38、hg19了  
GENOME_SIZE=前面处理3个注释文件中的chr size文件，提供全路径即可

然后就是限制酶切信息设置：

GENOME_FRAGMENT=前面处理的3个注释文件中的酶切片段bed文件，提供全路径即可  
LIGATION_SITE=酶切位点的相关序列，这个参考前面什么酶选什么切割序列，同样参考issue以及官方指南，找准自己用什么酶、位点如何、切割之后序列怎么样   

然后就是这里构建的互作map数据：

BIN_SIZE=实际上就是分辨率，以bp为单位，hic-pro有默认的bin size序列，juicer也有默认的bin size序列，当然以自己数据分析所需要的bin size为主进行构建，或者直接按照juicer序列构建所有bin size序列  
MATRIX_FORMAT=这里实际上是输出的矩阵是上三角还是全矩阵，注意这里按照默认设置即可，即upper，没有充分理由以及需求不要改为complete，就按照default即可！！！！！！！

总之以代码条内为准+白色图片，黑色图片仅供以前执行的经验，不作参考。
有问题看官方https://github.com/nservant/HiC-Pro/blob/master/doc/MANUAL.md。

例如：
我用的是arima，所有参数配置参考

参考https://github.com/nf-core/hic/issues/127；
其他无关参数建议空着默认即可https://github.com/nservant/HiC-Pro/issues/214。


注意：配置文件中不要加注释，会出现无法识别的报错，无论是中英文！
然后就是运行的时间，按照15左右的cpu以及如上配置，我大概处理7个样本，每个样本在300G左右即总共2T数据的情况下至少需要1周时间（确切的说是前6个样本花了7天，后1个样本花了2天左右，大概在1周+），所以建议nohup挂在后台即可。

在运行时还遇到样本数据量不够的情况，如上MB231，主要是数据深度不够，可以参考https://www.biostars.org/p/486204/，合并MB231的hic数据：



200million以上，应该是数据深度够了，可以看reads数目等。
总之，合并不合并，是否科学严谨，看自己数据分析需求，另外参考CNS经验，开发者经验等即可，上面仅供参考。