De novo组装#07 | 染色体挂载 （allhic，3d-
写在前面
初始组装经过基因组纠错（polish）以及去冗余（purge）之后，就可以将其挂载到染色体上，使其由contig/scaffold级别的基因组提升到chromosome级别的基因组。染色体挂载的方法有多种，我这里主要介绍基于HiC测序数据进行染色体挂载，用到了2套软件：allhic和3d-dna pipeline。
数据准备
- HiC双端二代测序数据:R1,R2
- 用于挂载的contig/scaffold级别基因组
allhic + juicebox
allhic为了组装多倍体甘蔗而开发的软件，适用于多倍体，基因组复杂，组装指标一般，序列条数较多的情况基因组，操作相对简单，把所有contig自动分组挂载到预设的30条染色体上。但挂载本身这个步骤相对其他步骤感觉还是要复杂一些，另外一般还需要用 juicebox手动调整一下自动挂载中的一些错误。
1.软件安装
allhic主页：https://github.com/tangerzhang/ALLHiC
Windows/Mac版juicebox下载地址：https://github.com/aidenlab/Juicebox/wiki/Download
## allhic安装
git clone https://github.com/tangerzhang/ALLHiC
cd ALLHiC
chmod +x bin/*
chmod +x scripts/*
export PATH=/your/path/to/ALLHiC/scripts/:/your/path/to/ALLHiC/bin/:$PATH ## 临时将这两个文件夹的脚本或程序添加环境变量
此外这软件还依赖其他软件：samtools v1.9+，bedtools，matplotlib v2.0+， 还有一些juicebox_scripts脚本用于格式转换
## samtools和bedtools用的比较多，应该都装了，这里装下matlock
## 自动安装
conda install -c bioconda matlock
## 手动安装
git clone --recursive https://github.com/phasegenomics/matlock.git
cd matlock
make
## juicebox_scripts下载即可使用，我们这里后续手动用juicebox调整时，会用到里面的agp2assembly.py脚本进行格式转换
git clone https://github.com/phasegenomics/juicebox_scripts.git
2.运行使用
allhic如果用于挂载多倍体基因组一般分为五步:pruning, partition, rescue, optimization, building。我这边用来挂二倍体基因组pruning和rescue步骤可不进行。
此外，allhic还有一个可选步骤：基于 hic 数据的比对情况，对基因组中潜在的组装错误进行打断，但该操作会明显降低基因组的连续性，可以先不做这步骤。如果最好挂载结果看到contig内部由有明显的hic信号错误，可以再来执行这一步骤。
0# 基因组打断（可选步骤）
## 把基因组和hic测序数据链接过来
ln -s /....../....../polished.purged.fasta seq.fasta
ln -s /....../....../clean.HiC_R1.fastq.gz HiC_R1.fastq.gz
ln -s /....../....../clean.HiC_R2.fastq.gz HiC_R2.fastq.gz
### 建立基因组index
bwa index seq.fasta
samtools faidx seq.fasta
### bwa将二代的hic数据比对到基因组上
bwa mem -SP5M -t 24 seq.fasta HiC_R1.fastq.gz HiC_R2.fastq.gz \
| samtools view -hF 256 - \
| samtools sort -@ 24 -o sorted.bam -T tmp.ali
samtools index sorted.bam
### 对contig的潜在组装错误进行打断
ALLHiC_corrector \
-m sorted.bam \
-r seq.fasta \
-o corrected.fasta \ ## 拿到一个打断纠错过的fasta
-t 24
如果不进行打断可从这步开始进行染色体挂载步骤，我是没有打算就直接开始的
1# index reference genome：建立基因组index
## 首先设置allhic的全局变量，重连服务器或者后台脚本运行需要重新运行一下这个
export PATH=/newlustre/home/jfgui/Wangtao/software/ALLHiC/scripts/:/newlustre/home/jfgui/Wangtao/software/ALLHiC/bin:$PATH
# 建立索引
ln -s ../../PacBio/polish/flye_polish/pilon.out/pilon3.out/racon3.pilon3.fasta assembly.fasta ## 还是把基因组链接过来
/newlustre/home/jfgui/Wangtao/software/bwa/bwa index assembly.fasta ## 建立基因组bwa比对的index
samtools faidx assembly.fasta ## 建立基因组序列提取的index，fai后缀
2# mapping：将二代的hic数据用bwa比对到基因组，并对并对文件过滤并排序
/newlustre/home/jfgui/Wangtao/software/bwa/bwa mem \ ## bwa mem比对
-SP5M -t 24 assembly.fasta ../Liver.clean.R1.fq.gz ../Liver.clean.R2.fq.gz \
| samtools view -hF 256 - \ ## 对比对文件bam过滤掉次比对数据
| samtools sort -@ 24 -o sorted.bam ## 对比对文件bam进行排序
3# filter bam：基于比对质量对bam文件进行过滤
samtools view -b -q 40 sorted.bam \ ## 过滤阈值为40
| samtools view -b -t assembly.fasta.fai > unique.bam ## -t 用上一步的fai索引文件生成header文件防止报错。
PreprocessSAMs.pl unique.bam assembly.fasta MboI ## 还是一个过滤过程，即对上一步bam再进行处理，保留酶切位点MboI 附近的比对数据。最后生成的bam文件：unique.REduced.paired_only.bam
4# partition：基于处理过后的bam文件，根据Hi-C建议的链接对链接的contigs进行聚类分组
ALLHiC_partition \
-r assembly.fasta \ ## reference genome
-e GATC \ # 酶切类型MboI
-k 30 \ ## 分组数，即染色体个数
-b unique.REduced.paired_only.bam ## 输入的bam文件，即上面几步处理过的
5# optimize：对组内的contigs进行定性排序
## for循环生成30个allhic optimize命令文件放在cmd.list里
rm cmd.list
for((K=1;K<=30;K++));do echo "allhic optimize unique.REduced.paired_only.counts_GATC.30g${K}.txt unique.REduced.paired_only.clm" >> cmd.list;done
## 用ParaFly对cmd.list里的30个命令进行并行运算，用conda安装ParaFly就行: conda install -c bioconda parafly
ParaFly -c cmd.list -CPU 24
6# build：最后连接contigs生成染色体级别的assembly（groups.asm.fasta ）
ALLHiC_build assembly.fasta
7# polt：画hic热图
samtools faidx groups.asm.fasta ## 对最终的fasta建立序列读取index，groups.asm.fasta.fai
cut -f1,2 groups.asm.fasta.fai > chrn.list ## 读取 groups.asm.fasta.fai里的第1/2列
ALLHiC_plot -b sorted.bam -a groups.agp -l chrn.list -s 50k -o heatmap-pdf ## 绘图
8# juicebox手动重新调整
因为第7步自动生成的hic热图很少能画得很完美，肯定多少会有点信号错误，这时需要手动调整，再导出图hic热图，这里用 juicebox。juicebox直接在windows系统里操作，需要两个输入文件：groups.assembly和groups.hic，
以下步骤用于生成这两个juicebox的输入文件：
## 基于read name 重新对bam排序
samtools sort -n -@ 24 -o aligned.sort_name.bam ../sorted.bam
## 用matlock将bam转换成merged_nodups.txt文件，即juicer软件产生的一种文件
matlock bam2 juicer aligned.sort_name.bam merged_nodups.txt
## 用agp2assembly.py脚本将agp 格式转为3d-dna的assembly 格式，即groups.assembly
/newlustre/home/jfgui/Wangtao/software/juicebox_scripts/juicebox_scripts/agp2assembly.py ../groups.agp groups.assembly
## 基于前两步生成的merged_nodups.txt和groups.assembly，生成用于juicebox输入的二进制hic文件，即groups.hic
bash /newlustre/home/jfgui/Wangtao/software/3d-dna/visualize/run-assembly-visualizer.sh -q 1 -p true groups.assembly merged_nodups.txt
用window版的juicebox输入groups.assembly和groups.hic，手动调整信号不对的地方，以下是我的简单调整过后的：
hic热图-简单调整即可划分出明显的30条染色体
该结果简单调整了下没有细调，也不准备调了，因为发现contig内部出现了明显的信号错误，这是刚开始用的基因组组装得就有问题，然后我又没有进行基因组打断步骤，所有后面准备用3d-dna执行打断这一步骤。
红色箭头的contig内部可以发现信号有问题（绿色框为contig，蓝色框为染色体）
juicer+ 3d-dna + juicebox
用3d-dna挂载这一套流程个人感觉与allhic差不多，可能还简单一点。大致就是把对hic数据的比对以及处理用另一个软件juicer来代替了， 3d-dna只用来挂载。最后也是用 juicebox手动再调整一下。
3d-dna挂载流程图
1.软件安装
juicer主页：https://github.com/aidenlab/juicer/
3d-dna主页：https://github.com/aidenlab/3d-dna
Windows/Mac版juicebox下载地址：https://github.com/aidenlab/Juicebox/wiki/Download
## juicer安装
git clone https://github.com/theaidenlab/juicer.git
## 后续要用这个软件需要将CPU目录下的程序拷到运行juicer的目录，<myJuicerDir> 为想要运行juicer的目录
## 这个步骤后面再用的时候也还会再讲
cp juicer/CPU/* <myJuicerDir>/scripts
cd scripts/common
wget 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
## 3d-dna安装
git clone https://github.com/aidenlab/3d-dna.git
3d-dna/run-asm-pipeline.sh –h
2.使用运行
如上面的流程图主要分三大步：juicer分析hic数据，3d-dna挂载染色体，juicebox手动调整&重新生成final文件。
- juicer分析hic数据
1# 构建基因组index，以及基因组内各contig长度文件
ln -s ../../../PacBio/polish/flye_polish/pilon.out/pilon3.out/racon3.pilon3.fasta genome.fa # 将polish 过的基因组链接过来
/newlustre/home/jfgui/Wangtao/software/bwa/bwa index genome.fa # 建立基因组index
seqkit fx2tab -l -n -i genome.fa > chr.size # 使用seqkit fx2tab -l -n -i 统计各contig的长度
2# 准备基因组内可能的酶切位点文件，用到的脚本为：juicer/misc/generate_site_positions.py
/newlustre/home/jfgui/Wangtao/software/juicer/misc/generate_site_positions.py MboI genome genome.fa # 后接酶的类型，输出前缀，以及基因组文件
3# 准备juicer script目录，上面安装的时候提到了，后面运行juicer就是调用这里面的程序。
mkdir scripts
cp -r /newlustre/home/jfgui/Wangtao/software/juicer/CPU/* ./scripts
cd scripts/common/
wget 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
cd ../../
4# 准备HiC数据文件目录
mkdir fastq
cd fastq
ln -s /newlustre/home/jfgui/Wangtao/XX_assembly/HiC/Liver.clean.R1.fq.gz HiC_R1.fastq.gz
ln -s /newlustre/home/jfgui/Wangtao/XX_assembly/HiC/Liver.clean.R2.fq.gz HiC_R2.fastq.gz
cd ..
5#上面所有的文件都准备好了之后，开始用运行juicer。
export PATH=/newlustre/home/jfgui/Wangtao/software/bwa: \ # juicer会调用bwa对hic测序reads进行比对，这里添加bwa的全局环境变量
/newlustre/home/jfgui/Wangtao/XX_assembly/HiC/3d-dna/juicer1/scripts: \
/newlustre/home/jfgui/Wangtao/XX_assembly/HiC/3d-dna/juicer1/scripts/common:$PATH
# 运行juicer
bash ./scripts/juicer.sh \
-y genome_MboI.txt \ # 酶切位点文件
-z genome.fa \ # 基因组文件
-s MboI \ # hic测序酶的类型
-p chr.size \ # contig长度信息文件
-D ./ \ # scripts 和 fastq文件所在目录
-t 24 \ # 线程数
--assembly \ # 不是很理解，可能就是生成merged_nodups.txt文件，帮助文档里为： For use before 3D-DNA; early exit and create old style merged_nodups
&> juicer.log # 运行时间长就还是加一个日志信息，之前因为bwa不在全局变量里导致没出结果，后面看日志文件才发现了问题。
最后会在aligned目录下生成merged_nodups.txt，用于后续3d-dna的输入。
- 3d-dna挂载染色体
## 同样地，先把基因组链接过来
ln -s ../../../PacBio/polish/flye_polish/pilon.out/pilon3.out/racon3.pilon3.fasta genome.fa
## 由于本人运行3d-dna时报错，提示存储临时文件的空间不够了，可能服务器用的人当时比较多，系统默认的存储临时文件目录满了，这里的命令把临时文夹放到当前目录。
export TMPDIR=/newlustre/home/jfgui/Wangtao/XX_assembly/HiC/3d-dna/3d-dna/
## 运行3d-dna挂载程序
/newlustre/home/jfgui/Wangtao/software/3d-dna/run-asm-pipeline.sh \ # 用3d-dna目录的run-asm-pipeline.sh脚本
-r 2 \ # 纠错次数，0表示不纠错，默认为2。对组装结果自信先用-r 0 试一下。
./genome.fa \ # 基因组文件
../juicer/aligned/merged_nodups.txt \ # juicer运行结果文件
&> 3d-dna.log #追加日志信息
输出的文件较多，其中genome.final.assembly 和genome.final.hic 这两个final文件用于juicebox 输入。而genome.FINAL.fasta 为3d-dna 的自动挂载的基因组序列文件，这个文件后续在用juicebox手动调整后可以重新生成覆盖。
- juicebox手动调整&重新生成final文件
1# juicebox手动调整
将genome.final.assembly 和genome.final.hic下载到本地，输入到juicebox进行手动简单调整了下，大致结果如下
3d-dna自动挂载后的结果，这里没加contig和scaffold边框线
图上没加contig和scaffold边框线，但我自己仔细看了下contig内部确实没有什么冲突的地方，因为我在运行3d-dna的时候设置了
-r 2参数进行了2轮纠错打断，所以把第一种策略allhic发现冲突的那个contig给修正回来了，代价是基因组更碎了更不连续😂。
3d-dna自动挂载后的结果，这里没加contig和scaffold边框线
此外，我们在右下角看到特别多没有挂载到染色体上的一些contig或者说scaffold，其出现这么多进行2轮基因组纠错打断是一个原因，更多的是因为用来的挂载的这个基因组
../../../PacBio/polish/flye_polish/pilon.out/pilon3.out/racon3.pilon3.fasta是没有用purge_dups软件去过冗余的版本，这些序列大部分应该都是一些重复冗余序列。因此在进行染色体挂载之前还是先进行基因组去冗余操作。
2# 重新生成final文件
juicebox 手动调整后，重新生成assemby文件review.assembly ，上传到服务器使用run-asm-pipeline-post-review.sh程序生成最终的 fasta 文件和hic 文件。注意，最终的genome.FINAL.fasta文件是按大到小排序过，且用500N填充过gap的fasta文件。
/newlustre/home/jfgui/Wangtao/software/3d-dna/run-asm-pipeline-post-review.sh \
-i # 去掉15k以下的scaffold
--build-gapped-map \ # 建一个互作图？这个map是用500个N填充完gap的
--sort-output \ # 按大小降序排列挂载好的染色体级别的scaffold
-r review.assembly \ # juicebox输出的assembly文件
../../data/genome.fa merged_nodups.txt # 后面接原始的基因组文件以及jucer输出的merged_nodups.txt文件
3.查看结果
用assembly-stats软件查看了下最终的挂载结果:
其他相关推荐
使用ALLHiC基于HiC数据辅助基因组组装 - 简书 (jianshu.com)
基于3D-DNA，ALLHiC挂载二倍体基因组 - 简书 (jianshu.com)
利用3D-DNA挂载基因组 - 简书 (jianshu.com)
3d-DNA的使用及juicebox调整挂载到染色体水平 ｜ HiC辅助基因组组装（二） | 生信技术 (lxz9.com)
