FermiKit is a de novo assembly based variant calling pipeline for deep
Illumina resequencing data. It assembles reads into unitigs, maps them to the
reference genome and then calls variants from the alignment to an accuracy
comparable to conventional mapping based pipelines (see evaluation in the
directory). The assembly does not only encode SNPs and short INDELs, but also
retains long deletions, novel sequence insertions, translocations and copy
numbers. It is a heavily reduced representation of raw data. Storing,
distributing and analyzing assemblies is much faster and cheaper at an
acceptable loss of information.
FermiKit is not a prototype. It is a practical pipeline targeting large-scale data and has been used to process hundreds of human samples. On a modern server with 16 CPU cores, FermiKit can assemble 30-fold human reads in one day with about 85GB RAM at the peak. The subsequent mapping and variant calling only take half an hour.
Installation and Usage
The only library dependency of FermiKit is zlib. To compile on Linux or Mac:
git clone --recursive https://github.com/lh3/fermikit.git cd fermikit make
This creates a
fermikit/fermi.kit directory containing all the executables.
You can copy the
fermi.kit directory anywhere and invoke the pipeline by
specifying absolute or relative path:
# assembly reads into unitigs (-s specifies the genome size and -l the read length) fermi.kit/fermi2.pl unitig -s3g -t16 -l150 -p prefix reads.fq.gz > prefix.mak make -f prefix.mak # call small variants and structural variations fermi.kit/run-calling -t16 bwa-indexed-ref.fa prefix.mag.gz | sh
prefix.mag.gz for the final assembly and
for filtered SNPs and short INDELs and
prefix.sv.vcf.gz for long deletions,
novel sequence insertions and complex structural variations. If you have
multiple FASTQ files and want to trim adapters before assembly:
fermi.kit/fermi2.pl unitig -s3g -t16 -l150 -p prefix \ "fermi.kit/seqtk mergepe r1.fq r2.fq | fermi.kit/trimadap-mt -p4" > prefix.mak
It is also possible to call SNPs and short INDELs from multiple BAMs at the same time and produce a multi-sample VCF:
fermi.kit/htsbox pileup -cuf ref.fa pre1.srt.bam pre2.srt.bam > out.raw.vcf fermi.kit/k8 fermi.kit/hapdip.js vcfsum -f out.raw.vcf > out.flt.vcf
FermiKit does not use paired-end information during assembly, which potentially leads to loss of power. In evaluations, the loss is minor for germline samples and even without pair information, FermiKit is more sensitive to short INDELs and long deletions. Furthermore, with longer upcoming Illumina reads, it is actually preferred to merge overlapping ends in a pair before assembly and treat the merged reads as regular single-end reads (see AllPaths-LG and DISCOVAR).
Another technical limitation of FermiKit is that the error correction phase may take excessive RAM when the error rate is unusually high. In practice, this concern is also minor. I have assembled ~270 human samples and none of them require more than ~90GB RAM.
Running FermiKit twice on the same dataset under the same setting is likely to result in two slightly different assemblies. Please see bfc/count.c for the cause in BFC. Unitig construction also has a random factor under the multi-threading mode. Nonetheless, FermiKit should call the same variants from the same assembly.