Preprinted in 2013 , bwa-mem is a popular mapper for 100+ bp short reads. It is often the first choice in variant calling pipelines for human data. However, bwa-mem has increasingly become the performance and the cost bottleneck when the upstream/downstream get optimized. This led to multiple forks, including bwa-mem2 , bwa-meme , Sentieon bwa and Parabricks bwa , which produce identical or near…
There have been multiple ongoing efforts (see Appendix below), and debates , regarding porting bioinformatics tools to Rust. These include two rewrites of minimap2. Here are some of my current reflections, which may change in the future. From a purely technical standpoint, I believe AI rewrites, especially in Rust, benefit the field. Rust will not only speed up bioinformatics tools in Python and R…
Those who inspect SVs often have probably noticed many SVs fall in low-complexity regions (LCRs) or tandem repeats and they are challenging to call. How many and how challenging? I could not find a good answer in the literature, so I decided to work with Qian (Alvin) Qin to measure by ourselves. To start with, we need to identify LCRs in the human genome. This turned out to be a non-trivial…
Sequencing by Expansion (SBX) is Roche’s new short-read sequencing technology. Roche released the first batch of public SBX data after a webinar yesterday. This blog post provides my impression with the data. If you are interested in SBX, I highly recommend watching the recording of the webinar, as it contains a lot more information than this blog post. Roche released duplex SBX (SBX-D) data…
TL;DR: the new preset splice:sr in minimap2 can align short RNA-seq reads. It is similar to STAR in resource usage, approaches STAR in junction accuracy, and is overall better at SNP calling. Minimap2 can align short or long genomic reads and long RNA-seq reads. Short RNA-seq read alignment is a notable omission. This is changed with the latest v2.29 release — minimap2 now has a new splice:sr…
BWA implements three algorithms: bwa-aln , bwa-sw and bwa-mem . Although bwa-aln is the oldest, it remains a popular choice for ancient DNA (aDNA) data. Is it because researchers working on aDNA are too lazy to change? Not quite. The following plot from Dolenz et al (2024) shows the mismatch rate as a function of read length for reads mapped with mapping quality (mapQ) 25 or higher. This plot…
My blog posts are often motivated by casual conversations with other scientists. This blog post is an example. Someone asked me how I thought of using a compression algorithm for alignment in BWA . The short answer is: I didn’t; I learned from Tak-Wah Lam . Here is the longer story. On a 1000G conference call in April 2008, as I remember, Tak-Wah gave a short talk on a prototype of SOAP2 and…
Historically, the term “pan-genome” (or more often “pangenome” in recent literature) was first coined by Sigaux in 2000 for a purpose distinct from current uses. Victor Tetz’s definition is fairly close to the modern definition except that it is applied to all living organisms. These definitions are rarely cited nowadays. In microbiome research, a pangenome refers to the non-redundant set of genes…
In the past couple of months, I have been asked several times about what language(s) to learn if someone wants to write high-performance programs. This is a sensitive topic that often triggers heated debate partly because many fast languages share similar features and are comparable in performance. My general take is that for small research projects only involving a few developers, the choice of…
A $k$-long sequence $P$ is a ($k$,$s$)-open- syncmer , $s\le k$, if $P[1,s]$ is the smallest among all $s$-mers in $P$. Suppose function $\phi$ is a bijective hash function of $k$-long sequences. $P$ is a random ($k$,$s$)-syncmer if $\phi(P)$ is an open syncmer. Because we often map $k$-mers to integers, $\phi$ can take the form of an invertible integer hash function . In practice, $\phi$ does not…
Command-line interface , or CLI in brief, specifies how a user interacts with a program on the command line. Torsten Seemann wrote a good article on creating CLI. This blog post adds a few more suggestions. 1. Keep the backward compatibility of CLI as much as possible Backward compatibility here means users can upgrade and run a tool without changing the command lines they used in the past. This…
Definition. The dual assembly of a diploid sample consists of two sets of contigs with each set representing one complete haploid genome. Similar to contigs in a primary assembly, contigs in a dual assembly may have occasional switches between parental haplotypes. I called such an assembly as partially phased assembly in an earlier post but decided to coin a new term in our new hifiasm preprint…
This is a short post on how to remap short reads in an aligned BAM using bwa-mem. My recommendation is (requiring bash) samtools collate -Oun128 in.bam | samtools fastq -OT RG,BC - \ | bwa mem -pt8 -CH <(samtools view -H in.bam|grep ^@RG) ref.fa - \ | samtools sort -@4 -m4g -o out.bam - Here, samtools collate groups the two reads in a read pair and outputs an uncompressed BAM stream. samtools…
This post is inspired by Vince’s tweet . It describes my thoughts on the design of command-line interface (CLI). Note that this article doesn’t necessarily represent the best practices; it just shows my personal preferences. First of all, I need to clarify the terminology. For an example, in the command line below rm -f file.txt rm is the the command . file.txt is the command-line argument , which…
Leonardo Martins tweeted that xz can compress a 1.4 million SARS-CoV-2 genomes in a 39GB FASTA down to 74MB. That is a very impressive compression ratio! This reminds me of my earlier work on FM-index construction . For an experiment, I downloaded ~400k SARS-CoV-2 genomes from EBI’s COVID-19 data portal ( GISAID has ~1.5M genomes but imposes restrictions) and generated an FM-index of these…
Concepts in phased assembly: Contig : a contiguous sequence in an assembly. A contig does not contain long stretches of unknown sequences (aka assembly gaps ). Scaffold: a sequence consists of one or multiple contigs connected by assembly gaps of typically inexact sizes. A scaffold is also called a supercontig , though this terminology is rarely used nowadays. Assembly: a set of contigs or…
Ian Holmes has a twitter poll right now on the use of “SNP” (single-nucleotide polymorphism) versus “SNV” (single-nucleotide variant). I have been bugged by the two terminologies for years, so I decided to write a blog post on it. Personally, I use “SNP” for germline events and “SNV” for somatic events , but I understand others think differently. Here are my thoughts. The wiki page for SNP defines…
Honestly, I didn’t know what minigraph would be good for when I was writing the code. When I was writing the paper , I pitched minigraph as a fast caller for structural variations (SVs). However, except performance and convenience, minigraph is not that special. In fact, in the paper, minigraph is not as good as read-based SV callers because it randomly misses one parental allele when most…
Why? It is usually easy to evaluate the contiguity of a de novo assembly – just compute N50. It is much harder to evaluate the correctness. We typically identify misassemblies by aligning contigs to a reference genome. However, it is tricky to interpret the results. In case of human, there are thousands of structural variations (SVs) between the reference and the sample being assembled.…
In an earlier post a few days ago, I said “discarding base quality dramatically reduces variant calling accuracy”. I didn’t provide evidence. This certainly doesn’t sound persuasive. In this post, I will show an experiement to support my claim. I downloaded high-coverage short reads for sample HG002 from GIAB ftp , converted to unsorted FASTQ with samtools collate , mapped them to hs37d5 (for…
This short post evaluates the effect of format and quality binning on file sizes. I am taking SRR2052362 as an example. It gives 4.3-fold coverage on the human genome. For 2-binning, I turned original quality 20 or above to 30 and turned original quality below 20 to 10. For 8-binning, I took the scheme from a white paper (PDF) published by Illumina. Illumina has been using quality binning for more…
Background Python and R are slow when they can’t rely on functionality or libraries backed by C/C++. They are inefficient not only for certain algorithm development but also for common tasks such as FASTQ parsing. Using these languages limits the reach of biologists. Sometimes you may have a brilliant idea but can’t deliver a fast implementation only because of the language in use. This can be…
Given a de novo assembly, we often measure the “average” contig length by N50. N50 is neither the real average nor median. It is the length of the contig such that this and longer contigs cover at least 50% of the assembly. A longer N50 indicates better contiguity. We can similarly define N x such that contigs no shorter than N x covers x % of the assembly. The N x curve plots N x as a function of…
I wrote a blog post on a potential reference pan-genome model. I had more thoughts in my mind. I didn’t write about them because they are immature . Nonetheless, a few readers raised questions related to my immature thoughts , so I decide to add this “Part II” as a response. Please note that this and the previous blog posts only represent my own limited view . A consortium will be formed to build…
In the last weekend, I made gfatools and minigraph open to the public. Both repos come with some documentations, but they haven’t explained the background and motivation behind. This blog post gives a more complete picture. The primary assembly of GRCh38, our current human reference genome, is largely the concatenation of individual haplotype segments. It aims to model a single human genome and…
I often hear developers saying “I do XYZ because it saves my time”. “XYZ” could be the selection of programming language, the use of 3rd-party libraries or other choices in programming. When I hear this, my immediate reaction is always: where to put users’ time into the equation? Here is how I think about it. In my view, the value of a feature is roughly measured by value(feature) =…
In 2006, Ruiqiang Li, the now CEO of Novogene, said to me in Chinese: “you can’t maintain the TreeFam database forever”. Considering TreeFam as my most significant work at the time, I said “I think I can”. Some of you know what happened next: one year later, I started to develop tools for next-gen sequencing data and forgot about TreeFam almost completely. Few can constantly and single-handedly…
I wrote a commentary on the SAM/BAM format a while ago. I now publish it as a blog post, 10 years after I released the first samtools to SourceForge. It gives an overview about what has happened in the past 10 years. Here “we” refers to the samtools dev team and the HTS file format commitee, as well as those who have contributed to the specification or samtools/htslib. I am only one of them. When…
Sequence identity is a way to measure the similarity between two sequences. For sequencing data, it is often thought as the opposite of sequencing error rate. When we say “the sequence divergence between two species is ABC” or “the sequencing error rate is XYZ”, we assume everyone knows how to compute identity. In fact, there are more than one ways to compute identity. This blog post discusses a…
On Twitter, Zach charlop-powers asked me to give a code walkthrough for seqtk . This post does that, to a certain extent. Seqtk is a fairly simple project. It uses two single-header libraries for hash table and FASTQ parsing, respectively. Its single .c file consists of mostly independent components, one for each seqtk command. I will start with the two single-header libraries. Buffered stream…
SAM is a text format that is typically used to store the alignment of high-throughput sequence reads against a reference genome. BAM is the first binary representation of SAM designed at the same time. BAM is smaller, faster to process and has additional features like random access. BAM is not optimal in terms of compression ratio. By reorganizing binary data and using more advanced compression…
My minimap2 paper has been accepted for publication in Bioinformatics. You can find the latest LaTeX source at OverLeaf or in the tex directory of minimap2 . I am intentionally delaying the publication process for personal reasons. It will take a while for you to see the published version at Bioinformatics. I thought to write this blog post when the paper comes out, but there have been a few…
In the SAM format, the “X” and “=” CIGAR operators were not part of the original spec. Nonetheless, they were among the first several features added after the initial release of the spec. I was resistant to this feature for several reasons. First, CIGAR describes alignment, but sequence matches and mismatches are not indispensable properties of alignment. Second, for this reason, most older…
By mathematical definition , a k -order (or k -dimensional) De Bruijn graph, or ${\rm DBG}(k)$ in brief, over the DNA alphabet uses k -mers at vertices. It has $4^k$ vertices and $4^{k+1}$ edges. DBG(k) has two interesting properties. First, DBG(k) is the line graph of DBG(k-1). Intuitively, this means an edge in DBG(k-1) uniquely corresponds to a vertex in DBG(k) and that the edge adjacency of…
TL;DR: If you map reads to GRCh37 or hg19, use hs37-1kg : ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/technical/reference/human_g1k_v37.fasta.gz If you map to GRCh37 and believe decoy sequences help with better variant calling, use hs37d5 : ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz If you map reads to GRCh38 or hg38, use the…
Introduction Illumina Inc. released NovaSeq earlier this year and provided sample data at BaseSpace several months later. Different from the HiSeq series, NovaSeq uses 2-color chemistry . It has been observed that the NextSeq series, which also uses 2-color chemistry, produced data of worse quality . One naturally wonders if NovaSeq has a similar problem. This post might give you some hints. Data…
Getting Started Firstly, a few basic concepts. Conda is a portable package manager primarily for Python and precompiled binaries. Miniconda is the base system of conda. It includes a standard python and a few required dependencies such as readline and sqlite. In conda, a channel contains a set of software typically managed by the same group. Bioconda is a channel of conda focusing on…
I have just reimplemented the symmetric DUST algorithm (SDUST) for masking low-complexity regions. The program depends on kdq.h (double-ended queue) and kvec.h (simple vector); the command line interface further requires kseq.h for FASTA/Q parsing. As I have tried on human chr11, the output is identical to the output by NCBI’s dustmasker except at assembly gaps. The speed is four times as fast. I…
GraphMap is a new long-read mapper initially tuned for error-prone ONT reads. There are quite a few interesting points methodologically. The following two comments are mostly about techical and practical aspects. Before you read the comments, please bear in mind that I am the developer of BWA-MEM. I could be biased. Comment 1: consensus quality The most striking point I found in the preprint (and…
Today, Google and Broad Institute (my employer) have announced that they are teaming up to tackle genomic data. One sentence caught my attention: “Broad Institute has … either sequenced or genotyped the equivalent of more than 1.4 million biological samples”. Can we get the data? Current data sharing model In my limited experience, the current data sharing model is largely trust-and-distribute.…
Installing docker on Mac With all the buzz around docker , I finally decided to give it try. I first asked Broad sysadmins if there are machines set up for testing docker applications. They declined my request for security concerns and suggested Kitematic for my MacBook. This means I can hardly run sequence analyses for human. Anyway, I followed their suggestion. Kitematic turns out to be easy to…
As is discussed in my previous post, a major but potentially fixable problem with the VCF model is that we allow to and sometimes have to squeeze multiple alleles in one VCF line. This post gives the solution, the unary representation. The representation was first conceived by Richard Durbin a couple of years ago, but in my view had a few practical issues initially and thus never openly presented.…
This and the next posts were mostly written on the plane when I felt tired and did not have Internet connections. The logical flow is not very clear. In addition, I have to admit that I have not thought through the topic when I was writting it up. Now I have a clearer picture after I finish the posts. I will still put them online for a historical record. Edit-based representation VCF represents a…
I enjoy writing 2-page Application Notes these days. It takes less time to write, giving me more time on solving other problems. More importantly, I don’t need to fight to claim significance and novelty which are subjective most of time. The downside of writing short manuscripts is the lack of extensive discussions. Here are something I have not said in my new error correction preprint . The role…
While I was looking for an ancient email on my old (first) macbook, I noticed the numerous email exchanges during the early days of the SAM/BAM format. Here is a brief summary. The ideas below were proposed by various people in the 1000 Genomes Project analysis group. 2008-10-21: SAM got its name. 2008-10-22: The first day: fixed columns and optional tags; extended CIGAR and binning index.…
A recent paper published by Phil Ashton et al has triggered some discussions which subsequently moved to my domain: read mapping. Phil then asked me to clarify how the upcoming bwa-mem works with Oxford Nanopore (ONT) reads. Here we go. Although the very first version of bwa-mem worked with PacBio reads (well, not crashing), the alignment it produced was too fragmented to be useful. I initially…
Illumina has recently released four lanes of NA12878 data from HiSeq X10. I was playing with this data set and found that my program had bad accuracy on two of them. I initially thought the data quality was different, so wrote some code to investigate the data quality. It turns out that my program was buggy, but the finding of the HiSeq X10 data quality might be of its own interest, which I am…
Introduction Ever since the advent of the so-called Next-Generation Sequencing (NGS), we have been thinking about encoding all the population variations in a graph. That was 2008. Now, six years later, the rapidly growing number of sequenced human individuals continueously presses for the necessity of a graphical representation of the existing sequences, which leads to many publications in this…
I was out of the town in the past few days, so have not been able to focus on GFA. Now I am back to work to give the first update on the format based on the comments from many people, which I appreciate a lot. In comparison to my initial proposal, the first and the major change is to name segments instead of the ends of segments. This seems the consensus so far. Secondly, I am thinking to move the…
PSMC is my program to infer the historical effective population size from a diploid genome. It was published in Nature three years ago and has been cited over 100 times so far. Whenever I see a PSMC plot in a paper, I feel a moment of joy both as a scientist and as a programmer. PSMC is okay, but now there are better models and implementations at least in theory. MSMC , which has recently been…