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 ...
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 ...
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 ...
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 ...
An epic question. Unfortunately, the short answer is: no, there are no widely used solutions.
For several thousand samples, BCF2, the binary representation of VCF, should work well. I don't see the need of new tools at this scale. For a larger sample size, ExAC people are using spark-based hail. It keeps all per-sample annotations (like GL, GQ and DP) in ...
It’s a matter of preference I guess but I recommend the Ensembl builds. Decide whether you want the toplevel or primary assembly, and whether you want soft-masked, repeat-masked or unmasked files. The naming schema is very straightforward; the combinations are described in the README file, and all files reside in one directory.
For example, if you want the ...
You can use the Jellyfish software to calculate the k-mer profiles up to length 31.
From the instructions in the user guide:
The basic command to count all k-mers is as follows:
jellyfish count -m 21 -s 100M -t 10 -C reads.fasta
To compute the histogram of the k-mer occurrences, use the histo subcommand (see section 3.1):
jellyfish histo mer_counts.jf
You can also use R. I give you an example of only chr1 and only kmer=4.
genome <- BSgenome.Hsapiens.UCSC.hg38
kmers <- oligonucleotideFrequency(genome$chr1, 4)
m <- as.matrix(kmers)
tl;dr: Just use the either the downloads on the Bowtie2 homepage or the Illumina iGenomes. Or just uncompress and concatenate the FASTA files found on UCSC goldenpath and then build the index.
A bit longer answer:
There are two components to "genome for a read mapper" such as Bowtie or BWA.
First, you need to choose the actual sequence (genome release ...
In my opinion, it is not very reliable. LiftOver is very limited in terms of transformations it can support. The LiftOver Chain format can capture only matching regions in the same order. It means that it can account for indels, but even simple structural variations become problematic.
For instance, when a newer assembly is available, it is usually a ...
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 ...
TLDR; For the most modern CHM13 assembly (v1.1; Complete T2T reconstruction of a human genome; https://github.com/marbl/CHM13), the best compression without any use of a reference would be between 624 MB and 641 MB.
Upper estimates of genome compressibility can be made using standard compressors (see, e.g., this comparison). In my experience, xz can achieve ...
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 ...
1. Adapter Trimming
One of the first things I do after encountering a set of reads is to remove the adapter sequences from the start and end of reads. Most basecalling software includes some amount of built-in adapter trimming, but it is almost always the case that some adapter sequence will remain. Removing adapters is helpful for mapping because it ...
At the moment, the standard reference genomes (e.g. hg19, hg38) are haploid genomes. We know that the human genome is diploid. Naturally, the latter would be the respectively correct representation of the human genome.
The premise of the OP's question is false. The natural reference representation of the human genome is not diploid.
Think of a reference ...
Here's an example for the mouse genome:
select(org.Mm.eg.db, c("GO:0048406"), c("GENENAME","SYMBOL"), c("GO"))
You get output like:
1 pregnancy zone protein Pzp
2 nerve growth factor receptor (TNFR superfamily, member 16) Ngfr
For simple variants like SNPs it would not really be a problem to use the current genome assembly for other ethnic groups. But for more complex variants this could be indeed problematic, however not only for ethnic groups but also for individuals within the same population. Think of very complex regions such as HLA or KIR.
In studies where they compare ...
You will probably be interested in the following UCSC wiki page, which explains how to go from most of the UCSC tables to GTF/GFF:
The basic gist is that UCSC doesn't store any data internally as GTF or GFF, and so you will need to use our genePredToGtf utility in order to convert from our ...
The FAQ offers an answer:
I'm setting up my own Blat server and would like to use the same parameter values that the UCSC web-based Blat server uses.
We almost always expect there to be some small differences between the hgBlat/gfServer and the stand-alone command-line blat. The best matches can be found using pslReps and pslCDnaFilter utilities. ...
For calling small variants, the standard way is to simply call diploid genotypes. You can already do a variety of research with unphased genotypes. You may further phase genotypes with imputation, pedigree or with long reads/linked reads, but not many are doing this because phasing is more difficult, may add cost and may not always give you new insight into ...
This can be done quite easily using Ensebl's BioMart.
Choose the Ensembl Regulation database:
Select the "Human Regulatory Features" dataset:
That's basically it right there, just click on "Results":
Export to file and click "Go":
This will download a file called mart_export.txt which looks like this (I chose TSV for tab separated values):
$ head ...
There's a benchmark that compared dozens of compressors on DNA data, including human genome: http://kirr.dyndns.org/sequence-compression-benchmark/
In particular, this table shows file size of compressed human genome, produced by the strongest setting of each compressor: http://kirr.dyndns.org/sequence-compression-benchmark/?d=Homo+sapiens+GCA_000001405.28+%...
I think that right now, the only human builds that are worth considering are hg19/GRCh37 as many data bases such as gnomAD still exclusively use this release. On the other hand, hg38/GRCh8 has many important fixes and the useful (but yet underused) feature of alternative loci.
Anything from older releases should be remapped to a more recent one.
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 ...
Grabbing chromosomes for hg38:
$ wget ftp://hgdownload.cse.ucsc.edu/goldenPath/hg38/chromosomes/*.fa.gz
$ for fn in `ls *.fa.gz`; do gunzip $fn; done
Via kmer-counter and Python, here's how to search for kmers of length 7 from chromosome chrY:
#!/usr/bin/env python ...
Here is one definition of the classical MHC region: https://www.ncbi.nlm.nih.gov/grc/human/regions/MHC?asm=GRCh38.p13
which defines it as the 5mb region chr6:28510120-33480577 in GRCh38 coordinates.
The "extended MHC" defined here  is 7.2 mb and encompasses an additional "extended Class I" region upstream of the region given above. The ...
You want FASTA files from here.
If you want the entirety of the genome, get the dna-toplevel file. Files with rm and sm have repeat sequences masked, the former by replacing the repeats with Ns (hard mask) the latter by designating repeats as lower case letters (soft masked).
See the README for more info.
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!