如何绘制简单漂亮的全基因组互作矩阵HiC matrix

使用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

产生的结果文件在

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"命令查看