VG RDF, the Ensembl bacteria E. coli genome hack attack

Some of us would like to have an VG of the Ensembl E.coli bacterial genomes that so that we can see what mapped annotations from Ensembl&UniProt would look like on a VG datamodel.

Getting started: making a VG from a bunch of E.coli genomes.

mkdir ecoli
cd ecoli
#k12
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
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.gz

As the Ensebml fasta is a bit odd we first need to change the header to have the assembly id be the first part of the header. e.g.

>ASM1942v1 C1:Chromosome dna:chromosome chromosome:ASM1942v1:Chromosome:1:4686137:1 REF

vg construct \
  --reference ecoli/Escherichia_coli_str_k_12_substr_dh10b.ASM1942v1.dna.chromosome.Chromosome.fa \
  -m 1000 \
> ecoli/Escherichia_coli_str_k_12_substr_dh10b.ASM1942v1.vg

We are making a reference graph from the first E.coli genome and limit the node size to 1000.

Then its time to index the E.coli VG which we started with as root path in the vg graph.

vg index -x ecoli/Escherichia_coli_str_k_12_substr_dh10b.ASM1942v1.xg \
  -g ecoli/Escherichia_coli_str_k_12_substr_dh10b.ASM1942v1.gcsa \
  -k 11 \
  ecoli/Escherichia_coli_str_k_12_substr_dh10b.ASM1942v1.vg

Then its time to align other ecoli genomes to the first one

vg msga \
  -g ecoli/Escherichia_coli_str_k_12_substr_dh10b.ASM1942v1.vg \
  -f ecoli/Escherichia_coli_o157_h7_str_sakai.ASM886v1.dna.chromosome.Chromosome.fa \
> ecoli/k12_and_o157.vg

This gives you a vg file with two paths. Quick sanity check is to validate that we have nodes of varying length. e.g.

for i in $(./vg view -t ../ecoli/k12_and_o157.vg | grep r:value | cut -f 2 -d '"');
do 
  echo ${#i};
done|sort -nu

Getting Ensembl gene's into

Get BED file from Ensembl Genomes (click export button)

Then the BED file first column needs to correspond with your VG paths, you can check the paths in your vg file with this command

vg paths -L k12_and_o157.vg 

We then use the VG annotate command to make a new file in VG speak a "GAM file".

vg annotate -p -x k12_and_o157.xg -b Escherichia_coli_str_k_12_substr_dh10b.ASM1942v1.bed.fixed

This is merged back into a new vg file using the command vg mod

vg mod -i Escherichia_coli_str_k_12_substr_dh10b.ASM1942v1.gam k12_and_o157.vg > k12_and_o157_v2.vg

Preparing Ensembl and VG RDF

We get the Ensembl RDF from

ensembl33bac9="ftp://ftp.ensemblgenomes.org/pub/release-33/bacteria/rdf/bacteria_9_collection"
wget "$ensembl33bac9/escherichia_coli_str_k_12_substr_dh10b/escherichia_coli_str_k_12_substr_dh10b.ttl.gz"
wget "$ensembl33bac9/escherichia_coli_o157_h7_str_sakai/escherichia_coli_o157_h7_str_sakai.ttl.gz"

We turn the merged vg into rdf with vg view

./vg view -t ./ecoli/k12_and_o157.vg > ./ecoli/k12_and_o157.ttl

Compress it if you want even if its tiny

xz -T 4 ./ecoli/k12_and_o157.ttl

Assuming recent virtuoso set up loading assume data is in ${HOME}/ecoli it needs to be an absolute path.

isql $port $user $pass \ 
    exec="ld_dir('${HOME}/ecoli','k12_and_o157.ttl.xz', 'http://example.org/vg_ecoli')"
isql $port $user $pass \
    exec="ld_dir('${HOME}/ecoli','escherichia_coli_o157_h7_str_sakai.ttl.gz','http://example.org/ensembl_ecoli')"
isql $port $user $pass \
    exec="ld_dir('${HOME}/ecoli','escherichia_coli_str_k_12_substr_dh10b.ttl.gz','http://example.org/ensembl_ecoli')"
isql $port $user $pass \
    exec="rdf_loader_run(log_enable=>2);"

Checking the VG and Ensembl data

Sanity check the loaded VG data

SELECT
  (COUNT(?forward) as ?forward_nodes)
  (COUNT(?reverse) as ?reverse_nodes) 
WHERE {
  {?forward a <http://example.org/vg/Forward>} 
    UNION
  {?reverse a <http://example.org/vg/Reverse>}
}

SELECT
  DISTINCT ?path
WHERE {
  ?node <http://example.org/vg/path> ?path
}

linking up the faldo:reference set in Ensembl RDF and the vg:path

PREFIX vg:<http://example.org/vg/>
PREFIX faldo:<http://biohackathon.org/resource/faldo#>
SELECT 
  DISTINCT ?path  
WHERE {
  ?node vg:path|faldo:reference ?path 
}

Gives the chromosomes, plasmids and paths as in the db.

We are now going to insert new triples to state that the path in the vg graph are the same logical thing as the related linear references as known in the Ensembl RDF.

PREFIX path:<http://example.org/vg//path/>
INSERT {
  GRAPH<http://example.org/vg_ecoli_edits> 
  {
    path:ASM1942v1 <http://www.w3.org/2002/07/owl#sameAs>  <http://rdf.ebi.ac.uk/resource/ensembl/33/escherichia_coli_str_k_12_substr_dh10b/ASM1942v1/Chromosome> .
    path:ASM886v1 <http://www.w3.org/2002/07/owl#sameAs>  <http://rdf.ebi.ac.uk/resource/ensembl/33/escherichia_coli_o157_h7_str_sakai/ASM886v1/Chromosome>
  }
} 

Get the start and end positions of all protein coding genes on the Ecoli sakai chromosome that is a path in our VG graph.

PREFIX owl:<http://www.w3.org/2002/07/owl#> 
PREFIX faldo:<http://biohackathon.org/resource/faldo#> 
SELECT 
  ?s ?bp ?ep
WHERE {
 ?s a <http://purl.obolibrary.org/obo/SO_0001217>; 
   faldo:location ?l . 
 ?l faldo:end ?e ; faldo:begin ?b .
 ?e faldo:position ?ep ; faldo:reference ?chromo . 
 ?b faldo:position ?bp ; faldo:reference ?chromo . 
 <http://example.org/vg//path/ASM886v1> owl:sameAs ?chromo 
}

Fixing up IRIs in the VG output.