基于3D-DNA，ALLHiC挂载二倍体基因组
2020-05-04 本文已影响0人
斩毛毛
关于3D-DNA，前面已经有简单介绍，本次主要介绍ALLHiC分析流程
ALLHiC的流程主要分为以下5个步骤：
- Prune: 修剪（二倍体基因组可以不用）
- Partition: 根据Hic对contig进行分组
- Rescue： 将未分组的contig继续进行分组
- 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
大概流程
若发现组装的scaffolds存在很多misjoin，则可以通过3D-DNA进行纠错，然后使用ALLHiC。
关于3D-DNA流程，可查看前面内容，3D-DNA可以得到
- 准备输入文件
perl software/ALLHiC/scripts/release3DDNA.pl 16 seq.FINAL.fasta
上述将3D-DNA组装好的结果根据片段长度从大到小排列，得到相应的chr，而后将chr中的NN将其分割成多个tig.correced.contig，并得到相应chr.tour 文件。最后用ALLHIC_build根据tig.correced.contig序列以及tour文件构建groups.agp文件，因为是根据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文件
酶切位点根据自己实验进行选择
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以内没有酶切位点的reads全部过滤；
- 保留双端均比对上的reads，得到*.paired_only.bam 文件
- 得到reads比对情况 ..flagstat文件
的过滤标准:
- 比对质量高于30(MQ)
- 只保留唯一比对(XT:A:U)
- 编辑距离(NM)低于5
- 错误匹配低于(XM)4
- 不能有超过2个的gap(XO,XG)
-
Partition
根据预先定的组（本项目选用16个组）对contig进行分组
ALLHiC_partition -b sample.clean.bam -r tig.HiCcorrected.fasta -e AAGCTT -k 16
#参数
-h 帮助信息
-b prunned bam（跳过prune，直接用sample.clean.bam)
-r 组装的基因组
-e 酶切位点 (HindIII: AAGCTT; MboI: GATC)
-k 组的数量 (根据物种染色体定)
-m 最少的酶切位点，默认25
-
Optimize
确定每个组contig的顺序和方向
# 提取CLM和每个contig的酶切数
allhic extract sample.clean.bam tig.HiCcorrected.fasta--RE AAGCTT
--RE: 默认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
# 构建chrn.list
grep 'merge.clean.counts_GATC' len.txt > chrn.list
# 画图（对染色体画，500K）
ALLHiC_plot sample.clean.bam groups.agp chrn.list 500k pdf
画图脚本即根据bam文件获取reads比对到的contig位置并进行计数，根据AGP文件将contig对应在染色体上进行作图
脚本戳这里getFaLen.pl
