Graph & annotations for two E. coli strains
This tutorial will go from linear assemblies and annotations of two E. coli bacterial genomes to an equivalent graph representation incorporating both strains at once. We'll also explore the resulting graph a bit to give you a taste of the vg toolkit. Let's get cracking!
Last update: July 2026. Run on vg v1.75.1 - Spike.
A pangenome reference can be made by smushing a bunch of linear references together. Thus, we start by acquiring some high-quality linear E. coli genomes.
Ensembl has some lovely strain genomes available for FTP.
mkdir ecoli
cd ecoli
# k12 strain
ensemblbac33fasta="ftp://ftp.ensemblgenomes.org/pub/bacteria/release-33/fasta/bacteria_9_collection"
wget "$ensemblbac33fasta/escherichia_coli_str_k_12_substr_dh10b/dna/Escherichia_coli_str_k_12_substr_dh10b.ASM1942v1.dna.chromosome.Chromosome.fa.gz"
# o157 strain
wget "$ensemblbac33fasta/escherichia_coli_o157_h7_str_sakai/dna/Escherichia_coli_o157_h7_str_sakai.ASM886v1.dna.chromosome.Chromosome.fa.gz"
# Some downstream tools do not like gzipped files
gunzip Escherichia_coli_str_k_12_substr_dh10b.ASM1942v1.dna.chromosome.Chromosome.fa.gz
gunzip Escherichia_coli_o157_h7_str_sakai.ASM886v1.dna.chromosome.Chromosome.fa.gzThe Ensebml FASTA is a bit odd. For example, Escherichia_coli_str_k_12_substr_dh10b.ASM1942v1.dna.chromosome.Chromosome.fa has the sequence identifier:
>Chromosome dna:chromosome chromosome:ASM1942v1:Chromosome:1:4686137:1 REF
But, strictly speaking, the sequence name is supposed to be everything between the > and the first space. Which means here it's "Chromosome". We first need to change the header to have the assembly ID be the first part of the header.
e.g.
>ASM1942v1 dna:chromosome chromosome:ASM1942v1:Chromosome:1:4686137:1 REF
We can do this with
sed -i "s/Chromosome/ASM1942v1/" Escherichia_coli_str_k_12_substr_dh10b.ASM1942v1.dna.chromosome.Chromosome.fa
sed -i "s/Chromosome/ASM886v1/" Escherichia_coli_o157_h7_str_sakai.ASM886v1.dna.chromosome.Chromosome.faFirst, we make a "stick graph" from the first E. coli genome: just the sequence turned into a bunch of nodes in a line. The o157/ASM886v1 strain has more base pairs (check with wc *.fa) so we'll use that as the base of our graph. We could make this graph "manually", but instead we'll take advantage of the automatic graph index construction pipeline. All we have to do is pass a FASTA:
vg autoindex \
--ref-fasta Escherichia_coli_o157_h7_str_sakai.ASM886v1.dna.chromosome.Chromosome.fa \
--workflow lr-giraffe \
--prefix Escherichia_coli_o157_h7_str_sakai.ASM886v1Now we will align the other E. coli genome to the stick graph using vg giraffe's long read mode. We try both long-read presets. Technically, neither is designed for the case where you misuse a long read aligner as an assembly aligner!
# Use parameter preset intended for PacBio HiFi reads
vg giraffe --parameter-preset hifi --threads 1 \
--gbz-name Escherichia_coli_o157_h7_str_sakai.ASM886v1.giraffe.gbz \
-f Escherichia_coli_str_k_12_substr_dh10b.ASM1942v1.dna.chromosome.Chromosome.fa \
> Escherichia_coli_str_k_12_substr_dh10b.ASM1942v1.hifi.gam
# Use parameter preset intended for R10 Nanopore reads
vg giraffe --parameter-preset r10 --threads 1 \
--gbz-name Escherichia_coli_o157_h7_str_sakai.ASM886v1.giraffe.gbz \
-f Escherichia_coli_str_k_12_substr_dh10b.ASM1942v1.dna.chromosome.Chromosome.fa \
> Escherichia_coli_str_k_12_substr_dh10b.ASM1942v1.r10.gamThen check the results:
Total alignments: 1
Total primary: 1
Total secondary: 0
Total aligned: 1
Total perfect: 0
Total gapless (softclips allowed): 0
Total paired: 0
Total properly paired: 0
Alignment score: mean 1.10037e+06, median 1.10037e+06, stdev 0, max 1.10037e+06 (1 reads)
Mapping quality: mean 60, median 60, stdev 0, max 60 (1 reads)
Insertions: 379495 bp in 25623 read events
Deletions: 1191808 bp in 67393 read events
Substitutions: 117199 bp in 108228 read events
Matches: 4189443 bp (4.18944e+06 bp/aligned)
Softclips: 0 bp (0% of bases, 0 bp/aligned) in 0 read events
Total time: 176.41 seconds
Speed: 0.0056686 reads/second
Total alignments: 1
Total primary: 1
Total secondary: 0
Total aligned: 1
Total perfect: 0
Total gapless (softclips allowed): 0
Total paired: 0
Total properly paired: 0
Alignment score: mean 1.00103e+06, median 1.00103e+06, stdev 0, max 1.00103e+06 (1 reads)
Mapping quality: mean 60, median 60, stdev 0, max 60 (1 reads)
Insertions: 391472 bp in 22905 read events
Deletions: 1203785 bp in 65956 read events
Substitutions: 111740 bp in 103601 read events
Matches: 4182925 bp (4.18292e+06 bp/aligned)
Softclips: 0 bp (0% of bases, 0 bp/aligned) in 0 read events
Total time: 159.884 seconds
Speed: 0.00625454 reads/second
Okay, these seem pretty similar. I'll go with the hifi one because it has slightly more matched bases and less variants, indicating maybe it managed to align some parts of the sequence better.
We want the "read" path to become a path in the graph. This is a job for vg augment --include-paths.
# GBZ is not editable, so first convert to an editable format
vg convert Escherichia_coli_o157_h7_str_sakai.ASM886v1.giraffe.gbz | \
vg augment --include-paths - Escherichia_coli_str_k_12_substr_dh10b.ASM1942v1.hifi.gam \
> two_strains.vgWe don't just want an alignment, though. We want to use this graph to understand genetic features in context. So, let's add some features!
Now we get GenBank annotations from NCBI for ASM1942v1 and ASM886v1,
- Select all genes (hit topmost checkbox in the table)
- Download with "Download" / "Download Table"
This will get TSVs. I saved as ASM1942v1.annot.tsv and ASM886v1.annot.tsv. We now need to convert to BED files; this is largely a matter of shifting around columns. Note that for later we're going to want all genes to have distinct names. Thus we suffix each gene with its location to provide a unique identifier.
# Chromosome, start, end, symbol, score (unimportant), strand
tail -n+2 ASM1942v1.annot.tsv | awk -v FS="\t" -v OFS="\t" '{print "ASM1942v1", $2, $3, $7 ":ASM1942v1:" $2 "-" $3, "0", $5}' > ASM1942v1.annot.bed
tail -n+2 ASM886v1.annot.tsv | awk -v FS="\t" -v OFS="\t" '{print "ASM886v1", $2, $3, $7 ":ASM886v1:" $2 "-" $3, "0", $5}' > ASM886v1.annot.bedWe use vg annotate to convert those BEDs into GAM files, indicating how the annotations align against the graph. We can then merge these annotations back into the original graph with another application of vg augment. That gives us a graph which has gene annotations within it.
# Generate BED vs graph annotations
vg annotate --xg-name two_strains.vg \
--bed-name ASM1942v1.annot.bed \
--bed-name ASM886v1.annot.bed \
> two_strains.annot.gam
# Merge annotations into the original graph
vg augment --include-paths two_strains.vg two_strains.annot.gam > two_strains.annot.vg(P.S. if you're noticing that -x/--xg-name can be passed a non-XG version graph, you are correct! Many such graph input options are flexible on input graph type and can take GBZ, VG, etc.)
We have this lovely graph now. What can we do with it?
For one, the graph is fundamentally a sequence alignment between the two strains. vg deconstruct will call variants based on paths in the graph. Since ASM886v1 was the backbone of our graph, we treat it as the reference (--path ASM886v1) and we will thus call variants relative to it.
vg deconstruct --path ASM886v1 --contig-only-ref two_strains.vg > two_strains.vcfThis will output a VCF treating every other path in the graph as a single-haplotype sample. We don't care about the path_cover_N paths; those are leftovers from augmenting. The variants for ASM1942v1 are those implicitly encoded within the graph structure. Now turned into a simple flat VCF for perusal. Here are the relevant columns from the first twenty lines (cut -f1-10 two_strains.vcf | head -20):
##fileformat=VCFv4.2
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##INFO=<ID=CONFLICT,Number=.,Type=String,Description="Sample names for which there are multiple paths in the graph with conflicting alleles">
##INFO=<ID=AC,Number=A,Type=Integer,Description="Total number of alternate alleles in called genotypes">
##INFO=<ID=AF,Number=A,Type=Float,Description="Estimated allele frequency in the range (0,1]">
##INFO=<ID=NS,Number=1,Type=Integer,Description="Number of samples with data">
##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">
##INFO=<ID=AT,Number=R,Type=String,Description="Allele Traversal as path in graph">
##contig=<ID=ASM886v1,length=5498450>
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT ASM1942v1
ASM886v1 58 >2>310698 C G 60 . AC=1;AF=0.0588235;AN=17;AT=>2>310697>310698,>2>460781>310698;NS=17 GT 1
ASM886v1 64 >310698>3 C T 60 . AC=1;AF=0.0588235;AN=17;AT=>310698>310699>3,>310698>460782>3;NS=17 GT 1
ASM886v1 225 >8>333602 CACCACCACCATCACCATT C 60 . AC=1;AF=0.0588235;AN=17;AT=>8>333601>333602,>8>333602;NS=17 GT 1
ASM886v1 327 >11>274837 T TT 60 . AC=1;AF=0.0588235;AN=17;AT=>11>274837,>11>460783>274837;NS=17 GT 1
ASM886v1 410 >13>270291 G T 60 . AC=1;AF=0.0588235;AN=17;AT=>13>270290>270291,>13>460784>270291;NS=17 GT 1
ASM886v1 576 >18>19 C T 60 . AC=1;AF=0.0588235;AN=17;AT=>18>286900>19,>18>460785>19;NS=17 GT 1
ASM886v1 605 >19>241330 A G 60 . AC=1;AF=0.0588235;AN=17;AT=>19>241329>241330,>19>460786>241330;NS=17 GT 1
ASM886v1 638 >20>335025 C T 60 . AC=1;AF=0.0588235;AN=17;AT=>20>335024>335025,>20>460787>335025;NS=17 GT 1
ASM886v1 686 >22>231812 A G 60 . AC=1;AF=0.0588235;AN=17;AT=>22>231811>231812,>22>460788>231812;NS=17 GT 1
ASM886v1 707 >23>327698 T C 60 . AC=1;AF=0.0588235;AN=17;AT=>23>327697>327698,>23>460789>327698;NS=17 GT 1
We have annotations in GAM format. Using vg filter --tsv-out, we can output information about those annotations. In particular, here we're going to output the nodes in the graph along which the annotations live.
vg filter --tsv-out "name;nodes" two_strains.annot.gam > two_strains.annot.nodes.tsvExporting to a plain-text format means we can run whatever bespoke analyses we like. For example, here's a brief Python script to find gene annotations with overlapping nodes:
# Set up a {node ID : {gene names}} dictionary
node_id_to_gene_name = dict()
with open('two_strains.annot.nodes.tsv') as file:
file.readline() # Skip header
for line in file:
parts = line.split('\t')
gene_name = parts[0]
# Sanitize node IDs from list like "476163+,476164+,"
nodes = [int(node_id.strip('+-')) for node_id in parts[1].split(',')[:-1]]
for node_id in nodes:
# Tie gene name to node ID
if not node_id in node_id_to_gene_name:
node_id_to_gene_name[node_id] = set()
node_id_to_gene_name[node_id].add(gene_name)
# Find genes with overlap
overlapping_genes = set()
for gene_names in node_id_to_gene_name.values():
if len(gene_names) > 1:
overlapping_genes.add(','.join(gene_names))
# Print out results
for overlap in overlapping_genes:
print(overlap)Put that code in a file called, say, find_overlaps.py, then run with python find_overlaps.py to see all the overlapping genes.
There is not always a specialized tool for the job you want to do. While the vg toolkit is vast (frankly, so vast that people frequently fail to find the function that does what they want, even if it exists), it is not infinite. In addition, many of the more specialized tools were developed for a particular person's use and may require preprocessing, data formats, etc. that are more of a headache than they're worth. If you have a simple question, sometimes the fastest way to answer it outside of vg's ecosystem.
Also see: vg view --gfa (plaintext graphs) & --json (plaintext alignments), vg paths --extract-fasta (plaintext sequences), vg pack --as-table (plaintext node depth), etc. There are plenty of ways to export raw graph-world data.
What if we want to check out the nodes around some of those genes with overlaps? Let's go for a short-looking pair in my output: tnaC:ASM1942v1:3984042-3984116 and ECs4644:ASM886v1:4677144-4677218. Visualize the paths within that chunk of graph:
# Grab the reference path & one node in any direction
vg chunk --path "ASM1942v1:3984042-3984116" --context-steps 1 \
--xg-name two_strains.annot.vg \
| vg view -dpn - | dot -Tsvg -o chunk.svg