使用purge_haplogs处理基因组杂合区域

  1. 云栖社区>
  2. 博客>
  3. 正文

使用purge_haplogs处理基因组杂合区域

徐洲更 2018-10-07 16:24:00 浏览1095
展开阅读全文

基因组某些区域可能有着比较高的杂合度,这会导致基因组该区域的两个单倍型被分别组装成primary contig, 而不是一个为primary contig, 另一个是associated haplotig. 如果下游分析主要关注于单倍型,这就会导致一些问题。

那么有没有解决方案呢?其实也很好办,就是找到相似度很高的contig,将他们拆分。 purge_haplogs根据minimap2的比对结果,通过分析比对read的覆盖度决定谁去谁留。该工具适用于单倍型组装软件,例如 Canu, FALCON或 FALCON-Unzip primary contigs, 或者是分相后的二倍体组装(Falcon-Unzip primary contigs + haplotigs 。

软件安装

purge_haplotigs依赖软件比较多,手动安装会很麻烦,但是他可以直接用bioconda装

conda create -n purge_haplotigs_env
conda activate purge_haplotigs_env
conda install purge_haplotigs

安装完成后需要一步测试

purge_haplotigs test

简明教程

数据准备。 需要下载的数据集分为两个部分,一个是FALCON-Unzip后的primary contig 和 halplotigs. 另一个则是已经比完后的BAM文件

mkdir purge_haplotigs_tutorial
cd purge_haplotigs_tutorial
wget https://zenodo.org/record/841398/files/cns_h_ctg.fasta
wget https://zenodo.org/record/841398/files/cns_p_ctg.aligned.sd.bam # 1.7G
wget https://zenodo.org/record/841398/files/cns_p_ctg.aligned.sd.bam.bai 
 wget https://zenodo.org/record/841398/files/cns_p_ctg.fasta
wget https://zenodo.org/record/841398/files/cns_p_ctg.fasta.fai

当然我们不可能直接就拿到比对好的BAM文件,我们一般是有组装后的基因组以及用于组装的subread,假设这两个文件命名为, genome.fa 和 subreads.fasta.gz.

minimap2 -ax map-pb genome.fa subreads.fasta.gz \
    | samtools view -hF 256 - \
    | samtools sort -@ 8 -m 1G -o aligned.bam -T tmp.ali

如果你有二代测序数据,也可以用BWA-MEM进行比对得到BAM文件。

第一步:使用purge_haplotigs readhist从BAM中统计read深度,绘制柱状图。

samtools mpileup -r "000005F|quiver" -f cns_p_ctg.fasta cns_p_ctg.aligned.sd.bam

也就是下图,你能明显的看到图中有两个峰,一个是单倍型的覆盖度,另一个二倍型的覆盖度,

img_33ea28b3f2acdf0bf04f2807deab56e8.png
高杂合基因组read-depth histogram

你可能还想知道高纯合基因组是什么样的效果,我也找了一个纯合的物种做了也做了read-depth 柱状图,

img_2e9837fe80c993532740cd0b64e912c3.png
纯合基因组read-depth histogram

之后你需要根据read-depth 柱状图 确定这两个峰的位置用于下一步。下面是两个例子。对于我们则是,20,65,190.

img_122689ec21fb1e1efefecf537e81fb13.png
两个例子

第二步: 根据read-depth信息选择阈值。

purge_haplotigs  contigcov  -i cns_p_ctg.aligned.sd.bam.gencov  -o coverage_stats.csv  -l 20  -m 75  -h 190

这一步生成的文件是"coverage_stats.csv"

第三步:区分haplotigs.

purge_haplotigs purge  -g cns_p_ctg.fasta  -c coverage_stats.csv  -b cns_p_ctg.aligned.sd.bam  -t 4  -a 60

这一步会得到如下文件

  • curated.artefacts.fasta:无用的contig,也就是没有足够覆盖度的contig.
  • curated.fasta:新的单倍型组装
  • curated.haplotigs.fasta:从原本组装分出来的haplotigs
  • curated.reassignments.tsv: 单倍型的分配信息
  • curated.contig_associations.log: 运行日志, 下面是其中一个记录,表示000004F_004和000004F_027是000004F_017的HAPLOTIG, 而000004F_017和000004F_013又是000004F,的HAPLOTIG。
000004F,PRIMARY -> 000004F_013,HAPLOTIG
                -> 000004F_017,HAPLOTIG 
                                        -> 000004F_004,HAPLOTIG
                                        -> 000004F_027,HAPLOTIG

由于我们用的是单倍型组装primary contigs而不是二倍体组装的parimary + haplotigs, 因此我们需要将FALCON_Unzip的haplotgi合并到重新分配的haplotigs中,这样子我们依旧拥有二倍体组装结果

cat cns_h_ctg.fasta >> curated.haplotigs.fasta

网友评论

登录后评论
0/500
评论
徐洲更
+ 关注