Hi-C染色体挂载（1）——juicer2处理Hi-C数据

Juicer是一款非常经典的Hi-C数据处理软件，但是配置运行起来稍微有点麻烦，花了小半天时间踩坑记录一下。需要非常注意的一点，直接从github官网拉取的juicer仓库是juicer2，在release版本中有一个稳定版juicer1.6，两者的配置和产生的结果文件是不一样的！！！

秉着软件用新不用旧的原则，本篇笔记全程用的是juicer2，如果有软件运行问题可以在3D Genomics - Google 网上论坛交流，或者在github官方提Issues（不建议）。

顺便提一下我使用Juicer2过程中碰到的版本问题：

• juicer2配置的juicer_tools版本要在2.0以上，否则会报错Exception in thread "main" java.lang.RuntimeException: Unknown command: statistics，该报错会直接导致无法生成inter.txt、inter_hists.m 、inter_30.txt和inter_30_hists.m文件，也就是Hi-C互作的统计数据和矩阵。

• juicer2结果文件中不会自动产生merged_nodups.txt文件，该文件原本作为3D-DNA的输入文件，在juicer2中被同名的merged_nodups.bam文件和merged*.txt代替。如果想要merged_nodups.txt文件，需要加上参数--assembly。

1. 下载和配置juice2

1
2
3
4
5
6
7
8
9
10
11
12
13

# 拉取最新版本的juicer2仓库（这里用的github镜像，集群可能有dns污染无法直接访问github）
git clone https://ghproxy.com/https://github.com/aidenlab/juicer.git
cd juicer

# juicer主目录下配置scripts（必需）
## 个人电脑的是建立软链接到CPU文件夹下，其他如slurm、LSF、AWS和Univa等集群或者云端跑juicer是建立软链接或者复制到对应文件夹的scripts文件夹下
ln -s CPU scripts

# scripts/common文件夹中下载juicer_tools（必需）
cd scripts/common
wget -c https://ghproxy.com/https://github.com/aidenlab/Juicebox/releases/download/v2.20.00/juicer_tools.2.20.00.jar
# 修改文件属性为可执行文件后，建立软链接
ln -s juicer_tools.2.20.00.jar juicer_tools.jar

本来想配置slurm跑的，因为塔大集群使用的是slurm作业调度系统。花了好大功夫配置环境发现运行SLURM/scripts底下的juicer.sh还需要配置CUDA，晕，学校的集群没有GPU（怨念 = =）……所以如果集群没有GPU的话就老老实实用CPU文件夹下的juicer.sh，按照CPU流程跑后面的程序，缺点是不能用集群提交多线程作业非常难受。

官方也是推荐小的Hi-C实验可以用CPU版本，如果数据量比较大，还是推荐带有GPU加速的集群（最好是SLURM）或者云端跑。如下是官方推荐的两种方式：

Running Juicer on a cluster · aidenlab/juicer Wiki (github.com)

ENCODE-DCC/hic-pipeline: HiC uniform processing pipeline (github.com)

juicer主目录下需要配置两个必需的文件夹reference和restriction_sites（第二个在juicer流程中非必须，但是做染色体挂载一定要）

1
2
3
4
5
6
7
8
9
10

# juicer目录下创建参考基因组文件夹，放入参考基因组（复制或者软链接）和构建参考基因组索引文件（必需用bwa）
mkdir reference && cd reference
cp /path/to/your/reference/genome.fa ./
bwa index genome.fa

# juicer目录下创建限制性酶切位点文件夹，MboI酶酶切参考基因组（根据自己测hic用的限制酶）
mkdir restriction_sites && cd restriction_sites
python ~/biosoft/juicer/misc/generate_site_positions.py MboI genome ~/biosoft/juicer/reference/genome.fa
## 生成的酶切图谱文件名为genome_MboI.txt，同个文件夹下生成contig长度文件
awk 'BEGIN{OFS="\t"}{print $1, $NF}' genome_MboI.txt > genome.chrom.sizes

官方的酶切参考基因组的脚本generate_site_positions.py支持HindIII、DpnII、MboI、Sau3AI和Arima4种限制性核酸内切酶，如果你Hi-C测序用的是其他酶，可以根据实际修改这个脚本中的字典类型数据patterns，根据实际情况加入键值对信息：

稍微解释一下官网给的awk 'BEGIN{OFS="\t"}{print $1, $NF}'这个命令，这里awk的程序部分由两部分组成：

• BEGIN{OFS="\t"}：BEGIN块是在处理文件之前执行的代码块。在这里，它设置输出字段分隔符（OFS）为制表符（\t）。这意味着输出的字段将使用制表符进行分隔。
• {print $1, $NF}：这是awk的主要部分，它定义了要执行的操作。在这里，它打印每行的第一个字段$1和最后一个字段（NF）。通过使用逗号分隔它们，它们将以制表符分隔的形式打印出来。

genome.chrom.sizes文件如下所示：

工作目录下必需配置一个文件夹fastq，其中存放Hi-C的双端测序数据，命名格式参考juicer.sh文件中的正则表达式*_R*.fastq*。

1
2
3
4

# 工作目录下创建存放hic测序数据的文件夹，我这里为了方便查看，工作目录就是juicer主目录，测序数据可以用软链接的形式存放
mkdir fastq && cd fastq
ln -s ~/hi-c/BM_R1.fq.gz Av_R1.fastq.gz
ln -s ~/hi-c/BM_R2.fq.gz Av_R2.fastq.gz

2. 运行juicer2

官网的juicer.sh参数非常之多，我这里只展示一些关键的，其他用默认参数即可。1.x版本可以参考Usage · aidenlab/juicer Wiki (github.com)，2.0版本具体可以用bash juicer.sh -h查看，两个版本指令稍有不同。CPU版本和其他集群版本也稍有不同，以CPU中的2.0版本为例：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

Usage: juicer.sh [-g genomeID] [-d topDir] [-s site]
                 [-a about] [-S stage] [-p chrom.sizes path]
                 [-y restriction site file] [-z reference genome file]
                 [-D Juicer scripts parent dir] [-b ligation] [-t threads]
                 [-T threadsHic] [-i sample] [-k library] [-w wobble]
                 [-e] [-h] [-f] [-j] [-u] [-m] [--assembly] [--cleanup] [--qc]
-g genomeID		必需项，指定基因组和版本，比如人类的"hg19" 或者小鼠的"mm10"，如果没有，可以用-z命令替代，指定参考基因组文件
-z reference genome file	指定参考基因组文件，此文件中必需有BWA索引文件
-d topDir		指定输出目录，默认是当前目录。[topDir]/fastq必需包含fastq文件
-s site			必需项，指定限制性核酸内切酶。如果不做片段级别的分析可以写none
-S stage		指定运行pipeline的阶段，固定是"chimeric", "merge", "dedup", "final", "postproc", "early"其中之一，报错调试用
	chimeric：	alignment结束或者从aligned文件开始
	merge:	alignment结束但是没有生成merged_sort文件
	dedup：	文件已经merge结束，但是没有生成merged_nodups文件
	final：	reads在merged_nodups文件中已经删除重复，但是尚未生成最终的state和hic文件
	postproc：	hic文件已经生成，只有特征注释尚未完成
	early：	在hic文件生成前提前结束程序
-p chrom.sizes path		指定染色体长度文件
-y restriction site file		指定参考基因组酶切文件
-D Juicer scripts parent dir		指定Juicer/scripts所在目录
-t threads		指定跑BWA的线程数
-T threadsHic	指定生成hic文件用的核数（2.0版本新增）
--assembly		接下游3D-DNA分析，会提前结束并生成老版本的merged_nodups文件（2.0新增）
--cleanup		如果pipeline成功运行，自动清除中间文件（2.0新增）

如果是在集群中跑juicer，以下两个参数也是必需要改的，否则用默认分区名会直接报错：
-q queue		指定跑alignments的队列分区（默认commons），slurm中的分区也就是partition，可以理解为LSF、PBS等作业调度系统中的队列
-l long queue	指定跑需要时间更长的job，比如生成hic文件的队列分区（默认long）

因为我不是做三维基因组，Hi-C只测了一个样（主要用于辅助基因组组装，染色体挂载），比如做三维基因组就会有大量的Hi-C数据，juicer也支持多个结果的合并（使用mega.sh），详细也可以参考上面的juicer Wiki。

上面的限速步骤主要是跑BWA，因为第一次跑juicer不知道要多久，就稍微申请多一点的核数（虽然集群没有GPU，但是CPU资源很充足平常没人用）。

1
2
3
4
5

# 创建一个slurm作业
vim juicer.slurm

# 提交slurm作业
sbatch juicer.slurm

juicer.slurm文件内容如下（再次提醒--assembly 参数非常重要）：

1
2
3
4
5
6
7
8
9
10
11
12
13

#!/bin/bash
#SBATCH -N 1
#SBATCH -n 30
#SBATCH -t 7200

/public/home/wlxie/biosoft/juicer/scripts/juicer.sh \
-z /public/home/wlxie/biosoft/juicer/reference/genome.fa \
-p /public/home/wlxie/biosoft/juicer/restriction_sites/genome.chrom.sizes \
-y /public/home/wlxie/biosoft/juicer/restriction_sites/genome_MboI.txt \
-s MboI \
-D /public/home/wlxie/biosoft/juicer \
-t 30 \
--assembly

运行时间可以参考一下我这里220 Mb的参考基因组，Hi-C测序数据39G（双端测序，gz文件），实际上运行了13小时才出结果。

运行pipeline成功后会在日志中提示(-: Pipeline successfully completed (-:

3. 结果文件

CPU模式下会在工作目录创建aligned和splits文件夹：

• splits 文件夹存放整个pipeline的临时文件，跑完流程并且确认结果没问题后可以运行cleanup.sh或者直接删除（按照帮助文档的说法，–cleanup参数应该可以自动删除，我没有试过）
• aligned 文件夹存放运行结果
• 如果你是用集群模式跑的，比如SLURM，还会生成一个debug文件夹（CPU模式下用SLURM跑的没有该文件夹），可以查看各个步骤的error报告和output报告

aligned文件夹有以下结果文件：

我这里实际上是多了一些文件（实际上跑了两遍），因为一开始没有加上--assembly参数，所以一直运行到出hic结果文件，而且没有merged_nodups.txt。

我做染色体挂载只需要最后一个文件，第二次运行加了上面的参数，并且加上-S final参数，这样就省去了前面比对和去重复的时间（90%以上的时间花在这里），一个小时可以重新跑完出结果。

顺便了解一下其他结果文件的用途：

• inter.hic / inter_30.hic: The .hic files for Hi-C contacts at MAPQ > 0 and at MAPQ >= 30, respectively
• merged_nodups.txt: The Hi-C contacts with duplicates removed. This file is also input to the assembly and diploid pipelines （后面跑3D-DNA的输入文件）
• inter.txt, inter_hists.m / inter_30.txt, inter_30_hists.m: The statistics and graphs files for Hi-C contacts at MAPQ > 0 and at MAPQ >= 30, respectively. These are also stored within the respective .hic files in the header. The .m files can be loaded into Matlab. The statistics and graphs are displayed under Dataset Metrics when loaded into Juicebox （可以后续导入juicer box）