使用二代數(shù)據(jù)或三代數(shù)據(jù)得到contig后菜枷，下一步就是將contig提升到染色體水平。對(duì)于二倍體物種而言叁丧，目前3D-DNA應(yīng)該是組裝效果最好的一個(gè)軟件啤誊。
一：工作流程
使用3D-DNA做基因組組裝的整體流程如下圖，分別為組裝拥娄，Juicer分析Hi-C數(shù)據(jù)蚊锹，3D-DNA進(jìn)行scaffolding，使用JBAT對(duì)組裝結(jié)果進(jìn)行手工糾正稚瘾，最終得到準(zhǔn)染色體水平的基因組牡昆。
二： 軟件安裝
在安裝之前，確保服務(wù)器上有了下面這些依賴軟件工具
- LastZ（僅在雜合基因組的二倍體模式下使用）
- Java >= 1.7
- GNU Awk >= 4.02
- GNU coreutils sort > 8.11
- Python >= 2.7
- scipy, numpy, matplotlib
- GNU Parallel >=20150322 (不必要摊欠，但是強(qiáng)力推薦)
- bwa
我們需要安裝兩個(gè)軟件丢烘，一個(gè)是3D-DNA柱宦，另一個(gè)是juicer。
CPU版本的juicer安裝
mkdir ~/soft
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
然后用~/soft/juicer/scripts/juicer.sh -h檢查是否有幫助信息輸出播瞳。
3D-DNA安裝
cd ~/soft
git clone https://github.com/theaidenlab/3d-dna.git
用sh ~/soft/3d-dna/run-asm-pipeline.sh -h查看是否有幫助文檔輸出捷沸。
參數(shù)詳解
以CPU版本的為例，juicer.sh的參數(shù)如下
Usage: juicer.sh [-g genomeID] [-d topDir] [-s site] [-a about] [-R end]
[-S stage] [-p chrom.sizes path] [-y restriction site file]
[-z reference genome file] [-D Juicer scripts directory]
[-b ligation] [-t threads] [-r] [-h] [-f] [-j]
參數(shù)說(shuō)明：
- -g: 定義一個(gè)物種名
- -s: 酶切類(lèi)型, HindIII(AAGCTAGCTT), MboI(GATCGATC) , DpnII(GATCGATC), NcoI(CCATGCATGG)
- -z : 參考基因組文件
- -y: 限制性酶切位點(diǎn)可能出現(xiàn)位置文件
- -p: 染色體大小文件
- -C: 將原來(lái)的文件進(jìn)行拆分狐史，必須是4的倍數(shù)痒给，默認(rèn)是90000000, 即22.5M reads
- -S: 和任務(wù)重運(yùn)行有關(guān)，從中途的某一步開(kāi)始,"merge", "dedup", "final", "postproc" 或 "early"
- -d: juicer的目錄 我們安裝在/soft/骏全，所以設(shè)置為/soft/juicer
- -D: juicer scripts的目錄苍柏，我們安裝在/soft/，所以設(shè)置為/soft/juicer/CPU
- -a: 實(shí)驗(yàn)的描述說(shuō)明姜贡，可以不用設(shè)置
- -t: 線程數(shù)
如果你的基因組不是復(fù)雜基因組试吁，比如說(shuō)高雜合，高重復(fù)序列楼咳，或者Hi-C數(shù)據(jù)測(cè)太少熄捍，那么3d-dna的流程更加簡(jiǎn)單, run-asm-pipeline.sh -h只有四個(gè)參數(shù)需要改:
- -i|--input: 過(guò)濾長(zhǎng)度低于給定閾值的contig/scaffold, 默認(rèn)是15000
- -r|--round: 基因組中misjoin的糾錯(cuò)輪數(shù)，默認(rèn)是2母怜，當(dāng)基因組比較準(zhǔn)確時(shí)余耽，設(shè)置為0，然后在JABT中調(diào)整會(huì)更好
- -m|--mode: 是否調(diào)用merge模塊苹熏，當(dāng)且僅當(dāng)在雜合度比較高的情況下使用碟贾，也就是組裝的單倍型基因組明顯偏大
- -s|--stage: 從polish, split, seal, merge 或finalize 的某一個(gè)階段開(kāi)始
但是，一旦基因組復(fù)雜起來(lái)轨域，那么需要調(diào)整的參數(shù)就非常多了, run-asm-pipeline.sh --help會(huì)輸出更多的信息袱耽，你需要根據(jù)當(dāng)前結(jié)果去確定每個(gè)階段的參數(shù)應(yīng)該如何調(diào)整。
最終的輸出文件最關(guān)鍵的是下面三類(lèi):
- .fasta: 以FINAL標(biāo)記的是最終結(jié)果
- .hic: 各個(gè)階段都會(huì)有輸出結(jié)果干发，用于在JABT中展示
- .assembly: 各個(gè)階段都會(huì)有輸出朱巨，一共兩列，存放contig的組裝順序
三：分析流程
準(zhǔn)備兩個(gè)數(shù)據(jù):
- reference:存放一個(gè)genome.fa, 為組裝的contigs枉长。
- fastq: 存放HiC二代雙端測(cè)序結(jié)果冀续，read_R1_fastq.gz, read_R2_fastq.gz
有了這兩個(gè)數(shù)據(jù)就可以開(kāi)始了。
1. 準(zhǔn)備一個(gè)新的基因組
第一步：為基因組建索引
bwa index genome.fa
第二步: 根據(jù)基因組構(gòu)建創(chuàng)建可能的酶切位點(diǎn)文件
python /data1/spider/ytbiosoft/soft/juicer/misc/generate_site_positions.py DpnII genome genome.fa
第三步: 運(yùn)行如下命令, 獲取每條contig的長(zhǎng)度
awk 'BEGIN{OFS="\t"}{print $1, $NF}' genome_DpnII.txt > genome.chrom.sizes
2.運(yùn)行juicer
保證當(dāng)前目錄下有fastq和reference文件夾搀暑，然后運(yùn)行如下命令沥阳，一定要設(shè)置-z,-p,-y這三個(gè)參數(shù)
bash /data1/spider/ytbiosoft/soft/juicer/scripts/juicer.sh -g Arg_aqu -d /data1/spider/ytbiosoft/soft/juicer -D /data1/spider/ytbiosoft/soft/juicer/CPU -z /data1/spider/ytbiosoft/test/hic_aqu/reference/genome.fa -y /data1/spider/ytbiosoft/test/hic_aqu/reference/genome_DpnII.txt -p /data1/spider/ytbiosoft/test/hic_aqu/reference/genome.chrom.size -s DpnII -t 36
輸出的結(jié)果文件都在aligned目錄下跨琳，其中"merged_nodups.txt"就是下一步3D-DNA的輸入文件之一自点。
3. 運(yùn)行3d-dna
3d-dna的運(yùn)行也沒(méi)有多少參數(shù)可以調(diào)整，如果對(duì)組裝基因組質(zhì)量的信心高脉让，就用-r 0, 否則用默認(rèn)的-r 2就行了桂敛。
~/soft/3d-dna/run-asm-pipeline.sh -r 2 reference/genome.fa aligned/merged_nodups.txt &> 3d.log &
然后在Juicer-Tools中對(duì)結(jié)果進(jìn)行可視化功炮，對(duì)可能的錯(cuò)誤進(jìn)行糾正。
最后輸出文件中术唬，包含F(xiàn)INAL就是我們需要的結(jié)果薪伏。
4.使用juicerbox進(jìn)行手工糾錯(cuò)
關(guān)于juicerbox的用法，可以看hoptop的https://www.bilibili.com/video/av65134634
最常見(jiàn)的幾種組裝錯(cuò)誤:
- misjoin: 切割
- translocations: 移動(dòng)
- inversions: 翻轉(zhuǎn)
- chromosome boundaries: 確定染色體的邊界
這些錯(cuò)誤的判斷依賴于經(jīng)驗(yàn)粗仓，所以只能靠自己多試試了嫁怀。
最后輸出genome.review.assembly用于下一步的分析
5. 再次運(yùn)行3d-dna
根據(jù)JABT手工糾正的結(jié)果, genome.review.assembly, 使用run-asm-pipeline-post-review.sh重新組裝基因組。
~/soft/3d-dna/run-asm-pipeline-post-review.sh \
-r genome.review.assembly genome.fa aligned/merged_nodups.txt &> 3d.log &
最后
假如你不小心設(shè)置了錯(cuò)誤的-p參數(shù)借浊，也不是特別的要緊塘淑，因?yàn)橹笤谧詈箅A段（final） 才會(huì)遇到了下面這個(gè)報(bào)錯(cuò)。
Could not find chromosome sizes file for: reference/genome.chrom.size
***! Can't find inter.hic in aligned/inter_30.hic
***! Error! Either inter.hic or inter_30.hic were not created
Either inter.hic or inter_30.hic were not created. Check aligned for results
即便遇到了這個(gè)報(bào)錯(cuò)也不要緊蚂斤，因?yàn)閕nter.hic 和 inter_30.hic在3d-dna流程中用不到存捺，所以不需要解決。
如果需要解決的話曙蒸，有兩個(gè)解決方案捌治，一種重新運(yùn)行命令，只不過(guò)多加一個(gè)參數(shù)-S final, 就會(huì)跳過(guò)之前的比對(duì)纽窟，合并和去重步驟肖油，直接到后面STATISTICS環(huán)節(jié)。但是這樣依舊會(huì)有一些不必要的計(jì)算工作臂港，所以另一種方法就是運(yùn)行原腳本必要的代碼
juiceDir=~/soft/juicer
outputdir=aligned
genomePath=reference/genome.chrom.size
site_file=reference/genome_DpnII.txt
ligation=GATCGATC
# output is inter.hic
${juiceDir}/scripts/common/juicer_tools pre -f $site_file -s $outputdir/inter.txt -g $outputdir/inter_hists.m -q 1 $outputdir/merged_nodups.txt $outputdir/inter.hic $genomePath
# output is inter_30.txt
${juiceDir}/scripts/common/statistics.pl -s $site_file -l $ligation -o $outputdir/inter_30.txt -q 30 $outputdir/merged_nodups.txt
# output is inter_30.hic
${juiceDir}/scripts/common/juicer_tools pre -f $site_file -s $outputdir/inter_30.txt -g $outputdir/inter_30_hists.m -q 30 $outputdir/merged_nodups.txt $outputdir/inter_30.hic $genomePath
參考
http://www.reibang.com/p/0969e1001d3b
https://github.com/theaidenlab/3d-dna
https://github.com/aidenlab/juicer
