From 495486f53c6b0f9b86cacae7a54f7070b853294f Mon Sep 17 00:00:00 2001 From: Sebastian Seung Date: Tue, 19 Nov 2024 13:53:22 -0500 Subject: [PATCH] figures for FormVisionPaper --- FormVisionPaper/DataS2.jl | 62 + FormVisionPaper/FiguresNonspatial.jl | 1011 +++++++++++++++ FormVisionPaper/FiguresSpatial.jl | 1738 ++++++++++++++++++++++++++ 3 files changed, 2811 insertions(+) create mode 100644 FormVisionPaper/DataS2.jl create mode 100644 FormVisionPaper/FiguresNonspatial.jl create mode 100644 FormVisionPaper/FiguresSpatial.jl diff --git a/FormVisionPaper/DataS2.jl b/FormVisionPaper/DataS2.jl new file mode 100644 index 0000000..1f91ebd --- /dev/null +++ b/FormVisionPaper/DataS2.jl @@ -0,0 +1,62 @@ +# --- +# jupyter: +# jupytext: +# formats: ipynb,jl:percent +# text_representation: +# extension: .jl +# format_name: percent +# format_version: '1.3' +# jupytext_version: 1.16.2 +# kernelspec: +# display_name: Julia 1.10.6 +# language: julia +# name: julia-1.10 +# --- + +# %% [markdown] +# # Visualizations of the mapping of cells to hexagonal lattice +# Visualize cells of a given type with constant `p` or constant `q`. +# +# Creates clickable neuroglancer links. +# +# Used to generate visualizations in Supplementary Data 2, which demonstrate mapping of types to hexagonal lattice. + +# %% +using OpticLobe + +# %% +using Missings, MissingsAsFalse + +# %% +pmax, qmax = size(pq2column) + +# %% [markdown] +# ## visualize L1 + +# %% +celltype = "L1" # substitute your desired cell type here + +# %% +ids = type2ids(celltype) +coords = id2pq.(ids) + +# %% +@mfalse for p = 1:pmax + ng_hyper(ids[passmissing(getindex).(coords, 1) .== p], anchor = "p=$p") |> display +end + +# %% +@mfalse for q = 1:qmax + ng_hyper(ids[passmissing(getindex).(coords, 2) .== q], anchor = "q=$q") |> display +end + +# %% [markdown] +# overall, Hungarian assignment looks really good for finding correct locations +# +# big gap in p=15 +# +# 720575940605423532 in two columns? p=26 not right. should be p=25 +# +# anterior three cells of p=30 questionable +# +# p=31 anterior cell also seems out of place diff --git a/FormVisionPaper/FiguresNonspatial.jl b/FormVisionPaper/FiguresNonspatial.jl new file mode 100644 index 0000000..a24e543 --- /dev/null +++ b/FormVisionPaper/FiguresNonspatial.jl @@ -0,0 +1,1011 @@ +# -*- coding: utf-8 -*- +# --- +# jupyter: +# jupytext: +# formats: ipynb,jl:percent +# text_representation: +# extension: .jl +# format_name: percent +# format_version: '1.3' +# jupytext_version: 1.16.2 +# kernelspec: +# display_name: Julia 1.10.4 +# language: julia +# name: julia-1.10 +# --- + +# %% [markdown] +# # Form vision paper (nonspatial figures) +# For reproducing figures in Seung, [Predicting visual function by interpreting a neuronal wiring diagram](https://doi.org/10.1038/s41586-024-07953-5), Nature 634:113-123 (2024). +# +# This notebook is restricted to figures not involving space. + +# %% +using OpticLobe + +# %% +using NamedArrays, SparseArrays + +# %% +using Plots, Measures + +# %% +using MissingsAsFalse, Missings + +# %% +using StatsBase + +# %% +using StatsPlots + +# %% +using Distances, Clustering + +# %% +# because labels requires a row vector +Base.transpose(s::String) = s + +# %% [markdown] +# ## canonical ordering of types + +# %% +Dm3types = ["Dm3v", "Dm3p", "Dm3q"] +TmYtypes = ["TmY4", "TmY9q", "TmY9q⊥"] +formtypes = vcat(Dm3types, TmYtypes) + +# %% [markdown] +# ## cell numbers +# - Dm3 is more numerous than TmY +# - vertical is least numerous for Dm3 +# - horizontal is most numerous for TmY + +# %% +[Dm3types length.(type2ids.(Dm3types))] + +# %% +[TmYtypes length.(type2ids.(TmYtypes))] + +# %% +length(type2ids("Tm1")) + +# %% +length(type2ids("T2a")) + +# %% [markdown] +# ## Figure 1b-d. Visualizations of all Dm3v cells +# These show all cells of each type, rather than the subsets in the paper. + +# %% +for Dm3type in Dm3types + ng_hyper(Dm3type) |> display +end + +# %% [markdown] +# ## Figure 2a-c. Visualizations of all TmY4 and TmY9 cells +# These show all cells of each type, rather than the subsets in the paper. + +# %% +for TmYtype in TmYtypes + ng_hyper(TmYtype) |> display +end + +# %% [markdown] +# ## population-to-population connectivity +# two populations: +# - Dm3 includes trio of Dm3 types +# - TmY includes trio of TmY types + +# %% +Wpp = [sum(Wtt[a, b]) for a in [Dm3types, TmYtypes], b in [Dm3types, TmYtypes]] # 2x2 matrix of synapse numbers +Wpout = [sum(Wtc[a, :]) for a in [Dm3types, TmYtypes]] # outgoing synapses from two populations +Wpin = [sum(Wct[:, a]) for a in [Dm3types, TmYtypes]] # incoming synapses to two populations + +# %% +plot(heatmap(Wpp, title = "number"), + heatmap(Wpp./Wpin', title = "in fraction"), + heatmap(Wpp./Wpout, title = "out fraction"), + size = (1200, 250), + layout = (1, 3), + ticks = strings2ticks(["Dm3", "TmY"]), + yflip = true +) + +# %% [markdown] +# ## Figure 3a, type-to-type connectivity + +# %% +# include LC in postsynaptic types +formtypespost = vcat(Dm3types, TmYtypes, ["LC15", "LC10e"]) + +# %% +NamedArray(collect(Wtt[formtypes, formtypespost]), names = (formtypes, formtypespost)) + +# %% +default(; # Plots defaults + fontfamily="Helvetica", + label="", # only explicit legend entries + guidefontsize = 12, + tickfontsize = 12, + legendfontsize = 12, + linewidth = 2, + left_margin = 10mm, + bottom_margin = 0mm, + c = :hot, + xticks = strings2ticks(formtypes), yticks = strings2ticks(formtypes), + yflip = true +# xrotation = 60, + ) + +# %% +# version with cell types as tick labels +heatmap(infraction[formtypes, formtypespost], + yticks = strings2ticks(formtypes), + xticks = strings2ticks(formtypespost), + xrot = 30, legend = :none +) + +# %% +# version without tick labels for inclusion in figure file +heatmap(infraction[formtypes, formtypespost], xticks = strings2ticks(formtypespost), axis = ([], false), legend = :none) + +# %% +savefig("Dm3TmYPopulationInteractions.pdf") + +# %% [markdown] +# ## Figure 3b. identify connections included in cartoon wiring diagram + +# %% +NamedArray(collect(infraction[formtypes, formtypespost]), names = (formtypes, formtypespost)) |> showall + +# %% +# connections included in cartoon wiring diagram Figure 3b +collect(infraction[formtypes, formtypespost] .> 0.03) + +# %% +collect(infraction[formtypes, formtypespost] .> 0.03 .&& infraction[formtypes, formtypespost] .< 0.04) + +# %% [markdown] +# ## functions for tracing pathways backwards and forwards + +# %% +""" +find cells of type `celltype` presynaptic to cell `id`. +connections containing more than `thres` synapses +""" +function backwardtrace(id, celltype, thres = 4) + strongpre = findall(W[:, id2ind(id)] .> thres) + return @mfalse ind2id[strongpre[ind2type[strongpre] .== celltype]] +end + +""" +find cells of type `celltype` postsynaptic to cell `id`. +connections containing more than `thres` synapses +""" +function forwardtrace(id, celltype, thres = 4) + strongpost = findall(W[id2ind(id), :] .> thres) + return @mfalse ind2id[strongpost[ind2type[strongpost] .== celltype]] +end + +# %% [markdown] +# ## Figure 3c. Dm3q-Dm3p-TmY9q visualization + +# %% +first(type2ids("Dm3p"), 10) + +# %% +id = 720575940605138860 # Dm3p cell + +# %% +cellids = vcat(backwardtrace(id, "Dm3q", 5), id, forwardtrace(id, "TmY9q", 5)) + +# %% +ind2type[id2ind.(cellids)] + +# %% +ng_hyper(cellids) + +# %% +W[Name.(cellids), Name.(cellids)] |> collect + +# %% +[backwardtrace(id, "Tm1") for id in cellids] # presynaptic Tm1 cells + +# %% [markdown] +# ## Figure 3d. TmY4 visualization +# visualize TmY4-TmY4 pair reciprocally connected by more than thres synapses + +# %% +ids = type2ids("TmY4") +WintraTmY4 = W[Name.(ids), Name.(ids)] +thres = 4 +strongreciprocalpairs = findall(WintraTmY4 .> thres .&& WintraTmY4' .> thres) +println(length(strongreciprocalpairs), " reciprocal pairs connected by >", thres, " synapses") + +# %% +# first pair in list +i, j = collect(Tuple(strongreciprocalpairs[1])) +println("$(ids[i])→$(ids[j]): ", WintraTmY4[i, j], " synapses") +println("$(ids[j])→$(ids[i]): ", WintraTmY4[j, i], " synapses") + +# %% +ng_hyper(ids[[i, j]]) + +# %% [markdown] +# ## Figure 6a. LC15 visualization + +# %% +id = 720575940604229152 # LC15 cell + +# %% +ng_hyper(vcat([id], [backwardtrace(id, TmYtype) for TmYtype in TmYtypes]...)) + +# %% [markdown] +# ## Figure 6b. LC10e visualization + +# %% +id = 720575940606274592 # LC10e cell + +# %% +ng_hyper(vcat([id], [backwardtrace(id, TmYtype, 3) for TmYtype in TmYtypes]...)) + +# %% [markdown] +# ## Dm3p-Dm3q-TmY9q⊥ visualization + +# %% +@mfalse ids = type2ids("Dm3q") + +# %% +id = 720575940600623020 # Dm3q cell + +# %% +backwardtrace(id, "Dm3p", 3) + +# %% +forwardtrace(id, "TmY9q⊥", 5) + +# %% +ng_hyper(vcat([id], backwardtrace(id, "Dm3p", 3), forwardtrace(id, "TmY9q⊥", 5))) + +# %% +cellids = [720575940605123701, 720575940600623020, 720575940617250233] + +# %% +ind2type[id2ind.(cellids)] + +# %% +W[Name.(cellids), Name.(cellids)] |> collect + +# %% [markdown] +# ## Figure S7. disynaptic pathways + +# %% +# defined as numerous and receive L input + R7-8 +hexel = ["L1", "L2", "L3", "L4", "L5", "Tm1", "Tm2", "Tm9", "Mi1", "Mi4", "Mi9"] +length.(type2ids.(hexel)) + +# %% +nothexel = setdiff(intrinsictypes, hexel) + +# %% +Plots.reset_defaults() + +# %% +default(; # Plots defaults + fontfamily="Helvetica", + label="", # only explicit legend entries + guidefontsize = 12, + tickfontsize = 10, + legendfontsize = 10, + linewidth = 2, + left_margin = 0mm, + bottom_margin = 0mm, +# xrotation = 60, + ) + +# %% +h = [] +n = 12 + +for target in vcat(Dm3types, TmYtypes, "LC15", "LC10e") + scores = [scorepath([c1, c2, target]) for c1 in hexel, c2 in nothexel] + scores = NamedArray(scores, names = (hexel, nothexel)) +# println(target) +# showall(sum(scores, dims = 2)) +# println("\n") + scoresums = NamedArray(sum(scores.array, dims = 1)[:], nothexel) + scoremaxs = NamedArray(maximum(scores.array, dims = 1)[:], nothexel) + p = sortperm(scoresums, rev = true) + + push!(h, plot(100*[scoresums[p[1:n]], + scoremaxs[p[1:n]], scores["Tm1", p[1:n]]], + xticks = strings2ticks(nothexel[p[1:n]]), + label = ["sum" "max" "Tm1"], + annotations = (:E, Plots.text(target, "Helvetica", :right)), + permute = (:y, :x), + xlabel = startswith(target, "Dm3") ? "intermediary type" : "", + leftmargin = startswith(target, "Dm3") ? 5mm : 0mm, + xflip = true) + ) +end + +# %% +plot(h[[1 4 7 2 5 8 3 6]]..., + layout = grid(3, 3), + size = (800, 1000), + legend = :bottomright, +) + +# %% +savefig("DisynapticInFractions.pdf") + +# %% +type2nt[["TmY10", "TmY11", "Tm20", "Tm25", "Tm27", "Tm8a", "Tm16", "Tm7"]] + +# %% [markdown] +# ## list top partners (collapse) +# These are modified versions of `toppre` and `toppost` utility functions defined in `OpticLobe`. +# +# number of synapses per cell + +# %% +Av = A[:, visualtypes] + +# %% +# function toppost2(ids::Vector{Int64}, nresults = 15, reduce=false) +# if reduce +# pre = sort(NamedArray(mean(W[id2ind.(ids), :]*Av*collapse, dims=1)[:], collapsetypes), rev=true) +# else +# pre = sort(NamedArray(mean(W[id2ind.(ids), :]*Av, dims=1)[:], visualtypes), rev=true) +# end +# (names(pre)[1][1:nresults], round.(Int64, pre.array[1:nresults])) +# end + +# %% +function toppost2(ids::Vector{Int64}, nresults = 15, reduce=false) + if reduce + pre = sort(NamedArray(mean(W[id2ind.(ids), :]*Av*collapse, dims=1)[:], collapsetypes), rev=true) + else + pre = sort(NamedArray(mean(W[id2ind.(ids), :]*Av, dims=1)[:], visualtypes), rev=true) + end + (disallowmissing(names(pre, 1))[1:nresults], round.(Int64, pre.array[1:nresults])) +end + +# %% +function toppre2(ids::Vector{Int64}, nresults = 15, reduce=false) + if reduce + post = sort(NamedArray(mean(collapse'*Av'*W[:, id2ind.(ids)], dims=2)[:], collapsetypes), rev=true) + else + post = sort(NamedArray(mean(Av'*W[:, id2ind.(ids)], dims=2)[:], visualtypes), rev=true) + end + (disallowmissing(names(post, 1))[1:nresults], round.(Int64, post.array[1:nresults])) +end + +# %% [markdown] +# ## create matrix that collapses all Dm3 types to "Dm3", and all TmY4/TmY9 types to "TmY" + +# %% +using LinearAlgebra + +# %% +ntype = length(visualtypes) + +# %% +mask = NamedArray(trues(ntype), visualtypes) + +mask["TmY9q"] = false +mask["TmY9q⊥"] = false +mask["Dm3p"] = false +mask["Dm3q"] = false + +# %% +collapse = NamedArray(Matrix(I, ntype, ntype), names = (visualtypes, visualtypes)) + +collapse = collapse[:, mask] + +# %% +collapse["Dm3p", "Dm3v"] = 1 +collapse["Dm3q", "Dm3v"] = 1 +collapse["TmY9q", "TmY4"] = 1 +collapse["TmY9q⊥", "TmY4"] = 1 + +# %% +collapsetypes = visualtypes[mask] + +# %% +collapsetypes[collapsetypes .== "TmY4"] .= "TmY4/9" + +# %% +collapsetypes[collapsetypes .== "Dm3v"] .= "Dm3" + +# %% +collapse = NamedArray(collapse.array, names = (visualtypes, collapsetypes)) + +# %% +visualtypes[collapse[:, "TmY4/9"] .> 0] + +# %% [markdown] +# ## Figure S1. Dm3 inputs and outputs + +# %% +@mfalse Dm3ids = vcat([ind2id[ind2type .== celltype] for celltype in ["Dm3v", "Dm3p", "Dm3q"]]...) + +# top presynaptic and postsynaptic partners of Dm3 +npartners = 20 +@mfalse Dm3prestrings = toppre2(Dm3ids, npartners)[1] +@mfalse Dm3poststrings = toppost2(Dm3ids, npartners)[1] + +Dm3strings = ["Dm3v", "Dm3p", "Dm3q"] +Dm3inputs = collect(A[:, Dm3prestrings]'*W*A[:, Dm3strings]) +Dm3inputs = Dm3inputs./collect(sum(A[:, Dm3strings], dims=1)) + +Dm3outputs = collect(A[:, Dm3strings]'*W*A[:, Dm3poststrings])' +Dm3outputs = Dm3outputs./collect(sum(A[:, Dm3strings], dims=1)) + +# %% +sum(Dm3outputs, dims = 2) + +# %% +@mfalse Dm3prestringsreduced = toppre2(Dm3ids, npartners, true)[1] +@mfalse Dm3poststringsreduced = toppost2(Dm3ids, npartners, true)[1] + +Dm3inputsreduced = collect((Av*collapse)[:, Dm3prestringsreduced]'*W*Av[:, Dm3strings]) +Dm3inputsreduced = Dm3inputsreduced./collect(sum(Av[:, Dm3strings], dims=1)) + +Dm3outputsreduced = collect(Av[:, Dm3strings]'*W*(Av*collapse)[:, Dm3poststringsreduced])' +Dm3outputsreduced = Dm3outputsreduced./collect(sum(Av[:, Dm3strings], dims=1)) + +# %% [markdown] +# ### Plots.jl defaults + +# %% +Plots.reset_defaults() + +# %% +default(; # Plots defaults + fontfamily="Helvetica", + label="", # only explicit legend entries + guidefontsize = 12, + tickfontsize = 12, + legendfontsize = 12, + linewidth = 2, + left_margin = 10mm, + bottom_margin = 0mm, +# xrotation = 60, + ) + +# %% +h = [plot(data, xticks = strings2ticks(partnertypes), + labels = transpose(Dm3strings), + xlabel = "$direction cell type", + ylabel = "# $direction synapses", + permute = (:x, :y), xflip = true, + ylim = (0, Inf), widen = true + ) +for (data, partnertypes, direction) in zip( + [Dm3inputs, Dm3outputs, Dm3inputsreduced, Dm3outputsreduced], + [Dm3prestrings, Dm3poststrings, Dm3prestringsreduced, Dm3poststringsreduced], + ["input", "output", "input", "output"] + ) + ] +plot(h..., + layout = grid(2, 2, heights=[0.5, 0.5]), + size = (800, 1000), + legend = :bottomright +) + +# %% +savefig("Dm3InOut.pdf") + +# %% [markdown] +# ## dendrogram for Dm3 and TmY types + +# %% +embeddings = vcat(infraction[visualtypes, formtypes], outfraction[formtypes, visualtypes]') + +# %% +dist = pairwise(jaccard, embeddings) +c = hclust(dist, linkage = :average) +plot(c, xticks = strings2ticks(formtypes[c.order]), size = (1200, 300), xrot = 30, bottom_margin = 7mm) + +# %% +clist = [] +for embeddings in [vcat(infraction[celltypes, formtypes], outfraction[formtypes, celltypes]') for celltypes in [intrinsictypes, visualtypes, alltypes]] + dist = pairwise(jaccard, embeddings) + push!(clist, hclust(dist, linkage = :average)) +end + +# %% +h = [plot(c, xticks = strings2ticks(formtypes[c.order]), size = (1200, 300), xrot = 30, bottom_margin = 7mm) for c in clist] + +# %% +plot(h..., layout = (1, 3)) + +# %% [markdown] +# ## dendrogram for Dm3, TmY, and other types connected with them + +# %% +involved = vcat(formtypes,[ "L3", "L5", "Tm1", "Tm2", "Tm9", "Mi2", "Mi4", "Mi9", "T2a", "T3", "Tm20", "Tm5f", "Dm15", "Dm18", "Tm21", "Tm25", "Tm7", "Tm8a", "Tm16", "TmY10", "TmY11", "Y3", "T5a", "T4a", "T5c", "T5d", "T4c", "T4d"]) + +# %% +embeddings = vcat(infraction[visualtypes, involved], outfraction[involved, visualtypes]') + +# %% +dist = pairwise(jaccard, embeddings) +c = hclust(dist, linkage = :average) +plot(c, xticks = strings2ticks(involved[c.order]), size = (1200, 300), xrot = 30, bottom_margin = 7mm) + +# %% [markdown] +# ## Figure S2. Tm1 inputs and outputs + +# %% +# top presynaptic and postsynaptic partners of Tm1 +#npartners = 25 +npartners = 26 +@mfalse Tm1prestrings = toppre2(ind2id[ind2type .== "Tm1"], npartners)[1] +@mfalse Tm1poststrings = toppost2(ind2id[ind2type .== "Tm1"], npartners)[1] + +Tm1inputs = collect(A[:, Tm1prestrings]'*W*A[:, ["Tm1"]]) +Tm1inputs = Tm1inputs./collect(sum(A[:, ["Tm1"]], dims=1)) +Tm1outputs = collect(A[:, ["Tm1"]]'*W*A[:, Tm1poststrings])' +Tm1outputs = Tm1outputs./collect(sum(A[:, ["Tm1"]], dims=1)) + +# %% +# top presynaptic and postsynaptic partners of Tm1 +@mfalse Tm1prestringsreduced = toppre2(ind2id[ind2type .== "Tm1"], npartners, true)[1] +@mfalse Tm1poststringsreduced = toppost2(ind2id[ind2type .== "Tm1"], npartners, true)[1] + +Tm1inputsreduced = collect((Av*collapse)[:, Tm1prestringsreduced]'*W*Av[:, ["Tm1"]]) +Tm1inputsreduced = Tm1inputsreduced./collect(sum(Av[:, ["Tm1"]], dims=1)) +Tm1outputsreduced = collect(Av[:, ["Tm1"]]'*W*Av*collapse[:, Tm1poststringsreduced])[:] +Tm1outputsreduced = Tm1outputsreduced./collect(sum(Av[:, ["Tm1"]], dims=1)) + +# %% +h = [plot(data, xticks = strings2ticks(partnertypes), + labels = "Tm1", + xlabel = "$direction cell type", + ylabel = "# $direction synapses", + permute = (:x, :y), xflip = true, + ylim = (0, Inf), widen = true) +for (data, partnertypes, direction) in zip( + [Tm1inputs, Tm1outputs, Tm1inputsreduced, Tm1outputsreduced], + [Tm1prestrings, Tm1poststrings, Tm1prestringsreduced, Tm1poststringsreduced], + ["input", "output", "input", "output"] + ) + ] +plot(h..., + layout = grid(2, 2, heights=[0.5, 0.5]), + size = (800, 1000), + legend = :bottomright +) + +# %% +savefig("Tm1InOut.pdf") + +# %% [markdown] +# ## Figure S4. TmY4 and TmY9 four panels + +# %% [markdown] +# ### TmY9 + +# %% +TmY9strings = ["TmY9q", "TmY9q⊥"] +@mfalse TmY9ids = vcat([ind2id[ind2type .== celltype] for celltype in TmY9strings]...) + +# top presynaptic and postsynaptic partners of TmY +npartners = 30 +@mfalse TmY9prestrings = toppre2(TmY9ids, npartners)[1] + +@mfalse TmY9poststrings = toppost2(TmY9ids, npartners)[1] + +TmY9inputs = collect(A[:, TmY9prestrings]'*W*A[:, TmY9strings]) +TmY9inputs = TmY9inputs./collect(sum(A[:, TmY9strings], dims=1)) +TmY9outputs = collect(A[:, TmY9strings]'*W*A[:, TmY9poststrings])' +TmY9outputs = TmY9outputs./collect(sum(A[:, TmY9strings], dims=1)) + +# %% [markdown] +# ### TmY4 + +# %% +# top presynaptic and postsynaptic partners of TmY4 +npartners = 20 +@mfalse TmY4prestrings = toppre2(ind2id[ind2type .== "TmY4"], npartners)[1] +@mfalse TmY4poststrings = toppost2(ind2id[ind2type .== "TmY4"], npartners)[1] + +TmY4inputs = collect(Av[:, TmY4prestrings]'*W*Av[:, ["TmY4"]]) +TmY4inputs = TmY4inputs./collect(sum(Av[:, ["TmY4"]], dims=1)) +TmY4outputs = collect(Av[:, ["TmY4"]]'*W*Av[:, TmY4poststrings])' +TmY4outputs = TmY4outputs./collect(sum(Av[:, ["TmY4"]], dims=1)) + +# %% +h = [plot(data, xticks = strings2ticks(partnertypes), + labels = transpose(visualtypes), + xlabel = "$direction cell type", + ylabel = "# $direction synapses", + permute = (:x, :y), xflip = true, + ylim = (0, Inf), widen = true) +for (data, partnertypes, visualtypes, direction) in zip( + [TmY4inputs, TmY4outputs, TmY9inputs, TmY9outputs], + [TmY4prestrings, TmY4poststrings, TmY9prestrings, TmY9poststrings], + [ "TmY4", "TmY4", TmY9strings, TmY9strings], + ["input", "output", "input", "output"] + ) + ] +plot(h..., + layout = grid(2, 2, heights=[0.4, 0.6]), + size = (800, 1000), + legend = :bottomright +) + +# %% +savefig("TmY4TmY9InOut.pdf") + +# %% [markdown] +# ## Figure S5. TmY4 and TmY9 four panels, reduced + +# %% [markdown] +# ### TmY9 + +# %% +TmY9strings = ["TmY9q", "TmY9q⊥"] +@mfalse TmY9ids = vcat([ind2id[ind2type .== celltype] for celltype in TmY9strings]...) + +# top presynaptic and postsynaptic partners of TmY +npartners = 30 +@mfalse TmY9prestringsreduced = toppre2(TmY9ids, npartners, true)[1] + +@mfalse TmY9poststringsreduced = toppost2(TmY9ids, npartners, true)[1] + +TmY9inputsreduced = collect((Av*collapse)[:, TmY9prestringsreduced]'*W*Av[:, TmY9strings]) +TmY9inputsreduced = TmY9inputsreduced./collect(sum(Av[:, TmY9strings], dims=1)) +TmY9outputsreduced = collect(Av[:, TmY9strings]'*W*(Av*collapse)[:, TmY9poststringsreduced])' +TmY9outputsreduced = TmY9outputsreduced./collect(sum(Av[:, TmY9strings], dims=1)) + +# %% [markdown] +# ### TmY4 + +# %% +# top presynaptic and postsynaptic partners of TmY4 +npartners = 20 +@mfalse TmY4prestringsreduced = toppre2(ind2id[ind2type .== "TmY4"], npartners, true)[1] +@mfalse TmY4poststringsreduced = toppost2(ind2id[ind2type .== "TmY4"], npartners, true)[1] + +TmY4inputsreduced = collect((Av*collapse)[:, TmY4prestringsreduced]'*W*Av[:, ["TmY4"]]) +TmY4inputsreduced = TmY4inputsreduced./collect(sum(Av[:, ["TmY4"]], dims=1)) +TmY4outputsreduced = collect(Av[:, ["TmY4"]]'*W*(Av*collapse)[:, TmY4poststringsreduced])' +TmY4outputsreduced = TmY4outputsreduced./collect(sum(Av[:, ["TmY4"]], dims=1)) + +# %% +h = [plot(data, xticks = strings2ticks(partnertypes), + labels = transpose(visualtypes), + xlabel = "$direction cell type", + ylabel = "# $direction synapses", + permute = (:x, :y), xflip = true, + ylim = (0, Inf), widen = true) +for (data, partnertypes, visualtypes, direction) in zip( + [TmY4inputsreduced, TmY4outputsreduced, TmY9inputsreduced, TmY9outputsreduced], + [TmY4prestringsreduced, TmY4poststringsreduced, TmY9prestringsreduced, TmY9poststringsreduced], + [ "TmY4", "TmY4", TmY9strings, TmY9strings], + ["input", "output", "input", "output"] + ) + ] +plot(h..., + layout = grid(2, 2, heights=[0.4, 0.6]), + size = (800, 1000), + legend = :bottomright +) + +# %% +savefig("TmY4TmY9InOutReduced.pdf") + +# %% [markdown] +# ## Y3 and TmY5a + +# %% +# top presynaptic and postsynaptic partners of Y3 +npartners = 25 +@mfalse Y3prestrings = toppre2(ind2id[ind2type .== "Y3"], npartners)[1] +@mfalse Y3poststrings = toppost2(ind2id[ind2type .== "Y3"], npartners)[1] + +Y3inputs = collect(Av[:, Y3prestrings]'*W*Av[:, ["Y3"]]) +Y3inputs = Y3inputs./collect(sum(Av[:, ["Y3"]], dims=1)) +Y3outputs = collect(Av[:, ["Y3"]]'*W*Av[:, Y3poststrings])' +Y3outputs = Y3outputs./collect(sum(Av[:, ["Y3"]], dims=1)) + +# %% +# top presynaptic and postsynaptic partners of TmY5a +npartners = 25 +@mfalse TmY5aprestrings = toppre2(ind2id[ind2type .== "TmY5a"], npartners)[1] +@mfalse TmY5apoststrings = toppost2(ind2id[ind2type .== "TmY5a"], npartners)[1] + +TmY5ainputs = collect(Av[:, TmY5aprestrings]'*W*Av[:, ["TmY5a"]]) +TmY5ainputs = TmY5ainputs./collect(sum(Av[:, ["TmY5a"]], dims=1)) +TmY5aoutputs = collect(Av[:, ["TmY5a"]]'*W*Av[:, TmY5apoststrings])' +TmY5aoutputs = TmY5aoutputs./collect(sum(Av[:, ["TmY5a"]], dims=1)) + +# %% +h = [plot(data, xticks = strings2ticks(partnertypes), + labels = transpose(visualtypes), + xlabel = "$direction cell type", + ylabel = "# $direction synapses", + permute = (:x, :y), xflip = true, + ylim = (0, Inf), widen = true) +for (data, partnertypes, visualtypes, direction) in zip( + [Y3inputs, Y3outputs, TmY5ainputs, TmY5aoutputs], + [Y3prestrings, Y3poststrings, TmY5aprestrings, TmY5apoststrings], + [ "Y3", "Y3", "TmY5a", "TmY5a"], + ["input", "output", "input", "output"] + ) + ] +plot(h..., + layout = grid(2, 2, heights=[0.5, 0.5]), + size = (800, 1000), + legend = :bottomright +) + +# %% +savefig("Y3TmY5aInOut.pdf") + +# %% [markdown] +# ## T2a + +# %% +T2astrings = ["T2a"] +@mfalse T2aids = vcat([ind2id[ind2type .== celltype] for celltype in T2astrings]...) + +# top presynaptic and postsynaptic partners of TmY +npartners = 20 +@mfalse T2aprestrings = toppre2(T2aids, npartners)[1] + +@mfalse T2apoststrings = toppost2(T2aids, npartners)[1] + +T2ainputs = collect(Av[:, T2aprestrings]'*W*Av[:, T2astrings]) +T2ainputs = T2ainputs./collect(sum(Av[:, T2astrings], dims=1)) +T2aoutputs = collect(Av[:, T2astrings]'*W*Av[:, T2apoststrings])' +T2aoutputs = T2aoutputs./collect(sum(Av[:, T2astrings], dims=1)) + +# %% +plot( + plot(T2ainputs, + xticks = (1:length(T2aprestrings), T2aprestrings), + labels = "T2a", + xlabel = "input cell type", + ylabel = "# of input synapses", + permute = (:x, :y), xflip = true, +# left_margin = 10mm, + ), +plot(T2aoutputs, + xticks = (1:length(T2apoststrings), T2apoststrings), + labels = "T2a", + xlabel = "output cell type", + ylabel = "# of output synapses", + permute = (:x, :y), xflip = true, +# left_margin = 5mm, +), + layout = grid(1, 2), + bottom_margin = 5mm, + size = (800, 500), + legend = :bottomright + ) + +# %% +savefig("T2aInOut.pdf") + +# %% [markdown] +# ## Figure S8. Y3 and T2a + +# %% +h = [plot(data, xticks = strings2ticks(partnertypes), + labels = transpose(visualtypes), + xlabel = "$direction cell type", + ylabel = "# $direction synapses", + permute = (:x, :y), xflip = true, + ylim = (0, Inf), widen = true) +for (data, partnertypes, visualtypes, direction) in zip( + [Y3inputs, Y3outputs, T2ainputs, T2aoutputs], + [Y3prestrings, Y3poststrings, T2aprestrings, T2apoststrings], + [ "Y3", "Y3", "T2a", "T2a"], + ["input", "output", "input", "output"] + ) + ] +plot(h..., + layout = grid(2, 2, heights=[0.5, 0.5]), + size = (800, 1000), + legend = :bottomright +) + +# %% +savefig("Y3T2aInOut.pdf") + +# %% [markdown] +# ## Figure S9. TmY10 TmY11 + +# %% +# top presynaptic and postsynaptic partners of TmY10 +npartners = 20 +@mfalse TmY10prestrings = toppre2(ind2id[ind2type .== "TmY10"], npartners)[1] +@mfalse TmY10poststrings = toppost2(ind2id[ind2type .== "TmY10"], npartners)[1] + +TmY10inputs = collect(Av[:, TmY10prestrings]'*W*Av[:, ["TmY10"]]) +TmY10inputs = TmY10inputs./collect(sum(Av[:, ["TmY10"]], dims=1)) +TmY10outputs = collect(Av[:, ["TmY10"]]'*W*Av[:, TmY10poststrings])' +TmY10outputs = TmY10outputs./collect(sum(Av[:, ["TmY10"]], dims=1)) + +# %% +# top presynaptic and postsynaptic partners of TmY10 +npartners = 20 +@mfalse TmY11prestrings = toppre2(ind2id[ind2type .== "TmY11"], npartners)[1] +@mfalse TmY11poststrings = toppost2(ind2id[ind2type .== "TmY11"], npartners)[1] + +TmY11inputs = collect(Av[:, TmY11prestrings]'*W*Av[:, ["TmY11"]]) +TmY11inputs = TmY11inputs./collect(sum(Av[:, ["TmY11"]], dims=1)) +TmY11outputs = collect(Av[:, ["TmY11"]]'*W*Av[:, TmY11poststrings])' +TmY11outputs = TmY11outputs./collect(sum(Av[:, ["TmY11"]], dims=1)) + +# %% +h = [plot(data, xticks = strings2ticks(partnertypes), + labels = transpose(visualtypes), + xlabel = "$direction cell type", + ylabel = "# $direction synapses", + permute = (:x, :y), xflip = true, + ylim = (0, Inf), widen = true) +for (data, partnertypes, visualtypes, direction) in zip( + [TmY10inputs, TmY10outputs, TmY11inputs, TmY11outputs], + [TmY10prestrings, TmY10poststrings, TmY11prestrings, TmY11poststrings], + [ "TmY10", "TmY10", "TmY11", "TmY11"], + ["input", "output", "input", "output"] + ) + ] +plot(h..., + layout = grid(2, 2, heights=[0.5, 0.5]), + size = (750, 1000), + legend = :bottomright +) + +# %% +savefig("TmY10TmY11InOut.pdf") + +# %% [markdown] +# ## Figure S10. LC10e and LC15 + +# %% +### new clustering +LC10ed = [720575940603823136, 720575940604935089, 720575940610075029, 720575940614496735, 720575940614746178, 720575940616179050, 720575940618360741, 720575940620416176, 720575940621175979, 720575940622744414, 720575940622805854, 720575940623979815, 720575940624787504, 720575940625143396, 720575940625981080, 720575940628023296, 720575940628359783, 720575940628615752, 720575940629385943, 720575940629398492, 720575940629958764, 720575940631576546, 720575940633027615, 720575940637401321, 720575940637706574, 720575940637731955, 720575940638701375, 720575940644393070] +LC10ev = [720575940604125920, 720575940606274592, 720575940608977988, 720575940609733387, 720575940609890402, 720575940612305250, 720575940612332401, 720575940614897662, 720575940621233157, 720575940621602047, 720575940622486057, 720575940622602992, 720575940623601193, 720575940624050407, 720575940624595893, 720575940628632064, 720575940629076042, 720575940629756556, 720575940630438651, 720575940631305467, 720575940632034385, 720575940634088602, 720575940640974555] + +# %% +# top presynaptic and postsynaptic partners of LC10e +npartners = 25 +@mfalse LC10eprestrings = toppre2(vcat(LC10ed, LC10ev), npartners)[1] +@mfalse LC10epoststrings = toppost2(vcat(LC10ed, LC10ev), npartners)[1] + +LC10evinputs = mean(Av[:, LC10eprestrings]'*W[:, id2ind.(LC10ev)], dims=2)[:] +LC10evoutputs = mean(W[id2ind.(LC10ev), :]*Av[:, LC10epoststrings], dims=1)[:] + +LC10edinputs = mean(Av[:, LC10eprestrings]'*W[:, id2ind.(LC10ed)], dims=2)[:] +LC10edoutputs = mean(W[id2ind.(LC10ed), :]*Av[:, LC10epoststrings], dims=1)[:] + +# %% +# top presynaptic and postsynaptic partners of LC15 +npartners = 25 +@mfalse LC15prestrings = toppre2(ind2id[ind2type .== "LC15"], npartners)[1] +@mfalse LC15poststrings = toppost2(ind2id[ind2type .== "LC15"], npartners)[1] + +LC15inputs = collect(Av[:, LC15prestrings]'*W*Av[:, ["LC15"]]) +LC15inputs = LC15inputs./collect(sum(Av[:, ["LC15"]], dims=1)) +LC15outputs = collect(Av[:, ["LC15"]]'*W*Av[:, LC15poststrings])' +LC15outputs = LC15outputs./collect(sum(Av[:, ["LC15"]], dims=1)) + +# %% +h = [plot(data, xticks = strings2ticks(partnertypes), + labels = transpose(visualtypes), + xlabel = "$direction cell type", + ylabel = "# $direction synapses", + permute = (:x, :y), xflip = true, + ylim = (0, Inf), widen = true) +for (data, partnertypes, visualtypes, direction) in zip( + [[LC10edinputs, LC10evinputs], [LC10edoutputs, LC10evoutputs], LC15inputs, LC15outputs], + [LC10eprestrings, LC10epoststrings, LC15prestrings, LC15poststrings], + [ ["LC10e dorsal", "LC10e ventral"], ["LC10e dorsal", "LC10e ventral"], "LC15", "LC15"], + ["input", "output", "input", "output"] + ) + ] +plot(h..., + layout = grid(2, 2, heights=[0.5, 0.5]), + size = (800, 1000), + legend = :bottomright +) + +# %% +savefig("LC10eLC15InOut.svg") + +# %% [markdown] +# ## search for VPNs with high infraction from TmY + +# %% +thres = 0.045 + +# %% +abovethresholdtypes = unique(visualtypes[sum(infraction[TmYtypes, visualtypes], dims = 1)[:] .> 0.08]) + +# %% +upstreamTmY = NamedArray(sum(infraction[TmYtypes, visualtypes], dims = 1)[:].array, visualtypes) + +# %% +vpntypes = vcat(filter(startswith("LC"), visualtypes), filter(startswith("LT"), visualtypes)) + +# %% +rankedlist = sort(upstreamTmY[vpntypes], rev=true) + +# %% +for celltype in names(rankedlist)[1][1:20] + @mfalse println(sum(ind2type .== celltype)) + toppre(celltype) +end + +# %% +toppre("LTe57") + +# %% +toppre("LTe37") + +# %% +toppre("LC45") + +# %% +toppre("LT84") + +# %% +toppre("LT86") + +# %% +toppre("LC46") + +# %% +toppre("LTe35") + +# %% +first(sort(infraction["TmY9q", :], rev = true), 20) + +# %% +toppre("LTe63") + +# %% +toppre("TmY9q", nresults = 30) + +# %% +toppre("LTe08") + +# %% +toppre("CB0472") + +# %% +toppre("LTe28") + +# %% +toppre("LTe19") + +# %% +toppre("LTe45") + +# %% +toppre("LC41") + +# %% +toppre("LC10e") + +# %% +toppre("LT80") + +# %% +toppre("LT72") + +# %% +first(sort(infraction["TmY9q⊥", :], rev = true), 20) + +# %% +toppre("LTe08") + +# %% +for id in type2ids("LCe08") + toppre([id]) +end diff --git a/FormVisionPaper/FiguresSpatial.jl b/FormVisionPaper/FiguresSpatial.jl new file mode 100644 index 0000000..07a28e3 --- /dev/null +++ b/FormVisionPaper/FiguresSpatial.jl @@ -0,0 +1,1738 @@ +# -*- coding: utf-8 -*- +# --- +# jupyter: +# jupytext: +# formats: ipynb,jl:percent +# text_representation: +# extension: .jl +# format_name: percent +# format_version: '1.3' +# jupytext_version: 1.16.2 +# kernelspec: +# display_name: Julia 1.10.6 +# language: julia +# name: julia-1.10 +# --- + +# %% [markdown] +# # Form vision paper (spatial figures) +# +# For reproducing figures in Seung, [Predicting visual function by interpreting a neuronal wiring diagram](https://doi.org/10.1038/s41586-024-07953-5), Nature 634:113-123 (2024). +# +# This notebook is restricted to figures about spatial analysis. + +# %% +using OpticLobe + +# %% +using Missings, MissingsAsFalse + +# %% +using StatsBase +using StatsPlots # boxplots +using Measures + +# %% +using NamedArrays, SparseArrays, LinearAlgebra + +# %% +using OrderedCollections + +# %% +#using LaTeXStrings, MathTeXEngine + +# %% +using ProgressMeter + +# %% +using Printf + +# %% +using Luxor + +# %% +using ColorSchemes + +# %% [markdown] +# ## cell types + +# %% +Dm3types = ["Dm3v", "Dm3p", "Dm3q"] +TmYtypes = ["TmY4", "TmY9q", "TmY9q⊥"] + +# %% +hexel = ["L1", "L2", "L3", "L4", "L5", "Tm1", "Tm2", "Tm9", "Mi1", "Mi4", "Mi9"] +nothexel = setdiff(intrinsictypes, hexel) + +# %% [markdown] +# ## helper functions + +# %% +""" +for orientation boxplots +""" +function wraparound(theta; center = 0) + if theta > center + pi/2 + return theta - pi + end + if theta < center - pi/2 + return theta + pi + end + return theta +end + +# %% +# Deprecated because LaTeX fonts don't play well with Illustrator +# had to give up on superscript ⊥ +#function type2label(s::String) +# return latexstring(replace(s, "⊥" => "^⊥")) +#end + +#function path2label(v::Vector{String}) +# return latexstring(replace(join(v, "→"), "⊥" => "^⊥")) +#end + +# %% +function type2label(s::String) + return s +end + +function path2label(v::Vector{String}) + return join(v, "→") +end + +# %% +# """ +# `im` - square matrix of size m = 2k + 1 +# output is always vector of length m +# locations are: +# -k:k for cardinal orientations +# (-k:k)*sqrt(3)/2 for orthogonal orientations +# """ +# function hexproject(im::Matrix, axis) +# @assert size(im, 1) == size(im, 2) +# m = size(im, 1) +# @assert isodd(m) +# w = [0.5, 1, 0.5] +# k = m ÷ 2 +# P = zeros(Int32, m, m) .+ collect(-k:k) +# Q = zeros(Int32, m, m) .+ collect(-k:k)' +# # cardinal orientations, lattice spacing +# if axis == :v # P + Q runs from -2k:2k, output will run from -k:k +# proj = counts(P + Q, weights(im)) +# proj = dropdims(conv(proj[:, newaxis, :], w[:, newaxis, newaxis], pad = 1, stride = 2), dims = 2) +# return proj +# elseif axis == :p # P + Q runs from -3k:3k. crop to -2k:2k +# proj = counts(2*P - Q, -2*k:2*k, weights(im)) +# proj = dropdims(conv(proj[:, newaxis, :], w[:, newaxis, newaxis], pad = 1, stride = 2), dims = 2) +# return proj +# elseif axis == :q +# proj = counts(2*Q - P, -2*k:2*k, weights(im)) +# proj = dropdims(conv(proj[:, newaxis, :], w[:, newaxis, newaxis], pad = 1, stride = 2), dims = 2) +# return proj +# # orthogonal orientations, sqrt(3)/2 times lattice spacing +# elseif axis == :h +# proj = counts(Q - P, -k:k, weights(im)) +# return proj +# elseif axis == :p⊥ +# proj = counts(Q, weights(im)) +# return proj +# elseif axis == :q⊥ +# proj = counts(P, weights(im)) +# return proj +# end +# end + +# %% [markdown] +# ## add Dm3 and TmY locations to `id2pq` + +# %% +# missing means that there is zero connectivity +for posttype in vcat(Dm3types, TmYtypes, "Y3") + rfs = preimage("Tm1", posttype); + @mfalse for (i, id) in enumerate(ind2id[ind2type .== posttype]) + center = findcenter(float(rfs[i]), seven) + if ~ismissing(center) + id2pq[id] = center + end + end +end + +# %% [markdown] +# ## Fig 1 Tm1-Dm3 connectivity maps + +# %% [markdown] +# ### Fig 1bcd average Tm1-Dm3 maps + +# %% +#Drawing(800, 800, "Tm1Dm3Averages.pdf") +Drawing(400, 450) +origin() +kernels = typetriad("Tm1", Dm3types, radius = 3, hexelsize = 16, normalize = :common) +finish() +preview() + +# %% [markdown] +# ### Fig 1f Tm1-Dm3v example cells on eye maps + +# %% +rfs = preimage("Tm1", "Dm3v"); +maximum.(rfs[6:8]) + +# %% +#Drawing(600, 700, "Tm1Dm3vExamples.pdf") +Drawing(600, 700) +origin() +eyetriad("Tm1", type2ids("Dm3v")[6:8]) +finish() +preview() + +# %% [markdown] +# ### Fig 1g Tm1-Dm3v example cells, aligned and cropped + +# %% +# Drawing(450, 500, "Tm1Dm3vExamplesAlignedCropped.pdf") +Drawing(450, 500) +origin() +celltriad("Tm1", type2ids("Dm3v")[6:8], radius = 5, hexelsize = 12, ellipse = true) +finish() +preview() + +# %% [markdown] +# ### Fig 1h Tm1-Dm3 ellipses + +# %% +#Drawing(600, 600, "Dm3Ellipses.pdf") +Drawing(600, 600) +origin() +scale(50) +setline(0.5) + +for (celltype, ellipsecolor) in zip(Dm3types, ["red", "green", "blue"]) + setcolor(ellipsecolor) + ellipses = ellipsesummary.(preimage("Tm1", celltype)) + for e in sample(ellipses, 50, replace=false) + if !ismissing(e) + Luxor.rotate(-pi/6 + e.theta) + Luxor.ellipse(0, 0, e.minor, e.major, :stroke) # sqrt(3) for conversion to units of latticeconstant + Luxor.rotate(pi/6 - e.theta) + end + end +end +sethue("black") +setline(4) +drawpqaxes(1) +finish() +preview() + +# %% [markdown] +# ### Fig 1ij Tm1-Dm3 angles and aspect ratio boxplots + +# %% +default(; # Plots defaults + fontfamily="Arial Unicode", + label="", # only explicit legend entries + guidefontsize = 12, + tickfontsize = 12, + legendfontsize = 10, +# linewidth = 2, +# ylim = (0, Inf), + widen = true, +# xrotation = 30 + ) + +# %% +ellipses = [ellipsesummary.(preimage("Tm1", Dm3type)) for Dm3type in Dm3types] +aspectratios = [[e.major/e.minor for e in skipmissing(ellipse)] for ellipse in ellipses] +angles = [[(wraparound(e.theta - pi/6)) for e in skipmissing(ellipse)] for ellipse in ellipses] +diamsmajor = [[e.major for e in skipmissing(ellipse)] for ellipse in ellipses] +diamsminor = [[e.minor for e in skipmissing(ellipse)] for ellipse in ellipses] +#diamsmajor = [[passmissing(x -> x.major)(e) for e in ellipse] for ellipse in ellipses] # need to skip because quantiles don't work with missing values +#diamsminor = [[passmissing(x -> x.major)(e) for e in ellipse] for ellipse in ellipses] +Dm3ellipses = OrderedDict(Dm3types .=> ellipses) +Dm3aspectratios = OrderedDict(Dm3types .=> aspectratios) +Dm3angles = OrderedDict(Dm3types .=> angles) +celltypes = vcat([fill(i, length(v)) for (i, v) in enumerate(values(Dm3aspectratios))]...) + +haspect = boxplot(celltypes, vcat(values(Dm3aspectratios)...), ylabel = "aspect ratio", + xtick = strings2ticks(Dm3types), ylim = (1, Inf) +) + +hangle = boxplot(celltypes, vcat(values(Dm3angles)...)*180/pi, ylabel = "orientation (deg)", + xtick = strings2ticks(Dm3types), yticks = [-60, 0, 60] +) + +plot(hangle, haspect, size = (800, 300), left_margin = 5Measures.mm, bottom_margin = 3Measures.mm) + +# %% +savefig("Tm1Dm3anglesaspects.pdf") + +# %% [markdown] +# ### Fig 1klm Dm3 projections + +# %% +paths = [["Tm1", Dm3type] for Dm3type in Dm3types] +rfs = inmaps.(tracetypes.(paths)) +radius = 4 +rfscrop = [passmissing(crop).(rf.array, id2pq.(type2ids(Dm3type)), radius, 1) for (rf, Dm3type) in zip(rfs, Dm3types)]; + +# %% +long = [hexproject.(skipmissing(ims), ax) for (ims, ax) in zip(rfscrop, [:v, :p, :q])] +trans = [hexproject.(skipmissing(ims), ax) for (ims, ax) in zip(rfscrop, [:h, :p⊥, :q⊥])]; + +# %% +Plots.reset_defaults() + +# %% +default(; # Plots defaults + fontfamily="Arial Unicode", + label="", # only explicit legend entries + guidefontsize = 12, + tickfontsize = 12, + legendfontsize = 10, +# linewidth = 2, +# ylim = (0, Inf), + widen = true, +# xrotation = 30 + ) + +# %% +h = [groupedbar(-radius:radius, [a c], yerror = [b d], + linewidth = 1, + label = ["longitudinal" "transverse"], + ylabel = "synapses", + xlabel = "displacement", + annotations = (-3, 25, Plots.text(e, "Helvetica")), + ylim = (0, 30) +) for (a, b, c, d, e) in zip(mean.(long), std.(long), mean.(trans), std.(trans), Dm3types)] + +# %% +plot(h..., layout = (1, 3), size = (1000, 200), bottom_margin = 10Measures.mm, left_margin = 5Measures.mm) + +# %% +savefig("Dm3projections.pdf") + +# %% [markdown] +# ## Fig 2 Tm1-TmY connectivity maps + +# %% [markdown] +# ### Fig 2def average Tm1-TmY maps + +# %% +#Drawing(800, 800, "Tm1TmYAverages.pdf") +Drawing(400, 450) +origin() +kernels = typetriad("Tm1", TmYtypes, radius = 3, hexelsize = 16, normalize = :common) +finish() +preview() + +# %% +#Drawing(800, 800, "Tm1TmYAverages.pdf") +Drawing(400, 450) +origin() +kernels = typetriad("Tm1", TmYtypes, radius = 3, hexelsize = 16, normalize = :separate) +finish() +preview() + +# %% [markdown] +# ### Fig 2g example TmY4 cells + +# %% +@mfalse ids = type2ids("TmY4") +# Drawing(800, 800, "Tm1TmY4ExamplesAligned.pdf") +@drawsvg rfs = celltriad("Tm1", ids[13:15], radius = 5, hexelsize = 12) 430 500 + +# %% [markdown] +# ### Fig 2h example TmY9q cells + +# %% +@mfalse ids = type2ids("TmY9q") +#Drawing(800, 800, "Tm1TmY9qExamplesAligned.pdf") +Drawing(800, 800) +@drawsvg celltriad("Tm1", ids[1:3], radius = 5, hexelsize = 12) 430 500 + +# %% [markdown] +# ### Fig 2i example TmY9q⊥ cells + +# %% +@mfalse ids = ind2id[ind2type .== "TmY9q⊥"] +#Drawing(800, 800, "Tm1TmY9qperpExamplesAligned.pdf") +Drawing(800, 800) +@drawsvg celltriad("Tm1", ids[1:3], radius = 5, hexelsize = 12) 430 500 + +# %% [markdown] +# ### Fig 2j Tm1-TmY ellipses + +# %% +#Drawing(600, 600, "TmYEllipses.pdf") +Drawing(600, 600) +origin() +scale(50) +setline(0.5) + +for (celltype, ellipsecolor) in zip(TmYtypes, ["red", "green", "blue"]) + setcolor(ellipsecolor) + ellipses = ellipsesummary.(preimage("Tm1", celltype)) + for e in sample(ellipses, 50, replace=false) + if ~ismissing(e) + Luxor.rotate(-pi/6 + e.theta) + Luxor.ellipse(0, 0, e.minor, e.major, :stroke) # sqrt(3) for conversion to units of latticeconstant + Luxor.rotate(pi/6 - e.theta) + end + end +end +sethue("black") +setline(4) +drawpqaxes(1) +finish() +preview() + +# %% [markdown] +# ### Fig 2kl Tm1-TmY angles and aspect ratios + +# %% +ellipses = [ellipsesummary.(preimage("Tm1", TmYtype)) for TmYtype in TmYtypes] +aspectratios = [[e.major/e.minor for e in skipmissing(ellipse)] for ellipse in ellipses] +angles = [[wraparound(e.theta - pi/6, center = pi/2) for e in skipmissing(ellipse)] for ellipse in ellipses] +diamsmajor = [[e.major for e in skipmissing(ellipse)] for ellipse in ellipses] +diamsminor = [[e.minor for e in skipmissing(ellipse)] for ellipse in ellipses] +TmYellipses = OrderedDict(TmYtypes .=> ellipses) +TmYaspectratios = OrderedDict(TmYtypes .=> aspectratios) +TmYangles = OrderedDict(TmYtypes .=> angles) +celltypes = vcat([fill(i, length(v)) for (i, v) in enumerate(values(TmYaspectratios))]...) + +# %% +default(; # Plots defaults + fontfamily="Arial Unicode", + label="", # only explicit legend entries + guidefontsize = 12, + tickfontsize = 12, + legendfontsize = 10, +# linewidth = 2, +# ylim = (0, Inf), + widen = true, +# xrotation = 30 + ) + +# %% +haspect = boxplot(celltypes, vcat(values(TmYaspectratios)...), ylabel = "aspect ratio", + xtick = strings2ticks(TmYtypes), ylim = (1, Inf) +) + +# %% +celltypes = vcat([fill(i, length(v)) for (i, v) in enumerate(values(TmYangles))]...) +hangle = boxplot(celltypes, vcat(values(TmYangles)...)*180/pi, ylabel = "orientation (deg)", + xtick = strings2ticks(TmYtypes), yticks = [60, 90, 150] +) + +# %% +plot(hangle, haspect, size = (800, 300), left_margin = 5Measures.mm, bottom_margin = 3Measures.mm) + +# %% +savefig("Tm1TmYanglesaspects.pdf") + +# %% [markdown] +# ### Fig 2mno TmY projections + +# %% +paths = [["Tm1", TmYtype] for TmYtype in TmYtypes] +rfs = inmaps.(tracetypes.(paths)) +radius = 4 +rfscrop = [passmissing(crop).(rf.array, id2pq.(type2ids(TmYtype)), radius, 1) for (rf, TmYtype) in zip(rfs, TmYtypes)]; + +# %% +trans = [hexproject.(skipmissing(ims), ax) for (ims, ax) in zip(rfscrop, [:v, :q⊥, :q])] +long = [hexproject.(skipmissing(ims), ax) for (ims, ax) in zip(rfscrop, [:h, :q, :q⊥])]; + +# %% +Plots.reset_defaults() + +# %% +default(; # Plots defaults + fontfamily="Arial Unicode", + label="", # only explicit legend entries + guidefontsize = 12, + tickfontsize = 12, + legendfontsize = 10, +# linewidth = 2, +# ylim = (0, Inf), + widen = true, +# xrotation = 30 + ) + +# %% +h = [groupedbar(-radius:radius, [a c], yerror = [b d], + linewidth = 1, + label = ["longitudinal" "transverse"], + ylabel = e == "TmY4" ? "synapses" : "", + xlabel = "displacement", + annotations = (-3, 15, Plots.text(e, "Arial Unicode")), + ylim = (0, 20) +) for (a, b, c, d, e) in zip(mean.(long), std.(long), mean.(trans), std.(trans), TmYtypes)] + +# %% +plot(h..., layout = (1, 3), size = (1000, 200), bottom_margin = 10Measures.mm, left_margin = 5Measures.mm) + +# %% +savefig("TmYprojections.pdf") + +# %% [markdown] +# ## Fig 4 + +# %% [markdown] +# ### Fig 4abcd Dm3 average ERFs + +# %% +# threshold eliminated. now just the six strongest disynaptic pathways for each target +hexelsize = 7 + +monosynaptic = [["Tm1", Dm3type] for Dm3type in Dm3types] +crfs = inmaps.(tracetypes.(monosynaptic)) +crfsellipses = [ellipsesummary.(z.array) for z in crfs] +crfscrop = [passmissing(crop).(ims.array, id2pq.(ids), 5) for (ims, ids) in zip(crfs, type2ids.(Dm3types))] +crfsavg = mean.(skipmissing.(crfscrop)) +crfsavgellipses = NamedArray(ellipsesummary.(crfsavg), Dm3types) # only the ellipse is used here + +corner = 260 +Drawing(1200, 1200, "Dm3ERF.pdf") +#Drawing(1200, 1200) + +target = "Dm3p" + +paths = [["Tm1", t, target] for t in intrinsictypes] +p = sortperm(scorepath.(paths), rev = true) +paths = paths[p[1:6]] +intermediaries = getindex.(paths, 2) + +rfs = inmaps.(tracebacktypes.(paths)) +rfscrop = [passmissing(crop).(rf.array, id2pq.(type2ids(target)), 7, 1) for rf in rfs] +rfsavg = [mean(skipmissing(rf)) for rf in rfscrop] +cseries = map(x -> x ? "green" : "magenta", type2nt[intermediaries].array .== "ACH") + +s = maximum.(rfsavg) +s[cseries .== "green"] ./= maximum(s[cseries .== "green"]) +s[cseries .!= "green"] ./= maximum(s[cseries .!= "green"]) + +origin() +Luxor.translate(-corner, -corner) +setline(2) +hexannulus(rfsavg, orbitscale = 1.5, text = type2label.(intermediaries), ellipses = ellipsesummary.(rfsavg), ellipsecolors = cseries, maxvals = true) +Luxor.scale(1.5) +Luxor.rotate(-pi/6) +Luxor.translate(-crfsavgellipses[target].center) +Luxor.rotate(pi/6) +setcolor("blue") +fontsize(12) +Luxor.text(target, Point(0, -35), halign = :center) +setline(3) +setdash("dotted") +drawellipse(crfsavgellipses[target]) +setdash("solid") +for (i, e) in enumerate(ellipsesummary.(rfsavg)) + setopacity(s[i]) + drawellipse(e, color = cseries[i]) +end +Luxor.rotate(-pi/6) +setline(1) +setcolor("black") +Luxor.poly(hexcenter.([HexagonEye(1, 0, hexelsize), HexagonEye(0, 0, hexelsize), HexagonEye(0, 1, hexelsize)]), :stroke) +Luxor.rotate(pi/6) + + +#################### + +target = "Dm3q" + +paths = [["Tm1", t, target] for t in intrinsictypes] +p = sortperm(scorepath.(paths), rev = true) +paths = paths[p[1:6]] +intermediaries = getindex.(paths, 2) + +rfs = inmaps.(tracebacktypes.(paths)) +rfscrop = [passmissing(crop).(rf.array, id2pq.(type2ids(target)), 7, 1) for rf in rfs] +rfsavg = [mean(skipmissing(rf)) for rf in rfscrop] +cseries = map(x -> x ? "green" : "magenta", type2nt[intermediaries].array .== "ACH") + +s = maximum.(rfsavg) +s[cseries .== "green"] ./= maximum(s[cseries .== "green"]) +s[cseries .!= "green"] ./= maximum(s[cseries .!= "green"]) + +origin() +Luxor.translate(corner, -corner) +setline(2) +hexannulus(rfsavg, orbitscale = 1.5, text = type2label.(intermediaries), ellipses = ellipsesummary.(rfsavg), ellipsecolors = cseries, maxvals = true) +Luxor.scale(1.5) +Luxor.rotate(-pi/6) +Luxor.translate(-crfsavgellipses[target].center) +Luxor.rotate(pi/6) +setcolor("blue") +fontsize(12) +Luxor.text(target, Point(0, -35), halign = :center) +setline(3) +setdash("dotted") +drawellipse(crfsavgellipses[target]) +setdash("solid") +for (i, e) in enumerate(ellipsesummary.(rfsavg)) + setopacity(s[i]) + drawellipse(e, color = cseries[i]) +end +Luxor.rotate(-pi/6) +setline(1) +setcolor("black") +Luxor.poly(hexcenter.([HexagonEye(1, 0, hexelsize), HexagonEye(0, 0, hexelsize), HexagonEye(0, 1, hexelsize)]), :stroke) +Luxor.rotate(pi/6) + + +################ +target = "Dm3v" + +paths = [["Tm1", t, target] for t in intrinsictypes] +p = sortperm(scorepath.(paths), rev = true) +paths = paths[p[1:6]] +intermediaries = getindex.(paths, 2) + +rfs = inmaps.(tracebacktypes.(paths)) +rfscrop = [passmissing(crop).(rf.array, id2pq.(type2ids(target)), 7, 1) for rf in rfs] +rfsavg = [mean(skipmissing(rf)) for rf in rfscrop] +cseries = map(x -> x ? "green" : "magenta", type2nt[intermediaries].array .== "ACH") + +s = maximum.(rfsavg) +s[cseries .== "green"] ./= maximum(s[cseries .== "green"]) +s[cseries .!= "green"] ./= maximum(s[cseries .!= "green"]) + +origin() +Luxor.translate(-corner, corner) +setline(2) +hexannulus(rfsavg, orbitscale = 1.5, text = type2label.(intermediaries), ellipses = ellipsesummary.(rfsavg), ellipsecolors = cseries, maxvals = true) +Luxor.scale(1.5) +Luxor.rotate(-pi/6) +Luxor.translate(-crfsavgellipses[target].center) +Luxor.rotate(pi/6) +setcolor("blue") +fontsize(12) +Luxor.text(target, Point(0, -35), halign = :center) +setline(3) +setdash("dotted") +drawellipse(crfsavgellipses[target]) +setdash("solid") +for (i, e) in enumerate(ellipsesummary.(rfsavg)) + setopacity(s[i]) + drawellipse(e, color = cseries[i]) +end +Luxor.rotate(-pi/6) +setline(1) +setcolor("black") +Luxor.poly(hexcenter.([HexagonEye(1, 0, hexelsize), HexagonEye(0, 0, hexelsize), HexagonEye(0, 1, hexelsize)]), :stroke) +Luxor.rotate(pi/6) + + +###### summary +origin() + +Luxor.translate(corner, corner) + +ellipses = [ellipsesummary.(rfs.array) for rfs in inmaps.(tracetypes.([["Tm1", "T2a", Dm3type] for Dm3type in Dm3types]))] +crfellipses = [ellipsesummary.(rfs.array) for rfs in inmaps.(tracetypes.([["Tm1", Dm3type] for Dm3type in Dm3types]))]; +println(length.(ellipses)) + +#Luxor.translate(0, corner/6) +scale(3) +setline(0.5) +setopacity(1) + +cseries = distinguishable_colors(3, [RGB(1,1,1), RGB(0,0,0), colorant"magenta", colorant"green"], dropseed=true) + +for (ellipse, crfellipse, ellipsecolor) in zip(ellipses, crfellipses, cseries) + sethue(ellipsecolor) + for icell in sample(1:length(ellipse), 50, replace=false) + if ~ismissing(ellipse[icell]) && ~ismissing(crfellipse[icell]) + Luxor.rotate(-pi/6) + Luxor.translate(-hexelsize.*crfellipse[icell].center) + Luxor.rotate(pi/6) + drawellipse(ellipse[icell], color = ellipsecolor, hexelsize = hexelsize) + Luxor.rotate(-pi/6) + Luxor.translate(hexelsize.*crfellipse[icell].center) + Luxor.rotate(pi/6) + end + end +end +sethue("black") +setline(6) +drawpqaxes(hexelsize) + + +# origin() + +# Luxor.translate(corner, corner) + +# paths = [["Tm1", "TmY9q⊥", "Dm3v"], ["Tm1", "TmY9q", "Dm3p"], ["Tm1", "TmY9q⊥", "Dm3q"]] +# ellipses = [ellipsesummary.(rfs.array) for rfs in inmaps.(tracetypes.(paths))]; + +# Luxor.translate(corner/2, -corner/2) +# scale(1.5) +# setline(0.5) +# setopacity(1) + +# hexelsize = 7 + +# for (ellipse, crfellipse, ellipsecolor) in zip(ellipses, crfellipses, ["red", "green", "blue"]) +# sethue(ellipsecolor) +# for icell in sample(1:length(ellipse), 50, replace=false) +# if ~ismissing(ellipse[icell]) && ~ismissing(crfellipse[icell]) +# Luxor.rotate(-pi/6) +# Luxor.translate(-hexelsize.*crfellipse[icell].center) +# Luxor.rotate(pi/6) +# drawellipse(ellipse[icell], color = ellipsecolor, hexelsize = hexelsize) +# Luxor.rotate(-pi/6) +# Luxor.translate(hexelsize.*crfellipse[icell].center) +# Luxor.rotate(pi/6) +# end +# end +# end +# sethue("black") +# setline(6) +# Luxor.line(hexelsize.*Point(-6, 0), hexelsize.*Point(-6, sqrt(3)), :stroke) + +finish() + +preview() + +# %% +paths = [["Tm1", "T2a", Dm3type] for Dm3type in Dm3types] +ellipses = [ellipsesummary.(rfs.array) for rfs in inmaps.(tracetypes.(paths))]; + +# %% +crfellipses = [ellipsesummary.(rfs.array) for rfs in inmaps.(tracetypes.([["Tm1", Dm3type] for Dm3type in Dm3types]))]; + +# %% +cseries = distinguishable_colors(3, [RGB(1,1,1), RGB(0,0,0), colorant"magenta", colorant"green"], dropseed=true) + +# %% +#Drawing(600, 600, "Dm3Ellipses.pdf") +Drawing(900, 900) +origin() +scale(5) +setline(0.5) + +hexelsize = 7 + +for (ellipse, crfellipse, ellipsecolor) in zip(ellipses, crfellipses, cseries) + sethue(ellipsecolor) + for icell in sample(1:length(ellipse), 50, replace=false) + if ~ismissing(ellipse[icell]) && ~ismissing(crfellipse[icell]) + Luxor.rotate(-pi/6) + Luxor.translate(-hexelsize.*crfellipse[icell].center) + Luxor.rotate(pi/6) + drawellipse(ellipse[icell], color = ellipsecolor, hexelsize = hexelsize) + Luxor.rotate(-pi/6) + Luxor.translate(hexelsize.*crfellipse[icell].center) + Luxor.rotate(pi/6) + end + end +end +sethue("black") +setline(6) +drawpqaxes(hexelsize) +finish() +preview() + +# %% [markdown] +# ### Fig 4efg T2a-Dm3 projections + +# %% +erfpaths = [["Tm1", "T2a", Dm3type] for Dm3type in Dm3types] +crfpaths = [["Tm1", Dm3type] for Dm3type in Dm3types] + +# %% +erfs = inmaps.(tracebacktypes.(erfpaths)) +radius = 7 +erfscrop = [passmissing(crop).(erf.array, id2pq.(type2ids(Dm3type)), radius, 1) for (erf, Dm3type) in zip(erfs, Dm3types)]; + +# %% +crfs = inmaps.(tracebacktypes.(crfpaths)) +radius = 7 +crfscrop = [passmissing(crop).(crf.array, id2pq.(type2ids(Dm3type)), radius, 1) for (crf, Dm3type) in zip(crfs, Dm3types)]; + +# %% +erfprojs = [hexproject.(skipmissing(ims), ax) for (ims, ax) in zip(erfscrop, [:v, :p, :q])] +crfprojs = [hexproject.(skipmissing(ims), ax) for (ims, ax) in zip(crfscrop, [:v, :p, :q])]; + +# %% +Plots.reset_defaults() + +# %% +default(; # Plots defaults + fontfamily="Arial Unicode", + label="", # only explicit legend entries + guidefontsize = 12, + tickfontsize = 12, + legendfontsize = 10, + linewidth = 2, +# ylim = (0, Inf), + widen = true, +# xrotation = 30 + ) + +# %% +h = [] +for (a, b, c, d, e) in zip(mean.(erfprojs), std.(erfprojs), mean.(crfprojs), std.(crfprojs), Dm3types) +plot(-radius:radius, 100*a, ribbon = 100*b, + xlabel = "displacement", + xflip = e == "Dm3p", + xticks = (-6:2:6), + ylim = (0, 0.3), + label = "T2a", + ylabel = e == "Dm3v" ? "input fraction (%)" : "", + legend = :topleft, + color = 1, + annotations = (:topcenter, Plots.text(e, "Helvetica", :center)), +) +push!(h, plot!(twinx(), -radius:radius, + 100*c, + ribbon = 100*d, + ylim = (0, 10), + yticks = (0:2:10), + label = "direct", + legend = :topright, + color = 2 +)) +end + +# %% +plot(h..., layout = (1, 3), size = (1000, 250), bottom_margin = 10Measures.mm, left_margin = 5Measures.mm) + +# %% +savefig("T2aDm3ERFprojections.pdf") + +# %% [markdown] +# ## Fig 5 + +# %% [markdown] +# ### Fig 5abcd TmY average ERFs + +# %% +thres = 0.0025 + +# %% +monosynaptic = [["Tm1", TmYtype] for TmYtype in TmYtypes] +crfs = inmaps.(tracetypes.(monosynaptic)) +crfsellipses = [ellipsesummary.(z.array) for z in crfs] +crfscrop = [passmissing(crop).(ims.array, id2pq.(ids), 5) for (ims, ids) in zip(crfs, type2ids.(TmYtypes))] +crfsavg = mean.(skipmissing.(crfscrop)) +crfsavgellipses = NamedArray(ellipsesummary.(crfsavg), TmYtypes) # only the ellipse is used here + +# %% +corner = 260 +hexelsize = 7 + +Drawing(1200, 1200, "TmYERF.pdf") +#Drawing(1200, 1200) + +target = "TmY4" + +intermediaries = intrinsictypes[[scorepath(["Tm1", celltype, target]) for celltype in intrinsictypes] .> thres] +paths = [["Tm1", t, target] for t in intermediaries] + +p = sortperm(scorepath.(paths), rev = true) +rfs = inmaps.(tracebacktypes.(paths[p])) + +rfscrop = [passmissing(crop).(rf.array, id2pq.(type2ids(target)), 7, 1) for rf in rfs] +rfsavg = [mean(skipmissing(rf)) for rf in rfscrop] +cseries = map(x -> x ? "green" : "magenta", type2nt[intermediaries[p]].array .== "ACH") + +s = maximum.(rfsavg) +s[cseries .== "green"] ./= maximum(s[cseries .== "green"]) +s[cseries .!= "green"] ./= maximum(s[cseries .!= "green"]) + +origin() +Luxor.translate(-corner, -corner) +setline(2) +hexannulus(rfsavg, orbitscale = 1.5, text = type2label.(intermediaries[p]), ellipses = ellipsesummary.(rfsavg), ellipsecolors = cseries, maxvals = true) +Luxor.scale(1.5) +Luxor.rotate(-pi/6) +Luxor.translate(-crfsavgellipses[target].center) +Luxor.rotate(pi/6) +setcolor("blue") +fontsize(12) +Luxor.text(type2label(target), Point(0, -35), halign = :center) +setline(3) +setdash("dotted") +drawellipse(crfsavgellipses[target]) +setdash("solid") +for (i, e) in enumerate(ellipsesummary.(rfsavg)) + setopacity(s[i]) + drawellipse(e, color = cseries[i]) +end +Luxor.rotate(-pi/6) +setline(1) +setcolor("black") +Luxor.poly(hexcenter.([HexagonEye(1, 0, hexelsize), HexagonEye(0, 0, hexelsize), HexagonEye(0, 1, hexelsize)]), :stroke) +Luxor.rotate(pi/6) + +#################### + +target = "TmY9q" +intermediaries = intrinsictypes[[scorepath(["Tm1", celltype, target]) for celltype in intrinsictypes] .> thres] +paths = [["Tm1", t, target] for t in intermediaries] + +p = sortperm(scorepath.(paths), rev = true) +rfs = inmaps.(tracebacktypes.(paths[p])) + +rfscrop = [passmissing(crop).(rf.array, id2pq.(type2ids(target)), 7, 1) for rf in rfs] +rfsavg = [mean(skipmissing(rf)) for rf in rfscrop] +cseries = map(x -> x ? "green" : "magenta", type2nt[intermediaries[p]].array .== "ACH") + +s = maximum.(rfsavg) +s[cseries .== "green"] ./= maximum(s[cseries .== "green"]) +s[cseries .!= "green"] ./= maximum(s[cseries .!= "green"]) + +origin() +Luxor.translate(corner, -corner) +setline(2) +hexannulus(rfsavg, orbitscale = 1.5, text = type2label.(intermediaries[p]), ellipses = ellipsesummary.(rfsavg), ellipsecolors = cseries, maxvals = true) +Luxor.scale(1.5) +Luxor.rotate(-pi/6) +Luxor.translate(-crfsavgellipses[target].center) +Luxor.rotate(pi/6) +setcolor("blue") +fontsize(12) +Luxor.text(type2label(target), Point(0, -35), halign = :center) +setline(3) +setdash("dotted") +drawellipse(crfsavgellipses[target]) +setdash("solid") +for (i, e) in enumerate(ellipsesummary.(rfsavg)) + setopacity(s[i]) + drawellipse(e, color = cseries[i]) +end +Luxor.rotate(-pi/6) +setline(1) +setcolor("black") +Luxor.poly(hexcenter.([HexagonEye(1, 0, hexelsize), HexagonEye(0, 0, hexelsize), HexagonEye(0, 1, hexelsize)]), :stroke) +Luxor.rotate(pi/6) + +################ +target = "TmY9q⊥" +intermediaries = intrinsictypes[[scorepath(["Tm1", celltype, target]) for celltype in intrinsictypes] .> thres] +paths = [["Tm1", t, target] for t in intermediaries] + +p = sortperm(scorepath.(paths), rev = true) +rfs = inmaps.(tracebacktypes.(paths[p])) + +rfscrop = [passmissing(crop).(rf.array, id2pq.(type2ids(target)), 7, 1) for rf in rfs] +rfsavg = [mean(skipmissing(rf)) for rf in rfscrop] +cseries = map(x -> x ? "green" : "magenta", type2nt[intermediaries[p]].array .== "ACH") + +s = maximum.(rfsavg) +s[cseries .== "green"] ./= maximum(s[cseries .== "green"]) +s[cseries .!= "green"] ./= maximum(s[cseries .!= "green"]) + +origin() +Luxor.translate(-corner, corner) +setline(2) +hexannulus(rfsavg, orbitscale = 1.5, text = type2label.(intermediaries[p]), ellipses = ellipsesummary.(rfsavg), ellipsecolors = cseries, maxvals = true) +Luxor.scale(1.5) +Luxor.rotate(-pi/6) +Luxor.translate(-crfsavgellipses[target].center) +Luxor.rotate(pi/6) +setcolor("blue") +fontsize(12) +Luxor.text(type2label(target), Point(0, -35), halign = :center) +setline(3) +setdash("dotted") +drawellipse(crfsavgellipses[target]) +setdash("solid") +for (i, e) in enumerate(ellipsesummary.(rfsavg)) + setopacity(s[i]) + drawellipse(e, color = cseries[i]) +end +Luxor.rotate(-pi/6) +setline(1) +setcolor("black") +Luxor.poly(hexcenter.([HexagonEye(1, 0, hexelsize), HexagonEye(0, 0, hexelsize), HexagonEye(0, 1, hexelsize)]), :stroke) +Luxor.rotate(pi/6) + +###### summary +origin() + +Luxor.translate(corner, corner) + +paths = [["Tm1", "TmY4", "TmY4"], ["Tm1", "TmY9q", "TmY9q"], ["Tm1", "TmY9q⊥", "TmY9q⊥"]] +ellipses = [ellipsesummary.(rfs.array) for rfs in inmaps.(tracetypes.(paths))]; + +crfellipses = [ellipsesummary.(rfs.array) for rfs in inmaps.(tracetypes.([["Tm1", TmYtype] for TmYtype in TmYtypes]))]; + +#Luxor.translate(-corner/2, -corner/2) +scale(3) +setline(0.5) +setopacity(1) + +cseries = distinguishable_colors(3, [RGB(1,1,1), RGB(0,0,0), colorant"magenta", colorant"green"], dropseed=true) + +hexelsize = 7 + +for (ellipse, crfellipse, ellipsecolor) in zip(ellipses, crfellipses, cseries) + sethue(ellipsecolor) + for icell in sample(1:length(ellipse), 50, replace=false) + if ~ismissing(ellipse[icell]) && ~ismissing(crfellipse[icell]) + Luxor.rotate(-pi/6) + Luxor.translate(-hexelsize.*crfellipse[icell].center) + Luxor.rotate(pi/6) + drawellipse(ellipse[icell], color = ellipsecolor, hexelsize = hexelsize) + Luxor.rotate(-pi/6) + Luxor.translate(hexelsize.*crfellipse[icell].center) + Luxor.rotate(pi/6) + end + end +end +sethue("black") +setline(6) +drawpqaxes(hexelsize) + +finish() +preview() + +# %% [markdown] +# ### Fig 5efg TmY-TmY projections + +# %% +erfpaths = [["Tm1", TmYtype, TmYtype] for TmYtype in TmYtypes] +crfpaths = [["Tm1", TmYtype] for TmYtype in TmYtypes] + +# %% +radius = 6 +erfs = inmaps.(tracebacktypes.(erfpaths)) +erfscrop = [passmissing(crop).(erf.array, id2pq.(type2ids(TmYtype)), radius, 1) for (erf, TmYtype) in zip(erfs, TmYtypes)]; + +# %% +crfs = inmaps.(tracebacktypes.(crfpaths)) +crfscrop = [passmissing(crop).(crf.array, id2pq.(type2ids(TmYtype)), radius, 1) for (crf, TmYtype) in zip(crfs, TmYtypes)]; + +# %% +erfprojs = [hexproject.(skipmissing(ims), ax) for (ims, ax) in zip(erfscrop, [:v, :p, :q])] +crfprojs = [hexproject.(skipmissing(ims), ax) for (ims, ax) in zip(crfscrop, [:v, :p, :q])]; + +# %% +Plots.reset_defaults() + +# %% +default(; # Plots defaults + fontfamily="Helvetica", + label="", # only explicit legend entries + guidefontsize = 12, + tickfontsize = 12, + legendfontsize = 10, + linewidth = 2, +# ylim = (0, Inf), + widen = true, +# xrotation = 30 + ) + +# %% +h = [] +for (a, b, c, d, e) in zip(mean.(erfprojs), std.(erfprojs), mean.(crfprojs), std.(crfprojs), TmYtypes) +plot(-radius:radius, 100*a, ribbon = 100*b, + xlabel = "displacement", + xticks = (-6:2:6), + yticks = (0:0.1:0.2), + ylim = (0, 0.2), + label = "TmY", + ylabel = e == "TmY4" ? "input fraction (%)" : "", + legend = :topleft, + annotations = (:topcenter, Plots.text(e, "Helvetica", :center)), + color = 1 +) +push!(h, + plot!(twinx(), -radius:radius, 100*c, ribbon = 100*d, + ylim = (0, 8), + yticks = (0:2:8), + label = "direct", + legend = :topright, + color = 2 +)) +end + +# %% +plot(h..., layout = (1, 3), size = (1000, 250), bottom_margin = 10Measures.mm, left_margin = 5Measures.mm) + +# %% +savefig("TmYTmYERFprojections.pdf") + +# %% [markdown] +# ## Fig 6 LC15 + +# %% [markdown] +# ### LC15 locations + +# %% +# anchor LC15 cells on Mi1→T3→LC15 maps +# missing means that there is zero connectivity +rfs = prepreimage("Mi1", "T3", "LC15") +for (i, id) in enumerate(type2ids("LC15")) + center = findcenter(float(rfs[i]), seven) + if ~ismissing(center) + id2pq[id] = center + end +end + +# %% [markdown] +# ### Fig 6b LC15 average disynaptic maps + +# %% +thres = 0.01 +scores = [scorepath([c1, c2, "LC15"]) for c1 in hexel, c2 in nothexel] +scores = NamedArray(scores, names = (hexel, nothexel)) +scoresums = NamedArray(sum(scores.array, dims = 1)[:], nothexel) +scoremaxs = NamedArray(maximum(scores.array, dims = 1)[:], nothexel) + +# no need to exclude inhibitory intermediaries here, as top six are predicted excitatory +# contrast with LC10ev +intermediaries = nothexel[scoresums.array .> thres] + +p = sortperm(scoresums[intermediaries], rev = true) + +scoresums[intermediaries[p]] + +# %% +# source hexel types +ignore, winners = findmax(scores[:, intermediaries], dims=1) +sources = hexel[getindex.(winners, 1)][:] +sources[p] + +# %% +paths = [[source, intermediary, "LC15"] for (source, intermediary) in zip(sources[p], intermediaries[p])] +rfs = inmaps.(tracebacktypes.(paths)); + +# %% +cseries = distinguishable_colors(6, ColorSchemes.hot[:], dropseed=true) + +# %% +rfscrop = [passmissing(crop).(rf.array, id2pq.(type2ids("LC15")), 8, 1) for rf in rfs] +rfsavg = [mean(skipmissing(rf)) for rf in rfscrop] + +#Drawing(700, 600) +Drawing(700, 600, "LC15.pdf") + origin() + ellipses = ellipsesummary.(rfsavg) + hexannulus(rfsavg, orbitscale = 1.5, text = type2label.(intermediaries[p]), ellipses = ellipses, ellipsecolors = cseries, maxvals = true) + Luxor.scale(1.5) + drawpqaxes(7) + Luxor.rotate(-pi/6) + Luxor.translate(-ellipses[1].center) # place T3 ellipse at center (doesn't matter so much for this plot) + Luxor.rotate(pi/6) + setline(3) + for (i, e) in enumerate(ellipses) + if in(i, [3 4 5]) # make TmY ellipses more visible + setopacity(1) + else + setopacity(0.3) + end + drawellipse(e, color = cseries[i]) + end +finish() +preview() + +# %% [markdown] +# ### Fig 6c LC15 ellipses + +# %% +ellipses = [ellipsesummary.(rf.array) for rf in rfs] # for individual cells, rfs computed above +hexelsize = 7 + +Drawing(500, 600, "Tm1LC15Ellipses.pdf") +#Drawing(500, 600) +origin() +scale(3) +setline(1) +for (e1, e2, e3, e4) in zip(ellipses[1], ellipses[3], ellipses[4], ellipses[5]) + Luxor.rotate(-pi/6) + Luxor.translate(-hexelsize * e1.center) # anchor on T3 ellipse, which isn't shown + Luxor.rotate(pi/6) + drawellipse(e2, color = cseries[3], hexelsize = hexelsize) + drawellipse(e3, color = cseries[4], hexelsize = hexelsize) + drawellipse(e4, color = cseries[5], hexelsize = hexelsize) + Luxor.rotate(-pi/6) + Luxor.translate(hexelsize*e1.center) + Luxor.rotate(pi/6) +end +setline(5) +setcolor("black") +drawpqaxes(hexelsize) + +finish() +preview() + +# %% [markdown] +# ## Fig 6 LC10e ventral + +# %% +### new clustering +LC10ed = [720575940603823136, 720575940604935089, 720575940610075029, 720575940614496735, 720575940614746178, 720575940616179050, 720575940618360741, 720575940620416176, 720575940621175979, 720575940622744414, 720575940622805854, 720575940623979815, 720575940624787504, 720575940625143396, 720575940625981080, 720575940628023296, 720575940628359783, 720575940628615752, 720575940629385943, 720575940629398492, 720575940629958764, 720575940631576546, 720575940633027615, 720575940637401321, 720575940637706574, 720575940637731955, 720575940638701375, 720575940644393070] +LC10ev = [720575940604125920, 720575940606274592, 720575940608977988, 720575940609733387, 720575940609890402, 720575940612305250, 720575940612332401, 720575940614897662, 720575940621233157, 720575940621602047, 720575940622486057, 720575940622602992, 720575940623601193, 720575940624050407, 720575940624595893, 720575940628632064, 720575940629076042, 720575940629756556, 720575940630438651, 720575940631305467, 720575940632034385, 720575940634088602, 720575940640974555] + +# %% +issetequal(type2ids("LC10e"), union(LC10ed, LC10ev)) + +# %% [markdown] +# ### ventral locations in pq coordinates + +# %% +# anchor LC10 cells on Tm1→TmY9q→LC10 maps +# missing means that there is zero connectivity +rfs = prepreimage("Tm1", "TmY9q", LC10ev); +# other possible anchors +# rfs = prepreimage("Tm1", "TmY9q⊥", LC10ev); +# rfs = prepreimage("Mi9", "Tm8a", LC10ev); +@mfalse for (i, id) in enumerate(LC10ev) + center = findcenter(float(rfs[i]), seven) + if !ismissing(center) + id2pq[id] = center + end +end + +# %% [markdown] +# ### Fig 6e LC10ev average disynaptic maps + +# %% +infractionLC10ev = sum(A'*W[:, Name.(LC10ev)], dims = 2)/sum(W[:, Name.(LC10ev)]) +infractionLC10ev = NamedArray(infractionLC10ev[:], alltypes) + +thres = 0.003 +scores = [infraction[c1, c2]*infractionLC10ev[c2] for c1 in hexel, c2 in nothexel] +scores = NamedArray(scores, names = (hexel, nothexel)) +scoresums = NamedArray(sum(scores.array, dims = 1)[:], nothexel) +scoremaxs = NamedArray(maximum(scores.array, dims = 1)[:], nothexel) + +intermediaries = nothexel[scoresums.array .> thres .&& type2nt[nothexel].array .== "ACH"] +#intermediaries = nothexel[scoresums.array .> thres] +p = sortperm(scoresums[intermediaries], rev = true) + +scoresums[intermediaries[p]] + +# %% +ignore, winners = findmax(scores[:, intermediaries], dims=1) +sources = hexel[getindex.(winners, 1)][:] +sources[p] + +# %% +paths = Vector{Union{String, Vector{<:Integer}}}[[source, intermediary, LC10ev] for (source, intermediary) in zip(sources[p], intermediaries[p])] +rfs = inmaps.(tracebacktypes.(paths)) +pathtexts = [[source, intermediary, "LC10ev"] for (source, intermediary) in zip(sources[p], intermediaries[p])] + +# %% +cseries = distinguishable_colors(6, ColorSchemes.hot[:], dropseed=true) + +# %% +rfscrop = [passmissing(crop).(rf.array, id2pq.(LC10ev), 8, 1) for rf in rfs] +rfsavg = [mean(skipmissing(rf)) for rf in rfscrop] + +Drawing(700, 600, "LC10ev.pdf") +#Drawing(700, 600) + origin() + ellipses = ellipsesummary.(rfsavg) + hexannulus(rfsavg, orbitscale = 1.5, text = type2label.(intermediaries[p]), ellipses = ellipses, ellipsecolors = cseries, maxvals = true) + Luxor.scale(1.5) + drawpqaxes(7) + Luxor.rotate(-pi/6) + Luxor.translate(-ellipses[1].center) + Luxor.rotate(pi/6) + setline(3) + for (i, e) in enumerate(ellipses) + if i in [1, 2] + setopacity(1) + else + setopacity(0.3) + end + drawellipse(e, color = cseries[i]) + end +finish() +preview() + +# %% [markdown] +# ### Fig 6f ventral ellipses + +# %% +# overload function (hexgraphics.jl) to handle LC10ev as list of IDs, rather than a string indicating a type +function OpticLobe.tracebacktypes(celltypes::Vector{Any}) + @mfalse inds = [isa(celltype, String) ? ind2type .== celltype : id2ind.(celltype) for celltype in celltypes] + return *([Wback[inds[i], inds[i+1]] for i = 1:length(inds)-1]...) +end + +# %% +hexelsize = 7 + +#paths = [["Mi9", "Tm8a", LC10ev], ["Tm1", "TmY9q", LC10ev], ["Tm1", "TmY9q⊥", LC10ev]] +paths = Vector{Union{String, Vector{<:Integer}}}[["Tm1", "TmY9q", LC10ev], ["Tm1", "TmY9q⊥", LC10ev]] + +rfs = inmaps.(tracebacktypes.(paths)) +ellipses = [ellipsesummary.(rf.array) for rf in rfs] + +Drawing(600, 800, "Tm1LC10evEllipses.pdf") +#Drawing(600, 800) +origin() +scale(3) +setline(1) +for (e1, e2) in zip(ellipses[1], ellipses[2]) + Luxor.rotate(-pi/6) + Luxor.translate(-hexelsize * e1.center) + Luxor.rotate(pi/6) + drawellipse(e1, color = cseries[1], hexelsize = hexelsize) + drawellipse(e2, color = cseries[2], hexelsize = hexelsize) + Luxor.rotate(-pi/6) + Luxor.translate(hexelsize*e1.center) + Luxor.rotate(pi/6) +end + +setline(5) +setcolor("black") +drawpqaxes(hexelsize) +finish() +preview() + +# %% [markdown] +# ### Fig 6f inset displacement vectors + +# %% +rot30 = [cos(pi/6) sin(pi/6); -sin(pi/6) cos(pi/6)] + +# %% +rfs = [prepreimage("Mi9", "Tm8a", LC10ev), prepreimage("Tm1", "TmY9q", LC10ev), prepreimage("Tm1", "TmY9q⊥", LC10ev)] +ellipses = [ellipsesummary.(ims) for ims in rfs] + +d = [collect(e.center)/sqrt(3) for e in ellipses[3]] - [collect(e.center)/sqrt(3) for e in ellipses[2]] + +# %% +hexelsize = 50 +Drawing(600, 600, "LC10eDisplacements.pdf") +#Drawing(600, 600) +origin() +setcolor("gray") + for x in d + Luxor.arrow(Point(0, 0), Point(hexelsize*rot30*x...), arrowheadlength=20, linewidth = 5) + end +Luxor.rotate(-pi/6) +setline(10) +setcolor("black") +Luxor.poly(hexcenter.([HexagonEye(1, 0, hexelsize), HexagonEye(0, 0, hexelsize), HexagonEye(0, 1, hexelsize)]), :stroke) +Luxor.rotate(pi/6) +finish() +preview() + +# %% [markdown] toc-hr-collapsed=true +# ## Fig S3 Tm1Mi4Tm2L3-Dm3 + +# %% +using Format + +# %% +function valuesandtext(values, textstring, fmtstring = "%3d") + setcolor("black") + fontsize(20) + fontface("Arial Unicode MS") + Luxor.text(cfmt(fmtstring, values[1]), Point(-170, 115)) + Luxor.text(cfmt(fmtstring, values[2]), Point(-190, 75)) + Luxor.text(cfmt(fmtstring, values[3]), Point(-150, 75)) + Luxor.text(textstring, Point(-150, 155), halign = :center) +end + +# %% +Drawing(900, 800, "Mi4Tm2L3Tm9-Dm3.pdf") +#Drawing(900, 800) +origin() +Luxor.translate(-200, -200) +rfs = typetriad("Mi4", Dm3types, radius = 3, hexelsize = 15, normalize = :separate) +valuesandtext(maximum.(rfs), "Mi4→Dm3", "%3.1f") + +origin() +Luxor.translate(200, -200) +rfs = typetriad("Tm2", Dm3types, radius = 3, hexelsize = 15, normalize = :separate) +valuesandtext(maximum.(rfs), "Tm2→Dm3", "%3.1f") + +origin() +Luxor.translate(-200, 200) +rfs = typetriad("L3", Dm3types, radius = 3, hexelsize = 15, normalize = :separate) +valuesandtext(maximum.(rfs), "L3→Dm3", "%3.1f") + +origin() +Luxor.translate(200, 200) +rfs = typetriad("Tm9", Dm3types, radius = 3, hexelsize = 15, normalize = :separate) +valuesandtext(maximum.(rfs), "Tm9→Dm3", "%3.1f") + +finish() +preview() + +# %% [markdown] toc-hr-collapsed=true +# ## Fig S6 lateral interactions Dm3-TmY circuit + +# %% +function maparray(pretypes, posttypes; radius = 6, spacing = 140) + rfscrop = [passmissing(crop).(preimage(pretype, posttype), id2pq.(type2ids(posttype)), radius, 1) for pretype in pretypes, posttype in posttypes] + rfscropavg = sum.(skipmissing.(rfscrop))./length.(rfscrop) + immax = maximum(hcat(rfscropavg[:]...)) + + for pre = 1:length(pretypes) + for post = 1:length(posttypes) + x, y = (post-(length(posttypes)+1)/2)*spacing, (pre-(length(pretypes)+1)/2)*spacing + Luxor.translate(x, y) + rect2hex(square2hex(get(ColorSchemes.hot, rfscropavg[pre, post]/immax))) + Luxor.translate(-x, -y) + end + end + return immax +end + +# %% +Drawing(1000, 1000, "Dm3TmYAllAverages.pdf") +spacing = 250 +#Drawing(1000, 1000) +origin() + +Luxor.translate(-spacing, -spacing) +maxvalue = maparray(Dm3types, Dm3types) +fontsize(16) +Luxor.text(@sprintf("maxvalue = %3.1f", maxvalue), Point(0, -220), halign = :center) + +origin() +Luxor.translate(spacing, -spacing) +maxvalue = maparray(Dm3types, TmYtypes) +fontsize(16) +Luxor.text(@sprintf("maxvalue = %3.1f", maxvalue), Point(0, -220), halign = :center) + +origin() +Luxor.translate(-spacing, spacing) +maxvalue = maparray(TmYtypes, Dm3types) +fontsize(16) +Luxor.text(@sprintf("maxvalue = %3.1f", maxvalue), Point(0, -220), halign = :center) + +origin() +Luxor.translate(spacing, spacing) +maxvalue = maparray(TmYtypes, TmYtypes) +fontsize(16) +Luxor.text(@sprintf("maxvalue = %3.1f", maxvalue), Point(0, -220), halign = :center) + +finish() +preview() + +# %% [markdown] +# ## Data S3 Dm3 individual cells +# TODO: update to `distinguishable_colors` + +# %% [markdown] +# ### Dm3p + +# %% +target = "Dm3p" +monosynaptic = names(first(sort(infraction[hexel, target], rev=true), 5), 1) + +scores = [infraction[c1, c2]*infraction[c2, target] for c1 in hexel, c2 in nothexel] +scores = NamedArray(scores, names = (hexel, nothexel)) +ignore, winners = findmax(scores, dims=1) +sources = hexel[getindex.(winners, 1)][:] + +scoresums = NamedArray(sum(scores.array, dims = 1)[:], nothexel) +scoremaxs = NamedArray(maximum(scores.array, dims = 1)[:], nothexel) + +#p = sortperm(scoremaxs, rev = true) +p = sortperm(scoresums, rev = true) + +disynaptic = [sources[p] nothexel[p]][1:10, :] + +paths = vcat( + [[t, target] for t in monosynaptic], + [[t..., target] for t in eachrow(disynaptic)] + ) + +rfs = inmaps.(tracebacktypes.(paths)); + +tocolor = length.(paths) .== 3 +tocolor[1] = true +ellipsecolors = Vector{Any}(nothing, length(paths)) +ellipsecolors[tocolor] .= distinguishable_colors(sum(tocolor)) + +@showprogress for id in type2ids(target) + ims = [rf[Name(id)] for rf in rfs] + montage(ims, fname = "$target/$id.pdf", labels = path2label.(paths), hexelsize = 6, ellipses = true, ellipsecolors = ellipsecolors, maxvals = true, centers = repeat([id2pq[id]], length(ims)), summary = 1) +end + +# %% [markdown] +# ### Dm3q + +# %% +target = "Dm3q" +monosynaptic = names(first(sort(infraction[hexel, target], rev=true), 5), 1) + +scores = [infraction[c1, c2]*infraction[c2, target] for c1 in hexel, c2 in nothexel] +scores = NamedArray(scores, names = (hexel, nothexel)) +ignore, winners = findmax(scores, dims=1) +sources = hexel[getindex.(winners, 1)][:] + +scoresums = NamedArray(sum(scores.array, dims = 1)[:], nothexel) +scoremaxs = NamedArray(maximum(scores.array, dims = 1)[:], nothexel) + +#p = sortperm(scoremaxs, rev = true) +p = sortperm(scoresums, rev = true) + +disynaptic = [sources[p] nothexel[p]][1:10, :] + +paths = vcat( + [[t, target] for t in monosynaptic], + [[t..., target] for t in eachrow(disynaptic)] + ) + +rfs = inmaps.(tracebacktypes.(paths)); + +tocolor = length.(paths) .== 3 +tocolor[1] = true +ellipsecolors = Vector{Any}(nothing, length(paths)) +ellipsecolors[tocolor] .= distinguishable_colors(sum(tocolor)) + +@showprogress for id in type2ids(target) + ims = [rf[Name(id)] for rf in rfs] + montage(ims, fname = "$target/$id.pdf", labels = path2label.(paths), hexelsize = 6, ellipses = true, ellipsecolors = ellipsecolors, maxvals = true, centers = repeat([id2pq[id]], length(ims)), summary = 1) +end + +# %% [markdown] +# ### Dm3v + +# %% +target = "Dm3v" +monosynaptic = names(first(sort(infraction[hexel, target], rev=true), 5), 1) + +scores = [infraction[c1, c2]*infraction[c2, target] for c1 in hexel, c2 in nothexel] +scores = NamedArray(scores, names = (hexel, nothexel)) +ignore, winners = findmax(scores, dims=1) +sources = hexel[getindex.(winners, 1)][:] + +scoresums = NamedArray(sum(scores.array, dims = 1)[:], nothexel) +scoremaxs = NamedArray(maximum(scores.array, dims = 1)[:], nothexel) + +#p = sortperm(scoremaxs, rev = true) +p = sortperm(scoresums, rev = true) + +disynaptic = [sources[p] nothexel[p]][1:10, :] + +paths = vcat( + [[t, target] for t in monosynaptic], + [[t..., target] for t in eachrow(disynaptic)] + ) + +rfs = inmaps.(tracebacktypes.(paths)); + +tocolor = length.(paths) .== 3 +tocolor[1] = true +ellipsecolors = Vector{Any}(nothing, length(paths)) +ellipsecolors[tocolor] .= distinguishable_colors(sum(tocolor)) + +@showprogress for id in type2ids(target) + ims = [rf[Name(id)] for rf in rfs] + montage(ims, fname = "$target/$id.pdf", labels = path2label.(paths), hexelsize = 6, ellipses = true, ellipsecolors = ellipsecolors, maxvals = true, centers = repeat([id2pq[id]], length(ims)), summary = 1) +end + +# %% [markdown] +# ## Data S4 TmY individual cells (supplementary data) + +# %% [markdown] +# ### TmY4 + +# %% +target = "TmY4" +monosynaptic = names(first(sort(infraction[hexel, target], rev=true), 4), 1) + +scores = [infraction[c1, c2]*infraction[c2, target] for c1 in hexel, c2 in nothexel] +scores = NamedArray(scores, names = (hexel, nothexel)) +ignore, winners = findmax(scores, dims=1) +sources = hexel[getindex.(winners, 1)][:] + +scoresums = NamedArray(sum(scores.array, dims = 1)[:], nothexel) +scoremaxs = NamedArray(maximum(scores.array, dims = 1)[:], nothexel) + +#p = sortperm(scoremaxs, rev = true) +p = sortperm(scoresums, rev = true) + +disynaptic = [sources[p] nothexel[p]][1:10, :] + +trisynaptic = [["Tm1", "TmY4", "Dm3v",]] + +paths = vcat( + [[t, target] for t in monosynaptic], + [[t..., target] for t in eachrow(disynaptic)], + [[t..., target] for t in trisynaptic] + ) + +rfs = inmaps.(tracebacktypes.(paths)); + +tocolor = length.(paths) .>= 3 +tocolor[1] = true +ellipsecolors = Vector{Any}(nothing, length(paths)) +ellipsecolors[tocolor] .= distinguishable_colors(sum(tocolor)) + +@showprogress for id in type2ids(target) + ims = [rf[Name(id)] for rf in rfs] + montage(ims, fname = "$target/$id.pdf", labels = path2label.(paths), hexelsize = 6, ellipses = true, ellipsecolors = ellipsecolors, maxvals = true, centers = repeat([id2pq[id]], length(ims)), summary = 1) +end + +# %% [markdown] +# ### TmY9q + +# %% +target = "TmY9q" +monosynaptic = names(first(sort(infraction[hexel, target], rev=true), 4), 1) + +scores = [infraction[c1, c2]*infraction[c2, target] for c1 in hexel, c2 in nothexel] +scores = NamedArray(scores, names = (hexel, nothexel)) +ignore, winners = findmax(scores, dims=1) +sources = hexel[getindex.(winners, 1)][:] + +scoresums = NamedArray(sum(scores.array, dims = 1)[:], nothexel) +scoremaxs = NamedArray(maximum(scores.array, dims = 1)[:], nothexel) + +#p = sortperm(scoremaxs, rev = true) +p = sortperm(scoresums, rev = true) + +disynaptic = [sources[p] nothexel[p]][1:10, :] + +trisynaptic = [["Tm1", "TmY9q", "Dm3p",]] + +paths = vcat( + [[t, target] for t in monosynaptic], + [[t..., target] for t in eachrow(disynaptic)], + [[t..., target] for t in trisynaptic] + ) + + +rfs = inmaps.(tracebacktypes.(paths)); + +tocolor = length.(paths) .>= 3 +tocolor[1] = true +ellipsecolors = Vector{Any}(nothing, length(paths)) +ellipsecolors[tocolor] .= distinguishable_colors(sum(tocolor)) + +@showprogress for id in type2ids(target) + ims = [rf[Name(id)] for rf in rfs] + montage(ims, fname = "$target/$id.pdf", labels = path2label.(paths), hexelsize = 6, ellipses = true, ellipsecolors = ellipsecolors, maxvals = true, centers = repeat([id2pq[id]], length(ims)), summary = 1) +end + +# %% [markdown] +# ### TmY9q⊥ + +# %% +target = "TmY9q⊥" +monosynaptic = names(first(sort(infraction[hexel, target], rev=true), 4), 1) + +scores = [infraction[c1, c2]*infraction[c2, target] for c1 in hexel, c2 in nothexel] +scores = NamedArray(scores, names = (hexel, nothexel)) +ignore, winners = findmax(scores, dims=1) +sources = hexel[getindex.(winners, 1)][:] + +scoresums = NamedArray(sum(scores.array, dims = 1)[:], nothexel) +scoremaxs = NamedArray(maximum(scores.array, dims = 1)[:], nothexel) + +#p = sortperm(scoremaxs, rev = true) +p = sortperm(scoresums, rev = true) + +disynaptic = [sources[p] nothexel[p]][1:10, :] + +trisynaptic = [["Tm1", "TmY9q⊥", "Dm3q",]] + +paths = vcat( + [[t, target] for t in monosynaptic], + [[t..., target] for t in eachrow(disynaptic)], + [[t..., target] for t in trisynaptic] + ) + +rfs = inmaps.(tracebacktypes.(paths)); + +tocolor = length.(paths) .>= 3 +tocolor[1] = true +ellipsecolors = Vector{Any}(nothing, length(paths)) +ellipsecolors[tocolor] .= distinguishable_colors(sum(tocolor)) + +@showprogress for id in type2ids(target) + ims = [rf[Name(id)] for rf in rfs] + montage(ims, fname = "$target/$id.pdf", labels = path2label.(paths), hexelsize = 6, ellipses = true, ellipsecolors = ellipsecolors, maxvals = true, centers = repeat([id2pq[id]], length(ims)), summary = 1) +end + +# %% [markdown] +# ## Data S5 LC10ev individual cells + +# %% +infractionLC10ev = sum(A'*W[:, Name.(LC10ev)], dims = 2)/sum(W[:, Name.(LC10ev)]) +infractionLC10ev = NamedArray(infractionLC10ev[:].array, alltypes) + +thres = 0.001 +scores = [infraction[c1, c2]*infractionLC10ev[c2] for c1 in hexel, c2 in nothexel] +scores = NamedArray(scores, names = (hexel, nothexel)) +scoresums = NamedArray(sum(scores.array, dims = 1)[:], nothexel) +scoremaxs = NamedArray(maximum(scores.array, dims = 1)[:], nothexel) + +intermediaries = nothexel[scoresums.array .> thres .&& type2nt[nothexel].array .== "ACH"] +p = sortperm(scoresums[intermediaries], rev = true) + +scoresums[intermediaries[p]] + +# %% +ignore, winners = findmax(scores[:, intermediaries[p]], dims=1) +sources = hexel[getindex.(winners, 1)][:] + +# %% +paths = Vector{Union{String, Vector{<:Integer}}}[[source, intermediary, LC10ev] for (source, intermediary) in zip(sources, intermediaries[p])] + +rfs = inmaps.(tracebacktypes.(paths)) + +pathtexts = [[source, intermediary, "LC10ev"] for (source, intermediary) in zip(sources, intermediaries[p])] + +# %% +cseries = distinguishable_colors(length(paths), ColorSchemes.hot[:], dropseed=true) + +# %% +hexelsize = 6 + +@showprogress for id in LC10ev + ims = [rf[Name(id)] for rf in rfs] + montage(ims, fname = "LC10ev/$id.pdf", + labels = path2label.(pathtexts), hexelsize = hexelsize, ellipses = true, ellipsecolors = cseries, + maxvals = true, centers = repeat([id2pq[id]], length(ims)), summary = 1) +end + +# %% [markdown] +# ## Data S5 LC15 individual cells + +# %% +#thres = 0.0025 +thres = 0.005 +target = "LC15" +scores = [scorepath([c1, c2, target]) for c1 in hexel, c2 in nothexel] +scores = NamedArray(scores, names = (hexel, nothexel)) +scoresums = NamedArray(sum(scores.array, dims = 1)[:], nothexel) +scoremaxs = NamedArray(maximum(scores.array, dims = 1)[:], nothexel) + +# %% +intermediaries = nothexel[scoresums.array .> thres .&& type2nt[nothexel].array .== "ACH"] +p = sortperm(scoresums[intermediaries], rev = true) + +# %% +ignore, winners = findmax(scores[:, intermediaries], dims=1) +sources = hexel[getindex.(winners, 1)][:] + +# %% +sources[p] + +# %% +paths = [[source, intermediary, "LC15"] for (source, intermediary) in zip(sources[p], intermediaries[p])] + +rfs = inmaps.(tracebacktypes.(paths)) + +# %% +cseries = distinguishable_colors(length(paths), ColorSchemes.hot[:], dropseed=true) + +# %% +hexelsize = 6 + +@showprogress for id in type2ids("LC15") + ims = [rf[Name(id)] for rf in rfs] + montage(ims, fname = "LC15/$id.pdf", + labels = path2label.(paths), hexelsize = hexelsize, ellipses = true, ellipsecolors = cseries, + maxvals = true, centers = repeat([id2pq[id]], length(ims)), summary = 1) +end + +# %%