前言
由于最近要开始做一些基因组的工作,然后老板安利了一个叫shovill的github仓库,主要作用就是加快SPAdes的拼接速度。这里记录一下它主要做的事情,以及benchmark一下加快的效果如何。
shovill简介
Shovill is a pipeline which uses SPAdes at its core, but alters the steps before and after the primary assembly step to get similar results in less time. Shovill also supports other assemblers like SKESA, Velvet and Megahit, so you can take advantage of the pre- and post-processing the Shovill provides with those too
以上是shovill自己的介绍,其中提到了其改变的步骤主要是primary assembly之前与之后的部分,并且也提到了其结果可能会发生一定的变化。
另外,除了以上的主要介绍以外,其中也有一些FAQ的回答,例如:
- 不能用于3代测序所产生的长序列测序结果
- 只适用于pair end 测序结果
- 不能用于metagenomic的结果
shovill主要过程
- 获取raw read统计信息(seqtk)
- 估计基因组的大小(kmc)
- 估计测序深度(1与2的结果的简单计算)
orig_depth = total_bp / genome_size
factor = depth / orig_depth
default depth: 100 - 通过seqtk同时对read1和read2进行subsample original data. **(其理由是Giving an assembler too much data is a bad thing) **
- 通过trimmomatic进行trimming(默认是不trim)
- 计算Kmer
default $KMER_READ_FRAC=0.75; maxK=127; minK=31
maxK = min(maxK, $KMER_READ_FRAC * length of read)
在127和3/4的read长的中选个最小的作为新的maxK
minK = 21 if avg_read_len < 75 else minK
如果平均read长小于75则minK等于21
number of K = 5
写死。。
step of K
通过min 和 max计算可得,保持全奇数的结果。 - 通过lighter对raw reads进行correct。 (校正测序数据,通过测序深度和预估基因组大小进行校正)
- 通过flash对Overlapping/stitcing PE reads进行校正
- 开始进行拼接。其中
--pe1-1 和--pe1-2 和--s2
均由flash产生。
spades.py --pe1-1 flash.notCombined_1.fastq.gz --pe1-2 flash.notCombined_2.fastq.gz --s2 flash.extendedFrags.fastq.gz --only-assembler --threads 16 --memory 100 -o . --tmp-dir /tmp -k 31,51,71,91,111 >> /dev/null 2>&1
- 通过Pilon进行拼接后的校正。首先会把r1和r2比对会拼接好的contig上。
_JAVA_OPTIONS=-Xmx100g pilon --genome contigs.fasta --frags shovill.bam --fix bases --output pilon --threads 16 --changes --mindepth 0.25 >> 80-pilon.log 2>&1
pilon的结果会被删除,只保留Count changes per contig
- 规整contigs.fa的header。
shovill的评估
对单菌株的全基因组测序数据(total bp:132M,预估genome size:4M,预估测序深度为300x)的测序数据总共耗时7min