vignettes/making-a-map-from-scratch.Rmd
making-a-map-from-scratch.Rmd
In this vignette we’re going to go through the process of making an antigenic map using Racmacs and viewing and saving the resulting map and images and interactive plots of the map.
The steps will be:
First of all we need to read in the titer data from a file, exactly how you do this depends on the format of your data but ultimately you need to end up with a matrix or dataframe of titers where rows represent antigens and columns represent sera.
If your titer data is in a simple format where you just have a single
header for serum names and a single first column of antigen names you
can try using the function read.titerTable()
, as in the
example below.
For this example we will be reading the titer data from the first antigenic map of H3N2 antigenic evolution, published in 2004 and included with the Racmacs package.
# Load the Racmacs package
library(Racmacs)
# Set an option for the number of computer cores to run in parallel when optimizing maps
# The default when running on CRAN is to use 2
options(RacOptimizer.num_cores = 1)
# However you can also set the number of cores to the maximum number like this
# options(RacOptimizer.num_cores = parallel::detectCores())
# Read in the titer table
path_to_titer_file <- system.file("extdata/h3map2004_hitable.csv", package = "Racmacs")
titer_table <- read.titerTable(path_to_titer_file)
The resulting titer table is a matrix of titers with column and row names:
print(titer_table[1:5,1:7])
#> HK/1/68 HK/107/71 BI/21793/72 EN/42/72 PC/1/73 BI/5930/74 SL/840/74
#> BI/15793/68 "354" "*" "*" "501" "109" "*" "*"
#> BI/16190/68 "640" "*" "640" "320" "80" "320" "40"
#> BI/16398/68 "1163" "*" "*" "1163" "291" "*" "73"
#> HK/1/68 "1280" "1280" "*" "2560" "53" "*" "<10"
#> BI/808/69 "320" "*" "*" "640" "120" "*" "80"
Now we have the titers in the correct format we can use it create an
antigenic map object or “acmap” object. For this we use the constructor
function acmap()
, akin to how data.frame()
is
used to make a data frame object or matrix()
for a
matrix.
# Create the acmap object, specifying the titer table
map <- acmap(
titer_table = titer_table
)
The resulting object, that we’ve given the name map
in
this case, only has properties associated with the table. If the table
had column names and row names, serum and antigen names are inferred
from these. Alternatively, you can set these explicitly when creating
the map with the ag_names
and sr_names
arguments.
Although we now have the acmap object there is no “map” yet as such,
since we need to invoke the optimiser function to try and determine the
antigenic map that best represents the patterns of reactivity seen in
the titer table. To do this we call the function
optimizeMap()
on the acmap object we’ve made, specifying
the number of optimization runs and other parameters like the minimum
column basis and number of dimensions.
In this case we’re going to try and find the best 2-dimensional solution for the map data, using 500 optimization runs.
# Perform some optimization runs on the map object to try and determine a best map
map <- optimizeMap(
map = map,
number_of_dimensions = 2,
number_of_optimizations = 500,
minimum_column_basis = "none"
)
#> a
#> Took 51.6 secs
#> Warning in optimizeMap(map = map, number_of_dimensions = 2,
#> number_of_optimizations = 500, : There is some variation (9.57 AU for one
#> point) in the top runs, this may be an indication that more optimization runs
#> could help achieve a better optimum. If this still fails to help see
#> ?unstableMaps for further possible causes.
By default the resulting map optimizations are sorted by their stress, so that the optimization that resulted in the lowest map stress (i.e. usually thought of as the best solution), will be first and will be the default optimization run selected and therefore will be the one used when querying anything like antigen and serum coordinates, map stress etc. unless otherwise specified.
The lowest stress optimization run will also be the one used by
default when viewing the map. To view the map, you can either do a
simple plot of the map using the plot()
function
plot(map)
or you can view the map interactively using the view()
function
view(map)
make.acmap()
Note that there is also a small utility function provided to allow
you to create the acmap object and run optimizations on it all in one
step. In the example below we use it to create a 3d version of the
antigenic map and then use view()
map to view the result
(plot()
will not work with 3 dimensional maps).
# Make the acmap object and run optimizations
map3d <- make.acmap(
titer_table = titer_table,
number_of_dimensions = 3,
number_of_optimizations = 500,
minimum_column_basis = "none"
)
#> a
#> Warning: The following SERA have do not have enough titrations to position in 3 dimensions. Coordinates were still optimized but positions will be unreliable
#>
#> 'HK/107/71'
#> Took 1.72 mins
#> Warning in optimizeMap(map = map, number_of_dimensions = number_of_dimensions,
#> : There is some variation (5.14 AU for one point) in the top runs, this may be
#> an indication that more optimization runs could help achieve a better optimum.
#> If this still fails to help see ?unstableMaps for further possible causes.
# View the result
view(map3d)