Important notice

Note that most of the R code included in the vignette is not executed during vignette build. Some of the rentrez (very good) package code triggers HTTP failure during the CRAN check and is thus not evaluated during the build. It should run smoothly on you desktop computer.

Summary

A phylogenetic placement corresponds to the position of a query sequence in a reference tree. Phylogenetic placement inference using pplacer requires the construction of “reference packages”, i.e. informations on the phylogeny and taxonomy of a given group of organisms. This vignette presents example of reference package construction using taxtastic, rppr and R.

The polerovirus example

Imagine you obtained some viral contigs that were assigned to the Polerovirus genus using a BLAST search. You want to obtain phylogenetic placements of these sequences using pplacer and now have to build a reference package for this genus.

A search on the ICTV website informs us that polerovirus belongs to the Luteovidiae family and that their “genome size is fairly uniform ranging from 5.6 kb to 6.0 kb”.

library("rentrez")
library("httr")
library("XML")
library("ape")
library("BoSSA")
set_config(config(http_version = 0))

r_search <- entrez_search(db="taxonomy", term="polerovirus")
r_search$ids

r_search <- entrez_search(db="nucleotide", term="txid119164[orgn] NOT srcdb_refseq[PROP]")
r_search$count

r_search <- entrez_search(db="nucleotide", term="txid119164[orgn] NOT srcdb_refseq[PROP] AND 5000:6500[slen]")
r_search$count

r_search <- entrez_search(db="nucleotide", term="txid119164[orgn] NOT srcdb_refseq[PROP] AND 5000:6500[slen]",retmax=1000)
gi <- r_search$ids
gi[1:10]

r_search <- entrez_search(db="nucleotide", term="txid119164[orgn] NOT srcdb_refseq[PROP] AND 5000:6500[slen]",use_history=TRUE)
all_recs <- NULL
for(seq_start in seq(1,ceiling(r_search$count/200)*200,200)){
     all_recs <- c(all_recs,entrez_fetch(db="nuccore", web_history=r_search$web_history,rettype="gbc",retmode="xml",retmax=200,retstart=seq_start))
}
rec_list <- do.call(c,lapply(all_recs,xmlToList))

write(sapply(rec_list,function(X){paste(">",X$INSDSeq_locus,"\n",X$INSDSeq_sequence,sep="")}),"polerovirus_from_genbank.fasta")

source_name <- t(sapply(rec_list,function(X){c(X$INSDSeq_locus,X$INSDSeq_organism)}))

Alignement and phylogeny

Now we have the sequences on the disk we can align the sequences (here using MAFFT) and construct a phylogenetic tree (here using FastTree)…

### in bash

mafft --adjustdirection --reorder polerovirus_from_genbank.fasta > polerovirus_from_genbank_MAFFT.fasta

FastTree -nt -gamma -gtr -log polerovirus_from_genbank_MAFFT.log polerovirus_from_genbank_MAFFT.fasta > polerovirus_from_genbank_MAFFT.tre

…and import the tree in R for control.

tree <- read.tree(paste(find.package("BoSSA"),"/extdata/polerovirus_from_genbank_MAFFT.tre",sep=""))
plot(tree,cex=0.4,no.margin=TRUE)
nodelabels(cex=0.8)

Note that using the “adjustdirection” during alignment, MAFFT generates reverse complement sequences and align them together with the remaining sequences. The best orientation is conserved for each sequence. When the reverse complement is aligned, the program add the “_R_” prefix in the name of the sequence. In this case and in order to keep the correspondance between informations files and sequences names, you’ll have to edit the sequences names.

After plotting the tree, and the node labels, it is apparent that the node 231 may be a good rooting point (may change depending on the sequence set you download i.e. as more poleroviruses sequences are upload to NCBI, the tree structure and node numbering will change).

root_tree <- root(tree,node=231,resolve=TRUE)
plot(root_tree,cex=0.4,no.margin=TRUE)

The rooted tree is then written to the disk.

write.tree(root_tree,"polerovirus_ROOTED.tre")

The Taxonomy

r_search <- entrez_search(db="taxonomy", term="txid119164[orgn]",retmax=1000)
tax_ids <-  r_search$ids
write(tax_ids,"taxonomy_polerovirus.id")

### in bash

taxit new_database

taxit taxtable -f taxonomy_polerovirus.id -o polerovirus_taxonomy.csv ncbi_taxonomy.db

taxo <- read.csv(paste(find.package("BoSSA"),"/extdata/polerovirus_taxonomy.csv",sep=""))
tax_id <- taxo$tax_id[match(source_name[,2],taxo$tax_name)]
info <- data.frame(seqname=source_name[,1],accession=source_name[,1],tax_id=tax_id,species=source_name[,2],is_type=rep("no",nrow(source_name)))
write.table(info,"polerovirus_info.csv",sep=",",row.names=FALSE)

Create the refpkg using taxit

The reference package is created using the taxit create command.

### in bash

taxit create -l polerovirus -P polerovirus.refpkg \
    --taxonomy polerovirus_taxonomy.csv \
    --aln-fasta polerovirus_from_genbank_MAFFT.fasta \
    --seq-info polerovirus_info.csv \
    --tree-stats polerovirus_from_genbank_MAFFT.log \
    --tree-file polerovirus_ROOTED.tre --no-reroot

### in bash

rppr check -c polerovirus.refpkg

rppr info -c polerovirus.refpk

refpkg_path <- paste(find.package("BoSSA"),"/extdata/polerovirus.refpkg",sep="")

refpkg(refpkg_path)

## ### Reference package summary
## 
## Path:/tmp/RtmpcxpYE8/Rinstaf84f3f26d2/BoSSA/extdata/polerovirus.refpkg
## 
## Tree with 204 tips 203 nodes
## 
## Classification:
## root 1 
## superkingdom 1 
## below_superkingdom 1 
## below_below_superkingdom 1 
## family 1 
## genus 1 
## below_genus 1 
## species 33 
## below_species 2

refpkg(refpkg_path,type="tree",cex.text=0.3,rank_tree="species")

refpkg(refpkg_path,type="pie",rank_pie="species",cex.text=0.6)

refpkg(refpkg_path,type="krona")

ImportText.pl for_krona.txt

Subsample the refpkg

It could be interesting to reduce the size of the reference package. Whereas, there is still a decent number of sequences in this polerovirus dataset, this number could be way higher and represent a computationnal burden for the rest of the analysis. Here is some code example to subset each species to a maximum of 10 sequences. T-Coffee could be usefull if you would like to subsample based on diversity.

d <- cophenetic.phylo(root_tree)

ids <- NULL
for(i in 1:length(unique(info$species))){
    usp <- unique(info$species)[i]
  sub <- info[as.character(info$species)==usp,]
  acc <- as.character(sub[,1])
  
  # the following code line prevent the code from crashing
  # i.e. new sequences not available in the example may be uploaded
  # when you will run the code
  acc <- acc[!is.na(match(acc,colnames(d)))]

  if(length(acc)<=10){
    ids <- c(ids,acc)
  }
  if(length(acc)>10){
    h <- hclust(as.dist(d[acc,acc]))
    grp <- cutree(h,k=10)
    ids <- c(ids,names(grp)[match(1:10,grp)])
  }
}

to_remove <- root_tree$tip.label[!root_tree$tip.label%in%ids]

Note that the highly recommanded practice is to subset the alignement, re-align and compute a new phylogenetic tree, one can also directly subset the already available phylogeny (in a “quick and dirty” way).

root_tree2 <- drop.tip(root_tree,to_remove)
write.tree(root_tree2,"polerovirus_ROOTED_SUBSAMPLED.tre")

Using these new files, you can create another smaller refpkg using the “taxit create” command as described above.

Citation

If you find BoSSA and/or its tutorials useful, you may cite:

citation("BoSSA")

## To cite package 'BoSSA' in publications use:
## 
##   Lefeuvre P (2020). _BoSSA: A Bunch of Structure and Sequence
##   Analysis_. R package version 3.7.
## 
## A BibTeX entry for LaTeX users is
## 
##   @Manual{,
##     title = {BoSSA: A Bunch of Structure and Sequence Analysis},
##     author = {Pierre Lefeuvre},
##     year = {2020},
##     note = {R package version 3.7},
##   }
## 
## ATTENTION: This citation information has been auto-generated from the
## package DESCRIPTION file and may need manual editing, see
## 'help("citation")'.

Other resources