juicer工具的使用

juicer工具是AidenLab推出的一款处理Hi-C数据的工具，来源于他们的大作Juicer Provides a One-Click System for Analyzing Loop-Resolution Hi-C Experiments

1. 安装

juicer最近的稳定版本可以在这里得到。(不用在这里下载，可以使用git clone一步到位)

具体的安装过程可以参考这里。

这篇博客里边写的建立链接的部分稍有问题，稍作修改之后，适合我们的版本是这样的:

# 创建juicer文件夹
mkdir /home/xxx/juicer
cd /home/xxx/juicer
# 创建juicer必须的文件夹
mkdir references
mkdir work
mkdir restriction_sites
# 从github clone最新版本的juicer
cd /home/xxx/package
git clone https://github.com/aidenlab/juicer.git
# 此时，juicer包被下载到package文件夹中，路径为/home/xxx/package/juicer
# 我们要建立，自己创建的juicer文件夹和juicer包中，scripts的链接关系
cd ../juicer
ln -s /home/xxx/package/juicer/CPU scripts
cd ./scripts/common
wegt https://hicfiles.tc4ga.com/public/juicer/juicer_tools.1.9.9_jcuda.0.8.jar
ln -s juicer_tools.1.9.9_jcuda.0.8.jar juicer_tools.jar

至此，juicer的安装就完成了。

2. 使用

在之前的步骤中，我们建立了references，work，restriction_sites三个文件夹，这三个文件夹的作用分别是：

/juicer/references # 存放参考基因组相关文件的文件夹
/juicer/work # 存放样本的序列文件，和分析结果的文件夹
/juicer/restriction_sites # 存放参考基因组酶切图谱的文件夹

在使用juicer进行数据处理之前，我们需要先准备相关的/references和restriction_sites文件夹。

准备references

我们要处理的数据是3dvi文中提到过的Li2019数据集，下载地址如图：

这是一个小鼠的数据集，因此，我们去ucsc上下载小鼠的参考基因组的序列文件，下载地址在这里。

ucsc提供了fasta文件，我们需要使用bwa为参考基因组文件建立索引。

bwa的安装和使用可以参考这篇博客。

服务器上的bwa已经安装好了，从fasta文件，我们需要得到bwa index的输出文件，即*.amb，*.ann，*.bwt，*.pac，*.sa这五个文件。

此时，mm9.fasta已经在references中准备好了。打开终端，执行以下命令：

bwa index mm9.fasta

准备restriction_sites

参考李国强等人的论文，他们在实验中所使用的酶是DpnII，juicer中支持直接输出这种内切酶的酶切图谱，因此，我们只需要运行以下命令：

按照博客中的说法，我们只需要执行以下命令：

cd ./references
python /home/penglab/package/juicer/misc/generate_site_positions.py DpnII mm9 mm9.fasta 
# 三个参数分别为 内切酶名称，参考基因组名称，参考基因组序列文件的路径

此时，mm9_DpnII.txt已经在references中准备好了。

另一个需要准备的文件是mm9.chrom.sizes，这个文件是参考基因组的染色体长度文件，可以从这里下载。

准备work

需要处理的文件从上图所示的路径中下载(此链接已经失效，如有需要可以找论文原作者进行询问)
软件运行时对样本文件的存放位置也有要求，必须位于work目录下，以样本名作为一个子目录，序列文件存放于fastq目录下。

这里，我仅从下载下来的数据中选择了一组进行验证性的操作。
将JL_392_1.ACTTGA.mm9.calmd.bam.mapQ10.noDup.bam.mapQ20.sam2juice.filter.gz移动到work目录下，新建JL_392_1文件夹，将上述文件解压至JL_392_1文件夹。

此时，全部的准备工作都已经完成了，我们可以开始运行juicer了。
此时，juicer文件夹下的文件结构如下图所示。

开始处理

juicer的用法如下：

cd ./scripts
bash juicer.sh \
-z /home/xxx/juicer/references/mm9.fasta 
-p /home/xxx/juicer/restriction_sites/mm9.chrom.sizes \
-y /home/xxx/juicer/restriction_sites/mm9_DpnII.txt \
-d /home/xxx/juicer/work/JL_392_1 \
-D /home/xxx/juicer \
# -z参数指定参考基因组fasta所在路径,在该路径下必须同时存在对应的bwa索引；
# -p参数指定染色体长度文件；
# -y指定基因组酶切图谱的路径；
# -d指定样本原始文件存放的路径；
# -D指定软件的安装路径

第一次跑踩了个坑，没有samtools。参考博客安装samtools。

这种做法不太适合我们的数据，使用下边的方法。

---

经过跟3dvi的作者进行讨论，发现可能是不需要使用juicer.sh对文件进行处理，直接使用juicer_tools进行处理即可。
参考这里的说明

我们使用下边的命令，可以得到全分辨率所有染色体的原始接触矩阵的.hic文件。

# "Create index for data at all resolutions, all chromosomes"
java -Xmx8g -jar /home/penglab/juicer/scripts/common/juicer_tools.jar pre -v -q 1 JL_392_1.ACTTGA.mm9.calmd.bam.mapQ10.noDup.bam.mapQ20.sam2juice.filter test.hic mm9

对每个文件都进行同样的操作即可，参考juicer作者给出的实例，可能在进行转化之前不解压也是可行的。

使用juicer将.hic文件转为.txt文件

此时，我们已经得到了.hic文件，使用juicer工具将其转化为我们熟悉的.txt文件。

java -jar /home/penglab/juicer/scripts/common/juicer_tools.jar dump observed NONE test.hic 1:0:200000000 1:0:200000000 BP 10000 ./test.txt

其他参数不用修改，每次修改BP后的参数，得到不同分辨率的稀疏矩阵。

在得到.txt文件之后，这时他的坐标值还需要除以我们选取的分辨率，得到真实的矩阵坐标，这个可以用任意语言完成。

对比实验

按照juicer作者的说法，long format格式的文件也可以直接bin以得到不同分辨率的接触矩阵。

<str1> <chr1> <pos1> <frag1> <str2> <chr2> <pos2> <frag2> <mapq1> <cigar1> <sequence1> <mapq2> <cigar2> <sequence2> <readname1> <readname2>

我们需要用到的是<chr1> <pos1> <frag1> <str2> <chr2> <pos2> <frag2>，将这六列提取出来，其他的操作与对其他数据集的操作大同小异。

将这样直接从long formatbin的结果和由juicer工具处理得到的结果进行对比。

目标

• 熟悉juicer的使用，熟悉hic矩阵的格式；
• 得到Li2019数据集不同分辨率的接触矩阵，感兴趣的主要为25k，40k和50k。

reference

• Juicer Provides a One-Click System for Analyzing Loop-Resolution Hi-C Experiments
• Juicer软件的安装详解
• 「Workshop」第二十五期：HiC 数据处理简介
• Normalization and de-noising of single-cell Hi-C data with BandNorm and scVI-3D
• Joint profiling of DNA methylation and chromatin architecture in single cells
• BWA使用详解
• Juicer Pre Wiki