|
| 1 | +# Analyzing output |
| 2 | + |
| 3 | +## Loading output |
| 4 | + |
| 5 | +### `PhysiCellSnapshot` |
| 6 | +The base unit of PhysiCell output is the [`PhysiCellSnapshot`](@ref). |
| 7 | +These are currently considered pcvct internals and so the API may change. |
| 8 | +Each snapshot records the path to the PhysiCell output folder, its index in the sequence of outputs, the time of the snapshot in the simulation, and optionally the cell, substrate, and mesh data at that snapshot. |
| 9 | + |
| 10 | +### `PhysiCellSequence` |
| 11 | +A [`PhysiCellSequence`](@ref) is the full sequence of snapshots corresponding to a single PhysiCell simulation. |
| 12 | +As with `PhysiCellSnapshot`'s, these are currently considered internals and their API may change. |
| 13 | +In addition to the path to the PhysiCell output folder and the vector of `PhysiCellSnapshot`'s, it holds metadata for the simulation. |
| 14 | + |
| 15 | +### `getCellDataSequence` |
| 16 | +The main function to get sequences of cell data is [`getCellDataSequence`](@ref). |
| 17 | +It accepts any of a simulation ID (`<:Integer`), a simulation (`::Simulation`), or a sequence (`::PhysiCellSequence`) and either a single label (`::String`) or a vector of labels (`::Vector{String}`). |
| 18 | +For each cell in the simulation (as determined by the cell ID), the output creates a dictionary entry (the key is the integer cell ID) whose value is a named tuple with the input labels as keys as well as `:time`. |
| 19 | +This means that if one sets |
| 20 | + |
| 21 | +```julia |
| 22 | +data = getCellDataSequence(1, "position") |
| 23 | +``` |
| 24 | +Then one can access the positions of the cell with ID 78 by |
| 25 | +```julia |
| 26 | +cell_78_positions = data[78].position # an Nx3 matrix for the N integer-indexed outputs (ignores the `initial_*` and `final_*` files) |
| 27 | +``` |
| 28 | +and plot the x-coordinates of this cell over time using |
| 29 | +```julia |
| 30 | +cell_78_times = data[78].time |
| 31 | + |
| 32 | +using Plots |
| 33 | +plot(cell_78_times, cell_78_positions[:,1]) |
| 34 | +``` |
| 35 | + |
| 36 | +**Note**: Each call to `getCellDataSequence` will load *all* the data unless a `PhysiCellSequence` is passed in. |
| 37 | +Plan your analyses accordingly as loading simulation data is not fast. |
| 38 | + |
| 39 | +## Population plots |
| 40 | + |
| 41 | +### Group by Monad |
| 42 | +Plotting population plots is one the most basic analysis tasks and pcvct makes it super easy! |
| 43 | +If you call `plot` on a `Simulation`, `Monad`, `Sampling`, or the return value of a call to `run` (though not for a sensitivity analysis), |
| 44 | +then a sequence of panels will be generated in a single figure. |
| 45 | +Each panel will correspond to a `Monad` (replicates using the same parameter values) and will plot mean +/- SD for each cell type. |
| 46 | + |
| 47 | +Finer-grained control of the output is possible, too! |
| 48 | +- to include dead cells in your counts: `plot(...; ..., include_dead=true, ...)` |
| 49 | +- select a subset of cell types to include: `plot(...; ..., include_cell_types="cancer", ...)` |
| 50 | +- select a subset of cell types to exclude: `plot(...; ..., exclude_cell_types="cancer", ...)` |
| 51 | + |
| 52 | +The `include_cell_types` and `exclude_cell_types` can also accept a `Vector{String}` to include or exclude certain cell types, respectively. |
| 53 | +Furthermore, if the value of `include_cell_types` is a `Vector` and one of its entries is a `Vector{String}`, pcvct will interpret this to sum up those cell types. |
| 54 | +In other words, to get the total tumor cell count in addition to the epithelial (`"epi"`) and mesenchymal (`"mes"`) components, you could use |
| 55 | +```julia |
| 56 | +using Plots |
| 57 | +plot(Monad(1); include_cell_types=["epi", "mes", ["epi", "mes"]]) |
| 58 | +``` |
| 59 | + |
| 60 | +Finally, this makes use of Julia's Plot Recipes (see [RecipesBase.jl](https://docs.juliaplots.org/stable/RecipesBase/)) so any standard plotting keywords can be passed in: |
| 61 | +```julia |
| 62 | +using Plots |
| 63 | +colors = [:blue :red] # Note the absence of a `,` or `;`. This is how Julia requires different series parameters to be passed in |
| 64 | +plot(Simulation(1); color=colors, include_cell_types=["cd8", "cancer"]) # will plot cd8s in blue and cancer in red. |
| 65 | +``` |
| 66 | + |
| 67 | +### Group by cell type |
| 68 | +Invert the above by including all data for a single cell type across all monads in a single panel with a call to `plotbycelltype`. |
| 69 | +This function works on any `T<:AbstractTrial` (`Simulation`, `Monad`, `Sampling`, or `Trial`) as well as any `PCVCTOutput` object (the return value to `run`). |
| 70 | +Everything above for `plot` applies here. |
| 71 | + |
| 72 | +```julia |
| 73 | +using Plots |
| 74 | +plotbycelltype(Sampling(1); include_cell_types=["epi", "mes", ["epi", "mes"]], color=[:blue :red :purple], labels=["epi" "mes" "both"], legend=true) |
| 75 | +``` |
| 76 | + |
| 77 | +## Substrate analysis |
| 78 | +pcvct supports two ways to summarize substrate information over time. |
| 79 | + |
| 80 | +### AverageSubstrateTimeSeries |
| 81 | +An [`AverageSubstrateTimeSeries`](@ref) gives the time series for the average substrate across the entire domain. |
| 82 | + |
| 83 | +```julia |
| 84 | +simulation_id = 1 |
| 85 | +asts = AverageSubstrateTimeSeries(simulation_id) |
| 86 | +using Plots |
| 87 | +plot(asts.time, asts["oxygen"]) |
| 88 | +``` |
| 89 | + |
| 90 | +### `ExtracellularSubstrateTimeSeries` |
| 91 | +An [`ExtracellularSubstrateTimeSeries`](@ref) gives the time series for the average substrate concentration in the extracellular space neighboring all cells of a given cell type. |
| 92 | +In a simulation with `cd8` cells and `IFNg` diffusible substrate, plot the average concentration of IFNg experienced by CD8+ T cells using the following: |
| 93 | + |
| 94 | +```julia |
| 95 | +simulation_id = 1 |
| 96 | +ests = ExtracellularSubstrateTimeSeries(simulation_id) |
| 97 | +using Plots |
| 98 | +plot(ests.time, ests["cd8"]["IFNg"]) |
| 99 | +``` |
| 100 | + |
| 101 | +## Motility analysis |
| 102 | +The `motilityStatistics` function returns the time alive, distance traveled, and mean speed for each cell in the simulation. |
| 103 | +For each cell, these values are split amongst the cell types the given cell assumed throughout (or at least at the save times). |
| 104 | +To calculate these values, the cell type at the start of the save interval is used and the net displacement is used to calculate the speed. |
| 105 | +Optionally, users can pass in a coordinate direction to only consider speed in a given axis. |
| 106 | + |
| 107 | +```julia |
| 108 | +simulation_id = 1 |
| 109 | +mss = motilityStatistics(simulation_id) |
| 110 | +all_mean_speeds_as_mes = [ms["mes"].speed for ms in mss if haskey(ms, "mes")] # concatenate all speeds as a "mes" cell type (if the given cell ever was a "mes") |
| 111 | +all_times_as_mes = [ms["mes"].time for ms in mss if haskey(ms, "mes")] # similarly, get the time spent in the "mes" state |
| 112 | +mean_mes_speed = all_mean_speeds_as_mes .* all_times_as_mes |> sum # start computing the weighted average of their speeds |
| 113 | +mean_mes_speed /= sum(all_times_as_mes) # finish computing weighted average |
| 114 | +``` |
| 115 | + |
| 116 | +```julia |
| 117 | +mss = motilityStatistics(simulation_id; direction=:x) # only consider the movement in the x direction |
| 118 | +``` |
0 commit comments