1 Python code to download annotation sizes from the BOSS

from intern.remote.boss import BossRemote
from intern.resource.boss.resource import *
from intern.utils.parallel import block_compute
import configparser
import requests
import numpy as np
from numpy import genfromtxt
import shutil
import blosc
import sys
import os
import itertools
from functools import partial
from multiprocessing import Pool
from multiprocessing.dummy import Pool as ThreadPool
from multiprocessing import cpu_count
import csv
import argparse


####
CONFIG_FILE    = 'neurodata.conf'

COLL_NAME      = 'collman' 
EXP_NAME       = 'M247514_Rorb_1_Site3Align2_EM' 
ANNO_NAME      = 'm247514_Site3Annotation_MN_global'
COORD_FRAME    = 'M247514_Rorb_1_Site3Align2_EM_m247514_Site3Annotation_MN_global'
OUTPUT         = 'annotation_sizes_pixels.csv'

config = configparser.ConfigParser()
config.read(CONFIG_FILE)
TOKEN = config['Default']['token']
boss_url = ''.join( ( config['Default']['protocol'],'://',config['Default']['host'],'/v1/' ) )
#print(boss_url)
#'https://api.boss.neurodata.io/v1/'

#intern
rem = BossRemote(CONFIG_FILE)

cf = CoordinateFrameResource(str(COLL_NAME + '_' + EXP_NAME))
cfr = rem.get_project(cf)
anno_res = ChannelResource(ANNO_NAME, COLL_NAME, EXP_NAME, 'annotation', datatype='uint64')

ex = {'x': cfr.x_stop, 'y': cfr.y_stop, 'z': cfr.z_stop}

blocks = block_compute(0,ex['x'],0,ex['y'],0,ex['z'],
           origin = (0,0,0), block_size = (512, 512, 16))

rid = []
for b in blocks:
    rid = rid + rem.get_ids_in_region(anno_res, 0, b[0], b[1], b[2], [0,1])

u = np.unique(np.asarray(rid))

def getSizeAnno(uniq):
    bb = rem.get_bounding_box(anno_res, 0, uniq, 'tight')
    out = [bb['id'], sizeAnno]
    return(out)

with ThreadTool(4) as tp:
    bbA  = np.asarray(tp.map(getSizeAnno, u))

with open(OUTPUT, 'w') as f1:
        wt = csv.writer(f1)
        wt.writerow(['id', 'size'])
        wt.writerows(bbA)

na <- read.csv("rorb_gaussianAvg_at.csv")
cnames <- colnames(na)
dat <- fread("rorb_gaussianAvg_at_GoogleDoc.csv")
colnames(dat) <- cnames
meta <- fread("GoogleDocData.csv")
loc <- fread("GoogleDocData.csv")[, 1:3]
ids <- as.data.frame(read.csv("GoogleDocData.csv")[, 4])
colnames(ids) <- 'id'

fdat <- cbind(dat, meta)
ccol <- c('blue', 'blue', 'blue', 'red', 'red', 'red', 'black', 'black', 'green', 'green', 'black', 'green')
ind <- order(ccol)
ccol <- sort(ccol)


## annotation ids in the BOSS are +1 from the ids in Forrest's gsheet.
annoSizes <- read.csv('annotation_sizes_pixels.csv')
annoSizes$id <- annoSizes$id - 1

tmp <- merge(ids, annoSizes, by = "id", all.x = TRUE)

dat <- merge(fdat, tmp, by = 'id')
dat <- dat[!is.na(size), ]
loc <- dat[, .(atx, aty, atz)]
sdat <- data.table(scale(dat[, cnames, with = FALSE], center = TRUE, scale = TRUE))
sdat$shaft <- dat$shaft
sdat$gaba <- dat$gabal
sdat$size <- dat$size
sdat$microm3 <- sdat$size * 3*3*50 / 1000^3

df1 <- sdat[, .(shaft =as.factor(shaft), size, microm3, PSD95, gaba)]

p <- ggplot(df1, aes(x = microm3, y = PSD95, color = shaft, shape = as.factor(gaba))) +
     geom_point(alpha = 0.4) + 
     scale_color_manual(values = c("red", "black")) +
     xlab(expression(paste(mu, "m^3"))) +
     ggtitle("PSD95 versus size colored by shaft")

p + geom_smooth(method = 'lm') + facet_grid(. ~ gaba, labeller = label_both)

2 KDE’s

p2 <- ggplot(sdat, aes(group = shaft))

p2 + geom_density(aes(x = PSD95, y = ..density.., col = as.factor(shaft))) + ggtitle("PSD95 | shaft") + facet_grid(. ~ gaba, labeller = label_both)

p2 + geom_density(aes(x = microm3, y = ..density.., col = as.factor(shaft))) + ggtitle("size | shaft") + facet_grid(. ~ gaba, labeller = label_both)

p2 + geom_density(aes(x = microm3, y = ..density.., col = as.factor(shaft))) + scale_x_log10() + ggtitle("size | shaft, x-axis on log10 scale") + facet_grid(. ~ gaba, labeller = label_both)

p2 + geom_violin(aes(x = as.factor(shaft), y = PSD95)) + geom_jitter(aes(x = as.factor(shaft), y = PSD95), alpha = 0.7) + facet_grid(. ~ gaba, labeller = label_both)

p2 + geom_violin(aes(x = as.factor(shaft), y = microm3)) + geom_jitter(aes(x = as.factor(shaft), y = microm3), alpha = 0.7) + facet_grid(. ~ gaba, labeller = label_both)

p3 <- ggplot(sdat[gaba == 1,], aes(group = shaft))
p3 + 
  geom_violin(aes(x = as.factor(shaft), y = microm3)) + 
  geom_jitter(aes(x = as.factor(shaft), y = microm3, color = microm3> 0.01), alpha = 0.7) + 
  facet_grid(. ~ gaba, labeller = label_both) + 
  scale_color_manual(values = c('black', 'red'))

#p2 + geom_boxplot(aes(x = as.factor(shaft), y = PSD95), notch = TRUE) + facet_grid(. ~ gaba, labeller = label_both)
#p2 + geom_boxplot(aes(x = as.factor(shaft), y = microm3), notch = TRUE) + facet_grid(. ~ gaba, labeller = label_both)

2.1 Synapses labeled as shaft gabergic and >0.01 micrometers^3

links
197 170 36
298 225 28
293 166 36
293 193 30
306 181 32