inst/doc/Measurable-Evolution.R

In dowser: B Cell Receptor Phylogenetics Toolkit

## ----eval=FALSE, warning=FALSE, message=FALSE---------------------------------
# library(dowser)
# 
# # load example AIRR tsv data
# data(ExampleAirr)
# 
# # Process example data using default settings
# clones = formatClones(ExampleAirr, traits="timepoint", minseq=3)
# 
# # Column shows timepoints in dataset
# print(table(ExampleAirr$timepoint))
# #0   7  14
# #62 102 225
# 
# # Calculate number of tissues sampled in tree
# timepoints = unlist(lapply(clones$data, function(x)
#   length(unique(x@data$timepoint))))
# 
# # Filter to multi-type trees
# clones = clones[timepoints > 1,]
# 
# # Build trees using maximum likelihood (can use alternative builds if desired)
# trees = getTrees(clones, build="pml")
# 

## ----eval=FALSE, warning=FALSE, message=FALSE---------------------------------
# 
# # correlation test with 10000 repetitions
# test = correlationTest(trees, permutations=10000, time="timepoint")
# print(test)
# 
# # A tibble: 6 × 12
# #  clone_id data       locus  seqs trees      slope     p corre…¹ random…²  min_p
# #     <dbl> <list>     <chr> <int> <list>     <dbl> <dbl>   <dbl>    <dbl>  <dbl>
# #1     3128 <airrClon> N        40 <phylo> -0.00205 0.859  -0.173 -0.0257  0.0667
# #2     3184 <airrClon> N        12 <phylo>  0.00111 0.497   0.649 -0.00429 0.5
# #3     3140 <airrClon> N         9 <phylo>  0.00156 0.335   0.630 -0.00835 0.333
# #4     3192 <airrClon> N         9 <phylo>  0.00739 0.498   0.956 -0.00306 0.5
# #5     3115 <airrClon> N         6 <phylo>  0.00159 0.244   0.565  0.00236 0.25
# #6     3139 <airrClon> N         6 <phylo>  0.00308 0.507   0.821  0.0112  0.5
# 
# 
# # use uniform correlaion test (more sensitive, but higher false positive rate)
# utest = correlationTest(trees, permutations=10000, time="timepoint", perm_type="uniform")
# print(utest)
# 
# # A tibble: 6 × 12
# #  clone_id data       locus  seqs trees      slope      p correlation random_c…¹
# #     <dbl> <list>     <chr> <int> <list>     <dbl>  <dbl>       <dbl>      <dbl>
# #1     3128 <airrClon> N        40 <phylo> -0.00205 0.856       -0.173   0.00146
# #2     3184 <airrClon> N        12 <phylo>  0.00111 0.0768       0.649  -0.00223
# #3     3140 <airrClon> N         9 <phylo>  0.00156 0.114        0.630   0.00205
# #4     3192 <airrClon> N         9 <phylo>  0.00739 0.110        0.956   0.000223
# #5     3115 <airrClon> N         6 <phylo>  0.00159 0.336        0.565   0.00409
# #6     3139 <airrClon> N         6 <phylo>  0.00308 0.172        0.821   0.00431
# 
# 

## ----eval=FALSE, warning=FALSE, message=FALSE---------------------------------
# library(ggtree)
# 
# # order trees by p value
# test = test[order(test$p),]
# 
# # Plot times on tree with lowest p value (not convincingly evolving)
# plotTrees(test)[[1]] +
#     geom_tippoint(aes(fill=timepoint), pch=21, size=2) +
#     scale_fill_distiller(palette="RdYlBu")

## ----eval=TRUE, echo=FALSE, warning=FALSE, message=FALSE----------------------
library(dowser)
data(TimeTrees)
TimeTrees = TimeTrees[order(TimeTrees$p),]
plotTrees(TimeTrees)[[1]] + 
    ggtree::geom_tippoint(aes(fill=timepoint), pch=21, size=2) +
    scale_fill_distiller(palette="RdYlBu")