基因组去重新方法解说 #04：purge_dups操作指南

写在前面

冗余序列的产生和多种因素有关，如 CLR 的测序错误，基因组自身的杂合性和重复序列的影响等等，purge_dups软件能根据read深度分析组装中haplotigs和overlaps（purge_dups is designed to remove haplotigs and contig overlaps in a de novo assembly based on read depth.）另外有时候我们又担心会过度 purge，所以purge的标准是什么呢？看了一个人这么说的：达到预估基因组大小，再进行BUSCUO评估，BUSCUO评估值又没有下降很多。软件主页的问答是这样的：
Q2: How can I validate the purged assembly? Is it clean enough or overpurged?
A2: There are many ways to validate the purged assembly. One way is to make a coverage plot for it which can also be hist_plot.py, the 2nd way is to run BUSCO and another way is to make a KAT plot with KAT (https://github.com/TGAC/KAT) or KMC (https://github.com/dfguan/KMC, use this if you only have a small memory machine) if short reads or some accurate reads are available.

工作流程

purge_dups软件的工作流程如上图分为三部分：1.首先将用于组装的三代测序数据比对到primary-assembly以计算基因组各部分覆盖度；2.primary-assembly自身的比对；3.将前两部分信息整合来推断出primary-assembly里的haplotigs和overlaps.

软件安装

purge_dups软件主页：https://github.com/dfguan/purge_dups

## purge_dups安装
git clone https://github.com/dfguan/purge_dups.git
cd purge_dups/src && make

Dependencies另外还有一个依赖软件需要安装，三代对比软件minimap2

## minimap2自动安装
conda install -c bioconda minimap2

## minimap2手动安装
wget https://github.com/lh3/minimap2/releases/download/v2.26/minimap2-2.26_x64-linux.tar.bz2
tar -jxvf minimap2-2.26_x64-linux.tar.bz2
./minimap2-2.26_x64-linux/minimap2

数据准备

这次的测试数据还是用的前面用4款软件分别组装出的4种primary-assembly，以及我们用于组装的CLR测序数据。
De novo组装#03 | 基因组拼接（flye, wtdbg2, mecat2, canu） - 简书 (jianshu.com)
下面以flye组装出的基因组为例，准备工作路径和测试数据

## 新建purge_dups/flye_purge目录作为purge工作目录
mkdir purge_dups
cd purge_dups
mkdir flye_purge
cd flye_purge
cp ../../flye.out/assembly.fasta ./  ## 将组装好的基因组拷贝到当前位置

软件运行

下面以flye组装出的基因组为例，其他软件都将以同样方式purge。

1#. minimap2将用于组装的三代测序数据比对到primary-assembly以计算基因组各部分覆盖度。（Run minimap2 to align pacbio data and generate paf files, then calculate read depth histogram and base-level read depth.）

## For PacBio CLR reads，用minimap2比对到基因组并压缩
/newlustre/home/jfgui/wy/anaconda3/envs/pengfang/bin/minimap2 \
      -x map-pb ./assembly.fasta ../../P01TYD20308306-1_r64030_20201110_065049_1_A02.subreads.fasta \ 
     | gzip -c - > pb_aln.paf.gz

## 统计paf, 输出PB.base.cov和PB.stat文件
/newlustre/home/jfgui/Wangtao/software/purge_dups/bin/pbcstat pb_aln.paf.gz   # produces PB.base.cov and PB.stat files
/newlustre/home/jfgui/Wangtao/software/purge_dups/bin/calcuts PB.stat > cutoffs 2>calcults.log  # bin/calcuts计算分界点

2#. 将初步组装基因组从N处进行打断，如果congtig中间没有N就不会被打断，然后同样是minimap2进行contig间的自身比对。( Split an assembly and do a self-self alignment. )

## 基因组拆分打断
/newlustre/home/jfgui/Wangtao/software/purge_dups/bin/split_fa ./assembly.fasta > assembly.fasta.split

## 基因组自身比对
/newlustre/home/jfgui/wy/anaconda3/envs/pengfang/bin/minimap2  \
     -xasm5 -DP assembly.fasta.split assembly.fasta.split  \
     | gzip -c - > assembly.fasta.split.self.paf.gz

3#.整合前两步得到的覆盖度及自我比对结果来对contig进行分类( Purge haplotigs and overlaps. )

## 分类信息在dups.bed文件里
/newlustre/home/jfgui/Wangtao/software/purge_dups/bin/purge_dups  \
-2 -T cutoffs -c PB.base.cov  assembly.fasta.split.self.paf.gz > dups.bed  2> purge_dups.log

4#.最后就是根据dups.bed文件里的分类信息，从原始的基因组种提取出Purged基因组( PGet purged primary and haplotig sequences from draft assembly. )

## 这一步输出两个结果 purged.fa 和 hap.fa，前者就是我们的最终结果
/newlustre/home/jfgui/Wangtao/software/purge_dups/bin/get_seqs  dups.bed  assembly.fasta

结果查看

我们用assembly-stats查看purge前后flye基因组的变化**
从下面的结果可以看出滤掉了336个contig共7Mb左右，似乎使得基因组的连续性变好了点。

##  purge前
$assembly-stats assembly.fasta
stats for assembly.fasta
sum = 771989320, n = 680, ave = 1135278.41, largest = 25186124
N50 = 9211426, n = 27
N60 = 7566198, n = 36
N70 = 6405554, n = 47
N80 = 5050486, n = 61
N90 = 2266136, n = 81
N100 = 493, n = 680
N_count = 0
Gaps = 0
*********************************************************************************************************************************************
##  purge后
$assembly-stats purged.fa
stats for purged.fa
sum = 764396646, n = 344, ave = 2222083.27, largest = 25186124
N50 = 9261417, n = 26
N60 = 7620965, n = 35
N70 = 6618835, n = 46
N80 = 5090803, n = 59
N90 = 3220038, n = 78
N100 = 611, n = 344
N_count = 0
Gaps = 0

下面我们也查看下其他版本的基因组purge前后变化:

wtdbg2版本基因组
过滤掉了近2000个cotig共45M左右

##  purge前
$assembly-stats wtdbg2_assembly.fa
stats for wtdbg2_assembly.fa
sum = 796796038, n = 2697, ave = 295437.91, largest = 24257430
N50 = 8405757, n = 33
N60 = 6518067, n = 44
N70 = 4827838, n = 58
N80 = 2584532, n = 80
N90 = 370995, n = 162
N100 = 2250, n = 2697
N_count = 0
Gaps = 0
*********************************************************************************************************************************************
##  purge后
$assembly-stats purged.fa
stats for purged.fa
sum = 751075917, n = 713, ave = 1053402.41, largest = 24257430
N50 = 8944165, n = 30
N60 = 6929004, n = 40
N70 = 5582188, n = 52
N80 = 3717695, n = 68
N90 = 1136788, n = 102
N100 = 2250, n = 713
N_count = 0
Gaps = 0

mecat2版本基因组
过滤掉了近500多个contig共24M左右

##  purge前
$assembly-stats contigs.fasta
stats for contigs.fasta
sum = 782440466, n = 1210, ave = 646645.01, largest = 20039767
N50 = 6788542, n = 36
N60 = 5033148, n = 50
N70 = 3470086, n = 68
N80 = 1490280, n = 101
N90 = 380864, n = 200
N100 = 641, n = 1210
N_count = 0
Gaps = 0
*********************************************************************************************************************************************
##  purge后
$assembly-stats purged.fa
stats for purged.fa
sum = 758769402, n = 672, ave = 1129121.13, largest = 20039767
N50 = 6804647, n = 35
N60 = 5417806, n = 47
N70 = 3860323, n = 63
N80 = 2086153, n = 91
N90 = 640976, n = 157
N100 = 669, n = 672
N_count = 0
Gaps = 0

最后

是否过度 purge呢？这个问题我打算最后再做，我打算先用未进行purge的基因组先进行染色体挂载，如果挂载率不高，很多contig如果挂不上，我在考虑再用pueged.fa去重新挂载，最后再做下两版基因组的评估工作，BUSCUO和KAT分析等

其他相关好文推荐：使用Purge_dups去冗余序列 - 简书 (jianshu.com)