HiC-Pro

2022-09-19 本文已影响0人 iBioinformatics

https://blog.csdn.net/qq_50637636/article/details/120827491

一、简介

HiC-Pro，可能是最好用的HiC数据标准化处理工具，源于“HiC data processing”。流程主要分为6个步骤，比对、过滤HiC比对结果、检测有效HiC序列、结果合并、构建HiC关联图谱以及关联图谱标准化。

HiC-Pro可以输出各个尺度的HiC标准化互作图谱，从每个窗口5Kb到1Mb，当然窗口越大，所生成的互作关系越少，计算时间就越少，反之则越多。此外HiC-Pro还可以进行HICUP、HiCdat等众多HiC数据处理软件所不能进行的等位基因特异性HiC分析。

HiC-Pro于2015年发表于Genome biology（https://genomebiology.biomedcentral.com/articles/10.1186/s13059-015-0831-x），有兴趣的可以去查看原文，更加深入了解。

完整的 pipeline 如下图所示

红色方框标记的是数据预处理部分，功能和HICUP软件类似，包括序列比对和筛选valid pairs；
预处理之后就是binning, 构建不同分辨率下的原始的交互矩阵contact map,
最后对原始的contact map进行归一化处理，得到校正后的contact map。

HiC-Pro的一个强大功能在于可以构建单倍型级别的Hi-C图谱，单倍型级别的Hi-C图谱有助于更加精细化理解基因组三维结构，进一步对基因调控等功能进行深入细致的研究。

整个处理过程分为以下几个步骤

1. 序列比对

HiC-Pro采用了两步比对的策略，如下所示

考虑到连接点在插入片段上的位置和测序读长的关系，第一步先将R1和R2端分别与基因组比对，对于没有比对上的reads, 可能是存在连接点的嵌合体reads, 也可能本身就是unmapping reads, 通过从3’端切除部分序列的方式，使得嵌合体序列也能够比对上基因组，两步策略保证了数据的利用率。

2. 筛选valid pairs

比对时将R1和R2端分开单独考虑，但是二者其实来自于同一个fragment, 这一步的筛选其实是能够代表染色质交互的有效fragment，这样的fragment肯定是一个嵌合体序列，有来自交互作用的两个染色质区域的序列构成，如下图所示

只有来自嵌合体fragment的reads才被定义为valid pairs, 然后进行后续分析。

3. 构建原始Hi-C图谱

根据指定的分辨率，统计两个bin区域内valid pairs的数目, 去除PCR重复之后，构建原始的交互矩阵。

4. 归一化

不同区域GC含量，mapping概率等系统误差都使得原始的交互矩阵不能够有效代表染色质交互信息，所以需要进行归一化。采用了一种迭代校正的归一化算法对原始的交互矩阵进行归一化，矫正系统误差。

HIC-Pro还提供了一系列的质控标准，如下图所示

一个高质量的文库绝大部分肯定都能够比对上基因组，如图A所示， R1和R2的比对率都很高。而比对上的reads中应该主要是unique mapping, 如图A第二张图所示，multiple hits和low quality也是文库质量的指标之一。

valid pairs的比例则是文库质量的最直接体现，valid pairs的比例至少要在50%以上。

将染色质交互作用进一步区分为染色质之间的inter-interaction. 对应B图中的trans contact, 和染色质内部的intra-interaction, 对应cis contact。对于cis contact, 根据距离阈值分成short和long两种。

一个高质量的文库首先intra-interaction的比例在40%以上，其次由于线性距离近的染色质更容易随机结合，引入系统误差，所以高质量文库的cis long contacts的比例在40%以上。

HiC-Pro所有的参数都放置在一个配置文件中，既可以一键化运行整个pipeline, 也可以分布运行，单独执行其中的某几步，灵活性很强，后续会介绍其详细用法。

二、准备三个文件：bed、table、bowtie2构建的索引

1. bowtie2索引构建，产生六个.bt2新文件

bowtie2-build [options]  <reference>  <bt2_index_base>

reference : 下载的参考基因组，genome.fa
bt2_index_base: 构建索引前缀

2. 使用digest_genome.py生成酶切片段文件

python HiC-Pro/bin/utils/digest_genome.py \
  -r [常用限制性内切酶序列]  \
  [-o OUT] fastafile

-r：常用限制性内切酶：

限制性内切酶	酶切位点，`^`为切割位点
MboI	^GATC
DpnII	^GATC
BglII	A^GATCT
HindIII	A^AGCTT

生成的bed文件

3. 生成基因组sizes文件,获得基因组每条染色体bases数bed文件

samtools faidx genome.fa
awk ‘{print $1 "t" $2}‘ genome.fa.fai > genome_sizes.bed

genome_sizes.bed

4. Hi-C数据准备

创建sample文件夹，一个文件夹放入一个样品的fastq文件（生物学重复可以放入）

5. 配置 Config 文件

vi ./config-install.txt

需要修改的参数有：
N_CPU:给定的CPU内存数，给的越多，运行的越快（根据服务器配置）；
LOGFILE:日志文件的名称；
JOB_MEM:内存的大小
PAIR1_EXT= _R1 :R1测序数据名称中有_R1
PAIR2_EXT = _R2:R2测序数据名称中有_R2
MIN_MAPQ: 最低的质量分数，用于筛选，表示低于该MAPQ值会被过滤
BOWTIE2_IDX_PATH: 基因组bowtie2索引路径，eg:/path/hg19
BOWTIE2_GLOBAL_OPTIONS: 默认GLOBAL比对设置
BOWTIE2_LOCAL_OPTIONS: 默认LOCAL比对设置
REFERENCE_GENOME： Bowtie2索引前缀
GENOME_SIZE：基因组sizes bed文件
GENOME_FRAGMENT: 基因组酶切文件,eg. /path/hg19_HindIII.bed
LIGATION_SITE: 酶切位点末端补平再次连接后形成的嵌合序列，eg. AAGCTAGCTT
MIN_FRAG_SIZE: 最小的理论酶切片段大小,eg. 100
MAX_FRAG_SIZE: 最大的理论酶切片段大小,eg. 100000
MIN_INSERT_SIZE: 最小的文库片段大小,eg.100
MAX_INSERT_SIZE: 最大的文库片段大小,eg.1000
BIN_SIZE:需要生成的矩阵分辨率（bp)
MATRIX_FORMAT：矩阵的形式，upper表示保留上半部分

1-12行不需要修改。

13-25行：
16行CPU数目，21行内存大小，越大计算越快，其余几行不用改。

26-32行：
原始的下机数据的前缀的最后要以_R1/_R2结尾，当然也可以在配置文件里修改成与raw_data前缀的最后一样。

37-42行：
39行改为构建索引的绝对路径，40、41不用改。

43-49行：
47行是索引的前缀，48行保存染色体/contig大小文件的绝对路径。

64-73行：
67行是基因组酶切的bed文件的绝对路径；
68行为酶切连接位点，Mobi酶切连接位点为GATCGATC；
69行酶切后最小的片段长度，70行酶切后最大的片段长度，这个范围越大，reads数越多；
71，72行插入片段的大小，HiC建库的插入片段长度一般在300-500bp。

85-91行：
89行是将bin的长度，我选择了10kb-1Mb之间，bin的值越小，分辨率越高，处理起来占用内存空间越大，后面做的图数据量能达到的话，结果会更精准吧。
90行是生成矩阵的形式，可选参数为：complete、asis、upper和lower，默认值为upper。如果后续还要做A/B compartment分析的话，需修改成complete，否则会影响PCA分析的结果。

6. HiC-Pro 运行

HiC-Pro  -i INPUT -o OUTPUT -c CONFIG [-s ANALYSIS_STEP]  [options]

-c: config文件路径
-o：结果生成路径
-i: 原始数据路径
-p: 集群运行

7. 结果解读

总目录

bowtie_results: 比对结果目录
hic_results: hic矩阵及分析结果目录
logs: 存放分析日志
rawdata：链接了原始数据
tmp：存放中间文件

7.1 Bowtie_result 目录

bwt2:存放合并后的bam文件和统计结果
bwt2_global:存放全局比对结果
bwt2_local:存放局部比对结果

7.2 hic_result 目录

data:存放validpair及其他无效数据文件
matrix:存放不同分辨率矩阵文件
pic:存放统计分析图片
stats:存放统计表

7.3 Data 文件

allVaildPairs:合并后的pairs数据
DEPairs:Dangling end pairs数据
DumpPairs:实际片段长度和理论片段长度

不同的数据
REPairs：酶切片段重新连接的pairs
FiltePairs:MAPQ过低的pairs
SCPairs：片段自连的pairs

Matrix文件

raw:原始矩阵

iced：ice标准化后的矩阵

Pic文件，出图