「三代组装」使用Pilon对基因组进行polish
软件安装
官方提供了编译好的jar包,方便使用
wget https://github.com/broadinstitute/pilon/releases/download/v1.23/pilon-1.23.jar
java -Xmx16G -jar pilon-1.23.jar
如果要顺利运行程序,要求JAVA > 1.7, 以及根据基因组大小而定的内存,一般而言是1M大小的基因对应1GB的内存。
总览
Pilon有如下作用
- 对初步组装进行polish
- 寻找同一物种不同株系间的变异,包括结构变异检测
他以FASTA和BAM文件作为输入,根据比对结果对输入的参考基因组进行提高,包括
- 单碱基差异
- 小的插入缺失(indels)
- 较大的插入缺失或者block替换时间
- 填充参考序列中的N
- 找到局部的错误组装
最后它输出polish后的FASTA文件, 以及包含变异信息的VCF文件(可选)
分析流程
推荐使用PCR-free建库测序得到的Illumina paired-end数据,这样子避免了PCR-duplication,有效数据更多,也不需要在分析过程中标记重复。
下面步骤,假设你的组装文件为draft.fa
, 质量控制后的illumina双端测序数据分别为read_1.fq.gz
和read_2.fq.gz
第一步:比对
bwa index -p index/draft draft.fa
bwa mem -t 20 index/draft read_1.fq.gz read_2.fq.gz | samtools sort [email protected] 10 -O bam -o align.bam
samtools index [email protected] 10 align.bam
第二步:标记重复(非PCR-free建库)
sambamba markdup -t 10 align.bam align_markdup.bam
第三步:过滤高质量比对的read
samtools view [email protected] 10 -q 30 align_markdup.bam > align_filter.bam
samtools index [email protected] 10 align_filter.bam
第三步:使用Pilon进行polish
MEMORY= #根据基因组大小而定
java -Xmx${MEMORY}G -jar pilon-1.23.jar --genome draft.fa --frags align_filer.bam \
--fix snps,indels \
--output pilon_polished --vcf &> pilon.log
一般Pilon迭代个2到3次就够了,所谓事不过三,过犹不及。
关于Pilon的一些参数说明:
-
--frags
表示输入的是1kb以内的paired-end文库,--jumps
表示 大于1k以上的mate pair文库,--bam
则是让软件自己猜测 -
-vcf
: 输出一个vcf文件,包含每个碱基的信息 -
--fix
: Pilon将会处理的内容,基本上选snps
和indels
就够了 -
--variant
: 启发式的变异检测,等价于--vcf --fix all,breaks
, 如果是polish不要使用该选项 -
minmq
: 用于Pilon堆叠的read最低比对质量,默认是0。
阅读日志输出
这个日志文件是标准输出而不是标准错误输出,不过保险起见用
&>
最开始,Pilon会输出他的版本信息(如下示例),以及将会对基因组做的调整,
Pilon version 1.14 Sat Oct 31 14:30:00 2015 -0400
Genome: genome.fasta
Fixing snps, indels
其中Fixing后面的含义为:
- "snps": 单碱基差异
- "indels":小的indel的差异
- "amb": 替换原有的N
- "gaps": 填充基因组的gap
- "local": 检测和修改错误组装
- "all": 上述所有
- "none": 不是上述的任何一种
接着Pilon会分染色体对BAM文件进行处理,根据BAM文件进行堆叠(pileup), 这个时候它会输出有效reads的深度,这里的有效reads包括未成对的read和正确成对的read。
Processing ctg1:1-5414473
frags align_mkdup.bam: coverage 19
Total Reads: 808985, Coverage: 19, minDepth: 5
从Pilon v1.14开始,Pilon还会输出基因组得到确认的碱基比例。
Confirmed 5403864 of 5414473 bases (99.80%)
后续是Pilon将会对原参考基因组做的一些调整的总体情况,如下表示纠正2个snp, 2个小的插入,4个缺失。
Corrected 2 snps; 0 ambiguous bases; corrected 2 small insertions totaling 12 bases, 4 small deletions totaling 6 bases
最后声明当前部分处理结束
Finished processing ctg1:1-5414473
如果,在--fix
中选了gaps
, 那么输出的内容还有如下内容。其中82048 -0 +276
解释为在坐标82428处移除0个碱基,插入276个碱基。
# Attempting to fill gaps
fix gap: scaffold00001:82428-93547 82428 -0 +276 ClosedGap
参考资料