--- title: "Tabular Model File" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Tabular Model File} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` ```{r setup} library(uqsa) library(errors) ``` The models we use in this package are written in a very loose tabular form: the model is a collection of TSV files (tab-separated-values). Tabs are the best separator in our opinion (as compared to commas or semi-colon separated files), because it _can_ be entered with a text editor into the file, by hand, but is quite hard to enter as the content of a field in a spreadsheet program. So, Tabs usually don't appear coincidentally in the files, but can be typed on purpose. Furthermore: - TSV files render nicely on sites such as GitHub - TSV files are text files and thus work well with git - TSV files do not require a library to read - Column-names can contain spaces, which is very human friendly We interpret the files partially using functions such as `grep` (in R), `tolower`, with settings such as `ignore.case` where possible, so that various different spellings are acceptable, e.g. "stdv", "sd", or "standard deviation" when describing a distribution. When reporting data values, we use parenthesized (concise) error notation, e.g.: `1.23(1)` which is intrepreted as `1.23` with an uncertainty of `0.01`, more examples: | string | value | standard-error | |---------------:|----------:|:---------------| | `1.23(1)` | 1.23 | 0.01 | | `1.23(1)e2` | 123 | 1 | | `-1.234(1)e-6` | -1.234e-6 | 1e-9 | | `1.23;0.01` | 1.23 | 0.01 | | `1.23±0.01` | 1.23 | 0.01 | The last two rows represent a fall-back notation with redundant zeros, which you can use if you really don't like parenthesized location, the parser automatically tries to find some separator if there are no parentheses in the string: ```{r} x <- parse_concise( c( "1.23(1)","1.23(1)e2","-1.234(1)e-6", # concise notation "1.23;0.01", "1.23±0.01" # fall-back ) ) print(class(x)) print(as.data.frame(x)) # this looks better as a data-frame ``` # Components of Models All quantities in the model have a unit. If your model is initially formulated using concentrations and reaction rates (kinetic laws), then you can still simulate it using the Gillespie solver. We convert all quantities with concentrations in the unit to particle counts automatically. So, feel free to formulate the model using concentrations. The role of each file is determined from the file's name: - **Constant.tsv**, list of constants that are not subject to optimization, calibration, or any other investigation - **Input.tsv**, known parameters that can change between different simulations - **Parameter.tsv**, possibly unknown parameters (e.g. ranges), subject to fitting/sampling/etc. - **Expression.tsv** - **Compound.tsv**, list of reacting compounds, e.g. `Ca` - **Reaction.tsv**, list of reactions - **Output.tsv**, algebraic functions that express a measurable value (or very close to measurable), + e.g. the total amount of something: `A + AB + AC`, the total amount of bound and free `A` - **Experiment.tsv**, list of experiments The first column in each TSV table contains some kind of unique _id_, the id of each _experiment_ (row) in the experiment table (Experiment.tsv) is also used as the name of the TSV file that contains the data for this experiment. # Data The TSV files correspond pretty directly to R `data.frame` objects, the import function `model_from_tsv` returns a list of data frames, named like the file. The model's name is taken from the directory (`dirname`) the models are stored in. This name is stored in the `comment` attribute. There are two main cases we want to distinguish in the same way that fungi and plants are distinct: 1. Experimental data that corresponds fairly well to an output function - there is an output function that can be compared to the data, perhaps up to some scaling constant - the data column is labeled exactly like the output function is named, e.g.: `AMPA_OUTPUT` 2. Experimental data that has a very complex relationship to the model - raw data, where several columns together can be used to make a comparison with something in the model ## One-to-One Correspondence The first case, with a one-to-one correspondence between output functions and data columns can be automatically parsed using the `experiments` function. This function reads several of the tables and determines the initial state, the input, the time-line for any given experiment and creates a data-matrix (`data`), with standard-errors using the `errors` package. Here, you can still create your own custom likelihood function that performs some kind of complex normalization to evaluate a simulation. It is not difficult to have a control experiment that contains a value we need to normalize with. The user supplied likelihood may perform such normalizations. ```{r data} m <- model_from_tsv(uqsa_example("AKAP79")) print(m$Reaction[,c(2,3,4)]) # an example table, the reactions print(head(m$Experiment)) # the list of experiments ``` This is how the output function corresponds to the data-label, **AKAR4pOUT**: ```{r label="data-and-output"} print(m$Output) rn <- rownames(m$Experiment) print(head(m[[rn[18]]])) # experiment 18, data found by name ``` So, the data-table (called like the row in the experiments table), has a column that corresponds to the output function. This is how we link the two together and know that the measured values have to somehow correspond to the output function `AKAR4pOUT`. ```{r, label="data-matrix"} x <- experiments(m) print(x[[18]]$data[,seq(8),drop=FALSE]) # a sub-set of the data-matrix ``` The default likelihood function will take the squared difference between the simulated values of `AKAR4pOUT` (the function) and the measured values in that column. ## Very Indirect Data In the second big case of very _cryptic_ data, you should not rely on the `data` matrix returned by the `experiments(m)` function and instead use the `measurements` field, which is a data frame, exactly as it was written in the TSV file, of just use the data data frame in `m`: ```{r measurements} print(head(m[[rn[18]]],12)) # these two should be the same print(head(x[[18]]$measurement,12)) # (a 12 row subset) ``` The only difference is that all numbers were parsed by `parse_concise` to resolve all parenthesized standard-errors (in `x`), in `m` the values are raw (either stings if they contain parentheses, or numbers if R could successfully coerce them. So, `m` can be used directly, if you have a very complex data-case and thus a very complex likelihood.