定义

之前的文章中有介绍过,HiC常用的几款软件的原理内容。可以点击链接访问了解一下
在这里不做赘述。

软件安装

3d-DNA

1
2
3
4
5
6
7
8
9
10
11

$ git clone https://hub.fastgit.org/aidenlab/3d-dna.git
$ cd 3d-dna
$ chmod 755 run-asm-pipeline.sh
$ chmod 755 run-asm-pipeline-post-review.sh

or
#github安装(2021年7月18日-目前的最新版本)
$ wget https://github.com/aidenlab/3d-dna/archive/refs/tags/201008.tar.gz
$ tar zxvf 201008.tar.gz

Juicer

1
2
3
4
5
6
git clone https://github.com/theaidenlab/juicer.git
cd juicer
ln -s CPU scripts
cd scripts/common
wget 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

要求环境

LastZ (version 1.03.73 released 20150708) – for diploid mode only
Java version >=1.7
Bash >=4
GNU Awk >=4.0.2
GNU coreutils sort >=8.11

3d-DNA使用

为基因组建索引

1
bwa index genome.fa

根据基因组构建创建可能的酶切位点文件

需要使用到juicer/misc/generate_site_positions.py

1
2
3
4
5
6
7
8
$ python /home/lixingze/software/juicer/misc/generate_site_positions.py 
Usage: /home/lixingze/software/juicer/misc/generate_site_positions.py <restriction enzyme> <genome> [location]

# <restriction enzyme> 选择对应的限制性内切酶
# <genome> 基因组序列文件
eg:
python /home/lixingze/software/juicer/misc/generate_site_positions.py DpnII genome genome.fa

运行如下命令, 获取每条contig的长度

1
awk 'BEGIN{OFS="\t"}{print $1, $NF}' genome_DpnII.txt > genome.chrom.sizes

运行juicer

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
$ /home/lixingze/software/juicer/scripts/juicer.sh -h
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 directory]
[-b ligation] [-t threads] [-h] [-f] [-j]
# -g: 自定义名称
# -s: 酶切类型, HindIII(AAGCTAGCTT), MboI(GATCGATC) , DpnII(GATCGATC), NcoI(CCATGCATGG),具体是什么酶切类型可以咨询测序公司
# -z : 参考基因组文件
# -y: 限制性酶切位点文件位置
# -p: 染色体大小文件
# -C: 将原来的文件进行拆分,必须是4的倍数,默认是90000000, 即22.5M reads
# -S: 和任务重运行有关,从中途的某一步开始,"merge", "dedup", "final", "postproc" 或 "early"
# -d: Hi-C 数据的存放目录
# -D: juicer的位置目录,我本人在/home/lixingze/software/juicer/
# -t: 线程数
# 必选项 -z -p -y

个人运行示例:

1
bash /home/lixingze/software/juicer/scripts/juicer.sh  -d /home/lixingze/data/HiC/05.3d-DNA-3cell/hic -D /home/lixingze/software/juicer/ -z ./genome.fa -y ./genome_DpnII.txt -p ./genome.chrom.sizes -s DpnII -t 70

输出的结果文件在aligned目录下,其中merged_nodups.txt就是下一步3D-DNA的输入文件之一。

运行3d-dna

在3d-dna目录下有个run-asm-pipeline.sh脚本,使用此脚本

1
2
3
4
5
6
7
USAGE: ./run-asm-pipeline.sh [options] <path_to_input_fasta> <path_to_input_mnd>
# <path_to_input_fasta> 输入参考基因组序列文件
# <path_to_input_mnd> 输入之前分析得到的merged_nodups.txt
# -m haploid/diploid 以特定模式运行,单倍体或二倍体(默认为单倍体)
# -i input_size 指定阈值输入Contigs/scaffolds大小(默认值为15000)。小于输入值的Contigs/scaffolds将被忽略。
# -r 2 指定错误联接更正的迭代轮数(默认值为2)。
# -s stage 快进到以后的组装步骤,可以进行polish, split, seal, merge, finalize

个人运行示例

1
nohup bash /home/lixingze/software/3d-dna/run-asm-pipeline.sh -r 2 genome.fa /home/lixingze/data/HiC/05.3d-DNA-3cell/hic/aligned/merged_nodups.txt &> 3d.log &

推荐使用 genome.0.hic文件以及 genome.0.assembly文件进行后续操作

juicebox调整3d-DNA输出的结果

这个调整过程需要细心耐心以及相关背景知识。

  1. aidenlab提供了在线的juicebox
    http://aidenlab.org/juicebox/

  2. 本地juicebox调整(推荐下载)
    https://github.com/aidenlab/juicebox/wiki/Download

网上有相关操作视频

调整完成之后将其保存为genome.review.assembly
如果是未发表的基因组,建议将染色体从大到小进行排列。

再次运行3d-DNA

这次使用run-asm-pipeline-post-review.sh脚本,用于在Juicebox Assembly Tools模块(由review.Assembly文件表示)中进行审阅,将程序集(由之前对齐的Hi-C reads和Juicer pipeline生成的)最终确定为染色体长度的fasta序列。该脚本将生成一个输出fasta文件、最终装配Hi-C map的assembly文件和一些补充注释文件,便于在Juicebox中查看结果。

1
2
3
4
5
6
7
8
9
USAGE: ./run-asm-pipeline-post-review.sh [options] -r <review.assembly> <path_to_input_fasta> <path_to_input_mnd>
# <review.assembly> 输入juicebox调整之后得到的genome.rawchrom.assembly
# <path_to_input_fasta> 输入参考基因组序列文件
# <path_to_input_mnd> 之前得到的merged_nodups.txt
# -r 查看“.assembly”文件的路径。
可选参数
# -i input_size 指定阈值输入Contigs/scaffolds(默认值为15000)。小于输入阈值的Contigs/scaffolds的将被忽略。应与运行原始脚本时使用的相同。
# -g gap_size 要在最终染色体长度支架中的支架序列之间添加的间隙大小(默认值为500)。
# --sort-output 选择按大小按降序排列染色体长度支架。

个人运行示例

1
nohup bash ~/software/3d-dna/run-asm-pipeline-post-review.sh -r genome.rawchrom.assembly genome.fa hic/aligned/merged_nodups.txt &> 3d.log 

得出最终的染色体水平文件 genome.FINAL.fasta

提升最后一步的速度

因为run-asm-pipeline-post-review.sh原始的速度太慢了。所以建议去修改一下源文件内容,大大提升最后一步的速度,可以参考链接