使用HiCPro绘制全基因组互作矩阵
为了图片的美观可以将挂载上的染色体先提取出来：genome.chr.fasta
①使用bowtie2-build为基因组建立索引
bowtie2-build genome.chr.fa genome
^基因组文件 ^产生的索引文件的前缀
②产生酶切位点信息.bed文件
digest_genome.py -r mboi -o genome_mboi.bed genome.chr.fa
-r:HiC测序使用的酶
-o:产生的bed文件名
genome.chr.fa 基因组文件
③产生基因组大小信息文件.sizes文件
seqkit fx2tab -l -n -i genome.chr.fasta > genome.sizes
cat genome.sizes
chromosome1 48776989
chromosome2 47893571
chromosome3 43383427
chromosome4 43078513
chromosome5 41810944
......
④修改HiCPro配置文件(请注意：为了避免报错，请都修改为绝对位置)
必须要修改的地方：PAIR1_EXT\PAIR1_EXT 确认HiC测序数据的特征符
BOWTIE2_IDX_PATH bowtie2建立索引的绝对位置
GENOME_SIZE 3步骤中产生的基因组大小信息文件.sizes文件的绝对位置
GENOME_FRAGMENT 2步骤中产生的酶切位点信息.bed文件的绝对位置
LIGATION_SITE 酶切位点
BIN_SIZE 生成矩阵的分辨率
vim config-hicpro.txt
#########################
# Please change the variable settings below if necessary
#########################################################################
## Paths and Settings - Do not edit !
#########################################################################
TMP_DIR = tmp
LOGS_DIR = logs
BOWTIE2_OUTPUT_DIR = bowtie_results
MAPC_OUTPUT = hic_results
RAW_DIR = rawdata
#######################################################################
## SYSTEM AND SCHEDULER - Start Editing Here !!
#######################################################################
N_CPU = 4 #运行的CPU数目
SORT_RAM = 500M #每个任务的运行内存
LOGFILE = hicpro.log
JOB_NAME =
JOB_MEM =
JOB_WALLTIME =
JOB_QUEUE =
JOB_MAIL =
#########################################################################
## Data
#########################################################################
PAIR1_EXT = _R1 #HiC测序数据的特征信息 请检查此处 你的测序数据应该包含“_R1”符号
PAIR2_EXT = _R2
#######################################################################
## Alignment options
#######################################################################
MIN_MAPQ = 10
BOWTIE2_IDX_PATH = /media/dell/04-HiC-matrix #bowtie2建立索引的目录 请修改此处
BOWTIE2_GLOBAL_OPTIONS = --very-sensitive -L 30 --score-min L,-0.6,-0.2 --end-to-end --reorder
BOWTIE2_LOCAL_OPTIONS = --very-sensitive -L 20 --score-min L,-0.6,-0.2 --end-to-end --reorder
#######################################################################
## Annotation files
#######################################################################
REFERENCE_GENOME = genome #产生的文件名
GENOME_SIZE = genome.sizes #3步骤中产生的基因组大小信息文件.sizes文件
#######################################################################
## Allele specific analysis
#######################################################################
ALLELE_SPECIFIC_SNP =
#######################################################################
## Capture Hi-C analysis
#######################################################################
CAPTURE_TARGET =
REPORT_CAPTURE_REPORTER = 1
#######################################################################
## Digestion Hi-C
#######################################################################
GENOME_FRAGMENT = genome_mboi.bed #2步骤中产生的酶切位点信息.bed文件
LIGATION_SITE = GATC #酶切位点HindIII(AAGCTAGCTT), mboi(GATC) , DpnII(GATCGATC), NcoI(CCATGCATGG)
MIN_FRAG_SIZE =
MAX_FRAG_SIZE =
MIN_INSERT_SIZE =
MAX_INSERT_SIZE =
#######################################################################
## Hi-C processing
#######################################################################
MIN_CIS_DIST =
GET_ALL_INTERACTION_CLASSES = 1
GET_PROCESS_SAM = 0
RM_SINGLETON = 1
RM_MULTI = 1
RM_DUP = 1
#######################################################################
## Contact Maps
#######################################################################
BIN_SIZE = 100000 200000 500000 1000000 #生成矩阵的分辨率
MATRIX_FORMAT = upper
#######################################################################
## Normalization
#######################################################################
MAX_ITER = 100
FILTER_LOW_COUNT_PERC = 0.02
FILTER_HIGH_COUNT_PERC = 0
EPS = 0.1
运行HiC-Pro
HiC-Pro -i INPUT -o OUTPUT -c CONFIG [-s ANALYSIS_STEP] [-p] [-h] [-v]
-i 存放HiC数据的目录位置
-o 输出的文件夹名
-c HiCPro配置文件
###########
如果运行结果中断，可以添加-s参数，并同时修改-i数据的位置，继续运行HiC-Pro，详见HiC-Pro -h。
产生的结果文件在
cd OUTPUT/hic_results/matrix/HiC
#将raw文件中相应分辨率下的文件拷贝到iced文件中相应的分辨率文件夹下
cp raw/1000000/* iced/1000000/
...
cd iced/1000000/
⑤矩阵可视化
#使用HiCPlotter.py进行可视化
python2.7 HiCPlotter.py -o genome \
-f genome_500000_iced.matrix \
-r 500000 -tri 1 \
-bed genome_500000_abs.bed \
-n genome \
-wg 1 -chr chromosome7
-o 输出的文件名
-f _500000_iced.matrix产生的矩阵文件
-r 矩阵的分辨率
-bed _500000_abs.bed产生的bed文件
-n 输出图片最上方的名字
-chr 最后一号染色体的名字 可使用"tail -n 1 *.bed"命令查看
结果
全基因组互作矩阵
