---
title: "DoOR analysis tools"
author: "Daniel Münch"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{DoOR analysis tools}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

The `DoOR.functions` package provides tools for the analysis of the data
provided by the `DoOR.data` package.

## Content
* [Loading data](#loading)
* [Identifying the sensillum we are recording from with `identify_sensillum()`](
#identify_sensillum)
* [Finding neuron-specific odorants with `private_odorant()`](#private_odorant)
* [Mapping response data from an unknown source with `map_receptor()`](
#map_receptor)
* [Changing the response unit with `back_project()`](#back_project)
* [Quantifying tuning curves with `sparse()`](#sparse)



## Loading data{#loading}
First we need to load packages and data:
```{r, echo = TRUE, message = TRUE, results='hide'}
#load data
library(DoOR.functions)
library(DoOR.data)
load_door_data(nointeraction = TRUE)
```

## Identifying the sensillum we are recording from with `identify_sensillum()` 
{#identify_sensillum}
Imagine we perform an electrophysiological recording from a _Drosophila_
sensillum (single sensillum recording, SSR) and we are not sure what sensillum
we are recording from. In order to identify the sensillum we used several
diagnostic odorants (maybe selected using
[`private_odorant()`](#private_odorant)) and got responses from the different
sensory neurons the sensillum houses. We can now pass our recorded data to
`identify_sensillum()`.

Let's make up some simple fake data. We pretend to have recorded with three
odorants (2,3-butanedione, ethanoic acid and carbon dioxide) and we could
separate the responses of two units. Unit1 responded strongly to 2,3-butanedione
only, unit2 only responded to carbon dioxide. We create a data.frame that
contains a column called `odorants` with the InChIKeys of our test odorants, and
one column for each unit (name the colnames as you like, e.g. unit1-n or Aneuro
if you are sure about the neuron).

```{r}
recording <- data.frame(
  odorants = c(trans_id(c("BEDN", "ETAS"), "Code"),
               trans_id("carbon dioxide", "Name")),
  unit1 = c(.9,.1,.1),
  unit2 = c(0, .1, 1)
)
```
Next we feed the recording to the function:

### using correlations
```{r, fig.width=7.1, fig.height=4.5}
identify_sensillum(recording, base_size = 8)
```
Note that the function tells us that it found hits for all units in ab1 and ab5,
meaning that e.g. within the four neurons housed in the ab1 sensillum both of
our units had good matches. You can set this correlation threshold with
`min.cor`. If we increase the threshold to 0.99 only ab1 is returned as a double
match:

```{r, fig.width=7.1, fig.height=4.5, fig.show='hide'}
identify_sensillum(recording, min.cor = .99)
```

We can define the number of best hits that we want to get returned (the default
is 10):
```{r, fig.width=7.1, fig.height=4.5}
identify_sensillum(recording, nshow = 5, base_size = 8)
```
And if we know e.g. that we are recording from a basiconic sensillum we can
restrict the search to one or a few sensilla types:
```{r, fig.width=7.1, fig.height=4.5}
identify_sensillum(recording, sub = "ab", nshow = 5, base_size = 8)
identify_sensillum(recording, sub = c("ac","at"), nshow = 5, base_size = 8)
```

### using Euclidean distances
Instead of correlations we can also use the Euclidean distance as a
(dis)similarity measure:

```{r, fig.width=7.1, fig.height=4.5}
identify_sensillum(recording, method = "dist", sub = "ab", nshow = 5, 
                   base_size = 8)
```

### returning data instead of plots
We can also return the correlation/distance data instead of the plot when
setting `plot =FALSE`:
```{r}
sensillumX <-
  identify_sensillum(recording,
  method = "dist",
  sub = "ab",
  plot = FALSE)
  head(sensillumX)
```

So apparently our fake recording came from the ab1 sensillum, which was
admittedly quite obvious as we had a strong carbon dioxide response and ab1
houses the carbon dioxide receptor :)



## Finding neuron-specific odorants with `private_odorant()` {#private_odorant}
There may be several cases where we might be interested in so called 
_private odorants_, odorants that specifically activate a given receptor or 
sensory neuron. Maybe we are looking for diagnostic odorants for sensillum
identification or we want to activate a specific neuronal pathway,
`private_odorant()` returns candidate odorants for that task.

Let's say we want to specifically activate Or22a neurons:
```{r}
private_odorant("Or22a")
```

We might want to return the odorant names instead of InChiKeys:
```{r}
private_odorant("Or22a", tag = "Name")
```
So according to the function sec-amyl acetate would be a good candidate. It
activates Or22a at 0.4 (DoOR response, max is 1) while the maximum activation in
all other tested responding units (receptors, neurons, glomeruli) is 0.016, a
difference of 0.40. Sounds good, but it was tested only in 4 other responding
units, so I would rather go for ethyl hexanoate with about the same difference
but being tested in 29 other responding units.

We can also restrict the search to the sensillum the responding units of
interest is related to:
```{r}
private_odorant("Or22a", tag = "Name", sensillum = T)
```
Ethyl 2-methylbutanoate sounds like a good hit, it has the same difference to
the other units as ethyl hexanoate but hardly elicits a response at all from the
other neuron. The n of 1 is fine as there are only 2 neurons housed in the ab3
sensillum.



## Mapping response data from an unknown source with `map_receptor()`
{#map_receptor}
Similar to `identify_sensillum()`, `map_receptor()` correlates a response vector
to all responding units of the existing DoOR consensus data. Let's grab a data
set from Or22a and see where it ends up:
```{r}
data <- data.frame(odorants  = Or22a$InChIKey, responses = Or22a$Hallem.2006.EN)
data <- na.omit(data)
head(data)
map_receptor(data = data, nshow = 5)
```
This example was a bit circular as the tested data contributed to the consensus
data...



## Changing the response unit (spikes, deltaF/F, ...) with `back_project()`
{#back_project}
The DoOR consensus data is normalized to values between 0 and 1. If we want to
compare the DoOR data to one of our own recordings, it would be great to have
the DoOR data in the same unit as our owne data (e.g. spikerate). This is what
we can do with `back_project()`,  we can rescale the DoOR data to fit a given
response template.

As an example, let's take the data Hallem et al. recorded from Or22a _via_
calcium imaging and rescale the DoOR responses accordingly. The template has to
have 2 columns named `odorants` and `responses`:
```{r, fig.width = 4.5, fig.height = 4.5}
template <- data.frame(odorants  = Or22a$InChIKey,
                       responses = Or22a$Hallem.2006.EN)

bp <- back_project(template, responding.unit = "Or22a")
plot(bp$back_projected$original.data,
     bp$back_projected$back_projected.data,
     xlab = "DoOR consensus response",
     ylab = "back_projected data [spikes, Hallem.2006.EN]"
)

head(bp$back_projected)

```
All the yellow lines in the first plot represent odorant responses that were not
available in the original data set but were projected onto the fitted function
and rescaled to the units in "Hallem.2006.EN". The second plot shows the
relationship between the rescaled data and the original consensus responses.


## Quantifying tuning curves with `sparse()`{#sparse}
The width of a tuning curve, i.e. for example to how many odorants a receptor
shows strong responses, can be quantified using different sparseness measures.
We implemented kurtosis[^1] and sparseness[^2] in `sparse()`. A high kurtosis or
sparseness value indicates a narrow tuning curve (see also `dplot_tuningCurve()`
in the [DoOR visualization vignette](DoOR_visualizations.html#tuningCurve)).
While kurtosis is able to deal with negative values, sparseness can't and thus
all values need to be be transformed to absolute values first. Sparseness scales
between 0 and 1, kurtosis between -∞ and ∞, a kurtisis of 0 corresponds to the
Gaussian distribution.

```{r}
rm.SFRreset <- reset_sfr(door_response_matrix, "SFR")

sparse(x = rm.SFRreset[,"Or69a"], method = "ltk")
sparse(x = rm.SFRreset[,"Or69a"], method = "lts")

sparse(x = rm.SFRreset[,"Gr21a.Gr63a"], method = "ltk")
sparse(x = abs(rm.SFRreset[,"Gr21a.Gr63a"]), method = "lts")

```



[^1]: Willmore, B., Tolhurst, D.J., 2001. Characterizing the sparseness of
neural codes. Network 12, 255–270. dx.doi.org/10.1080/net.12.3.255.270
[^2]: Bhandawat, V., Olsen, S.R., Gouwens, N.W., Schlief, M.L., Wilson, R.I.,
2007. Sensory processing in the Drosophila antennal lobe increases reliability
and separability of ensemble odor representations. Nature neuroscience 10,
1474–82. dx.doi.org/10.1038/nn1976
