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) # 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 explicitely when creating the map with the
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 1.43 mins
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
or you can view the map interactively using the
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 #> #> 'SERUM 2' #> Took 2.93 mins # View the result view(map3d)