Bioinformatics & Computational Biology: Chapter 2 — Databases and file formats

2.1 Primary archives and layers built on them #

2.1.1 Three nucleotide archives, one collection #

GenBank, the European Nucleotide Archive, and the DNA Data Bank of Japan form the International Nucleotide Sequence Database Collaboration. A submission to one is exchanged with the others, so the three archives are redundant on purpose. They check format. Acceptance does not mean the gene exists, the species is right, or the function was assayed.

A database record is not a peer-reviewed fact. Review, when it happens, belongs to a paper, which may cover only part of the record, and many records have no paper. Later curation removes some errors and does not replace an experiment. A browser can be more useful than the raw submission, and it can also copy one bad label widely.

2.1.2 RefSeq #

RefSeq is NCBI’s curated layer, built from the archive and not identical to it. The prefix is part of the evidence.

An `XP_` hit and an `NP_` hit are not the same support for a function. RefSeq can still be wrong, and it lags new assemblies. “In RefSeq” means “in this release”.

2.1.3 UniProt #

UniProt is the protein archive used when the question is function. Swiss-Prot entries are reviewed: a curator has compared the sequence with the literature and tied function to evidence. The text can still be incomplete, or later overturned.

TrEMBL entries are automatic translations. The product name is usually transferred, often from a name that was itself transferred. TrEMBL is larger, so the top hit is often unreviewed. The accession is stable. The mnemonic entry name can change.

2.1.4 The Protein Data Bank #

The Protein Data Bank archives structure depositions. A four-character code does not encode the biology, and one protein can have many entries. The file is a coordinate model plus the method, the construct, and, for diffraction and cryo-electron microscopy, a resolution.

Resolution is detail in the map, not a grade for biological truth. A 3 Å map can show the fold and still be the wrong object for a bond length. NMR entries are ensembles and have no crystallographic resolution. The chain may be a fragment. Cite the code, the chain, the method, and whether the coordinates were measured. Building the map is Structural Biology.

2.1.5 Gene Ontology #

The Gene Ontology is a controlled vocabulary, not a sequence database, with three namespaces: molecular function, biological process, and cellular component. A child term implies its parents. Mixing the namespaces produces a sentence that sounds annotated and cannot be tested. Every annotation should carry an evidence code.

Most annotations in a large database are electronic. An IEA or ISS term was copied, sometimes from a term that was itself copied, so the chain may contain no experiment. A term with no code cannot be read. A direct assay still has a scope: binding in a tube is not a map of where the protein acts in a cell.

2.2 Files a student must be able to read #

The study of Files a student must be able to read provides fundamental principles governing cellular architecture, molecular interactions, and physiological regulation.

🔬 Key Illustrated Scientific Figure

Full-Screen Zoom

Click to enlarge on this page

2.2.1 FASTA #

A record is a `>` header and a sequence. The first token is the identifier; the rest is free text and may be false. Programs often split on whitespace, so a space inside the identifier is lost downstream.

>toyA Hypothetical example, not a real accession
ATGACTGAGACCTCTGAACAGTAA

Case is usually ignored, and one file may hold many records. FASTA carries no quality and no coordinates. Hyphens in an alignment stored this way are a private convention, not part of the format.

2.2.2 FASTQ and the Phred score #

FASTQ adds a quality string of the same length as the sequence. The four lines are `@` header, sequence, `+`, and qualities.

@read1 toy
ATGC
+
IIII

Phred Q = −10 log₁₀(p), where p is the caller’s probability that the base is wrong. Q20 is 0.01, Q30 is 0.001, and Q40 is 0.0001. In Phred+33, used by Sanger and by modern short-read files, Q is the ASCII code minus 33: `I` (code 73) is Q40, `?` is Q30, and `5` is Q20. Older Illumina files used Phred+64. The character means nothing until the encoding is known.

Quality lines may themselves contain `@` and `+`, so a parser that starts a new record at every `@` will break the file. Q is only as honest as the caller. A miscalibrated Q30 is not one error in a thousand.

2.2.3 The GenBank flat file #

A `CDS` gives coordinates, often a join of exons, plus `/translation` and `/product`. Translate the interval yourself: disagreement with `/translation` means the record is inconsistent. The product string may rest on an assay, a weak similarity, or nothing.

LOCUS       TOY_EXAMPLE    24 bp    DNA     linear
DEFINITION  Toy record. Not a real submission.
ACCESSION   TOY000001
VERSION     TOY000001.2
FEATURES             Location/Qualifiers
     CDS             1..24
                     /product="toy protein"
                     /translation="MTETSEQ"
ORIGIN
        1 atgactgagacctctgaacagtaa

`atg act gag acc tct gaa cag taa` is MTETSEQ and a stop. The translation omits the stop. The format did not check that. This record is a toy.

2.2.4 SAM and BAM #

SAM is a text table of reads placed on a reference. BAM is the same table compressed. These notes do not specify that binary; convert it when a person must look. The columns that matter are name, flag, reference, start, mapping quality, CIGAR, sequence, and base qualities.

r1  0  chrToy  1  40  8M  *  0  0  ATGACTGA  IIIIIIII

Read `r1` is placed at position 1 of `chrToy`, mapping quality 40, CIGAR `8M`. The flag is bits: unmapped 4, reverse strand 16, secondary 256, duplicate 1024, supplementary 2048. Counting a duplicate invents depth. Counting a secondary line counts one molecule twice.

Mapping quality is the aligner’s doubt, on a Phred-like scale, that the read belongs at this locus. It is not a base quality. In CIGAR, `M` consumes one read base and one reference base and does not promise they match. `I` is extra in the read, `D` is missing from the read, `N` skips reference (an intron in an RNA-seq alignment, not a deletion), and `S` is soft-clipped sequence kept but not aligned. Operations that consume the read must sum to its length.

2.2.5 VCF #

A VCF line is a variant call against a stated reference, not a fact about the species. The columns are chromosome, position, identifier, reference allele, alternate, quality, filter, information, a format string, and sample columns.

chrToy  10  .  A  G  99  PASS  .  GT:DP  0/1:28

Position 10, reference A, alternate G, caller quality 99, filters passed, heterozygous, depth 28. `1/1` is homozygous alternate. A missing genotype is not homozygous reference. `PASS` means this caller’s filters accepted the site, not that the variant is real or pathogenic. Structural variants fit the same columns awkwardly, and their absence from a file of single-base calls is not evidence they are absent from the genome.

2.3 Names, versions, and trust #

2.3.1 Accessions #

An accession names the record. A version names the sequence at one moment. `TOY000001.1` and `TOY000001.2` may differ by one base or by the whole string, so the accession alone does not name a sequence. Older papers cite NCBI GI numbers; those were retired in favour of accession.version.

Coordinates need an assembly accession. “Chromosome 1, base 100000” is not a site until the build is named. Gene symbols are reused across species, and nucleotide and protein accessions are different namespaces.

2.3.2 How a false label spreads #

A weak similarity is deposited as `putative kinase`. The next pipeline copies `kinase`, and a third database imports that word as if it had been reviewed. No assay occurred. This transitive annotation error is ordinary.

Read the evidence code, and prefer a reviewed entry for a functional claim without treating the review as the assay. Check that the alignment covers the residues the function needs, and prefer a family model to the product name on the best hit. Identity to the wrong organism is a warning: a bacterial contig identical to a human sequence is contamination until shown otherwise. “Hypothetical protein” means the archive will not defend a function. It does not mean the chain does nothing.

2.4 Key terms, summary, and self-check #

Key terms. INSDC; GenBank; ENA; DDBJ; RefSeq; Swiss-Prot; TrEMBL; PDB; Gene Ontology; evidence code; FASTA; FASTQ; Phred score; SAM; BAM; CIGAR; VCF; accession; version; transitive annotation.

Summary. GenBank, ENA, and DDBJ are one exchanged nucleotide collection. RefSeq and Swiss-Prot are curated; TrEMBL is automatic; the PDB holds experimental coordinate models; GO terms mean nothing without an evidence code. FASTA is sequence, FASTQ adds a per-base error probability, a GenBank feature is an annotation, SAM is an alignment table and BAM its binary form, and VCF is a list of variant claims. An accession without a version does not name a sequence, and a copied function can have no experiment behind it.

Self-check.

1. Why are GenBank and ENA redundant, and why does that not make either a reviewer?

2. An `XP_` and an `NP_` align equally well. Why is the functional claim still different?

3. In Phred+33, what error probability does the character `5` claim?

4. In CIGAR, how does `N` differ from `D`?

5. The VCF filter says `PASS` and the gene product says `kinase`. What has been shown?