Repository navigation
Tutorial
This tutorial runs BACoN on public data: whole-genome Nanopore reads of 28 potato (Solanum tuberosum) cultivars from the Ural region, sequenced to compare their plastomes (NCBI BioProject PRJNA807056; MinION, reads longer than 6 kb, 100–700 Mb per sample). Potato cultivars carry a few types of plastome ("cytoplasm types"); T-type plastomes carry a 241 bp deletion between the ndhC and trnV-UAC genes that the other types lack (Kawagoe & Kikuta 1991; Hosaka 2002). Let's see what BACoN finds.
About 7 GB. The script downloads the 28 runs from ENA, named after their cultivar, and checks their MD5:
mkdir potato && cd potato
curl -s "https://www.ebi.ac.uk/ena/portal/api/filereport?accession=PRJNA807056&result=read_run&fields=sample_alias,fastq_ftp,fastq_md5&format=tsv" \
| tail -n +2 > runs.tsv
mkdir -p reads
while IFS=$'\t' read -r run alias url md5; do
name=$(echo "$alias" | tr ' ' '_')
curl -s -C - --retry 5 -o "reads/$name.fastq.gz" "https://$url"
echo "$md5 reads/$name.fastq.gz"
done < runs.tsv > md5.txt
md5sum -c md5.txt
# The reference: the plastome of cultivar Désirée, as fasta and as the annotated GenBank record
curl -s "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&id=NC_008096.2&rettype=fasta" \
> NC_008096.2.fasta
curl -s "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&id=NC_008096.2&rettype=gbwithparts&retmode=text" \
> NC_008096.2.gbIf an MD5 check fails, run the loop again: curl -C - resumes interrupted downloads.
bacon -r NC_008096.2.fasta -i reads/ -o bacon_potato -t 32 -p 8 --annotation NC_008096.2.gbThe defaults: baiting with minimap2, Filtlong capping each sample at 100x, templated assembly with samtools,
core SNPs with SKA2, and a FastTree tree. It takes about 3 minutes. --annotation gives the report the genes
and regions of the plastome and the effect of each SNP; it changes nothing else, so it can be added to a
finished run (nothing is redone). -r NC_008096.2.gb alone does the same: the sequence of the GenBank record is
the reference and its features the annotation, and reference.fasta is identical to the one made from the
fasta.
Open bacon_potato/report.html in a browser for everything below on one page. The report of this run is
online:
view it
or download it
(local paths removed). The next sections look at the files behind it.
Sample Status Raw_reads Raw_bases Baited_reads Baited_pct Filtered_reads Est_depth Assembly_length N_bases
12_22_134 ok 18388 202618009 3553 17.437 1433 100.0 155386 133
14_4_1 ok 15023 182354649 3018 18.349 1313 100.1 155383 153
14_6_3 ok 20396 234625686 4122 17.970 1292 100.2 155166 0
...
- 10–23% of the bases are plastid reads, typical of leaf DNA; three runs (Alaska, Argo, Shah) were already filtered to plastid reads (nearly 100%).
- Every sample reaches 64–100x after filtering, far more than needed.
- The consensus lengths fall into two groups: about 155,170 bp and about 155,390 bp. Samples of the second group
have 119–153
Nbases, nearly all in one region, next to position 52,580.
The SNP alignment (ska.snps.fasta) has 118 SNP sites between the 28 plastomes and the reference. The cultivars fall into
three groups:
| Group | Cultivars | SNPs to the reference | Within the group |
|---|---|---|---|
| T-type | 14_6_3, 16-35-5, 16_1_2, Bagira, Bankir, Iskra, Luks, Shah, Terra, Zdraven | 1 | identical |
| Lineage A | 15-27-1, Legenda | 67 | identical |
| Lineage B | 12_22_134, 14_4_1, 15_22_4, 16_4_3, Alaska, Amur, Argo, Baron, Bravo, Gornyak, Irbitskiy, Kamenskiy, Mishka, Otrada, Start, Utro_ranneye | 67–76 | 0–13 |
Lineage A is 66 SNPs from the T-type group, lineage B 66–75, and the two lineages are 76–85 SNPs apart. Lineage B holds three distinct plastomes: 11 identical cultivars, a group of four (14_4_1, 16_4_3, Baron, Start) 5 SNPs away, and 15_22_4, 12–13 SNPs from both.
snps.vcf gives the genotype of each cultivar at 135 SNP positions of the reference, more than the 118 SNP sites
of the alignment: the alignment keeps the SNPs found in every genome and counts a SNP of the inverted repeat once,
while the VCF lists every SNP position of the reference (Outputs).
With the annotation, the report's genome map shows the LSC, IRb, SSC and IRa regions and the genes, and the SNP table under it says where each SNP falls: 91 of the 135 are in the LSC, 31 in the SSC and 13 in the inverted repeats; 57 are in coding sequences (30 synonymous, 26 missense, and one that changes the stop codon of ndhF into another), 14 in introns, 2 in rRNA genes, 3 in pseudogenes, 1 in sprA (a small plastid RNA, annotated as a gene only) and 59 are intergenic (one SNP, where the end of ndhF overlaps the ycf1 pseudogene, is counted for both). ndhA and ycf1 carry 7 SNPs each, ndhF 6. For example, the SNP at position 2,669 (C>T, carried by the 16 cultivars of lineage B) changes codon 333 of matK from GAC to AAC (D333N), and the one at 883 (T>C, in 15-27-1; Legenda has no call there) codon 243 of psbA from GAA to GGA (E243G).
The T-type group shares the reference's plastome but for one SNP. Lineages A and B are not T-type (next section); which of the other cytoplasm types they are is not determined here.
Sample metadata colours the report by whatever is known about the samples (Usage). As an
illustration, the published report uses the three lineages above as metadata, in a TSV given with
--metadata lineages.tsv to the same command (nothing is redone; only the report is rebuilt):
# Lineages derived in the tutorial (section 4) from the SNP distances: an illustration of --metadata
sample lineage
14_6_3 T-type
15-27-1 lineage A
12_22_134 lineage B
...
The tree then shows the lineage after each cultivar's name with its marker (a colour and a shape for each lineage), the heatmap gets a band of lineage markers (the groups of identical genomes turn grey, so that colour means the lineage only), and a table under the identical genomes counts the cultivars of each lineage in each set of identical plastomes (lineage A is one set of 2; lineage B three sets of 11, 4 and 1; the T-type cultivars one set of 10). Real metadata (origin, breeding programme, year) is used the same way.
The N bases of lineages A and B sit in the ndhC–trnV-UAC spacer (reference positions 51,834 to
52,696). Aligning an assembly of each group to the reference shows why:
minimap2 -cx asm5 --cs NC_008096.2.fasta bacon_potato/3_assembled/all_assemblies/Alaska.fastathe plastomes of lineages A and B have an insertion of about 240 bp at position 52,578 relative to the T-type reference:
the known T-type deletion, seen from the other side. The templated assembly places it, but the reads do not
agree on all of its bases, hence the N. A de novo assembly resolves it completely:
bacon -r NC_008096.2.fasta -i reads/ -o bacon_potato -a flye -t 48 -p 12reuses the baited and filtered reads and assembles each sample with Flye (about 10 minutes). The Flye
assembly of Alaska has an insertion of exactly 241 bp at position 52,578. Flye also reports 17 of the 28
plastomes as one circular contig. Its SNP distances (4_compared/ska/snp_distances.tsv, replaced) are
identical to those of the templated assembly for all 406 pairs of genomes (the 28 cultivars and the
reference).
- Genome skimming data from a few hundred megabases per sample is enough for complete plastomes.
- The templated assembly gives the SNPs and small indels in minutes; a de novo assembly resolves the insertions and the structure, and confirms the SNPs independently.
- Identical plastomes (distance 0) are common: these 28 cultivars carry only five distinct plastomes.
- Kawagoe Y., Kikuta Y. (1991) Chloroplast DNA evolution in potato (Solanum tuberosum L.). Theoretical and Applied Genetics 81:13–20. https://doi.org/10.1007/BF00226106
- Hosaka K. (2002) Distribution of the 241 bp deletion of chloroplast DNA in wild potato species. American Journal of Potato Research 79:119–123. https://doi.org/10.1007/BF02881520