This is the second half of practical 10, done here instead of in Galaxy. Two people, a mother and her child, each sequenced from genomic DNA enriched for mitochondria by long-range PCR — so the files still carry a good deal of nuclear DNA, which in this experiment is contamination. The job is to map what is there against the mitochondrial genome, throw out the reads that landed badly, call what differs, and then work out what fraction of the copies in each person carries it.

First the reference. GRCh38 is 3.26 GB and a browser cannot hold it, but it does not need to: the FASTA has a .fai beside it, the index turns a locus into a byte offset, and the host answers a range request. And chrM is 16,569 bases, so this is not a window onto a genome — it is the whole mitochondrial genome, fetched entire, for about seventeen kilobytes.

Now the reads, which are the practical's own. Download the four files from https://zenodo.org/record/1251112/files — raw_child-ds-1.fq, raw_child-ds-2.fq, raw_mother-ds-1.fq and raw_mother-ds-2.fq, about 14 MB each — and bring them onto the page with the cells below.

Look at what came off the sequencer before mapping any of it. Quality falls away at the ends of Illumina reads, and a tail of bad bases is how a mapper is made to place a read somewhere it does not belong — which matters more than usual here, because a variant present in three copies out of a hundred has to be told apart from an error present in three reads out of a hundred.

Before fifty thousand reads, one. The cell below pulls the first record out of what fastp just wrote, and the one after it puts that single read beside the reference and names every base where the two differ. Nothing in the practical gives you this look, and it is the one that makes a CIGAR string mean something.

If that read did not place at all, that is a result and not a failure: these samples are full of nuclear DNA by design, and a read from chromosome 4 has no business on chrM. Change RECORD in the cell above to look at another one, or map the lot below and see what fraction of them lands.

Now the mapping and the filtering, as commands. The practical does these as four Galaxy forms — BWA-MEM, merge, MarkDuplicates, Filter BAM — and the two that carry the argument are the mapping and the filter, so those are the ones written out here. Paste them into the terminal below, one line at a time.

Map the child's trimmed pair against the mitochondrial genome. BWA-MEM in the class; minimap2's short-read preset here, under its own name.

minimap2 -ax sr chrM.fa raw_child-ds-1.trimmed.fq raw_child-ds-2.trimmed.fq > child.sam

The practical's Filter BAM step, as the two numbers it really is: mapping quality 20 or better, and properly paired only. Then sort what survived into position order, and ask samtools how much of it there is.

samtools view -h -q 20 -f 2 child.sam > child.filtered.sam samtools sort child.filtered.sam -o child.sorted.sam samtools flagstat child.sorted.sam

Now the mother's sample, the same way, from the files you brought in.

fastp -i raw_mother-ds-1.fq -I raw_mother-ds-2.fq -o m1.fq -O m2.fq -q 20 -l 50 --cut_tail minimap2 -ax sr chrM.fa m1.fq m2.fq > mother.sam samtools view -h -q 20 -f 2 mother.sam > mother.filtered.sam

And call hers. FreeBayes in the class; bcftools here. This is exactly what the caller cell below does to the child's sample, which is the point of doing it by hand once.

samtools sort mother.filtered.sam -o mother.sorted.bam samtools index mother.sorted.bam bcftools mpileup -Q 20 -f chrM.fa mother.sorted.bam > mother.pileup.vcf bcftools call -mv mother.pileup.vcf > mother.vcf

How much of it mapped, how deep it is, how much is duplicate. With mitochondria-enriched DNA the mapped fraction is the interesting number: it says how well the enrichment worked, and it is the first thing anybody will ask you about this dataset.

A table saying there is a difference at a position is a claim; the reads stacked under it, with a minority carrying another base, is the evidence. Set the region to whatever position the caller flags below and come back to look at it.

Now let the caller decide. bcftools mpileup walks the alignment and reports what it saw at each position; bcftools call decides which of those are variants. The cell sorts and indexes first, because no caller will look at an unsorted alignment.

This is the practical's last step, and where it is worth going one further. Counting how many reads carried each allele gives a percentage; what it does not give is how well that percentage is known. Ninety reads out of a hundred and nine out of ten are both ninety per cent, and only one of them would survive somebody asking. So the cell below counts, and then puts an interval round the count, and where both samples were called it asks whether the two really differ.

Has anybody seen these positions before? dbSNP is 29.5 GB and the answer costs a few hundred kilobytes, because its tabix index turns a locus into a byte range: the database stays at NCBI and the question goes to it. The mitochondrial genome is one of the best described stretches of DNA there is, so most of what a healthy person carries here already has a name.

The practical ends this section by installing IGV and opening the files in it, and its first half ends by opening JBrowse. This is IGV — igv.js, the same code base, under its own name — over the reference, the reads and the calls already on this page. The cell sorts and indexes the alignment first, and shows you those commands and an IGV batch script that reproduces the same view on the desktop.

One step past where the practical stops. Thirteen of the mitochondrial genome's genes are proteins, all of them parts of the respiratory chain, and a base that changes one of them changes a residue in a machine you can look at. The cell below folds MT-CO1, cytochrome c oxidase subunit 1, by its UniProt accession; change the accession to whichever gene your own call fell in.

Everything above is compute, and compute is the cheap half. None of these tools decided whether the enrichment worked well enough to trust a three per cent allele, or whether a difference between mother and child is a real germ-line bottleneck or the two libraries having been made on different days. That judgement is the whole job, and it is the part the form in Galaxy hides by having a button for everything else.