44
$\begingroup$

I'd like to learn which format is most commonly used for storing the full human genome sequence (4 letters without a quality score) and why.

I assume that storing it in plain-text format would be very inefficient. I expect a binary format would be more appropriate (e.g. 2 bits per nucleotide).

Which format is most common in terms of space efficiency?

$\endgroup$
  • $\begingroup$ see: biostars.org/p/75178 "Why Don't We Use Binary Format?" $\endgroup$ – Pierre May 16 '17 at 18:11
  • 4
    $\begingroup$ Also important question: is the goal to create the smallest footprint on disk for an isolated single genome, or a multiple genomes? $\endgroup$ – woemler May 16 '17 at 18:13
  • $\begingroup$ @GWW Even if you had 5 letters (i.e. with N) you could get away with 3 bits per nucleotide and still have room for 3 more nucleotide encodings, maybe for like U, mC, hmC. $\endgroup$ – Jon Deaton Jul 27 '17 at 16:55
37
$\begingroup$

Genomes are commonly stored as either fasta files (.fa) or twoBit (.2bit) files. Fasta files store the entire sequence as text and are thus not particularly compressed.

twoBit files store each nucleotide in two bits and contain additional metadata that indicates where there's regions containing N (unknown) bases.

For more information, see the documentation on the twoBit format at the UCSC genome browser.

You can convert between twoBit and fasta format using the faToTwoBit and twoBitToFa utilities.

For the human genome, you can download it in either fasta or twoBit format here: http://hgdownload.cse.ucsc.edu/goldenPath/hg38/bigZips/

$\endgroup$
22
$\begingroup$

The standard formats for storing sequence data are fasta and fastq. Fasta is used if you only need the raw sequence data, fastq is used if you want to store the sequence data along with the quality information from base calling. Each of these can be compressed using gzip or another standard compression algorithm.

Typically we want to keep the quality information along with the raw sequence data, but the quality information accounts for half the storage space required. Some people have developed algorithms for lossy compression of the quality data that allow us to reduce the storage requirements.

If you are interested in storing variant calling data, the standard format for that is VCF. VCF is useful if you want to store quality information of the variant calls, genomic positions, and any annotations you might have about the position. VCFs can be compressed and indexed using bgzip and tabix. Many tools require variant data to be compressed and indexed using these tools.

$\endgroup$
  • 4
    $\begingroup$ fastq is generally used for storing read information garnered from sequencing. The information in a fastq is certainly not in chromosome or positional order and can be duplicated with multiple reads representing a position. $\endgroup$ – morgantaschuk May 16 '17 at 19:23
  • $\begingroup$ I can confirm with FASTA. Even Matlab's Bioinformatics Toolbox has a function (fastaread) for the import of such data. Plus, all of the genome data (at least that I have used) on NCBI is available in fasta. $\endgroup$ – Ashton Hellwig May 17 '17 at 3:31
  • 1
    $\begingroup$ I read the question as being about storing the reference sequence, not sequencing data. $\endgroup$ – Ian Sudbery May 18 '17 at 9:01
  • $\begingroup$ Yes, this answer is about raw sequence data from sequencing projects which is a very inefficient format for storing stable, large sequences. FASTA is indeed very popular, but not FASTQ, not for things like genomes. $\endgroup$ – terdon May 18 '17 at 10:57
18
$\begingroup$

The standard and the most common sequence format is FASTA for sure. You can compress it with a compressor. For the ~3GB human genome, gzip reduces the size to ~900MB, depending on the option in use.

Another often used format is UCSC's 2-bit format. This format keeps each A/C/G/T with 2 bits. As I remember, it keeps non-A/C/G/T bases and lowercases in two separate lists. These lists basically tell you that bases between offset x and y are all "N"/lowercase. The 2-bit format loses IUB codes that GRCh37 has. UCSC's hg19 differs from GRCh37 at a few bases.

BWA also produces its own 2-bit format with indexing. You can generate it separately with:

bwa fa2pac -f hg19.fa

Unlike UCSC, BWA keeps all IUB codes but loses letter cases. BWA does not provide utilities to convert its 2-bit representation to FASTA, either.

The 2-bit format typically reduces the file size down to 1/4 of its original size, unless there are too many scattered ambiguous bases. For human genome, you get a file ~784MB in size. You can compress it further with gzip, but that actually doesn't work well. A gzip'd 2-bit file is only ~5-10% smaller.

If you want to achieve an even smaller file size, you can compress the BWT of 2-bit file. This gives you a ~633MB file:

bwa pac2bwtgen hg19.fa.pac tmp.bwt && gzip tmp.bwt

A bit-aware compression algorithm may achieve an even higher compression ratio. However, such BWT-based compression prevents you from extracting subsequences. It is probably of little use in practice.

$\endgroup$
  • $\begingroup$ BWA replaces ambiguous bases by random nucleotides. See the original BWA paper: Section 2.7.1 ... "Non-A/C/G/T bases on the reference genome are converted to random nucleotides. Doing so may lead to false hits to regions full of ambiguous bases. Fortunately, the chance that this may happen is very small given relatively long reads. We tried 2 million 32 bp reads and did not see any reads mapped to poly-N regions by chance." $\endgroup$ – Karel Brinda May 30 '17 at 15:47
  • 3
    $\begingroup$ @Karel All the ambiguous bases are kept in the .amb file. BWA can reconstruct every base, well, at least in principle. $\endgroup$ – user172818 May 30 '17 at 16:21
  • $\begingroup$ Thank you very much for this information!! I have never used .amb files and they seem to be very useful for me. I wish I were aware of them earlier. Btw. I think we have some code to reconstruct the original sequences from BWA .bwt files. During our work on ProPhyle, we played a little bit with this type of compression. Maybe we will create a separate program for bwt2fa. $\endgroup$ – Karel Brinda May 30 '17 at 16:50
12
$\begingroup$

There are several things to consider when asking for "the most efficient" way to store data, it all depends on your use case. Do you just need ACGT, or are there also IUPAC codings for combinations? Do you need additional data (like quality values)? What kind of application are you using the data for (does it need to load all at once or in in chunks? Once or multiple times? Sequential or random access? etc.pp)?

E.g., most efficient for:

  1. Lowest footprint on disk, without a lot of hassle: use either FASTA or 2bit, but run through standard compressor (gzip, bzip2, others). The literature you want to consult here is that of standard text compression I think. Also of interest Large Text Compression Benchmark
  2. Keeping the file on disk, but ultra-fast loading small subsets into memory, being able to work in memory with character sized entities: a simple dump of the DNA as characters to disk, maybe combined with an index file to know which chromosome starts where. Then use mmap
  3. Storing quality values: See papers like Compression of FASTQ and SAM Format Sequencing Data or Sequence squeeze: an open contest for sequence compression
  4. Any combination of the above use cases + a lot more
$\endgroup$
  • $\begingroup$ I would classify 2bit format as 1b. Lowest footprint on disk, permitting some hassle. In order to use it for anything other than storage, it will need to be converted back into plaintext (fasta), or compressed plaintext (e.g. fasta.gz). $\endgroup$ – morgantaschuk May 16 '17 at 20:03
  • 1
    $\begingroup$ Your case 2 doesn’t require storing in uncompressed character data format. In fact, A compressed index (for instance) can make access faster by avoiding cache thrashing. $\endgroup$ – Konrad Rudolph May 17 '17 at 9:51
7
$\begingroup$

It is not yet standardized, but graph format has the potential for being the most space-efficient method for storing genomes. The idea is this: rather than store a genome as a linear string of sequenced nucleotides, genomes are stored as overlapping graphs, where sequence variants branch off from the reference genome, and then rejoin when the alignment continues. Basically, you start with a reference genome, and for every subsequent genome added to the graph, only the differences are stored. This could allow for an enormous gain in space efficiency.

$\endgroup$
7
$\begingroup$

I think the question is a bit ambiguous so please excuse this answer that's a bit redundant from the rest of the ones provided.

As others have mentioned, if you want to store a full genome, FASTA and 2bit formats are appropriate. For some context, hg19 is about 900Mb compressed for the FASTA file and about 780Mb compressed for the 2bit file . hg19 is a reference and is haploid so doesn't represent a "full" human genome that would normally have two alleles for the autosome (non-sex chromosomes).

A common format for representing variant information is Variant Call Format (VCF). The VCF format represent differences from a reference (hg19, say) that can be used to recover the original full sequence by using the reference and the differences encoded in the VCF file. I've seen VCF files in the range of 100Mb, but a reference file is still needed to recover the full genome sequence which is the range of 800Mb+, as mentioned above.

If you're considering just one "whole genome" in isolation, then the answer is pretty clear: 2bit format is probably approaching the entropy limit of the human genome and you probably won't be able to do much better. The reason why your question is a bit ambiguous is that as soon as you start encoding more than one genome, a population of genomes, say, then you can start exploiting the redundancy of the genome as shared by the population.

For example, say you want to store two "whole genomes". You could download the hg19 reference and download two VCF files which would give around 1Gb worth of data (around 800Mb for the 2bit file and around 200Mb for both of the VCF files). Now you've been able to represent a "whole genome" in 500Mb instead of the 800Mb. You can see a similar argument for downloading 3 VCF files and more.

The minimum amount of information needed to represent a population of genomes is, as far as I know, unknown, but I would guess in the 2.5Mb-5Mb range. For example, see "Human genomes as email attachments" by Christley, Lu, Li and Xie which claims a 4Mb encoding of a genome.

Things get tricky because you have to ask what you're claiming as a "whole genome". VCF files are notoriously bad because older versions of the specification only store high quality differences from reference, throwing away high quality called sections. If you want to store low quality information, the encoding is now going to depend on the sequencing technology in weird ways.

Insertions, deletions, mobile insertion elements, copy number variants, other structural variants, etc. all complicate this matter further. Genome Graphs are trying to tackle at least some of these problems but the focus is on variant calling rather than efficient individual whole genome representation, though perhaps can be adapted in the future.

$\endgroup$
  • $\begingroup$ How can there be 2 orders of magnitude differences between the 2bit compression (hundreds of Mb) and the "e-mail attachment" (a few Mb)? Is the "e-mail attachment" case a really impractical storage format so that no one actually uses it? The abstract doesn't say it, but it seems that what they store is actually variations than the full standalone info. $\endgroup$ – bli Jun 8 '17 at 20:50
  • 1
    $\begingroup$ @bli, the order of magnitudes difference comes from exploiting the redundancy in a population. Storing one persons (encoded) DNA requires ~800Mb. Storing thousands of peoples takes (probably) a couple megabytes (each). When dealing with a population of DNA data that you want to encode, there are many ways of doing it. One way is to store variants from a reference. Another is to store a library of short "reads" and then reference that library. The paper is meant as a proof of concept to answer the underlying question of "what is the theoretical minimum space required to store a whole genome". $\endgroup$ – abetusk Jun 13 '17 at 0:59
4
$\begingroup$

In terms of raw storage capacity 2 bits per nucleotide, and then further compressed with standard compression techniques would be the most efficient. However, you'd still have other storage considerations. Like what to do about non-standard bases: like if you want to indicate a gap or ambiguity.

I'd also query if it is really necessary to store them as binary since it reduces the readability of the data. It is quite convenient having a whole bunch of unix and programming tools that can operate on the string level in text files.

$\endgroup$
  • $\begingroup$ In terms of storage space '2 bits per nucleotide' and compressed would be a way to save resources. $\endgroup$ – Severian May 16 '17 at 18:16
3
$\begingroup$

In all seriousness, the most efficient way to store DNA sequence data is...you guessed it...in DNA. (Church, Gao, and Kasuri, 2012) and others have used DNA synthesis and sequencing as an information writing/reading mechanism.

Practical? Not yet.

Storage efficiency? Unparalleled!

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.