How to Read Your Own DNA Sequence From Raw Files

By the end of this you will be able to take a gene name, pull the exact nucleotide sequence for that region out of your own genome, apply your own variants to it, translate it in the right frame on the right strand, and then go back to the aligned reads to check whether the sequence you just read is supported by real evidence or by two reads and an aligner’s guess. You need a machine with 16 GB of RAM and a few hundred GB of disk, a Unix shell, and your data: either a BAM/CRAM aligned to a named reference plus a VCF, or FASTQs if you are starting further back. Install samtools, bcftools, bedtools, seqkit, and gffread (conda-forge and bioconda have all of them: mamba install -c conda-forge -c bioconda samtools bcftools bedtools seqkit gffread). The web translators that rank for “DNA sequence reader” work fine on a 600 bp insert you pasted from a plasmid map. They are the wrong tool for a human genome, and the reason is that a human gene is spliced, strand-dependent, and diploid, none of which a text box knows about.

1. Identify what files you have

Four formats cover almost everything you will be handed.

.fastq.gz is unaligned reads with per-base quality. Each record is four lines: header, sequence, +, quality string in Phred+33 ASCII. Q30 means a 1-in-1000 chance the base call is wrong. A 30x human WGS FASTQ pair runs 40 to 60 GB compressed.

.bam / .cram is those reads aligned to a reference. CRAM is reference-compressed and roughly a third the size, which matters, and it is useless without the exact reference it was made against.

.vcf.gz is the list of positions where you differ from the reference, plus genotypes. A single-sample WGS VCF has 4 to 5 million variants, the large majority of them common SNPs.

.fa / .fa.gz is the reference genome itself, not your sequence.

Check what you have before anything else:

samtools quickcheck -v *.bam *.cram
samtools view -H sample.cram | grep -E '^@SQ' | head -3
samtools view -H sample.cram | grep -E '^@PG|^@RG'
bcftools view -h sample.vcf.gz | grep -E 'reference|contig=<ID=chr1,'

The @SQ lines tell you the reference. If contigs are named chr1, you are on a UCSC-style GRCh38 or hg19. If they are 1, you are on Ensembl/GRC naming. If SN:chr1 has LN:248956422, that is GRCh38. LN:249250621 is hg19. LN:248387328 is T2T-CHM13v2.0. Getting this wrong is the single most common way to read a sequence that is off by a few hundred kilobases and looks plausible anyway.

2. Get and index the matching reference

wget https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa
samtools faidx GRCh38_full_analysis_set_plus_decoy_hla.fa

faidx writes a .fai sidecar: contig name, length, byte offset, bases per line, bytes per line. That index is what makes random access to any coordinate an O(1) seek instead of a 3 GB scan. Confirm the md5s match your BAM header if it has M5: tags on the @SQ lines. If they do not match, stop and find the right reference. A CRAM decoded against the wrong reference will silently produce wrong bases at every position where the two references differ.

We use GRCh38 for anything where you want to compare against public annotation, ClinVar, gnomAD, or an existing VEP cache, because that is where the annotation lives. T2T-CHM13 is the better assembly and closes the segmental duplications and centromeres that GRCh38 leaves as gaps or collapsed copies, but the annotation ecosystem around it is still thinner. Pick one per project and do not mix.

3. Pull the sequence for a region

Coordinates are the second trap. SAM, VCF, GFF/GTF, and samtools region strings are 1-based and inclusive on both ends. BED files and the BAM binary format are 0-based half-open. A BED interval chr7 117559590 117559620 is the same 30 bases as the region string chr7:117559591-117559620. Off-by-one here shifts your reading frame and hands you a completely different protein.

Get an annotation so you are not typing coordinates by hand:

wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_45/gencode.v45.annotation.gtf.gz
zcat gencode.v45.annotation.gtf.gz \
  | awk -F'\t' '$3=="gene"' \
  | grep 'gene_name "CFTR"' \
  | cut -f1,4,5,7
# chr7  117287120  117715971  +

Then:

samtools faidx GRCh38_full_analysis_set_plus_decoy_hla.fa chr7:117287120-117715971 > CFTR_locus.fa
seqkit stats CFTR_locus.fa

That is 428,852 bases of genomic sequence. Almost all of it is intron. The coding sequence is about 4,443 bases across 27 exons. If you paste the genomic sequence into a web ORF finder, you will get back nonsense, because introns interrupt every frame. This is the central reason the generic tools fail on human data.

To get the actual transcript, use gffread, which does the splicing for you and handles strand:

zcat gencode.v45.annotation.gtf.gz | grep 'gene_name "CFTR"' > cftr.gtf
gffread cftr.gtf -g GRCh38_full_analysis_set_plus_decoy_hla.fa \
  -w cftr_transcripts.fa -x cftr_cds.fa -y cftr_protein.fa
seqkit fx2tab -nl cftr_cds.fa | head

-w is spliced exonic sequence, -x is CDS only, -y is the translated protein. For a gene on the minus strand, gffread reverse-complements for you. Verify by eye: every CDS should start with ATG and end with TAA, TAG, or TGA.

4. Strand and direction

DNA sequence is written 5’ to 3’, left to right, always, by convention. A reference FASTA gives you the plus strand (also called the forward or top strand) of each chromosome. Genes on the plus strand have an mRNA that reads the same as the reference with T in place of U. Genes on the minus strand, roughly half of them, have an mRNA that is the reverse complement of the reference at that locus.

The polymerase synthesizes 5’ to 3’ and reads the template 3’ to 5’. Both statements are true and they describe different molecules. What matters operationally: if the GTF says - in column 7, reverse complement before you translate.

seqkit seq -r -p region.fa > region_revcomp.fa   # -r reverse, -p complement

Note the two flags are separate. seqkit seq -r alone gives you a reversed sequence that is biologically meaningless.

5. Frames, codons, and the genetic code

Three bases per codon, 64 codons, 61 amino acid codons plus three stops (TAA, TAG, TGA). Reading a codon chart is just a base-3 lookup: first base picks the row block, second the column, third the sub-row. ATG is methionine and is also the canonical start. The code is degenerate mostly in the third position, which is why synonymous changes cluster there.

Any linear DNA has six reading frames: three on the plus strand (start at offset 0, 1, 2) and three on the minus strand.

seqkit translate -f 1 cftr_cds.fa | head -2      # frame 1 only
seqkit translate -f 6 region.fa                  # all six frames
seqkit translate --trim --clean -f 1 cftr_cds.fa # trim trailing stop, X out ambiguous

Use -T 2 for the vertebrate mitochondrial code if you are translating anything from chrM. The mitochondrial code differs in ways that will burn you: TGA is tryptophan rather than a stop, ATA is methionine, and AGA/AGG are stops. Translating chrM with the standard table gives truncated proteins that look like nonsense mutations.

For open reading frames in sequence you have not annotated (a plasmid, a viral contig, a de novo assembly), EMBOSS getorf is the tool we reach for:

getorf -sequence contig.fa -outseq orfs.fa -minsize 300 -find 1

-find 1 returns translations between stops, -find 3 returns nucleotide sequence, and -minsize is in nucleotides for -find 0/2/3 and amino acids for the protein modes, which is documented but easy to misread. In human genomic sequence, getorf output is not meaningful as gene prediction. Use the GENCODE annotation.

6. Make it your sequence, not the reference

Everything so far reads the reference. To read your sequence, apply your variants:

bcftools index sample.vcf.gz
samtools faidx GRCh38.fa chr7:117287120-117715971 > region_ref.fa

bcftools consensus -f region_ref.fa -s SAMPLE01 -H 1 \
  -M N --mark-del '-' sample.vcf.gz > region_hap1.fa
bcftools consensus -f region_ref.fa -s SAMPLE01 -H 2 \
  -M N --mark-del '-' sample.vcf.gz > region_hap2.fa

-H 1 and -H 2 take the first and second allele of each genotype. Unless your VCF is phased (check for | instead of / in the GT field), these are not real haplotypes. They are an arbitrary per-site assignment, and a compound scenario across two variants in the same gene will be represented incorrectly. Short-read WGS gives you phasing only locally, from read-backed phasing tools like WhatsHap, or across longer blocks from linked-read or long-read data. If phase matters for what you are looking at, say so out loud and get the right data.

-M N replaces positions with missing genotypes with N rather than silently emitting the reference base. We always set it. The default behavior makes low-coverage regions look like confident reference matches, which is the exact failure you are trying to avoid.

bcftools consensus writes a report to stderr of how many variants it applied. Read it. If you asked for a 400 kb region and it applied 12 variants, either your sample is unusually homozygous-reference there or your VCF and FASTA disagree on contig naming.

7. Go back to the reads

A VCF is a set of claims. The reads are the evidence. Short-read sequencing has systematic failure modes in repeats, segmental duplications, homopolymers, and GC extremes, and variant calls in those regions carry much lower confidence than the file format suggests ¹.

samtools tview -p chr7:117559590 -d T sample.cram GRCh38.fa
samtools mpileup -f GRCh38.fa -r chr7:117559590-117559600 sample.cram
samtools depth -a -r chr7:117559590-117559600 sample.cram

For a heterozygous SNP in a clean region you expect an allele fraction near 0.5 and depth near your genome-wide average. An allele fraction of 0.25 at a site with 90x coverage in a 30x genome is a collapsed duplication, meaning reads from two paralogous loci piled onto one, and the “variant” is a paralog difference. Check the MQ and MQ0 tags in the VCF INFO field. Mapping quality zero means the aligner could not place the read uniquely.

For anything you plan to act on, IGV is worth the GUI. Load your CRAM and the reference, jump to the position, group reads by strand, and look for strand bias (all support on one strand) and for variants clustering at read ends.

Interpretation is where most of the difficulty lives, not sequencing. The generation and annotation of variant lists is largely mechanical now. Determining which of your millions of variants has any functional consequence is the open problem ². Most of your variation sits outside protein-coding exons, in regions whose regulatory function has to be inferred from assay data rather than read off the sequence ³.

8. Annotate before you interpret

vep -i sample.vcf.gz -o sample.vep.vcf --vcf --cache --offline \
    --assembly GRCh38 --everything --fork 8 \
    --fasta GRCh38.fa --hgvs --canonical --pick_allele_gene

--hgvs gives you the standard nomenclature: NM_000492.4:c.1521_1523delCTT and NP_000483.3:p.Phe508del for the common CFTR deletion. Always carry the transcript accession with a version number. c. numbering is transcript-relative, so the same genomic position produces different c. coordinates on different transcripts of the same gene, and this is a frequent source of confusion when comparing to a published paper.

--pick_allele_gene chooses one consequence per gene per allele. Without a pick flag, VEP emits every transcript, which for a gene with 20 isoforms means 20 rows and a spreadsheet you will misread.

Anything you find that looks clinically significant belongs with a clinician and a CLIA-certified confirmatory test. Research-grade WGS is not a diagnostic result, and the difference is not pedantry: per-base error rates, coverage gaps, and the absence of orthogonal confirmation all bear on whether a call is real.

Common problems

Contig naming mismatch. chr1 vs 1. bcftools consensus will apply zero variants and exit 0. Fix with bcftools annotate --rename-chrs chr_map.txt.

CRAM cannot find its reference. Symptoms are a hang or an error mentioning M5. samtools falls back to downloading reference slices from EBI over the network. Set export REF_PATH=/path/to/cache/%2s/%2s/%s and export REF_CACHE= the same, or pass -T reference.fa explicitly to every samtools invocation.

Off-by-one from BED vs region string. If your translated protein starts with a plausible amino acid but the whole thing is gibberish after residue 1, you are one or two bases off. Shift and retranslate.

Translating a minus-strand gene without reverse complementing. You get a protein full of stops in every frame. Check column 7 of the GTF.

Using the genomic sequence instead of the CDS. Introns. Use gffread -x.

Assuming unphased heterozygous calls are haplotypes. They are not. bcftools consensus -H 1 on unphased data produces a chimera.

Missing genotypes silently become reference. Set -M N.

Treating an indel position as exact. Indels in homopolymer runs have ambiguous left-alignment. Run bcftools norm -f ref.fa -m -any before comparing VCFs from different pipelines, or identical variants will look different.

Gaps you cannot see. GRCh38 still has regions where short reads map poorly or not at all, and coverage-zero regions are absent from your VCF rather than flagged. Generate a callable-regions BED with mosdepth or GATK CallableLoci and check whether the locus you care about is in it. The functional density of the human genome is high enough that these gaps are not safely ignorable ³, and the assembly and annotation of even a single well-studied chromosome took years of dedicated effort to resolve ⁴. The throughput that made personal sequencing affordable came from massively parallel chemistry, not from any reduction in the difficulty of reading hard regions ⁵.

Oak builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

1. Daniel C. Koboldt, Li Ding, Elaine R. Mardis, et al. Challenges of sequencing human genomes. Briefings in Bioinformatics, 2010. https://doi.org/10.1093/bib/bbq016 ↩

2. Tuuli Lappalainen, Alexandra J. Scott, Margot Brandt, et al. Genomic Analysis in the Age of Human Genome Sequencing. Cell, 2019. https://doi.org/10.1016/j.cell.2019.02.032 ↩

3. Ewan Birney, Paul Flicek, Damian Keefe, et al. Identification and analysis of functional elements in 1% of the human genome by the ENCODE pilot project. Nature, 2007. https://doi.org/10.1038/nature05874 ↩ ↩²

4. Andrew J. Mungall, Sophie Palmer, Sarah Sims, et al. The DNA sequence and analysis of human chromosome 6. Nature, 2003. https://doi.org/10.1038/nature02055 ↩

5. Stephen P. A. Fodor. DNA SEQUENCING: Massively Parallel Genomics. Science, 1997. https://doi.org/10.1126/science.277.5324.393 ↩