關(guān)于3D-DNA，前面已經(jīng)有簡單介紹，本次主要介紹ALLHiC分析流程
ALLHiC的流程主要分為以下5個步驟：
- Prune: 修剪（二倍體基因組可以不用）
- Partition: 根據(jù)Hic對contig進行分組
- Rescue： 將未分組的contig繼續(xù)進行分組
- Optimize: 確定每組的順序和方向
- Build：得到fasta序列和AGP文件
軟件安裝
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+
conda create -y -n allhic python=3.7 samtools bedtools matplotlib
大概流程
若發(fā)現(xiàn)組裝的scaffolds存在很多misjoin，則可以通過3D-DNA進行糾錯，然后使用ALLHiC。
關(guān)于3D-DNA流程，可查看前面內(nèi)容，3D-DNA可以得到
- 準(zhǔn)備輸入文件
perl software/ALLHiC/scripts/release3DDNA.pl 16 seq.FINAL.fasta
上述將3D-DNA組裝好的結(jié)果根據(jù)片段長度從大到小排列，得到相應(yīng)的chr，而后將chr中的NN將其分割成多個tig.correced.contig，并得到相應(yīng)chr.tour 文件。最后用ALLHIC_build根據(jù)tig.correced.contig序列以及tour文件構(gòu)建groups.agp文件，因為是根據(jù)N進行分割的，所以agpW和U一次排列
得到作為ALLHiC輸入文件
- 建索引
samtools faidx tig.HiCcorrected.fasta
bwa index -a bwtsw tig.HiCcorrected.fasta
- HiC reads回帖基因組
bwa aln -t 24 tig.HiCcorrected.fasta reads_R1.fastq.gz > sample_R1.sai
bwa aln -t 24 tig.HiCcorrected.fasta reads_R2.fastq.gz > sample_R2.sai
bwa sampe tig.HiCcorrected.fasta sample_R1.sai sample_R2.sai reads_R1.fastq.gz reads_R2.fastq.gz > sample.bwa_aln.sam
- 過濾SAM文件
酶切位點根據(jù)自己實驗進行選擇
PreprocessSAMs.pl sample.bwa_aln.sam tig.HiCcorrected.fasta HindIII
#如果有BAM文件
#PreprocessSAMs.pl sample.bwa_aln.bam tig.HiCcorrected.fasta HindIII
filterBAM_forHiC.pl sample.bwa_aln.REduced.paired_only.bam sample.clean.sam
samtools view -bt tig.HiCcorrected.fasta.fai sample.clean.sam > sample.clean.bam
其中
- 將500bp以內(nèi)沒有酶切位點的reads全部過濾；
- 保留雙端均比對上的reads，得到*.paired_only.bam 文件
- 得到reads比對情況 ..flagstat文件
的過濾標(biāo)準(zhǔn):
- 比對質(zhì)量高于30(MQ)
- 只保留唯一比對(XT:A:U)
- 編輯距離(NM)低于5
- 錯誤匹配低于(XM)4
- 不能有超過2個的gap(XO,XG)
-
Partition
根據(jù)預(yù)先定的組（本項目選用16個組）對contig進行分組
ALLHiC_partition -b sample.clean.bam -r tig.HiCcorrected.fasta -e AAGCTT -k 16
#參數(shù)
-h 幫助信息
-b prunned bam（跳過prune，直接用sample.clean.bam)
-r 組裝的基因組
-e 酶切位點 (HindIII: AAGCTT; MboI: GATC)
-k 組的數(shù)量 (根據(jù)物種染色體定)
-m 最少的酶切位點，默認(rèn)25
-
Optimize
確定每個組contig的順序和方向
# 提取CLM和每個contig的酶切數(shù)
allhic extract sample.clean.bam tig.HiCcorrected.fasta--RE AAGCTT
--RE: 默認(rèn)GATC
# 確定contig順序和方向 (可并行)
allhic optimize sample.clean.counts_AAGCTT.16g1.txt sample.clean.clm
allhic optimize sample.clean.counts_AAGCTT.16g2.txt sample.clean.clm
...
allhic optimize sample.clean.counts_AAGCTT.16g16.txt sample.clean.clm
- 獲取染色體級別組裝
ALLHiC_build tig.HiCcorrected.fasta
- 繪制熱圖
# 獲取染色體長度
perl getFaLen.pl -i groups.asm.fasta -o len.txt
# 構(gòu)建chrn.list
grep 'merge.clean.counts_GATC' len.txt > chrn.list
# 畫圖（對染色體畫，500K）
ALLHiC_plot sample.clean.bam groups.agp chrn.list 500k pdf
畫圖腳本即根據(jù)bam文件獲取reads比對到的contig位置并進行計數(shù)，根據(jù)AGP文件將contig對應(yīng)在染色體上進行作圖
腳本戳這里getFaLen.pl
