之前使用的是3D-DNA流程做Hi-C的輔助組裝矮台，它的最大優(yōu)勢(shì)就是輸出結(jié)果可以對(duì)接下游的JBAT(juicerbox with Assembly Tools)進(jìn)行手動(dòng)矯正乏屯。然而它點(diǎn)缺陷也很明顯，處理速度不夠快瘦赫，且對(duì)植物的優(yōu)化不行辰晕，同時(shí)目前許久不更新了。
最近我發(fā)現(xiàn)了一套新的組合确虱，chromap + yahs 完全替代之前3D-DNA流程含友。它的依賴工具如下
- chromap: 高效Hi-C數(shù)據(jù)比對(duì)
- samtools: sam轉(zhuǎn)bam
- yahs: 另一個(gè)Hi-C scaffolding工具。糾錯(cuò)上準(zhǔn)確性高校辩，排序上略強(qiáng)3d-dna窘问，遠(yuǎn)超SALSA2。
- juicer_tools: 用于輸出導(dǎo)入JuiceBox
chrompa, samtools, yahs可以直接用conda進(jìn)行安裝宜咒，juicer_tools依賴Java環(huán)境惠赫，并需要單獨(dú)下載
conda create -n hic-scaffolding -c bioconda -c conda-forge chromap samtools yahs samtools assembly-stats openjdk
# 1.19.02版本就行了, 最新的3.0不向下兼容
wget https://s3.amazonaws.com/hicfiles.tc4ga.com/public/juicer/juicer_tools_1.19.02.jar
具體分析步驟如下，我們需要提供前期組裝結(jié)果故黑，以及Hi-C數(shù)據(jù)
contigsFasta=/到/你的/contig.fa的路徑
r1Reads=/到/你的/Hi-C R1測(cè)序的路徑
r2Reads=/到/你的/Hi-C R2測(cè)序的路徑
第一步儿咱，數(shù)據(jù)比對(duì)
# 建立索引
samtools faidx $contigsFasta
chromap -i -r $contigsFasta -o contigs.index
# alignment
chromap \
--preset hic \
-r $contigsFasta \
-x contigs.index \
--remove-pcr-duplicates \
-1 $r1Reads \
-2 $r2Reads \
--SAM \
-o aligned.sam \
-t 50
# 排序
samtools view -bh aligned.sam | samtools sort -@ 50 -n > aligned.bam
rm aligned.sam
按照read的名字進(jìn)行排序和按照位置排序或未排序的結(jié)果會(huì)有一些不同庭砍。
第二步，又快又好的scaffolding混埠。默認(rèn)只需要兩個(gè)輸入逗威，組轉(zhuǎn)的contig.fa和比對(duì)的bam，和C語言一樣簡(jiǎn)潔岔冀。
yahs $contigsFasta aligned.bam
在輸出結(jié)果中
- inital_break 表示糾錯(cuò)的中間輸出
- _scaffolds_final.agp和_scaffolds_final.fa則是最終結(jié)果
對(duì)于輸出結(jié)果凯旭，我們希望進(jìn)行可視化，此時(shí)可以使用yahs提供的jucier工具
第三步使套，為juicer_tools準(zhǔn)備輸入
juicer pre -a -o out_JBAT \
yahs.out.bin \
yahs.out_scaffolds_final.agp \
$contigsFasta.fai
# -o out_JBAT表示輸出文件名的前綴
一共包括如下幾個(gè)文件
- out_JBAT.assembly
- out_JBAT.assembly.agp
- out_JBAT.hic
- out_JBAT.liftover.agp
- out_JBAT.txt
out_JBAT.txt就作為下游的輸入
JUICER=/路徑/到/juicer_tools_1.19.02.jar
asm_size=$(awk '{s+=$2} END{print s}' $contigsFasta.fai)
java -Xmx36G -jar $JUICER \
pre out_JBAT.txt out_JBAT.hic <(echo "assembly ${asm_size}")
輸出的out_JBAT.hic就可以導(dǎo)入到JBAT進(jìn)行組裝,導(dǎo)出為out_JBAT.review.assembly
將手動(dòng)修改的結(jié)果傳遞給juicer罐呼，進(jìn)行scaffolding。
juicer post -o out_JBAT out_JBAT.review.assembly out_JBAT.liftover.agp contigs.fa
輸出結(jié)果為 out_JBAT.FINAL.agp, out_JBAT.FINAL.fa
