---
title: "Plotting alongside phylogeny"
author: "Fernando Roa"
date: "02 10 2019"
pkgdown:
  as_is: true
output:
  prettydoc::html_pretty:
    theme: leonids
    highlight: github
    includes:
      in_header: myheader.html
    toc: true
    toc_depth: 2
    number_sections: true
    pandoc_args: ["--lua-filter=refs/multiple-bibliographies.lua", "--number-offset=2"]
    
#bibliography: [refs/allrefs.bib,refs/iqtree.bib,refs/revbayes.bib,refs/packages.bib, refs/idiogramFISH.bib, refs/Rrefs.bib]

bibliography_normal: [refs/allrefs.bib,refs/iqtree.bib,refs/revbayes.bib]
bibliography_software: [refs/packages.bib, refs/idiogramFISH.bib, refs/Rrefs.bib]
# bibliography_docs: [refs/docPackages.bib]

csl: refs/clima5.csl
link-citations: yes
resource_files:
  - refs/multiple-bibliographies.lua
  - refs/allrefs.bib
  - refs/iqtree.bib
  - refs/revbayes.bib
  - refs/packages.bib
  - refs/idiogramFISH.bib
  - refs/Rrefs.bib
  - refs/docPackages.bib
  - refs/clima5.csl
  - firstplot2.png
  - secondplot2.png
  - thirdplot2.png
  - firstplot.png
  - secondplot.png
  - thirdplot.png
  
# params:
  # packageCheck: FALSE
vignette: >
  %\VignetteIndexEntry{3. Plotting alongside phylogeny}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---
<!--<a name="top"></a> -->
```{r setup, include=FALSE}
#Create myheader.html
if(Sys.info()['sysname']=="Windows"){
res<-!as.logical(system(paste("ping", "www.google.com")) )
  if(res){
  fileConn <- file("myheader.html")
  writeLines('<script src="https://kit.fontawesome.com/af0a13599b.js" crossorigin="anonymous"></script>', fileConn)
  close(fileConn)
  }
} else {
  fileConn <- file("myheader.html")
  writeLines('<script src="https://kit.fontawesome.com/af0a13599b.js" crossorigin="anonymous"></script>', fileConn)
  close(fileConn)
}
```

``` {r, results="asis", echo=FALSE, message=FALSE}
# <!-- pkgdown --> 
# <!-- jquery --><script src="js/jquery.min.js" crossorigin="anonymous"></script>
myfile<-"js/jquery.min.js"
if(file.exists(myfile)){
cat(paste0('<script src="',myfile,'" crossorigin="anonymous"></script> <!-- # -->'))
}
# <!-- clipboard.js --><script src="js/clipboard.min.js"  crossorigin="anonymous"></script>
myfile<-"js/clipboard.min.js"
if(file.exists(myfile)){
cat(paste0('<script src="',myfile,'"crossorigin="anonymous"></script>'))
}
# <!-- Font Awesome icons --><link rel="stylesheet" href="css/all.minMod.css"  crossorigin="anonymous">
myfile<-"css/all.minMod.css"
if(file.exists(myfile)){
cat(paste0('<link rel="stylesheet" href="',myfile,'"  crossorigin="anonymous">'))
}
# <!-- Bootstrap --><link rel="stylesheet" href="css/bootstrap.minO.css" crossorigin="anonymous">
myfile<-"css/bootstrap.minO.css"
if(file.exists(myfile)){
cat(paste0('<link rel="stylesheet" href="',myfile,'"  crossorigin="anonymous">'))
}
# <!-- # <script src="js/bootstrap.min.js"  crossorigin="anonymous"></script> -->
myfile<-"js/bootstrap.min.js"
if(file.exists(myfile)){
cat(paste0('<script src="',myfile,'" crossorigin="anonymous"></script> <!-- # -->'))
}
myfile<-"js/pkgdown2.js"
if(file.exists(myfile)){
cat(paste0('<script src="',myfile,'"></script> <!-- # -->'))
}
```
<link rel="stylesheet" href="css/mystyle.css">

<div class="goright" id="top">
Vignette [1. Plotting chromosomes](AplotIdiogramsVig.html)  
Vignette [2. Circular Plots](BcircularVig.html)  
Vignette [4. Human karyotype](EhumanVig.html)  
Vignette [5. Credits](index.html)  

<!-- # there are 2 spaces at the end of those lines -->
</div>

``` {r, echo=F, message=FALSE, fig.show = "hold", fig.align = "default", results="asis"}
if (requireNamespace("RCurl", quietly = TRUE)  ) {
# version of manual
v<-sub("Version: ","",readLines("../DESCRIPTION")[3])
# v<-tryCatch(suppressWarnings(rvcheck:::check_github_gitlab("ferroao/idiogramFISH", "gitlab")$latest_version), error=function(e) NA )
pkg<-"idiogramFISH"
link<-tryCatch(suppressWarnings(badger::badge_custom("Documentation", paste(pkg,v), "cornflowerblue") ), error=function(e) NA )
  if(!is.na(link)){
  svglink<-gsub("\\[|\\]|!|\\(|\\)","", link)
  manual_cont <- tryCatch(suppressWarnings(RCurl::getURLContent(svglink) ), error=function(e) NA )
    if (!is.na(manual_cont)){
    manual_contFile <- "../man/figures/manual.svg"
    writeLines(manual_cont, con = manual_contFile)
    manual_contFile <- normalizePath(manual_contFile)
    knitr::include_graphics(manual_contFile)
    # cat(paste0("&nbsp;![''](",knitr::include_graphics(manual_contFile),")" ) )
    }
  }
} # rcurl
```

# Plotting alongside phylogeny

This guide shows the files to plot idiograms alongside a phylogeny 

## Load package

visit gitlab for installation instructions https://gitlab.com/ferroao/idiogramFISH

```{r, results="hide", message=FALSE, warning=FALSE, eval=TRUE}

#load package
library(idiogramFISH) 
```

## ggtree of iqtree and monocentrics

idiogramFISH comes with two trees and data.frames with chr. and marks' data for the correspondent OTUs, first we will plot a tree produced with iqtree [@Nguyen2015c]

### Load the iqtree: {-}

We will use phytools for that [@phytools2012]

```{r , message=FALSE, echo=FALSE}
packageCheck<-all(unlist(invisible(lapply(c("ggtree","grid","ggpubr","phytools","treeio","plyr"), function(pkg) requireNamespace(pkg, quietly=TRUE)  ) ) ) )
```

``` {r , message=FALSE, eval=packageCheck}
require(ggplot2)
require(phytools)
require(ggpubr)
require(grid)   #pushViewport
require(ggtree)
# list.files(system.file('extdata', package = 'my_package') )

# find path of iqtree file
iqtreeFile    <- system.file("extdata", "eightSpIqtree.treefile", package = "idiogramFISH")

# load file as phylo object
iqtreephylo   <- read.newick(iqtreeFile) # phytools

# transform tree
iqtreephyloUM <- force.ultrametric(iqtreephylo, method= "extend") # phytools

```

### Make a ggtree [@ggtree2018] {-}

```{r, message=FALSE, eval=packageCheck}
ggtreeOf8 <- ggtree(iqtreephyloUM) + geom_tiplab(size=6)
```

Modify optionally graphical parameters with ggplot and ggpubr: [@ggplot22016; @R-ggpubr]

```{r, message=FALSE,eval=packageCheck}
gbuil2      <-  ggplot_build(ggtreeOf8)       # get ggplot_built
gtgbuild    <-  ggplot_gtable(gbuil2)         # get gtable from ggplot_built
gtgbuild$layout$clip[gtgbuild$layout$name == "panel"] <- "off"                # modify gtable
ggtreeOf8b   <- as_ggplot(gtgbuild)            # back to ggplot
gtgbuildgg2 <- ggtreeOf8b +  theme(plot.margin = unit(c(1,9.5,3,1.5), "cm") ) # top right bottom left - modify margins
```

Our plot of the tree is ready, let's get the order of species in tree

```{r, message=FALSE,eval=packageCheck}
ggtreeOf8TIPS<-ggtreeOf8$data[which(ggtreeOf8$data$isTip),]
desiredOrder  <- rev(ggtreeOf8TIPS[order(ggtreeOf8TIPS$y),]$label)
```
 
### Order OTUs of data.frame of chr. data {-}

Apply order of phylogeny to data.frame

```{r, message=FALSE,eval=packageCheck}
# make a vector without missing OTUs
desiredFiltered <- intersect(desiredOrder,allChrSizeSample$OTU)

# establish desired order
allChrSizeSample$OTU <- factor(allChrSizeSample$OTU, levels = desiredFiltered)

# order
allChrSizeSample     <- allChrSizeSample[order(allChrSizeSample$OTU),]
```

Now we have to establish where are the OTUs in the tree, that don't have chr. data

```{r, message=FALSE,eval=packageCheck}
# Establish position of OTUs before missing data OTUs
matchres <- match(desiredOrder,desiredFiltered)
matchres[is.na(matchres)] <- "R"
reps     <- rle(matchres)
posOTUsBeforeMissing      <- as.numeric(matchres[which(matchres=="R")-1][which(matchres[which(matchres=="R")-1]!="R")] )

# This are the OTUs that come before missing chr. data OTUs
BeforeMissing             <- desiredFiltered[posOTUsBeforeMissing]

# This is the amount of missing OTUs, spaces to add (ghost karyotypes)
valuesOfMissRepsBefore    <- reps$lengths[which(reps$values=="R")]
```

### Plotting {-}

Now we are ready to plot adding those arguments for `addMissingOTUAfter` and `missOTUspacings`

```{r, message=FALSE, comment=NA, results="hide",eval=packageCheck}
# plot to png file
png(file="firstplot.png" ,width=962,height=962 )

par(omi=rep(0,4) , mar=c(0,1,2,1), 
    mfrow=c(1,2) )   # one row two columns
par(fig=c(0,.3,0,1)) # location of left ghost plot
plot.new()           # ghost plot to the left
par(fig=c(.3,1,0,1)) # location of right plot

plotIdiograms(dfChrSize = allChrSizeSample,    # data.frame of Chr. Sizes
              dfMarkPos = allMarksSample,      # d.f. of Marks (inc. cen. marks) 
              dfMarkColor =  mydfMaColor,      # d.f. of mark characteristics
              
              roundness = 4,                   # roundness of vertices
              lwd.chr=.5,                      # width of lines
              orderBySize = FALSE,             # don't order chr. by size
              centromereSize = 1.3,            # apparent cen. size
              chrWidth =.75,                   # width of chr.
              chrSpacing = .25,                # horizontal spacing of chr.
              indexIdTextSize=.4,              # font size of indices and chr. names
              
              karHeight = 4.8,                 # karyotype vertical relative size without spacing
              karHeiSpace = 6.5,               # karyotype vertical relative size with spacing
              
              nameChrIndexPos=4,               # move the name of chr. indexes to left
              morpho="both",                   # add chr. morphology
              chrIndex = "both",               # add chr. indices
              karIndex = TRUE,                 # add karyotype indices
              yTitle= "",
              
              markLabelSpacer = 0              # spaces from rightmost chr. to legend
              
              ,ylimTopMod = -.1                # modify ylim top margin
              ,ylimBotMod=1.6                  # modify ylim bottom margin
              
              ,rulerPos = -0.5                 # position of rulers
              ,rulerNumberSize = .35           # font size of ruler number
              ,rulerNumberPos = .4             # position of ruler numbers
              ,ruler.tck = -.004               # tick size and orient.
              # ,asp=1                         # y x aspect
              
              ,addMissingOTUAfter = BeforeMissing          # OTUs after which there are ghost karyotypes - empty spaces
              ,missOTUspacings    = valuesOfMissRepsBefore # number of ghost karyotypes
)

# plot to the left the ggtree
pushViewport(viewport(layout = grid.layout(1, 2)))
pushViewport(viewport(layout.pos.col = 1, layout.pos.row = 1))
print(gtgbuildgg2,newpage=F) 

# close png
dev.off()
```

```{r, results="asis", comment=NA, echo=FALSE}
if(packageCheck){
# cat(paste0("![](",myt,")"))
cat(paste0("![](firstplot.png)" ) )
# cat(paste0("![](firstplot.svg)" ) )
} else {
cat(paste0("![](../man/figures/firstplot2.png)" ) )
# img1_path <- "../man/firstplot2.svg"
# if(file.exists(img1_path)) {
# cat(paste0("<img src=",img1_path," width=\"100%\">") )
# }
}
```
<div class="goright95">
<a href="#TOC">top <i class="fa fa-caret-square-o-up" aria-hidden="true"></i></a>
</div>
## plot of revBayes tree and holocentrics

Function `plotIdiogramsHolo` deprecated after ver. 1.5.1 

Now we are going to plot a tree from revBayes [@Hohna2017]

### First, load the revBayes tree: {-}

``` {r , message=FALSE, eval=packageCheck}
require(ggplot2)
require(phytools)
require(ggpubr)
require(grid)   #pushViewport
require(ggtree)
require(treeio)

# find path of iqtree file
revBayesFile    <- system.file("extdata", "revBayesTutorial.tree", package = "idiogramFISH")

# load file as phylo object
revBayesPhylo   <- read.beast(revBayesFile) # ggtree or treeio

# transform tree
revBayesPhyloUM <- force.ultrametric(revBayesPhylo@phylo, method= "extend") # phytools

```

Then, get order of OTUs in tree

```{r, message=FALSE,eval=packageCheck}
is_tip           <- revBayesPhyloUM$edge[,2] <= length(revBayesPhyloUM$tip.label)
ordered_tips     <- revBayesPhyloUM$edge[is_tip, 2]
desiredorderRevB <- rev(revBayesPhyloUM$tip.label[ordered_tips])
```
 
### Order OTUs of data.frame of chr. data {-}

First, create some data for holocentrics

```{r, message=FALSE,eval=packageCheck}
allChrSizeSampleHolo <- allChrSizeSample
allChrSizeSampleHolo <- allChrSizeSampleHolo[,c("OTU","chrName","longArmSize")]
colnames(allChrSizeSampleHolo)[which(names(allChrSizeSampleHolo)=="longArmSize")]<-"chrSize"

allMarksSampleHolo   <- allMarksSample
allMarksSampleHolo   <- allMarksSampleHolo[which(allMarksSampleHolo$chrRegion!="cen"),]
allMarksSampleHolo   <- allMarksSampleHolo[c("OTU","chrName","markName","markDistCen","markSize")]
colnames(allMarksSampleHolo)[which(names(allMarksSampleHolo)=="markDistCen")] <- "markPos"
allMarksSampleHolo[which(allMarksSampleHolo$markName=="5S"),]$markSize <- .5
```

Apply order of phylogeny to data.frame

```{r, message=FALSE,eval=packageCheck}
# make a vector without missing OTUs
desiredFiltered <- intersect(desiredorderRevB,allChrSizeSampleHolo$OTU)

# establish desired order
allChrSizeSampleHolo$OTU <- factor(allChrSizeSampleHolo$OTU, levels = desiredFiltered)

# order
allChrSizeSampleHolo <- allChrSizeSampleHolo[order(allChrSizeSampleHolo$OTU),]
```

Now we have to establish where are the OTUs in the tree, that don't have chr. data

```{r, message=FALSE,eval=packageCheck}
# Establish position of OTUs before missing data OTUs
matchres <- match(desiredorderRevB,desiredFiltered)
matchres[is.na(matchres)]   <- "R"
reps     <- rle(matchres)
posOTUsBeforeMissing        <- as.numeric(matchres[which(matchres=="R")-1][which(matchres[which(matchres=="R")-1]!="R")] )

# This are the OTUs that come before missing chr. data OTUs
BeforeMissingPlot2          <- desiredFiltered[posOTUsBeforeMissing]

# This is the amount of missing OTUs, spaces to add (ghost karyotypes)
valuesOfMissRepsBeforePlot2 <- reps$lengths[which(reps$values=="R")]
```

### Plotting {-}

Now we are ready to plot adding those arguments for `addMissingOTUAfter` and `missOTUspacings`

```{r, message=FALSE, comment=NA, results="hide",eval=packageCheck}
# plot to png file
png(file=paste0("secondplot.png" ),width=962,height=700)

{
par(omi=rep(0,4) , mar=c(0,0,0,0), mfrow=c(1,2))
par(fig=c(0,.27,0,1)) 
par(mar=c(2,0,2,0)) # b l t r

plot(revBayesPhyloUM)
par(fig=c(0.27,1,0,1), new=TRUE)
par(mar=c(0,0,0,0)) # b l t r

# Function plotIdiogramsHolo deprecated after ver. 1.5.1 

plotIdiograms(allChrSizeSampleHolo,               # chr. size data.frame
              dfMarkPos = allMarksSampleHolo,     # data.frame of marks' positions
              dfMarkColor =  mydfMaColor,         # d.f. of mark characteristics
              
              roundness = 4,                      # vertices roundness
              karHeight = 2.8,                    # karyotype height
              karHeiSpace = 4.5,                  # vertical size of kar. including spacing
              yTitle="",

              karIndex = TRUE,                    # add karyotype index
              indexIdTextSize=.4                  # font size of indices and chr. names
              
              ,addMissingOTUAfter = BeforeMissingPlot2           # add ghost OTUs after these names
              ,missOTUspacings    = valuesOfMissRepsBeforePlot2  # how many ghosts, respectively
              ,lwd.chr=.5                         # line width

              ,markLabelSpacer = 0                # dist. of legend to rightmost chr.
              ,legendWidth = 2.3                  # width of square or dots of legend
              
              ,rulerPos = -0.5                    # position of ruler
              ,rulerNumberSize = .35              # font size of number of ruler
              ,rulerNumberPos = .4                # position of ruler number
              ,ruler.tck=-.004                    # tick of ruler size and orient.

              ,ylimTopMod = -4                    # modify ylim of top
              ,ylimBotMod = -4                    # modify ylim of bottom
              ,xlimRightMod = 3                   # modify xlim right argument
              # ,asp=1                            # y x aspect
              )
}

# close png
dev.off()
```

```{r, results="asis", comment=NA, echo=FALSE}
if(packageCheck){
# myt2<-normalizePath(myt2)
# cat(paste0("![](",myt2,")"))
cat(paste0("![](secondplot.png)" ) )
} else {
cat(paste0("![](../man/figures/secondplot2.png)" ) )
# img1_path <- "../man/secondplot2.svg"
# if(file.exists(img1_path)) {
# cat(paste0("<img src=",img1_path," width=\"100%\">") )
# }
}
```
<div class="goright95">
<a href="#TOC">top <i class="fa fa-caret-square-o-up" aria-hidden="true"></i></a>
</div>
## plot of revBayes tree and holocentrics and monocentrics

Available for ver. > 1.5.1

Create data.frames with both types of karyotypes [@plyr2011]

```{r, message=FALSE,eval=packageCheck}
require(plyr)
# Select this OTU from the monocen.
monosel<-c("Species_F","Species_C","Species_A")
# chr.
allChrSizeSampleSel  <- allChrSizeSample [which(allChrSizeSample$OTU  %in% monosel  ),]
# marks
allMarksSampleSel    <- allMarksSample   [which(allMarksSample$OTU    %in% monosel  ),]

# Select the others from the holocen.
holosel    <- setdiff(unique(allChrSizeSampleHolo$OTU),monosel)
# chr.
allChrSizeSampleHoloSel <- allChrSizeSampleHolo[which(allChrSizeSampleHolo$OTU %in% holosel  ),]
# marks
allMarksSampleHoloSel   <- allMarksSampleHolo  [which(allMarksSampleHolo$OTU   %in% holosel  ),]

# merge chr d.fs
mixChrSize <- plyr::rbind.fill(allChrSizeSampleSel,allChrSizeSampleHoloSel)

# merge marks' d.fs
mixMarks   <- plyr::rbind.fill(allMarksSampleSel,allMarksSampleHoloSel)
```

Get the `desiredorderRevB` object from above to continue:

```{r, message=FALSE,eval=packageCheck}
# make a vector without missing OTUs
desiredFiltered <- intersect(desiredorderRevB, mixChrSize$OTU)

# establish desired order
mixChrSize$OTU <- factor(mixChrSize$OTU, levels = desiredFiltered)

# order data.frame
mixChrSize <- mixChrSize[order(mixChrSize$OTU),]

# Establish position of OTUs before missing data OTUs
matchres <- match(desiredorderRevB,desiredFiltered)
matchres[is.na(matchres)]   <- "R"
reps     <- rle(matchres)
posOTUsBeforeMissing        <- as.numeric(matchres[which(matchres=="R")-1][which(matchres[which(matchres=="R")-1]!="R")] )

# This are the OTUs that come before missing chr. data OTUs
BeforeMissingPlot2          <- desiredFiltered[posOTUsBeforeMissing]

# This is the amount of missing OTUs, spaces to add (ghost karyotypes)
valuesOfMissRepsBeforePlot2 <- reps$lengths[which(reps$values=="R")]
```

### Plotting {-}

Now we are ready to plot adding those arguments for `addMissingOTUAfter` and `missOTUspacings`

```{r, message=FALSE, comment=NA, results="hide",eval=packageCheck}
# plot to png file
png(file=paste0("thirdplot.png" ),width=1100,height=1000)
{
  par(omi=rep(0,4) , mar=c(0,0,0,0), mfrow=c(1,2))
  par(fig=c(0,.25,0,1)) 
  par(mar=c(1,0,0,0))
  
  plot(revBayesPhyloUM)
  par(fig=c(0.25,1,0,1), new=TRUE)
  par(mar=c(0,0,0,0))
  
plotIdiograms(mixChrSize,                         # chr. size data.frame
              dfMarkPos = mixMarks,               # data.frame of marks' positions (inc. cen. marks)
              dfMarkColor = mydfMaColor,          # d.f. of mark characteristics
              
              origin="b",                         # position measured from bottom of chr.
              
              karHeight = 2.8,                    # vertical size of kar. including spacing
              karHeiSpace = 4.5,                  # vertical size of kar. including spacing
              roundness = 5,                      # vertices roundness
              chrSpacing = .25,                   # horizontal spacing among chr.
              yTitle = "",
              
              karIndex = TRUE                     # add karyotype index
              ,indexIdTextSize=.4                 # font size of indices and chr. names
              ,distTextChr=0.7                    # dist. among chr. and chr. name
              
              ,addMissingOTUAfter = BeforeMissingPlot2           # add ghost OTUs after these names
              ,missOTUspacings    = valuesOfMissRepsBeforePlot2  # how many ghosts, respectively
              ,lwd.chr=.5                         # line width
              
              ,markLabelSpacer = 0                # dist. of legend to rightmost chr.
              ,legendWidth = 2                    # width of square or dots of legend
              
              ,ylimTopMod = -2                    # modify ylim of top
              ,ylimBotMod = -2                    # modify ylim of bottom
              
              ,rulerPos = -0.5                    # position of ruler
              ,rulerNumberSize = .35              # font size of number of ruler
              ,rulerNumberPos = .4                # position of ruler number
              ,ruler.tck=-.004                    # ruler tick size and orient.
              ,OTUfont=3                          # italics
              ,OTUfamily="Courier New"            # for OTU name
)
}
# close png
dev.off()
```

```{r, results="asis", comment=NA, echo=FALSE, eval=TRUE}
if(packageCheck){
cat(paste0("![](thirdplot.png)" ) )
} else {
cat(paste0("![](../man/figures/thirdplot2.png)" ) )
# img1_path <- "../man/logo.svg"
# if(file.exists(img1_path)) {
# cat(paste0("<img src=",img1_path," width=\"100%\">") )
# }
}
```
<div class="goright95">
<a href="#TOC">top <i class="fa fa-caret-square-o-up" aria-hidden="true"></i></a>
</div>
`r if (knitr:::is_html_output()) '
## References {-}
'`
::: {#refs_normal}
:::

### R-packages references {-}

::: {#refs_software}
:::

<!-- # Documentation references {-} -->

<!-- ::: {#refs_docs} -->
<!-- ::: -->



