De novo组装#07

2024-07-04 17:05| 来源: 网络整理| 查看: 265

写在前面

初始组装经过基因组纠错（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> 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的目录，为想要运行juicer的目录 ## 这个步骤后面再用的时候也还会再讲 cp 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 ## 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)

【本文地址】

De novo组装#07

De novo组装#07

今日新闻

推荐新闻