DISCLAIMER: I’ve created this demo with an early version of gggenomes. Not all of the code will function with current releases. However, the general workflow both with respect to the external tools and gggenomes code is still valid and this demo should therefore be understood as an hopefully inspiring guide - not as an exactly reproducible code example (Those kind of examples you can find in the examples sections of the documentation).
This is a real-life example demonstrating the use of gggenomes
to explore viral genomes. We start with just a bunch of viral contigs and ask: Is there anything interesting going on here?
We use a few bioinformatics commandline tools to run some analyzes, and visualize the results using gggenomes
. By successively adding new data as new tracks we build a rich plot that ultimatly reveals novel insights into an exciting system.
This example is bundled data(package="gggenomes")
To rerun the this example including the bioinformatics analyzes download the raw data from github.
We start with a fasta file of 33 viral genomes. We read sequence length and some metadata from the header lines using `readfai`…
library(gggenomes)
# parse sequence length and some metadata from fasta file
emale_seqs <- read_fai("emales.fna") %>%
tidyr::extract(seq_desc, into = c("emale_type", "is_typespecies"), "=(\\S+) \\S+=(\\S+)",
remove=F, convert=T) %>%
dplyr::arrange(emale_type, length)
# plot the genomes - first six only to keep it simple for this example
emale_seqs_6 <- emale_seqs[1:6,]
p1 <- gggenomes(emale_seqs_6) +
geom_seq() + geom_bin_label()
p1
# https://github.com/thackl/seq-scripts
seq-join -n emales-concat < emales.fna > emales-concat.fna
# Annotate genes | https://github.com/hyattpd/Prodigal
prodigal -n -t emales-prodigal.train -i emales-concat.fna
prodigal -t emales-prodigal.train -i emales.fna -o emales-prodigal.gff -f gff
# A little help to clean up the prodigal gff | https://github.com/thackl/seq-scripts
gff-clean emales-prodigal.gff > emales.gff
emale_genes <- read_gff("emales.gff") %>%
dplyr::rename(feature_id=ID) %>% # we'll need this later
dplyr::mutate(gc_cont=as.numeric(gc_cont)) # per gene GC-content
p2 <- gggenomes(emale_seqs_6, emale_genes) +
geom_seq() + geom_bin_label() +
geom_gene(aes(fill=gc_cont)) +
scale_fill_distiller(palette="Spectral")
p2
It is known that these type of viruses often have linear genomes with terminal inverted repeats (TIRs). So let’s look for those next.
# split into one genome per file | https://bioinf.shenwei.me/seqkit/
seqkit split -i emales.fna
# self-align opposite strands
for fna in `ls emales.fna.split/*.fna`; do
minimap2 -c -B5 -O6 -E3 --rev-only $fna $fna > $fna.paf;
done;
cat emales.fna.split/*.paf > emales-tirs.paf
# prefilter hits by minimum length and maximum divergence
emale_tirs_paf <- read_paf("emales-tirs.paf") %>%
dplyr::filter(seq_id1 == seq_id2 & start1 < start2 & map_length > 99 & de < 0.1)
emale_tirs <- bind_rows(
dplyr::select(emale_tirs_paf, seq_id=seq_id1, start=start1, end=end1, de),
dplyr::select(emale_tirs_paf, seq_id=seq_id2, start=start2, end=end2, de))
p3 <- gggenomes(emale_seqs_6, emale_genes, emale_tirs) +
geom_seq() + geom_bin_label() +
geom_feature(size=5) +
geom_gene(aes(fill=gc_cont)) +
scale_fill_distiller(palette="Spectral")
p3