Hi-C染色体挂载（2）——3D-DNA配置运行和手动调整Hi-C互作图

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

# 拉取3D-DNA仓库
git clone https://ghproxy.com/https://github.com/aidenlab/3d-dna.git
## 将run-asm-pipeline-post-review.sh和run-asm-pipeline.sh属性分别改成可执行，对应chmod的0755权限
chmod 755 run-asm-pipeline-post-review.sh
chmod 755 run-asm-pipeline.sh

# 安装LastZ（跑二倍体模式需要用，一个DNA序列比对工具）
conda install LastZ

# 安装GNU Parallel
## 官方推荐用该shell下实现并行运算的工具，加快程序运行速度。集群用户如果没有预装可以在~/bin中以如下方式编译安装
wget http://ftp.gnu.org/gnu/parallel/parallel-latest.tar.bz2
tar -xvjf parallel-latest.tar.bz2
cd parallel-20230622
./configure --prefix=$HOME && make && make install

3D-DNA有两个pipeline脚本，12个独立的功能模块（每个模块一个文件夹），官方给出的pipeline如下：

三个橘色框是需要手动调整的地方，稍微了解下这几个独立模块分别在做什么（具体参考manual_180319 (aidenlab.org)）：

主要的8个独立模块（对应上面的流程）：

• scaffold：对给定的scaffolds进行排序和定向，转换基于3D proximity的一系列单个片段序列成为单个的megascaffold，这个megascaffold用于检测错误的连接或者保留用于后续处理（split模块还会用到）

• visualize：与后续的Juicebox Assembly Tools处理对接，也就是可视化组装结果。在3D-DNA中用于连接各个处理阶段，便于审查和调整参数。

• edit：包含检测3D-DNA错误连接的脚本和一些校正算法，这个模块在做特殊处理的时候最可能需要调整参数

• polish：校正3D-DNA scaffolding算法可能引起的错误，类型基因组的polish过程

• split：分割megascaffold成为染色体

• seal：消除上一步结果的假阳性（通过引入碎片的方式，详细看手册说明）

• merge：融合代表替代单倍型的组装片段（仅在二倍体模式下）起作用

• finalize：生成最终的fasta格式输出，在最终fasta生成过程中，考虑到可能的顺序和方向，将被识别为同一染色体内部的scaffold连接起来，同时在每个scaffold之间添加固定大小（500bp）的间隙

附加的4个独立模块：

• utils：一些其他模块用到的核心脚本，包括从.fasta文件生成.assembly文件
• lift：一些从基因组草图到最后组装结果的核心脚本
• supp：几个附加脚本，包括生成色谱图的脚本
• data：一些数据表

如果要研究各个步骤的算法就可以在各个模块中找源码，现在的我只是个调包侠，程序不出问题就不深入研究了（以后有时间一定= =）

两个pipeline脚本也是对应不同阶段使用的：

run-asm-pipeline.sh 将几个独立模块串联，组装基因组（也就是上面流程图用的脚本）

run-asm-pipeline-post-review.sh 用于Juicebox Assembly Tools (JBAT)手动调整之后执行，只会执行上个脚本的后几个阶段（finalize或者seal和finalize），生成最终的hic文件（用于质控）和fasta序列

了解不同脚本和模块作用后，接下来就是很简单的设置脚本参数运行了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

# 创建run.slurm脚本
vim run.slurm

# 脚本内容如下
---------------------------------------------------------------
#!/bin/bash
#SBATCH -n 30
#SBATCH -N 1
#SBATCH -t 7200

/public/home/wlxie/biosoft/3d-dna/run-asm-pipeline.sh \
		/public/home/wlxie/biosoft/juicer/reference/genome.fa \
		-r 0 \
		/public/home/wlxie/biosoft/juicer/aligned/merged_nodups.txt
---------------------------------------------------------------

# 运行脚本
sbatch run.slurm

220Mb参考基因组，50G的merged_nodups.txt文件，实际运行结束时长为10小时。

run-asm-pipeline.sh脚本用法和参数如下：

1
2
3
4
5
6
7
8
9
10

USAGE: ./run-asm-pipeline.sh [options] <path_to_input_fasta> <path_to_input_mnd>

path_to_input_fasta：组装的基因组fasta文件位置
path_to_input_mnd：juicer处理后得到的去重复的Hi-C reads比对基因组文件位置，也就是merged_nodups.txt

OPTIONS:
-m|--mode haploid/diploid	单倍体/二倍体模式（默认单倍体）
-i|--input input_size	指定输入的 contig/scaffold size阈值，默认15000，小于15000的将被忽略
-r|--rounds number_of_edit_rounds 	指定纠错轮数（默认为2轮）
-s|--stage stage	指定运行的阶段，可以是polish, split, seal, merge和finalize

对于运行模式，作者团队建议是先运行默认的haploid模式，检查拼接后的图谱中是否存在与杂合度不足相关的信号，diploid模式会尝试合并一些单倍型。

The diploid mode is to be used when one suspects large amounts of under collapsed heterozygosity to be present in the draft genome assembly (this was indeed the case for AaegL2). Running in the diploid mode will attempt to analyze the final assembly to remove the interspersed haplotigs and merge overlaps not recognized by the contigger because of variant differences. In most cases the default haploid mode should be used.

可以看官方的这篇文章帮助理解下两种模式的区别，以及软件LastZ在diploid模式中起的作用：Extended_Hi-C_methods.docx (dropbox.com)

至于纠错轮数，最好是按照默认的参数跑2轮，然后选择效果最好的一次结果导入juicerbox进行调整。

那么问题来了，怎么确认用哪个结果呢？首先要确认下结果文件有些啥。

这个软件产生的结果文件非常多，中间步骤产生的文件是不会自动删除的，看一下主要的几个文件：

• .fasta 基因组序列文件，最终组装的基因组序列名称为.FINAL.fasta

• .hic 高度压缩的二进制文件，与下面的.assembly文件都可以载入Juicebox Assembly Tools，不同阶段产生不同的hic文件，0表示未校正

• .assembly 空格分隔的文本，记录对基因组草图执行的指令，包括拆分、更改顺序、方向和锚定到染色体中（不同阶段都会产生，同时产生同名的.cprops和.asm，这两种文件已弃用）

• .scaffold_track.txt & .superscaf_track.txt scaffold和super scaffold的染色体边界文件，Juicebox 2D 注释格式，使用juicerbox的时候可以用到，但是使用Juicebox Assembly Tools不需要用（因为边界注释会根据. assembly 文件自动生成）
• .bed & .wig 对pipeline进行故障审查的时候用到，可以与.hic文件导入juicebox，暂时与Juicebox Assembly Tools不兼容

可以看出很多中间文件还是为了导入juicebox做分析产生的，而我们仅仅是为了辅助基因组组装，确认染色体的挂载情况，用到的是Aiden团队开发的juicebox的桌面应用拓展模块，也就是Juicebox Assembly Tools(JBAT)，很多中间文件已经用不到了，这里我们需要用到的是.hic和.assembly两种类型文件，且这两种文件名称是一一对应的。

如果是默认参数跑的话还会有genome.1.hic和genome.2.hic两个文件，分别代表校正了一次和两次的结果，都可以导入JBAT看结果。

看很多教程其实陷入了误区，需要校正多少轮，用哪一组数据其实没有固定说法，有的人用genome.0.hic，有的人用genome.1.hic，都是可以的，这几个只是polish前和polish后产生的hic文件，别忘了我们还有其他阶段产生的hic文件，都是可以用的。这里我用genome.rawchrom.hic和对应的genome.rawchrom.assembly这组文件，因为这组文件实际上区分染色体的情况最好。

下一步就是导入JBAT手动调整染色体边界和修改组装过程中产生的错误。

这一步在个人电脑上完成，最好保证电脑有4G以上的运存（2023年了个人电脑都有这个配置了吧），否则会巨卡无比。

我用的官方1.11.08版本的windows工具：

https://s3.amazonaws.com/hicfiles.tc4ga.com/public/Juicebox/Juicebox_1.11.08.exe

不要安装直接可以运行，首先导入.hic文件：

导入hic文件之后Assembly菜单就可以点了，导入对应名字的.assembly文件：

完整导入后界面如图所示，主要是用到以下几个菜单栏：

• ①：标准化方式，选择balance看起来舒服点= =
• ②：分辨率，调整显示的分表率大小
• ③：色彩范围，调整表示互作强度的色彩范围
• ④右下角几个颜色块代表图上线条框的区域处于什么状态。默认黄色线框的是待编辑区（这里没有，对区块进行操作的时候才有），紫色线框是染色体区块，绿色线框是scaffold区域

我这个互作图没有明显的需要染色体易位、倒位的调整，只需要确定染色体边界（也就是蓝色的线）就行。

染色体易位和倒位的调整在官方视频里就有，b站有人做了翻译，看一个视频足够了翻译 | Juicebox Assembly Tools教程_哔哩哔哩_bilibili

调高分辨率，比如下面我要做的就是将紫色框分成两个：

放大到一定分辨率后，scaffold的交界处鼠标会变成L型，点一下就可以加上染色体边界，这里截图没截出来……

还有一种方法，长按shift键，单击拖动鼠标框选住你要操作的scaffold（不需要完全覆盖scaffold，只要框的面积里有你要的scaffold就行），松开shift键，这个时候选中的scaffold会处于待编辑状态（黄色的框线）。右键，选择Add chr boundaries即可为选中的scaffold区域添加染色体边界。

如果想要撤销操作，右键后点击Undo即可。

整体调整完以后是这样：

染色体条数不是无缘无故添加的，既要以图为准，又要由实验确定，或者通过已经发表的文献（该物种或者近缘物种）佐证，纠正算法产生染色体边界的误差。我是通过查阅文献得到这个物种染色体条数为11，调整后的Hi-C互作图也能明显分为11个互作区块，证实组装是没有问题的。

官方的操作手册上也有详尽的手动调整Hi-C互作图教程，可以参考manual_180322.pdf (dropbox.com)。

手动调整完后，导出调整后的.assembly文件，我们需要用这个文件重新回到3D-DNA生成最终的染色体水平基因组fasta文件。

简单来说，我们在JBAT导出的最终.assembly文件记录了各个scaffold的坐标位置和互作信息，只需要在3D-DNA中跑最后finalize流程就可以了，软件提供了脚本run-asm-pipeline-post-review.sh：

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

USAGE: ./run-asm-pipeline-post-review.sh [options] -r <review.assembly> <path_to_input_fasta> <path_to_input_mnd> 

path_to_input_fasta：指定组装的草图基因组fasta文件
path_to_input_mnd：指定juicer产生的merged_nodups.txt

OPTIONS:
-r|--review path_to_review_assembly	指定上一步JBAT导出的.assembly文件
-i|--input input_size	指定输入的 contig/scaffold size阈值，默认15000，小于15000的将被忽略
-s|--stage stage	指定运行开始的阶段，可以是seal或者finalize，默认finalize
-q|--mapq mapq	指定最终map的可视化MAPQ阈值，默认1
-g|--gap_size gap_size	指定最终scaffold之间的间隔，默认500
--sort-output	对最终输出的染色体级别的scaffold长度排序，降序
--build-gapped-map	Option to output an additional contact map corresponding to the assembly after the gaps have been added between scaffolded sequences.（我也不知道干啥的）

比较重要的就三个文件位置参数：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

# 创建run_postview.slurm脚本
vim run_postview.slurm

# 脚本内容如下
---------------------------------------------------------------
#!/bin/bash
#SBATCH -n 30
#SBATCH -N 1
#SBATCH -t 7200

/public/home/wlxie/biosoft/3d-dna/run-asm-pipeline-post-review.sh \
         -r /public/home/wlxie/biosoft/3d-dna/baima_diploid/genome_rawchrom.review.assembly \
         /public/home/wlxie/biosoft/juicer/reference/genome.fa \
         /public/home/wlxie/biosoft/juicer/aligned/merged_nodups.txt
---------------------------------------------------------------

# 运行脚本
sbatch run_postview.slurm

静静等待生成最终的fasta序列就可以啦，如果前面做过contig基因组草图注释的话，还需要将注释的基因信息转坐标到最终的染色体水平基因组上；如果没做过注释的话，跑一遍注释流程吧！

最终生成的基因组序列要看genome.FINAL.fasta，发现染色体数目和前面手动标定染色体边界的数目不一致，这很正常。可以用seqkit工具统计下各个染色体的长度，如下：

可以看到前11条是染色体级别的contig，这和手动标记的染色体数目一致，后面的contig长度占比非常小，可以标记区分一下，这些都要放在组装的基因组文件里。

顺便放一下自己写的python代码，同样的处理效果，顺手就当python练习了（运行前提：fasta文件中没有空行）。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

fasta_file = 'genome.FINAL.fasta'	# 换成自己的fasta文件路径

def fasta_length(file_path):
    sequences = {}
    with open(file_path, 'r') as file:
        content = file.read()   # 读整个文件

    blocks = content.split('>')     # 根据关键字“>”分块

    for block in blocks[1:]:    # 第一个块为空，跳过
        index = block.find('\n')    # 查找第一个换行符索引
        newline_number = block.count('\n')    # 统计换行符数量
        header = block[:index]  # 序列名称
        sequence_length = len(block) - index - newline_number    # 序列长度
        sequences[header] = sequence_length # 序列名称和长度存入字典
    return sequences

sequence_lengths = fasta_length(fasta_file)

for header, length in sequence_lengths.items():
    print(f'{header}:{length}')

考虑到juicebox调整的这个图不能直接用在论文里（看起来有点不直观），这里再推荐两个项目：

Atvar2/plotHicGenome: The tools are used for displaying hic signal of genome (github.com)

biotools/Hic/juicerbox2matrix at main · yukaiquan/biotools (github.com)

第一个工具可以将Juicebox调整的结果（.review.assembly文件）和juicer2的结果（merged_nodups.txt）结合，做全基因组层面的Hi-C互作图。第二个工具是另一个作者优化了生成bin matrix矩阵的过程，测试了下吃个饭的功夫就跑完了。其实也可以用.hic以及.review.assembly文件做matrix，然后用HiCPlotter作图（怎么用HiCPlotter可以看之前的这篇博客 三维基因组学——如何用Hi-C数据分析拓扑关联结构域），效果基本上是差不多的。

这里记录下使用上面项目的过程和最后跑出来的效果：

1
2
3
4
5
6
7
8
9
10
11

#!/bin/bash
#SBATCH -n 8
#SBATCH -t 7200

assembly=/public/home/wlxie/biosoft/3d-dna/baima_diploid/genome_rawchrom.review.assembly
merged_nodups=/public/home/wlxie/biosoft/juicer/aligned/merged_nodups.txt

./juicerbox2matrix -a ${assembly} -c 11 -i ${merged_nodups} -o run

plotHicGenome juicer ${merged_nodups} ${assembly} -H run/Hicmatrix.txt -W whole -n 11 -s False -l t -F 4 -r 500000 -X 2 -w 0.5 -d 3 -S 'dashed' -i 300 -z 6,6 -C 'black' -L 0.8 -A 0.8 -B '1%' -D 0.2 -o Juicerbox.pdf -R ./run