## ----eval=FALSE, warning=FALSE, message=FALSE---------------------------------
# library(dowser)
#
# # load example AIRR tsv data
# data(ExampleAirr)
#
# # trait value of interest
# trait="biopsy"
#
# # Process example data using default settings
# clones = formatClones(ExampleAirr,
# traits=trait, num_fields="duplicate_count", minseq=3)
#
# # Column shows which biopsy the B cell was obtained from
# print(table(ExampleAirr[[trait]]))
# #Lung Nose
# # 145 244
#
# # Calculate number of tissues sampled in tree
# tissue_types = unlist(lapply(clones$data, function(x)
# length(unique(x@data[[trait]]))))
#
# # Filter to multi-type trees
# clones = clones[tissue_types > 1,]
#
# # Build trees using maximum likelihood (can use alternative builds if desired)
# trees = getTrees(clones, build="pml")
#
## ----eval=FALSE, warning=FALSE, message=FALSE---------------------------------
# # the location of the igphyml executable
# # this is location in Docker image, will likely be different if you've set it up yourself
# # note this is the location of the compiled executable file, not just the source folder
# igphyml_location = "/usr/local/share/igphyml/src/igphyml"
#
# # build trees as before, but use IgPhyML to reconstruct the states of internal
# # nodes using maximum parsimony
# trees = getTrees(clones, build="pml", trait=trait, igphyml=igphyml_location)
#
# # show internal node (edge) predictions based on maximum parsimony
# plotTrees(trees, tips=trait, nodes=TRUE, palette="Set1")[[1]]
## ----eval=TRUE, echo=FALSE, warning=FALSE, message=FALSE----------------------
library(dowser)
# Load data instead of running phylip
data(BiopsyTrees)
plotTrees(BiopsyTrees, tips="biopsy", nodes=TRUE, palette="Set1")[[1]]
## ----eval=FALSE, warning=FALSE, message=FALSE---------------------------------
# # the location of the igphyml executable
# # this is location in Docker image, will likely be different if you've set it up yourself
# # note this is the location of the compiled executable file, not just the source folder
# igphyml_location = "/usr/local/share/igphyml/src/igphyml"
#
# # calculate switches along trees compared to 100 random permutations
# # this may take a while, and can be parallelized using nproc
# switches = findSwitches(trees, permutations=100, trait=trait,
# igphyml=igphyml_location, fixtrees=TRUE)
#
# # Perform PS test on switches
# ps = testPS(switches$switches)
# print(ps$means)
# # A tibble: 1 x 6
# # RECON PERMUTE PLT DELTA STAT REPS
# # <dbl> <dbl> <dbl> <dbl> <chr> <int>
# #1 8 8.6 0.4 -0.6 PS 100
#
#
# sp = testSP(switches$switches, alternative="greater")
# print(sp$means)
# # A tibble: 2 x 8
# # Groups: FROM [2]
# # FROM TO RECON PERMUTE PGT DELTA STAT REPS
# # <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr> <int>
# #1 Lung Nose 0.131 0.323 1 -0.192 SP 100
# #2 Nose Lung 0.869 0.677 0 0.192 SP 100
#
## ----eval=FALSE, warning=FALSE, message=FALSE---------------------------------
# # calculate switches along bootstrap distribution of trees
# # build using the 'pml' maximum likelihood option
# # in a real analysis it's important to use at least 100 permutations
# switches = findSwitches(trees, permutations=10, trait=trait,
# igphyml=igphyml_location, fixtrees=FALSE, build="pml")
#
# sp = testSP(switches$switches, alternative="greater")
# print(sp$means)
# # A tibble: 2 x 8
# # Groups: FROM [2]
# # FROM TO RECON PERMUTE PGT DELTA STAT REPS
# # <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr> <int>
# #1 Lung Nose 0.168 0.358 1 -0.190 SP 10
# #2 Nose Lung 0.832 0.642 0.1 0.190 SP 10
#
## ----eval=FALSE, warning=FALSE, message=FALSE---------------------------------
# sp = testSP(switches$switches, alternative="greater", permuteAll=TRUE)
# print(sp$means)
# # A tibble: 2 x 8
# # Groups: FROM [2]
# # FROM TO RECON PERMUTE PGT DELTA STAT REPS
# # <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr> <int>
# #1 Lung Nose 0.168 0.241 0.6 -0.0736 SP 10
# #2 Nose Lung 0.832 0.759 0.4 0.0736 SP 10
## ----eval=FALSE, warning=FALSE, message=FALSE---------------------------------
# # Downsample each tree to a tip-to-switch ratio of 10 instead of 20
# # this will reduce the false positive rate but also (likely) power
# switches = findSwitches(trees, permutations=100, trait=trait,
# igphyml=igphyml_location, fixtrees=TRUE, tip_switch=10)
#
# # didn't have much effect for this dataset
# sp = testSP(switches$switches, alternative="greater")
# print(sp$means)
# # A tibble: 2 x 8
# # Groups: FROM [2]
# # FROM TO RECON PERMUTE PGT DELTA STAT REPS
# # <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr> <int>
# #1 Lung Nose 0.168 0.358 1 -0.190 SP 10
# #2 Nose Lung 0.832 0.642 0.1 0.190 SP 10
#
## ----eval=FALSE, warning=FALSE, message=FALSE---------------------------------
# # the location of the igphyml executable
# # this is location in Docker image, will likely be different if you've set it up yourself
# # note this is the location of the compiled executable file, not just the source folder
# igphyml_location = "/usr/local/share/igphyml/src/igphyml"
#
# # constant region column name
# trait = "c_call"
#
# # Vector of human isotypes in the proper order. Isotype switching
# # can only go from left to right (e.g. IGHM to IGHA1). One exception
# # is IGHD, which can switch to IGHM.
# isotypes = c("IGHM","IGHD","IGHG3","IGHG1","IGHA1","IGHG2",
# "IGHG4","IGHE","IGHA2")
#
# # Process example data using default settings with "c_call" as a trait value
# clones = formatClones(ExampleAirr, traits=trait, minseq=3)
#
# # Column shows which constant region associated with a BCR
# print(table(ExampleAirr[[trait]]))
# # IGHA1 IGHA2 IGHD IGHG1 IGHG2 IGHG3 IGHG4 IGHM
# # 55 56 11 58 64 60 63 22
#
# # Calculate number of istoypes sampled in each tree
# isotype_counts = unlist(lapply(clones$data, function(x)
# length(unique(x@data[[trait]]))))
#
# # make model file with irreversibility constraints
# # Will prohibit switches from right to left in the "states" vector
# # IGHD to IGHM switching listed as an exception, since this can occur
# makeModelFile(file="isotype_model.txt", states=isotypes,
# constraints="irrev", exceptions=c("IGHD,IGHM"))
#
# # Build trees and predict states at internal nodes using maximum parsimony
# trees = getTrees(clones[isotype_counts > 1,], trait=trait, igphyml=igphyml_location,
# modelfile="isotype_model.txt")
#
# # show internal node (edge) predictions based on maximum parsimony
# plotTrees(trees, tips=trait, nodes=TRUE, palette="Paired", ambig="grey")[[1]]
## ----eval=TRUE, echo=FALSE, warning=FALSE, message=FALSE----------------------
data(IsotypeTrees)
plotTrees(IsotypeTrees, tips="c_call", nodes=TRUE, palette="Paired", ambig="grey")[[1]]
## ----eval=FALSE, warning=FALSE, message=FALSE---------------------------------
# # Downsample each tree to a tip-to-switch ratio of 10 instead of 20
# # this will reduce the false positive rate but also (likely) power
# switches = findSwitches(trees, permutations=100, trait=trait,
# igphyml=igphyml_location, fixtrees=TRUE, tip_switch=10,
# modelfile="isotype_model.txt")
#
# # didn't have much effect for this dataset
# sp = testSP(switches$switches, alternative="greater")
# print(sp$means,n=42)
# # A tibble: 42 × 8
# # Groups: FROM [7]
# # FROM TO RECON PERMUTE PGT DELTA STAT REPS
# # <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr> <int>
# # 1 IGHA1 IGHA2 0.0972 0.0994 0.56 -0.00219 SP 100
# # 2 IGHA1 IGHD 0 0 1 0 SP 100
# # 3 IGHA1 IGHG1 0 0 1 0 SP 100
# # 4 IGHA1 IGHG2 0.00815 0.00833 0.58 -0.000184 SP 100
# # 5 IGHA1 IGHG3 0 0 1 0 SP 100
# # 6 IGHA1 IGHG4 0.00352 0.00409 0.59 -0.000575 SP 100
# # 7 IGHA2 IGHA1 0 0 1 0 SP 100
# # 8 IGHA2 IGHD 0 0 1 0 SP 100
# # 9 IGHA2 IGHG1 0 0 1 0 SP 100
# #10 IGHA2 IGHG2 0 0 1 0 SP 100
# #11 IGHA2 IGHG3 0 0 1 0 SP 100
# #12 IGHA2 IGHG4 0 0 1 0 SP 100
# #13 IGHD IGHA1 0 0 1 0 SP 100
# #14 IGHD IGHA2 0 0 1 0 SP 100
# #15 IGHD IGHG1 0.0144 0.0127 0.42 0.00173 SP 100
# #16 IGHD IGHG2 0.00926 0.00981 0.54 -0.000545 SP 100
# #17 IGHD IGHG3 0.0157 0.0113 0.3 0.00440 SP 100
# #18 IGHD IGHG4 0.00203 0.00271 0.58 -0.000681 SP 100
# #19 IGHG1 IGHA1 0.00632 0.00757 0.59 -0.00126 SP 100
# #20 IGHG1 IGHA2 0.00326 0.00419 0.51 -0.000932 SP 100
# #21 IGHG1 IGHD 0 0 1 0 SP 100
# #22 IGHG1 IGHG2 0.0417 0.0545 0.86 -0.0127 SP 100
# #23 IGHG1 IGHG3 0 0 1 0 SP 100
# #24 IGHG1 IGHG4 0.0451 0.0493 0.66 -0.00427 SP 100
# #25 IGHG2 IGHA1 0 0 1 0 SP 100
# #26 IGHG2 IGHA2 0.00542 0.00511 0.45 0.000316 SP 100
# #27 IGHG2 IGHD 0 0 1 0 SP 100
# #28 IGHG2 IGHG1 0 0 1 0 SP 100
# #29 IGHG2 IGHG3 0 0 1 0 SP 100
# #30 IGHG2 IGHG4 0.0715 0.0581 0.08 0.0133 SP 100
# #31 IGHG3 IGHA1 0.0716 0.0600 0.16 0.0116 SP 100
# #32 IGHG3 IGHA2 0.0508 0.0492 0.37 0.00154 SP 100
# #33 IGHG3 IGHD 0 0 1 0 SP 100
# #34 IGHG3 IGHG1 0.159 0.181 0.96 -0.0225 SP 100
# #35 IGHG3 IGHG2 0.228 0.203 0.05 0.0250 SP 100
# #36 IGHG3 IGHG4 0.161 0.175 0.8 -0.0137 SP 100
# #37 IGHG4 IGHA1 0 0 1 0 SP 100
# #38 IGHG4 IGHA2 0.00596 0.00437 0.27 0.00159 SP 100
# #39 IGHG4 IGHD 0 0 1 0 SP 100
# #40 IGHG4 IGHG1 0 0 1 0 SP 100
# #41 IGHG4 IGHG2 0 0 1 0 SP 100
# #42 IGHG4 IGHG3 0 0 1 0 SP 100
## ----eval=FALSE, warning=FALSE, message=FALSE---------------------------------
# sp = testSP(switches$switches, alternative="greater", to="IGHA2")
# print(sp$means)
# # A tibble: 8 × 8
# # Groups: FROM [8]
# # FROM TO RECON PERMUTE PGT DELTA STAT REPS
# # <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr> <int>
# #1 IGHA1 IGHA2 0.620 0.634 0.59 -0.0144 SP 100
# #2 IGHD IGHA2 0 0 1 0 SP 100
# #3 IGHE IGHA2 0 0 1 0 SP 100
# #4 IGHG1 IGHA2 0.0116 0.0292 0.85 -0.0176 SP 100
# #5 IGHG2 IGHA2 0.0371 0.0252 0.32 0.0119 SP 100
# #6 IGHG3 IGHA2 0.297 0.283 0.45 0.0144 SP 100
# #7 IGHG4 IGHA2 0.0346 0.0290 0.36 0.00563 SP 100
# #8 IGHM IGHA2 0 0 1 0 SP 100