---
title: "Reference package construction from scratch"
output:
  prettydoc::html_pretty:
    theme: architect
    highlight: github
vignette: >
  %\VignetteIndexEntry{Reference package construction from scratch}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}  
---

```{r style, echo = FALSE, results = 'asis'}
```

Please report comments or bugs to Pierre Lefeuvre - <pierre.lefeuvre@cirad.fr>

[BoSSA CRAN page](https://cran.r-project.org/package=BoSSA)

# 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](https://matsen.fhcrc.org/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](https://github.com/fhcrc/taxtastic), [rppr](http://matsen.github.io/pplacer/generated_rst/rppr.html) 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](https://talk.ictvonline.org/ictv-reports/ictv_9th_report/positive-sense-rna-viruses-2011/w/posrna_viruses/265/luteoviridae) 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".

### Loading the required packages

We will use the [rentrez package](https://cran.r-project.org/package=rentrez) to query NCBI from R. Along with a vignette available at the [rentrez CRAN webpage](https://cran.r-project.org/package=rentrez), examples of rentrez use are available [here](https://docs.ropensci.org/rentrez/).

```{r }
library("rentrez")
library("httr")
library("XML")
library("ape")
library("BoSSA")
set_config(config(http_version = 0))
```

### How many polerovirus sequences are available on NCBI ?

We first search for the taxonomic id associated to the *Polerovirus* genus in NCBI...

```{r, eval=F}
r_search <- entrez_search(db="taxonomy", term="polerovirus")
r_search$ids
```

... and use this taxonomic id (i.e. 119164) to query the nucleotide database. Note that we remove the reference sequence set from the search (using "NOT srcdb_refseq[PROP]") as these are duplicates from records already available in the nucleotide database.

```{r, eval=F}
r_search <- entrez_search(db="nucleotide", term="txid119164[orgn] NOT srcdb_refseq[PROP]")
r_search$count
```

This number represents the total number of sequences from the genus polerovirus. We will now search for full genomes using the "AND 5000:6500[slen]" query.

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

To recover all the GenBank identifiers (gi), we have to increase the retmax parameter (default is 20) to something superior to r_search$count.

```{r, eval=F}
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]
```

### Download all these sequences

#### Get the genbank informations in XML format

We now use these gi to download the information relative to the sequences in XML format. XML is a handy format as it allow to directly extract specific fields form the genbank file. Because the number of sequence may be a bit large, we will use the history option in our rentrez query. Sequences will then be obtained by batch of 200.

```{r, eval=F}
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))
```
#### Extract the fasta from the XML

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

#### Extract the taxonomy from the XML

Also, we need to extract the polerovirus classification of each sequence.

```{r, eval=F}
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.

```{r }
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).

```{r }
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.

```{r}
write.tree(root_tree,"polerovirus_ROOTED.tre")
```

# The Taxonomy

### Download the polerovirus taxonomy IDs

We need to obtain the complete taxonomy information for the polerovirus. Taxastic has a function for that. We first download the full taxonomy from NCBI using taxit and subset it for the poleroviruses tax ids.

```{r, eval=F}
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
```

### The information file

We can now create and write the information file for the sequences from our phylogenetic tree.

```{r, eval=F}
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
```

### Check the refpkg using rppr

Rppr offers the possibility to check the reference package using the following commands.

```
### in bash

rppr check -c polerovirus.refpkg

rppr info -c polerovirus.refpk
```

### Refpkg stats using BoSSA

Using BoSSA, we can extract some summary statistics and draw some plots.

- A reference package summary 

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

refpkg(refpkg_path)
```

- A tree with tips colored according to a taxonomic level, here the species

```{r tree1, fig.width=5, fig.height=5}
refpkg(refpkg_path,type="tree",cex.text=0.3,rank_tree="species")
```

- A pie chart summarizing the taxonomy with here again, the species level. Note there is a slight decay between the text labels and slices...

```{r pie1, fig.width=5, fig.height=5}
refpkg(refpkg_path,type="pie",rank_pie="species",cex.text=0.6)
```

- Alternatively, a input file for KronaTools to generate a krona chart could be generated.

```{r}
refpkg(refpkg_path,type="krona")
```

A file (default is for_krona.txt) is generated and can be used to generate a krona chart with the ImportText.pl scrit from the [Krona software](https://github.com/marbl/Krona/wiki/KronaTools):

```
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](http://tcoffee.readthedocs.io/en/latest/tcoffee_main_documentation.html#extracting-removing-sequences-with-the-identity) could be usefull if you would like to subsample based on diversity.

```{r, eval=F}
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).

```{r, eval=F}
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:

```{r}
citation("BoSSA")
```

# Other resources

### On phylogenetic placements
[pplacer website](https://matsen.fhcrc.org/pplacer/) and [documentation](http://matsen.github.io/pplacer/)

[taxtastic](https://github.com/fhcrc/taxtastic)

### Other R package with a related topic
[ggtree](https://bioconductor.org/packages/release/bioc/html/ggtree.html) and [clstutils](https://bioconductor.org/packages/release/bioc/html/clstutils.html)

### On krona charts
[krona](https://github.com/marbl/Krona/wiki/KronaTools)

### R packages used by BoSSA
[RSQLite](https://cran.r-project.org/package=RSQLite) and [jsonlite](https://cran.r-project.org/package=jsonlite) to read files, [ape](https://cran.r-project.org/package=ape) and [phangorn](https://cran.r-project.org/package=phangorn) to manipulate phylogenetic trees and [plotrix](https://cran.r-project.org/package=plotrix) for pie charts.