Putting It All Together and Where to Go Next

You have now met every stage of a genomics workflow one lesson at a time, and in this final lesson we line them all up so you can see the whole journey from raw data to a meaningful answer. A workflow is simply a fixed series of steps, and a genomics workflow turns messy sequencing output into a short, readable list of the exact places where a sample's DNA differs from the reference. Those differences are called variants, and finding and interpreting them is the heartbeat of modern genomics.

The Two Ingredients

Every analysis begins with two files. The reference genome is stored in FASTA format, the plain-text file that lists a species' DNA as one long string of A, C, G, and T letters, and it acts as the agreed-upon map everyone compares against. Your reads live in FASTQ format, which stores the short DNA fragments your sequencer produced from a real sample together with a quality score that says how sure the machine was of every letter.

# Prepare the reference so the tools can find any position instantly
bwa index reference.fasta
samtools faidx reference.fasta

# 1. MAPPING  (FASTQ -> BAM): line the reads up against the reference
#    R1 and R2 are the two ends of each paired-end DNA fragment
bwa mem reference.fasta reads_R1.fastq reads_R2.fastq > aligned.sam
samtools sort aligned.sam -o aligned.sorted.bam
samtools index aligned.sorted.bam

# 2. VARIANT CALLING  (BAM -> VCF): find where the sample differs
bcftools mpileup -f reference.fasta aligned.sorted.bam \
  | bcftools call -mv -Oz -o variants.vcf.gz
bcftools index variants.vcf.gz

# 3. ANNOTATION: label each variant with its gene and predicted effect
snpEff GRCh38.105 variants.vcf.gz > variants.annotated.vcf

The first two lines index the reference, which builds a lookup table so the tools can jump straight to any position instead of scanning three billion letters each time. The mapping step then runs bwa (short for Burrows-Wheeler Aligner), which finds where each short read fits best along the reference, and samtools, which sorts those alignments by position and saves them as a BAM file, a compressed table listing every read and exactly where it landed.

The variant-calling step feeds that BAM to bcftools: bcftools mpileup reads the pile of letters stacked at every position, and bcftools call decides where the sample truly differs from the reference. The -m flag selects its standard calling model, -v keeps only the sites that are real variants, and -Oz writes the result as a compressed VCF (Variant Call Format), a text file with one line per difference. SnpEff then reads that VCF and adds an annotation to each variant naming the gene it lands in and predicting whether it changes a protein.

bcftools view -H variants.annotated.vcf | head -1

20	32434638	.	A	G	67.2	.	DP=28;MQ=60;ANN=G|missense_variant|MODERATE|ASXL1|ENSG00000171456|transcript|ENST00000375687|protein_coding|c.1652A>G|p.Gln551Arg	GT:PL	1/1:95,18,0

Here the -H flag told bcftools to hide the long header block and print only the data line. Reading it left to right, the variant sits on chromosome 20 at position 32,434,638 where the reference letter A is replaced by G. Its INFO column is shown trimmed to a few key tags: DP=28 means 28 reads covered the site, MQ=60 is their average mapping quality, and SnpEff's ANN field, shortened here to its most useful parts, flags the change as a missense_variant, meaning it swaps one amino acid for another (glutamine for arginine, written p.Gln551Arg) in the ASXL1 gene. To see this with your own eyes, open a genome browser such as IGV (the Integrative Genomics Viewer), load the reference, the BAM, and the VCF, then jump to that position to watch the stacked reads and the highlighted difference.

Where To Go Next

From here you can branch into RNA-seq, which measures how strongly each gene is switched on rather than which letters differ, or into population genomics, which compares variants across thousands of people, and clinical genomics, which links variants to disease and treatment. Keep a file-format cheat sheet within reach: FASTQ holds raw reads, BAM holds alignments, VCF holds variants, and BED and GFF describe regions and gene features along the genome. For trustworthy help, lean on the samtools and bcftools manuals, the GATK best-practices documentation, the Ensembl and UCSC genome browsers, and the Biostars question-and-answer community.