Reads to variants: does this person carry CYP2C192*?

You will take raw sequencing reads from a human exome and turn them into a list of genetic variants – checking the reads, aligning them to the reference genome, sorting and indexing the alignments, calling and filtering variants – and finish by finding one variant that matters for prescribing: CYP2C19*2. You will run the pipeline twice: by typing commands, and in a Galaxy-style workflow engine.

About 2½ hoursNo software to installChrome, Edge, Firefox or Safari on a laptop or desktop

The clinical question

Clopidogrel is an antiplatelet drug given after heart attacks, stents and strokes. It is a prodrug: liver cytochrome P450 enzymes – mainly CYP2C19 – convert it to the active metabolite that blocks the platelet P2Y₁₂ receptor. People who carry a no-function CYP2C19 allele make less active drug and are at higher risk of clots on standard treatment. The commonest no-function allele, CYP2C19*2 (rs4244285), is carried on about 15% of chromosomes in Europeans and 30–36% in East and South Asians (1000 Genomes). Guidelines now recommend genotype-guided prescribing – for example CPIC for acute coronary syndromes and NICE (DG59) after ischaemic stroke or TIA.

Your job: start from the raw reads of one person and find out whether she carries CYP2C19*2 – and how confident you can be.

The data

The person is NA12878, a woman from a Utah family of northern and western European ancestry (CEPH collection) whose DNA has been sequenced more times than anyone else's. Her genome is the reference sample that laboratories use to test their pipelines (HapMap, 1000 Genomes, Genome in a Bottle). You have exome reads – DNA fragments captured from the exons – sequenced on an Illumina machine as paired-end 76-base reads. To keep things fast we only give you the read pairs from the region of chromosome 10 that holds the CYP2C gene cluster (CYP2C18, CYP2C19, CYP2C9, CYP2C8) – about 150,000 pairs instead of tens of millions.

Everything happens on this page

Instructions are on the left. The workbench on the right has four tabs: a Terminal (a Linux command line), a Report viewer for FastQC, the IGV genome browser and a Workflow engine that works like Galaxy.

samtools and bcftools are the real programs, compiled to run inside your browser (WebAssembly). FastQC and BWA are too big for a browser, so they replay real runs on exactly these data. Your files never leave your computer.
buttons type a command into the terminal for you – press Enter to run it. Try typing it yourself first: typing is how you learn it.
Activity steps tick themselves when the page sees you do them; you can also tick them yourself. Answers are saved in this browser – My answers (top right) downloads them.

Your name (optional – it goes on your downloaded answers):

The pipeline

FASTQ reads→quality control→align to hg19→sort & index (BAM)→call variants (VCF)→filter→inspect in IGV

Chapter	What you will do	Time
	Find your way around; look inside a FASTQ file; count reads	20 min
	Run FastQC and read the report	15 min
	Align the reads with BWA-MEM; read a SAM file	20 min
	Make a sorted, indexed BAM; where did the reads go?	15 min
	bcftools mpileup and call; read and filter a VCF; find CYP2C192*	25 min
	Check the call in the reads; read the codon	15 min
	Run the same pipeline in a Galaxy-style engine; finish and run a workflow	30 min
	Command line vs workflow engine; quiz	10 min

Sections marked Extension are optional and are not counted in your progress bar. The next practical, Variants to meaning, starts from the variants you make today, and your coursework uses the same steps on a patient's exome.

By the end you should be able to…

describe the FASTQ, SAM/BAM and VCF formats and read the key fields of each;
judge read quality from a FastQC report;
run an alignment and variant-calling pipeline on the command line, and explain what each step and option does;
check a variant call against the reads in a genome browser;
run the same pipeline in a workflow engine, find the command behind any result, and weigh up the two approaches.

The command line

≈ 20 minWorkbench: Terminal

Almost every sequence-analysis program is a command-line program: you type its name, options and files, and it writes files back. Genomics servers and clusters have no windows and mice, only a terminal like this one. The same few commands carry you a long way.

1.1 Where am I?

Activity 1.1 · Look around

Click in the terminal, type pwd (print working directory) and press Enter.
List what is here with ls.
Read the notes: cat README.txt
List the data and reference folders with sizes: ls -lh data ref (-l long format, -h human-readable sizes).

Command	What it does
`pwd`	Where am I? (the working directory)
`ls` · `ls -lh FOLDER`	List files · with sizes and dates
`cd FOLDER` · `cd ..` · `cd ~`	Change directory · go up one level · go home
`cat FILE` · `head -n 20 FILE` · `tail FILE`	Print a whole file · its first 20 lines · its last 10 lines
`zcat FILE.gz`	Print a gzip-compressed file, uncompressed
`wc -l` · `grep PATTERN`	Count lines · keep only lines containing PATTERN
Tab · ↑ / ↓	Complete a file name (saves typing and typos) · recall earlier commands

Paths such as data/NA12878_chr10_R1.fastq.gz are relative to where you are. The same file has the absolute path /home/student/practical/data/NA12878_chr10_R1.fastq.gz, and ~ is short for your home folder, /home/student.

The reference genome file ref/hg19.fa is about 3 GB. Why so big?

It contains every read ever sequenced from the human genome It is the whole human genome sequence: about 3.1 billion letters at one byte each It is compressed binary data used only by BWA

1.2 Inside a FASTQ file

Sequencers write reads in FASTQ format. The files are big, so they are almost always gzip-compressed (.gz). zcat uncompresses a file on the fly, and the pipe | sends its output straight into the next program – here head, which keeps the first lines.

Activity 1.2 · The first two reads

Show the first 8 lines of read file 1: zcat data/NA12878_chr10_R1.fastq.gz | head -n 8
Now the first read of file 2: zcat data/NA12878_chr10_R2.fastq.gz | head -n 4. Compare its name with the first read of file 1.

Line	Example	Meaning
1	`@SRR098401.94186032/1`	`@` then the read name. `/1` and `/2` mark the two reads of a pair.
2	`TTCTATTTTTATAAGG…`	The bases called by the sequencer (N = no call)
3	`+`	Separator (sometimes repeats the name)
4	`CAEFCCDDDDDCDGIH…`	One quality character per base

Each quality character encodes a Phred score Q = −10 log₁₀(p), where p is the probability that the base is wrong. The character's ASCII code minus 33 gives Q (“Phred+33”, also called Sanger or Illumina 1.8+).

Character	`#`	`+`	`5`	`?`	`@`	`I`
Phred score	2	10	20	30	31	40
Chance the base is wrong	63%	1 in 10	1 in 100	1 in 1,000	≈ 1 in 1,260	1 in 10,000

How long are the reads (in bases)?

76 bases – every read in these files is 76 long (FastQC's Basic Statistics agrees).

The second read in file 1 ends in #######. What quality is that, and what should an aligner or variant caller do with those bases?

# is ASCII 35 → Q2: a 63% chance of being wrong, i.e. no real information. Illumina's software marks the unreliable end of a read this way. The bases should count for (almost) nothing: BWA-MEM can soft-clip poor ends that do not match, and variant callers weight every base by its quality – bcftools ignores bases below Q13 by default. Trimming the tails off before alignment is another option.

The first reads of the two files have the same name except for /1 and /2. What does that tell you, and why is sequencing both ends of each fragment useful?

They are the two ends of the same DNA fragment (a read pair): read 1 from one end, read 2 from the other end on the opposite strand. Pairs are useful because the two reads must map close together, facing each other: that places reads in repeats and paralogous genes that a single short read could not place, flags poorly mapped reads, estimates the fragment (insert) size, and reveals structural variants when pairs map too far apart or in the wrong orientation. The files must stay in the same order – read n of file 1 pairs with read n of file 2.

1.3 How many reads?

Activity 1.3 · Count the reads

Count the lines in file 1: zcat data/NA12878_chr10_R1.fastq.gz | wc -l. Each read takes 4 lines.
Another way to count reads is to count the name lines, which start with @: zcat data/NA12878_chr10_R1.fastq.gz | grep -c '^@' (^ means “at the start of the line”, -c counts matches).
Try again with a longer pattern: zcat data/NA12878_chr10_R1.fastq.gz | grep -c '^@SRR'

How many reads are in file 1?

612,720 lines ÷ 4 = 153,180 reads (and the same number in file 2: 153,180 pairs, 306,360 reads in all).

grep -c '^@' gave a bigger number than the number of reads. Why? What does this teach you about checking your commands?

@ is also a quality character (Q31), so some quality lines (line 4) start with @ too – here 24,186 of them in file 1. '^@SRR' only matches read names and gives the right answer, 153,180. Lesson: a command that runs without an error can still give a wrong answer – check results against something you know (here, lines ÷ 4). Programs that understand the format (FastQC, seqkit stats) are safer.

Quality control with FastQC

≈ 15 minWorkbench: Terminal, Report

Before spending hours on alignment, check that the reads are worth aligning. FastQC summarises millions of reads in a few plots, each marked ✓Pass, !Warning or ✗Fail. The flags are rules of thumb designed for whole-genome shotgun data – you still have to think about what they mean for your data.

Activity 2.1 · Run FastQC

Run FastQC on both files, writing the reports into results/: fastqc -o results data/NA12878_chr10_R1.fastq.gz data/NA12878_chr10_R2.fastq.gz
See what it made: ls results
Open the first report: open results/NA12878_chr10_R1_fastqc.html – it appears in the Report tab.
Open the report for file 2 as well (↑ recalls the last command – edit it).

Work down each report (click a module in the summary box to jump to it; each has a short explanation). Hover over the plots to read the values.

Basic Statistics – number of reads, read length, GC content and the quality encoding.
Per base sequence quality – box plots of quality at each position along the reads. Green, orange and red bands mark very good (≥ 28), reasonable and poor (< 20) quality.
Sequence Duplication Levels and Adapter Content – signs of PCR over-amplification and of fragments shorter than the reads.

Which file has a Fail, in which module – and what does the plot show?

File 1 – Adapter Content: adapters in most reads File 2 – Per base sequence quality: quality falls towards the end of the reads Both files – Sequence Duplication Levels

Why does quality fall along the reads, and why is read 2 usually worse than read 1?

Illumina sequencing reads one base per chemistry cycle. In each cycle a few molecules in a cluster fail to extend or extend twice, so the copies in a cluster drift out of step (phasing and pre-phasing); the signal gets noisier with every cycle and the quality scores fall. Read 2 is sequenced after read 1 and after the fragments have been re-synthesised on the flow cell, so its clusters have been through many more cycles and chemistry steps – it is typically worse.

Would you trim these reads before aligning them? Give one reason for and one against.

For: about a fifth of the read 1s and nearly a third of the read 2s end in a run of Q2 (#) bases; trimming them removes noise and helps some aligners. Against / not essential here: BWA-MEM does local alignment and soft-clips ends that do not match; bcftools discards bases below Q13 and weights the rest by quality, so the poor tails hardly affect the calls. Adapter content is below 1%, so adapter trimming is not needed. Many pipelines trim anyway (e.g. fastp, Trimmomatic, Cutadapt); what matters is to decide, and to record the decision.

FastQC estimates that about 79% of the reads would remain if exact duplicates were removed. Why do exome data usually show some duplication, and what do pipelines do about it?

Exome libraries are amplified by PCR and the captured targets are sequenced very deeply, so many reads are copies of the same original fragment (and deep sequencing of a small target also produces some identical fragments by chance). Duplicates are not independent evidence, and a PCR error copied many times can look like a variant. Pipelines therefore mark duplicates after alignment (e.g. samtools markdup, Picard MarkDuplicates), and callers such as bcftools ignore reads flagged as duplicates. We skip that step today to keep the pipeline short – it is one of the “clean-up” steps a production pipeline adds.

Aligning the reads to the genome

≈ 20 minWorkbench: Terminal

A read on its own is 76 letters with no address. Alignment (or mapping) finds where in the reference genome each read comes from, allowing for sequencing errors and for real differences – the variants we are looking for.

3.1 The reference genome

Activity 3.1 · Meet hg19

Look at the start of the reference: head -n 3 ref/hg19.fa – a FASTA file: a >name line, then the sequence.
The .fai index lists every sequence and its length: cut -f 1,2 ref/hg19.fa.fai

How long is chromosome 10 in hg19?

135,534,747 bases.

Genome builds

The human reference is updated from time to time, and positions change between versions. hg19 (UCSC name) = GRCh37 (Genome Reference Consortium name), released in 2009; the current build is GRCh38/hg38 (2013). CYP2C19*2 is at chr10:96,541,616 on GRCh37 but chr10:94,781,859 on GRCh38. Use one build for everything – reads, annotation and databases. Your coursework data are aligned to hg19, so this practical uses hg19 too (and so should VEP and gnomAD when you look variants up: choose GRCh37).

3.2 The index

Scanning 3.1 billion bases for each of 300,000 reads would take far too long, so aligners first build an index of the genome. BWA builds an FM-index, a compressed index based on the Burrows–Wheeler transform (the same transform the bzip2 compression program uses): with it BWA can find where any short sequence occurs in the genome almost instantly. The index files sit next to the reference: hg19.fa.amb, .ann, .bwt, .pac and .sa.

Activity 3.2 · Index (already done)

Type the indexing command to see what happens: bwa index ref/hg19.fa

3.3 Align with BWA-MEM

bwa mem takes the reference (it finds the index files from the name) and the two read files, and writes alignments in SAM format to the screen. The redirect > sends that output into a file instead. Messages about progress go to a separate channel (standard error), so they still appear on screen.

Activity 3.3 · Run the alignment

Align the read pairs: bwa mem ref/hg19.fa data/NA12878_chr10_R1.fastq.gz data/NA12878_chr10_R2.fastq.gz > results/aligned.sam
While it runs, read the messages. BWA reads the pairs in batches and, from the pairs that map uniquely, estimates the size of the DNA fragments (the insert size).

What is the mean insert size (fragment length) BWA estimated? Give the number in bp.

About 169 bp (standard deviation about 64 bp; median 158 bp). Two 76-base reads cover 152 bases, so nearly half of the fragments are short enough for the two reads to overlap in the middle – the same bases are read twice, once from each end. Adapter sequence appears only when a fragment is shorter than a single read (76 bases), which is rare.

3.4 Reading a SAM file

Activity 3.4 · Inside aligned.sam

The file starts with a header of lines beginning with @: head -n 5 results/aligned.sam
Skip the header to see the first read pair: grep -v '^@' results/aligned.sam | head -n 2 (-v keeps the lines that do not match).

Each alignment line has 11 tab-separated columns, then optional tags:

#	Column	Meaning
1	QNAME	Read name (the same for both reads of a pair)
2	FLAG	A number whose bits answer yes/no questions about the read (table below)
3–4	RNAME, POS	Chromosome and leftmost position of the alignment
5	MAPQ	Mapping quality: Phred-scaled confidence that this is the right place. 60 is BWA's maximum; 0 means the read fits equally well elsewhere
6	CIGAR	How the read aligns: `76M` = 76 aligned bases; `S` soft-clipped, `I` insertion, `D` deletion
7–9	RNEXT, PNEXT, TLEN	Where the mate aligned (`=` same chromosome) and the fragment length
10–11	SEQ, QUAL	The bases and qualities (reverse-complemented if the read aligned to the reverse strand)
12+	tags	`NM` number of differences, `MD` where they are, `AS` alignment score, `XS` score of the next-best place

FLAG bit	Meaning	FLAG bit	Meaning
1	read is paired	64	first read of the pair
2	pair mapped properly (right distance and orientation)	128	second read of the pair
4	read unmapped	256	secondary alignment
8	mate unmapped	1024	PCR or optical duplicate
16	read on the reverse strand	2048	supplementary (part of a split read)
32	mate on the reverse strand

The first alignment line has FLAG 99. Break it into bits: which bits are set?

99 = 64 + 32 + 2 + 1: the read is paired (1), in a proper pair (2), its mate is on the reverse strand (32) and it is the first read of the pair (64). Its mate, on the next line, has FLAG 147 = 128 + 16 + 2 + 1: second read, on the reverse strand, proper pair. The two reads face each other, as a properly paired fragment should.

3.5 How well did it map?

Activity 3.5 · samtools flagstat

Summarise the alignments by their flags: samtools flagstat results/aligned.sam

What percentage of the reads mapped?

98.85%. 97.38% were properly paired: both reads mapped to the same chromosome, facing each other, a sensible distance apart. The total (306,496) is a little more than the 306,360 reads because a read split across two places counts once more as a supplementary alignment (136 here).

Sort, index and summarise

≈ 15 minWorkbench: Terminal

The SAM file lists the reads in the order they came off the sequencer. To look at – or call variants in – one region, programs need the reads sorted by position, stored compactly as BAM (binary, compressed SAM), and indexed (a .bai file that says where in the BAM each region starts).

Activity 4.1 · Make a sorted, indexed BAM

Sort by position and write BAM: samtools sort -o results/aligned.sorted.bam results/aligned.sam
Index it: samtools index results/aligned.sorted.bam – this runs the real samtools in your browser.
Compare the sizes: ls -lh results

How big are the SAM and the BAM files? Why do we keep BAM rather than SAM?

SAM ≈ 77 MB, sorted BAM ≈ 22 MB – about 3.5 times smaller. BAM is compressed binary (blocks of gzip, “BGZF”), so it is smaller and faster to read, and – once sorted and indexed – programs can jump straight to any region instead of reading the whole file. A full human exome BAM is typically 5–15 GB; the SAM would be several times larger.

Try this Extension

What happens if you do the steps in the wrong order? Try samtools index results/aligned.sam and read the error. Error messages are clues, not failures.

4.2 Where did the reads go?

Activity 4.2 · samtools idxstats

Count the reads on each chromosome (this uses the index): samtools idxstats results/aligned.sorted.bam. The columns are: name, length, mapped reads, unmapped reads.
Too long to read? Sort by the third column, biggest first, and keep the top 5: samtools idxstats results/aligned.sorted.bam | sort -k3,3nr | head -n 5

How many reads mapped to chr10?

300,965 – 99.3% of the mapped reads.

These reads were chosen because they came from chromosome 10. Why have some mapped to other chromosomes?

About 2,000 reads (0.7% of the mapped reads) landed elsewhere – a few dozen to a couple of hundred per chromosome. 3,516 reads did not map at all: most are listed beside their mapped mate (idxstats' last column), and 390 from pairs where neither read mapped are on the * line. The reads were chosen because the 1000 Genomes pipeline had placed them in chr10:95–98 Mb, but a short read that matches several places equally or almost equally well (repeats such as Alu and LINE elements, pseudogenes, related genes) can be placed differently by another aligner version or reference build. Most of these reads have low mapping quality, and very few pass the depth filter when variants are called. With a whole exome, the same effect spreads a few reads from every gene across the genome.

4.3 One position

With an index, samtools can pull out just the reads covering a region, written chromosome:start-end.

Activity 4.3 · Reads at the CYP2C19*2 position

Count the reads that cover chr10:96,541,616: samtools view -c results/aligned.sorted.bam chr10:96541616-96541616
Look at two of them: samtools view results/aligned.sorted.bam chr10:96541616-96541616 | head -n 2

How many reads cover chr10:96,541,616?

109 alignments overlap that base. That is deep coverage – plenty to call a genotype confidently. (The variant caller will count a few less – chapter 6 explains why.)

Calling variants

≈ 25 minWorkbench: Terminal

At every position, the reads either agree with the reference or show another base. A variant caller weighs that evidence – how many reads, how good their bases and mappings are – and decides the most likely genotype: two copies of the reference allele (0/0), one of each (0/1, heterozygous) or two copies of the other allele (1/1, homozygous).

Two steps in bcftools

bcftools mpileup “piles up” the reads at every covered position and computes genotype likelihoods: how probable the reads would be if the genotype were 0/0, 0/1 or 1/1. bcftools call combines those likelihoods with a prior and reports the variant sites and their most likely genotypes. Because the pile-up has a line for every covered base (millions), it is piped straight into call rather than saved.

5.1 Run the caller

Activity 5.1 · mpileup | call

First, peek at the pile-up for five bases around the CYP2C19*2 position: bcftools mpileup -a AD,DP -f ref/hg19.fa -r chr10:96541614-96541618 results/aligned.sorted.bam | grep -v '^##'. Look at the line for 96541616: the ALT column offers A and the AD tag counts the reads for each allele.
Now call variants across all the reads (about a minute in the browser): bcftools mpileup -Ou -a AD,DP -f ref/hg19.fa results/aligned.sorted.bam | bcftools call -mv -Ov -o results/calls.vcf

Option	Meaning
`-f ref/hg19.fa`	the reference, to know what base is expected at each position
`-a AD,DP`	also report, for each sample, the number of reads supporting each allele (AD) and the depth (DP)
`-Ou`	output uncompressed BCF (binary VCF) – fastest for passing to the next program
`-m`	the multiallelic calling model (recommended; you must choose it or the older `-c` model)
`-v`	output variant sites only (not every position)
`-Ov -o FILE`	write uncompressed VCF to FILE

5.2 Reading a VCF file

Activity 5.2 · Inside calls.vcf

The header lines start with ## and describe every field; skip them to see the column names and the first records: grep -v '^##' results/calls.vcf | head -n 4
Count the variants (every line that does not start with #): grep -c -v '^#' results/calls.vcf
What do the FORMAT fields mean? grep '^##FORMAT' results/calls.vcf

Column	Meaning
CHROM, POS	Where (1-based position of the first REF base)
ID	A database identifier such as an rs number, if added (`.` = none)
REF, ALT	Reference allele and alternative allele(s)
QUAL	Phred-scaled confidence that the site is truly variant (QUAL 30 = 1 in 1,000 chance it is not)
FILTER	`PASS`, `.` (not filtered), or the name of a filter the site failed
INFO	Site information, e.g. `DP` (read depth), `MQ` (mapping quality)
FORMAT + sample	Per-sample fields named by FORMAT: `GT` genotype, `PL` genotype likelihoods, `DP` depth, `AD` reads per allele

How many variants did bcftools call?

2,744. That is far more than the region's exons could hold – many are unreliable calls from a handful of reads. Next we filter.

5.3 Filter the calls

A raw call set always contains false positives: calls made from one or two reads, calls in reads that map badly, sequencing or PCR errors. The simplest filters require a minimum call quality and read depth.

Activity 5.3 · Keep the confident calls

Keep calls with QUAL ≥ 30 and at least 10 reads: bcftools view -i 'QUAL>=30 && INFO/DP>=10' results/calls.vcf > results/calls.filtered.vcf
Count what is left (↑ and edit the earlier grep).

How many variants pass the filter?

218 of 2,744 (8%). Most of the removed calls had fewer than 10 reads: exome capture pulls down the exons strongly, but some reads come from off-target DNA around them, where coverage is thin.

Filtering is a trade-off. What do you lose by making the thresholds stricter (say QUAL ≥ 100 and DP ≥ 30)? When might you accept that?

Stricter filters remove more false positives but also more true variants (lower sensitivity) – especially in regions with low coverage, which in an exome include the ends of exons and GC-rich exons. For a diagnostic search you usually want high sensitivity and then check candidates by eye (IGV) or by another method (Sanger sequencing); for a large research set you may prefer high specificity. The filter thresholds are part of your methods and must be reported.

5.4 Find CYP2C192*

Activity 5.4 · The answer to the question

Search the filtered calls for the position: grep -w 96541616 results/calls.filtered.vcf (-w: whole word)

What are the REF and ALT alleles, and the genotype (GT)? How many reads support each allele (AD)?

REF G, ALT A, GT 0/1 – heterozygous. AD = 44,48: 44 reads carry G and 48 carry A, close to the 50:50 expected for a heterozygote. QUAL is 222 – as high as bcftools usually goes (confident calls top out around 220–230). NA12878 carries one copy of CYP2C19*2.

One copy of a no-function allele: what would you predict for NA12878's clopidogrel response? What else would you want to check before saying so in a report?

With one no-function allele she is most likely a CYP2C19 intermediate metaboliser (e.g. *1/*2): less active metabolite, weaker platelet inhibition; CPIC advises an alternative (prasugrel or ticagrelor) after ACS/PCI if there is no contraindication. Before reporting you would check: the call in the reads (IGV); other CYP2C19 alleles – e.g. *3, and the increased-function *17 in the promoter, which an exome may not cover (no call is not the same as a reference genotype); and ideally confirm with an independent method. The next practical works through this.

Extension: a neat table, and a second caller Extension

Activity 5.5 · bcftools query and FreeBayes

bcftools query prints chosen fields: bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%QUAL\t[%GT\t%AD]\n' results/calls.filtered.vcf | head
The original version of this practical used FreeBayes, a haplotype-based caller. Run it: freebayes -f ref/hg19.fa results/aligned.sorted.bam > results/freebayes.vcf, count the calls, and look for 96541616.

Do the two callers agree about CYP2C19*2? Do they find the same number of variants? Why might they differ?

Both call chr10:96541616 G>A as heterozygous. FreeBayes with default settings reports 4,142 sites – it reports almost any candidate and leaves filtering to you; with --standard-filters --min-coverage 10 (Galaxy's “simple diploid calling with filtering and coverage”) it reports 268. The callers use different models (bcftools: per-position genotype likelihoods; FreeBayes: short haplotypes, so neighbouring variants and indels are described differently) and different default read and base filters, so call sets differ mostly at low-coverage and messy sites. Pipelines are usually compared on a truth set – for NA12878, the Genome in a Bottle benchmark calls.

See it in the genome browser

≈ 15 minWorkbench: IGV

Never trust a variant call you have not looked at. The Integrative Genomics Viewer (IGV) shows the reads stacked against the reference so you can judge the evidence yourself – this is the same igv.js used by the web version of IGV at igv.org.

Activity 6.1 · Load your files

Load your BAM and filtered VCF and go to the variant – from the terminal: igv results/aligned.sorted.bam results/calls.filtered.vcf chr10:96,541,616 or with this button:
If you are not there yet, type chr10:96,541,616 in IGV's search box and press Enter.

From top to bottom you see: the chromosome ideogram, your VCF (one mark per variant), the coverage track (read depth – a bar is coloured where many reads disagree with the reference), the reads (coloured by strand; mismatching bases are coloured letters), the reference sequence and the RefSeq genes.

Activity 6.2 · Weigh the evidence

Click the coverage bar at the centre line. The pop-up gives the total depth and the count of each base.
Click the variant in the VCF track to see its QUAL, genotype and AD.
Are the A bases on reads from both strands (both colours)? Are there reads with MAPQ 0 (drawn faint)? Anything else suspicious – nearby clusters of mismatches, soft-clipped ends?

What counts does IGV give for G and A? Why might they differ slightly from the AD in your VCF?

IGV shows about 51 G and 58 A (all reads). bcftools counted 44 and 48 because it ignores bases with quality below 13, judged after two adjustments: where both reads of a fragment cover the position (9 fragments here) it counts that DNA molecule only once, and its BAQ step lowers the quality of bases that might be misaligned. (Add -Q 0 to bcftools mpileup and it also reports 51 and 58.) The two views agree on what matters: roughly half the reads carry A, on both strands, in well-mapped reads – a clean heterozygous call.

Activity 6.3 · Which codon?

Zoom out a little (−) to see where you are in CYP2C19. The small arrows on the gene show the direction of transcription. Which exon is it?
Zoom right in until the amino acids appear under the exon (letters in alternating shades). Hover over the amino acid at the centre line to see its number.

Which amino acid is encoded at the variant, and what is its number?

Proline 227 (P227), near the start of exon 5. CYP2C19 is on the plus strand here (transcribed left to right), so you read the codon directly: CCG → CCA. Both codons encode proline – the change is synonymous.

The change does not alter the amino acid – yet CYP2C19*2 makes no working enzyme. Which is the most likely explanation?

Proline is translated more slowly from CCA than from CCG It changes how the mRNA is spliced The call is a sequencing error

Activity 6.4 · A figure for your report

Frame a clear view of the variant (zoomed in, variant centred) and press Save image above the browser. Figures like this, with a clear legend, are expected in your coursework.

The same analysis in a workflow engine

≈ 30 minWorkbench: Workflow engine

A workflow engine runs a chain of tools for you and records exactly what was done. Galaxy (usegalaxy.eu, usegalaxy.org) does this through a web page: you pick tools from menus instead of typing, every result is kept in a history, and a sequence of steps can be saved as a workflow and re-run on new data with one click. Research groups also use code-based workflow managers such as Nextflow and Snakemake, where the workflow is a script. The engine on this page is a small Galaxy look-alike that runs the same programs as your terminal.

7.1 A tour

Activity 7.1 · Find your way around

Open the Workflow engine tab. (The tutorial panel narrows to make room – drag the divider to change it.)
The history (right) already holds the two FASTQ files, imported from the course data library. Click a dataset's name to expand it, and its eye icon to view it.
The tools (left) are grouped by job; type in the search box to find one.

7.2 Run one tool

Activity 7.2 · FastQC, the Galaxy way

Click FastQC in the tool panel.
Choose dataset 1 (the R1 reads) and press Run tool. Watch the new datasets in the history: grey = queued, yellow = running, green = done.
Click the eye icon of the Webpage output – it opens in the Report tab, the same report you made in the terminal.
Click the ⓘ (dataset details) icon of the same output. Scroll to Command line.

Compare Galaxy's command line with the one you typed in chapter 2. What is the same? What did Galaxy add, and why?

The core is the same program: fastqc on the same file. Galaxy first links the dataset into a job folder under a readable name (ln -s …dataset_1.dat NA12878_chr10_R1.fastq.gz), then adds options for running unattended on a server: --outdir ., --threads (the processors it was given), --quiet, --extract (unzip the results so the raw data can become a second output) and -f fastq (don't guess the format from the name). Finally it copies the outputs to the dataset files shown in the history. Every result in Galaxy carries this record – its provenance.

7.3 Finish the workflow

Your lecturer has shared a workflow that chains the steps of chapters 2–5 – but it stops at bcftools call: the filtering step is missing.

Activity 7.3 · Edit the workflow

Click Workflows in the top bar of the engine.
Press Edit on “Reads to variants (bcftools)”. Each box is a step; each wire carries a dataset from an output (right) to an input (left). The small labels show the data types. Click a step to see its settings on the right.
With the editor open, click bcftools view in the tool panel to add it.
Try connecting the wrong thing: drag from FastQC's Webpage (html) output to bcftools view's input. What happens?
Now drag from bcftools call's vcf output to bcftools view's input.
Click the bcftools view step and check its Include expression: QUAL>=30 && INFO/DP>=10 – the same filter you used on the command line.

Why did the engine refuse the first connection? How does this help you?

Every dataset has a type (fastqsanger.gz, bam, bcf, vcf, html…) and every tool input accepts only certain types. An HTML report cannot be filtered as a VCF, so the connection is refused. Type checking catches mistakes when you build the workflow, not hours later when a job fails – one of the safety nets a workflow engine gives you (on the command line nothing stops you from giving a program the wrong file).

7.4 Run it

Activity 7.4 · One click, eight jobs

Press Run (in the editor, or on the Workflows list). Check that the forward reads input is the R1 dataset and the reverse reads is R2, then press Run workflow.
Watch the history fill up. The mpileup step reads the whole BAM file, so give it a minute or two.
Open (eye) the final VCF from bcftools view. How many variants does its history entry report?
Click display with IGV on it, and go to chr10:96,541,616.
Open the ⓘ details of the final VCF – the command line, parameters and input datasets are all recorded.

How many variants are in the workflow's filtered VCF – the same as your command-line result? Why?

218 – identical, because the workflow ran the same programs, with the same versions, on the same data and reference, with the same parameters. The only difference is that Galaxy's BWA-MEM step pipes straight into samtools sort (no SAM file is kept), which does not change the alignments. (One subtlety: Galaxy runs BWA-MEM with two threads. BWA-MEM reads its input in batches whose size depends on the number of threads and estimates the insert size batch by batch, so on a real server a handful of alignments can differ between thread counts unless the batch size is fixed with -K. Here the engine re-uses the terminal's alignments, so they match exactly.) Same inputs + same tools + same parameters = same result: that is reproducibility – and a workflow engine records all three for you.

Extension: add a second caller Extension

Add FreeBayes to the workflow, connect it to the BWA-MEM BAM output, run the workflow again and compare the two VCFs in IGV. Which differences do you see near CYP2C19?

Command line or workflow engine?

≈ 10 min

You have now run the same pipeline twice. Both are used in real laboratories – often side by side.

	Command line	Workflow engine (Galaxy)
Getting started	Steeper: commands, options and file names must be exact	Gentle: forms, menus and sensible defaults
Record of what was done	Only if you keep it: shell history, scripts, notes	Automatic: every dataset keeps its tool, version, parameters and inputs
Re-running on new data	Write a script (or a Nextflow/Snakemake pipeline)	Run the saved workflow on a new history
Flexibility	Any tool, any option, any combination; easy to automate and loop	The tools and options installed on that server
Computing power	Your own machine or cluster – you install and manage the software	The server's – shared, with queues and storage quotas
Data	Stays where you put it	Uploaded to the server – check the rules for patient data
Checking your work	Nothing stops a wrong file or option	Typed connections catch some mistakes

Your history is a script

Type history in the terminal. The commands that worked, in order, are your pipeline. Saved in a file they become a script that anyone can re-run:

#!/bin/bash
# reads-to-variants.sh – NA12878, CYP2C region, hg19
set -euo pipefail
fastqc -o results data/NA12878_chr10_R1.fastq.gz data/NA12878_chr10_R2.fastq.gz
bwa mem ref/hg19.fa data/NA12878_chr10_R1.fastq.gz data/NA12878_chr10_R2.fastq.gz > results/aligned.sam
samtools sort -o results/aligned.sorted.bam results/aligned.sam
samtools index results/aligned.sorted.bam
bcftools mpileup -Ou -a AD,DP -f ref/hg19.fa results/aligned.sorted.bam | bcftools call -mv -Ov -o results/calls.vcf
bcftools view -i 'QUAL>=30 && INFO/DP>=10' results/calls.vcf > results/calls.filtered.vcf

Activity 8.1 · Reflect

Type history and compare it with the script above.

A clinical lab must analyse 200 exomes a month and show an inspector exactly how every result was produced. A PhD student wants to try five aligners with different settings on one sample. Which approach suits each, and why?

Clinical lab: a fixed, validated, version-controlled pipeline run by a workflow system (Galaxy, or Nextflow/Snakemake on a cluster) – automatic provenance, the same steps every time, easy to audit, scales to many samples. PhD student: the command line (or a notebook/script) – maximum flexibility to change tools and options quickly; once the best combination is found, it can be captured in a script or workflow. Both need records: a lab's standard operating procedure, a student's scripts and notes.

Check your understanding

A base has Phred quality 30. The chance it is wrong is…

1 in 10 1 in 100 1 in 1,000 1 in 10,000

Which order is right?

align → index → sort → call align → sort → index → call call → align → sort → index

You email a colleague aligned.sorted.bam but IGV on their computer refuses to show it. What did you forget?

The index file, aligned.sorted.bam.bai The SAM file The reference genome

A VCF line shows GT:PL:DP:AD 1/1:230,60,0:20:0,20. The person is…

heterozygous homozygous for the ALT allele homozygous for the reference allele

Why did some reads from chromosome 10 map to other chromosomes?

BWA was told to use every chromosome equally Their sequence also occurs, or nearly occurs, elsewhere in the genome The sample was contaminated with bacterial DNA

Six months later a reviewer asks which filter you used on your variants. Where is the answer if you used the workflow engine?

In the dataset details (ⓘ) of the filtered VCF Nowhere – you should have written it down

What next?

You know where NA12878 differs from the reference. In Variants to meaning you will annotate the same variants – which gene, which consequence, how common, how damaging – and look at the affected proteins in 3D. Your coursework then applies the whole pipeline to a patient's exome.

Cheat sheet

The commands

Step	Command
Look at reads	`zcat R1.fastq.gz \| head -n 8`
Count reads	`zcat R1.fastq.gz \| wc -l` (÷ 4)
Quality control	`fastqc -o results R1.fastq.gz R2.fastq.gz`
Index a genome (once)	`bwa index ref.fa` · `samtools faidx ref.fa`
Align	`bwa mem ref.fa R1.fastq.gz R2.fastq.gz > aligned.sam`
Sort + index	`samtools sort -o aligned.sorted.bam aligned.sam` · `samtools index aligned.sorted.bam`
Summaries	`samtools flagstat FILE.bam` · `samtools idxstats FILE.bam`
Reads in a region	`samtools view -c FILE.bam chr10:96541616-96541616`
Call variants	`bcftools mpileup -Ou -a AD,DP -f ref.fa FILE.bam \| bcftools call -mv -Ov -o calls.vcf`
Filter	`bcftools view -i 'QUAL>=30 && INFO/DP>=10' calls.vcf > filtered.vcf`
Table of fields	`bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t[%GT]\n' filtered.vcf`

File formats

Format	Holds	Text or binary
FASTA (.fa)	Reference sequences	text
FASTQ (.fastq.gz)	Reads with qualities	text, usually gzipped
SAM (.sam)	Alignments	text
BAM (.bam) + BAI (.bai)	Alignments (compressed) + index	binary
VCF (.vcf, .vcf.gz) / BCF	Variant calls	text / binary

Glossary

Exome	The protein-coding exons of all genes (about 1–2% of the genome), captured by hybridisation to probes before sequencing
Coverage / depth	How many reads cover a base
Paired-end	Both ends of each DNA fragment are sequenced (read 1 and read 2)
Insert size	Length of the sequenced DNA fragment
MAPQ	Mapping quality – confidence the read is placed correctly
Genotype 0/0, 0/1, 1/1	Homozygous reference, heterozygous, homozygous alternative
*Star allele (2)**	Pharmacogene haplotype names defined by PharmVar
Provenance	The record of how a result was produced: tools, versions, parameters, inputs

Data, software and credits

Reads: 1000 Genomes Project exome sequencing of NA12878 (run SRR098401; 1000 Genomes Project Consortium (2015) Nature 526:68–74); read pairs from chr10:95–98 Mb of the project's alignment, re-aligned here to hg19. Reference: UCSC hg19 / GRCh37 (Genome Reference Consortium).
samtools 1.17, bcftools 1.10 and htslib (Danecek P et al. (2021) GigaScience 10:giab008; MIT licence) compiled to WebAssembly by the biowasm project and run with Aioli (Aboukhalil R, MIT licence).
FastQC 0.12.1 (Andrews S, Babraham Bioinformatics), BWA-MEM 0.7.19 (Li H (2013) arXiv:1303.3997) and FreeBayes 1.3.10 (Garrison E & Marth G (2012) arXiv:1207.3907) were run on the same data to record their output.
Genome browser: igv.js (Robinson JT et al. (2023) Bioinformatics 39:btac830; MIT licence). Gene models: NCBI RefSeq via UCSC.
The workflow engine imitates Galaxy (The Galaxy Community (2024) Nucleic Acids Res 52:W83–W94); tool names and versions follow the Galaxy tool shed.
Pharmacogenetics: CPIC guideline for CYP2C19 and clopidogrel (Lee CR et al. (2022) Clin Pharmacol Ther 112:959–967); NICE diagnostics guidance DG59 (2024).
Adapted from the module's Galaxy variant-calling practical and IGV tutorials.