4_Connectomics.Rmd

---
title: 4_Connectomics
output: html_notebook
---

## Setup
Basic installation/update:
```{r, eval=FALSE}
natmanager::install('natverse')
natmanager::install(pkgs = 'jefferis/vfbconnectr')
install.packages('forcats')
# update Python package
vfbconnectr::vc_install()
```

Reinstall catmaid package since we have updated it and add a new package, forcats
```{r}
natmanager::install(pkgs = 'rcatmaid')
if(!requireNamespace('forcats'))
  install.packages('forcats')
```

```{r, message=FALSE}
# load packages
library(vfbconnectr)
library(natverse)
library(dplyr)
#library(plyr)
library(ggplot2)

# if Running on Rstudio server (==cloud) use plotly as visualisation default
#if(isTRUE(Sys.getenv('RSTUDIO_PROGRAM_MODE')=="server"))
options(nat.plotengine = 'plotly')
```

Set up connections to VFB APIs and VFB CATMAID server.

```{r}
vc = VfbConnect()
# Connect to the VFB CATMAID server hosting the FAFB data
catmaid::vfbcatmaid('fafb')
```


### VFB provides a simple way in to exploring connectomics data

Get connections between defined upstream and downstream neuron types:

```{r}
vc$get_connected_neurons_by_type(upstream_type='GABAergic neuron', downstream_type='adult descending neuron', weight=10)
```

Get neurons downstream of a neuron

```{r}
vc$get_neurons_downstream_of('DA3_adPN_R - 1703683361', weight = 20)
```

Get neurons downstream of a neuron, filtered by a classification.

```{r}
vc$get_neurons_downstream_of('DA3_adPN_R - 1703683361', classification="'Kenyon cell'", weight = 20)
```

### More sophisticated connectomics queries require direct queries of connectomics DB APIs

#### CATMAID connectivity queries

Broadly speaking, CATMAID lets you fetch connectivity data as either lists of up- and downstream partners or as whole adjacency matrices.

These examples use the VFB FAFB server.

```{r}
# Using VFB to get neurons by type
DA1_tab = vc$get_instances("'adult antennal lobe projection neuron DA1'", summary=TRUE) %>%
  vc_df()
DA1_tab
```

Where do these data come from?

```{r}
table(DA1_tab$data_source)
```

So a mix of fafb and hemibrain (and in fact it's not an accident that there are
approximately 2x as many FAFB as hemibrain). However there are 6 without a data
source. How about them?

```{r}
DA1_tab %>%
  filter(data_source=="") %>%
  select(label, id, dataset, templates)
```
Ok so those are from light level datasets, although there is something funny to me about the labels, since they are all identical and the fact that some do not
have templates.

Now let's map neurons to CATMAID Skeleton IDs (skids)
```{r}
# note that db='catmaid_fafb' filters the returned values
# while reverse_return 
da1_id_list = vc$neo_query_wrapper$vfb_id_2_xrefs(
  DA1_tab[,'id'], 
  db='catmaid_fafb', 
  reverse_return=TRUE)
# the names of each element in the returned list are the actual ids
da1_skids = names(da1_id_list)
da1_skids
```

We can now fetch the neurons like this:
```{r}
da1 = read.neurons.catmaid(da1_skids)
```

Let's generate a connectivity table for neurons downstream of DA1 neurons

```{r}
da1_ds = catmaid_get_connector_table(da1_skids, direction = 'outgoing')
da1_ds
```
That's a lot of rows. Let's do some quick summaries:

```{r}
da1_ds %>%
  count(skid)
```
```{r}
library(ggplot2)
da1_ds %>%
  count(skid) %>%
  # note that we turn skid into a factor to make it categorical
  mutate(skid=factor(skid)) %>%
  # and sort the levels by frequency
  mutate(skid=forcats::fct_reorder(skid, desc(n))) %>% 
  # xlim insists that the X axis starts at 0
  qplot(n, skid, data=.) 
```

I wonder why there are such big differences in synapse count. I have a hunch
that they might not have been uniformly reconstructed. Let's add in side
information by using the metadata we fetched earlier. The *label* field includes
the side of the brain `Uniglomerular mALT DA1 lPN#L2 (FAFB:2319457)`.

```{r}
da1_ds %>%
  count(skid) %>%
  merge(y=DA1_tab, by.y="accession", by.x='skid') %>% 
  mutate(side=stringr::str_match(label, '#([LR])')[,2]) %>%
  # note that we turn skid into a factor to make it categorical
  mutate(skid=factor(skid)) %>%
  # and sort the levels by frequency
  mutate(skid=forcats::fct_reorder(skid, desc(n))) ->
  da1_ds.summary


  # xlim insists that the X axis starts at 0
  qplot(n, skid, xlim=c(0, NA), col=side, data=da1_ds.summary) 
```

So a cautionary tale!

We can also visualise the individual cell types, like so:

```{r}
qplot(n, skid, xlim=c(0, NA), label=side, col=parents_label, data=da1_ds.summary, geom='text') 
```

By the way you can mix python and R analysis

```{python}
import pymaid
rm=pymaid.connect_catmaid(server="https://fafb.catmaid.virtualflybrain.org/", api_token=None, max_threads=10)
bates = pymaid.find_neurons(annotations='Paper: Bates and Schlegel et al 2020')
adj = pymaid.adjacency_matrix(bates)

import matplotlib.pyplot as plt
import seaborn as sns 
ax = sns.clustermap(adj, vmax=10, cmap='Greys')
plt.show()
```



```{r}
# Get an adjacency matrix between all Bates, Schlegel et al. neurons
adj = catmaid_adjacency_matrix('Paper: Bates and Schlegel et al 2020')
adj[1:5,1:5]
```

```{r}
heatmap(adj, col=grey.colors(20, start = 1, end=.3))
```

Let's look at axo-axonic connections between two different types of DA1 PNs

```{r}
cn = catmaid_get_connectors_between(2863104, 1811442)
cn
```

```{r}
nclear3d()
plot3d(da1[[1]], WithConnectors=TRUE)
```

#### neuprint connectivity queries

To programmatically interface with neuprint, we will use neuprintr (link). It requires an API token which you can get via the website and is bound to the Google account that you use to log into neuprint. For this workshop we provide such a token as environment variable but you will need to start using your own token after the workshop is over.

These examples use the hemibrain v1.1 dataset.

### neuprintr
First we have to initialize the connection.

```{r}
library('neuprintr')
neuprint_login()
```

Most functions in `neuprintr` allow queries that can be used to look for specific body IDs, types, etc:

```{r}
?neuprintr
```

### Fetching neurons

Let's say we want to find all antennnal lobe projection neurons (PNs). Their type nomenclature adheres to `{glomerulus}_{lineage}PN` (e.g. `DA1_lPN`)for uniglomerular PNs and a `M_{lineage}PN{tract}{type}` (e.g. `M_vPNml50` = "**m**ultiglomerular **v**entral lineage **PN** **m**edio**l**ateral tract type **50**) for multiglomerular PNs. 

To get them all, we need to use regex patterns (see this [cheatsheet](https://cheatography.com/davechild/cheat-sheets/regular-expressions/)):

```{r}
# get neurons that match a filter criteria (using regexp like here)
pns = neuprint_search('.*?_.*?PN.*?')

print(paste0('The number of PNs found were :',nrow(pns)))
head(pns)

```

```{r}
# Check that the regex did not have any accidental by-catch
unique(pns$type)
```

Find neuron types with VFB instead:

```{r}
ALPNs = vc$get_instances("'adult antennal lobe projection neuron'", summary=TRUE)
ALPNs_df = do.call(rbind, ALPNs)
```

### Fetching synaptic partners

Looks good! Next: What's downstream of those PNs?

```{r}
ds = neuprint_connection_table(pns$bodyid, prepost = "POST", details = TRUE)
colnames(ds)[colnames(ds) == 'name'] <- 'instance_post'
colnames(ds)[colnames(ds) == 'type'] <- 'type_post'
ds_merged = merge(ds, pns, by.x='bodyid', by.y='bodyid')[, c(1:6,15)]
colnames(ds_merged)[colnames(ds_merged) == 'name'] <- 'type_pre'
ds_merged = ds_merged[,c("bodyid","partner","prepost","weight","type_pre","type_post")]
ds_merged = ds_merged[complete.cases(ds_merged), ]
head(ds_merged)
```

Each row is now a connections from a single up- to a single downstream neuron. The "weight" is the number of synapses between the pre- and the postsynaptic neuron. Let's simplify by grouping by type:

```{r}
groupColumns = c("type_pre","type_post")
dataColumns = c("weight")
by_type = plyr::ddply(ds_merged, groupColumns, function(x) colSums(x[dataColumns]))
head(by_type)
```

The strongest connections are between PNs and Kenyon Cells (KCs). That's little surprising since there are thousands of KCs. For the sake of the argument let's say we want to know _where_ these connections occur:

```{r}
kcs = neuprint_search('KC.*?')
pns_kcs = c(pns$bodyid, kcs$bodyid)
adj = neuprint_connection_table(bodyids=pns_kcs, by.roi = T)

adj_select = adj[adj$bodyid %in% pns_kcs & adj$partner %in% pns_kcs,]

```

```{r}
# Group by region of interest (ROI)
groupColumns = c("roi")
dataColumns = c("ROIweight")
by_roi = plyr::ddply(adj_select, groupColumns, function(x) colSums(x[dataColumns]))
by_roi = by_roi[complete.cases(by_roi), ]

```

```{r}
g <- ggplot(by_roi, aes(roi))
# Number of cars in each class:
g <- g + geom_bar(aes(weight = ROIweight))

g

```

#### Querying paths
Let's say we want to find out how to go from a PN (second order olfactory neurons) all the way to a descending neuron (presumably leading to motor neurons in the VNC).

```{r}
# First fetch the DNs
dns = neuprint_search('.*DN[^1]{0,}.*|Giant Fiber',type = "name")
head(dns)
```

Neuprint lets you query paths from a single source to a single target. For multi-source or -target queries, your best bet is to download the entire graph and run the queries locally using networkx or igraph.

```{r}
# Find all paths from A PN to A DNs 
#@Greg: This function simply hangs,so I just used one bodyid..
paths = neuprint_get_shortest_paths(body_pre= pns$bodyid, body_post = dns$bodyid[1])
```

So it looks like there are three separate 7-hop paths to go from `M_vPNml53` to `DN1a`. Let's visualize the neurons involved!

### Fetching meshes & skeletons
You can fetch skeletons as SWCs directly via `neuprintr`. For visualization however it's easiest to load neuron morphologies via `nat`. For that `nat` wraps `neuprintr` and adds some convenience functions:


```{r}
neuprint_datasets()
```

```{r}
nl = neuprint_read_neurons(bodyids = paths$from[1:2])
```

```{r}
al_r.mesh = neuprint_ROI_mesh(roi = "AL(R)")
lh_r.mesh = neuprint_ROI_mesh(roi = "LH(R)")
ca_r.mesh = neuprint_ROI_mesh(roi = "CA(R)")
```

```{r}
options(nat.plotengine = 'plotly')
nclear3d()
plot3d(nl,plotengine = 'plotly')
wire3d(al_r.mesh, col = 'red', alpha = 0.1)
wire3d(lh_r.mesh,col = 'magenta', alpha = 0.1)
wire3d(ca_r.mesh,col = 'blue', alpha = 0.1)
```