Init

.gitattributes

@@ -0,0 +1 @@
+*.csv filter=lfs diff=lfs merge=lfs -text

Main.m

@@ -0,0 +1,38 @@
+% This script takes calcium imaging data as input and returns the underlying neuronal connectome 04/06/2020
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%                             MAIN SCRIPT                             %%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+warning ('off','all'); clear; clc; close all; 
+path = genpath('subfunctions'); addpath(path); clear path;
+
+%% INPUTS
+INPUTS.ops = sp_deconv_inputs(); % Taking imaging rate and fluorescence decay rate as inputs
+
+%% DATA
+INPUTS.calcium = import_data(); % NN by NT (NN: number of neurons, NT: number of time frames)
+
+%% Spike Deconvolution  
+RESULTS.spikes_raw = sp_deconv(INPUTS.ops,INPUTS.calcium); % Non-negative deconvolution using Suite2P and OASIS packages
+
+%% Spike Denoising
+RESULTS.spikes_denoised = sp_denoise(RESULTS.spikes_raw); % Removing the noise present in deconvolved spikes using an optimal threshold
+
+%% Spike Smoothing
+RESULTS.spikes_smooth = sp_smooth(RESULTS.spikes_denoised); % Smoothing the spikes to improve the accuracy of connectome inference
+
+%% Partial Correlation
+RESULTS.partial_correlation = partial_corr(RESULTS.spikes_smooth); % Computing partial correlation using matrix inversion 
+
+%% Network
+[RESULTS.network,RESULTS.edges] = connectome_inf(RESULTS.partial_correlation); % Connectome Inference based on partial correlation coeff.
+
+%% Circular Graph
+circular_graph = circ_graph(RESULTS.network); % Circular graph of the inferred neuronal connectome
+
+%% Graph with positions
+pos_graph = graph_w_positions(RESULTS.network); % inferred neuronal connectome based on actual spatial arrangement of neurons 
+                                                %(if distance/position information is available)
+
+%% Network Dynamics
+net_dyn(RESULTS.spikes_raw,RESULTS.partial_correlation);

README.md

@@ -0,0 +1,7 @@
+# FARCI
+
+
+![alt text](https://github.com/Sabermea/FARCI/blob/master/docs/logo.jpg?raw=true)
+
+
+FARCI provides an efficient implementation of neuronal connectome inference from Calcium imaging neuronal activity data. FARCI outputs a partial correlation matrix, which is depicted as a circular graph, representing the neuronal connectome. When neuronal spatial information is available, the neuronal connectome also accounts for the position of neurons.

subfunctions/circ_graph.m

@@ -0,0 +1,16 @@
+function circular_graph = circ_graph(network)
+
+if ~iscell(network)
+    figure;
+    myColorMap = lines(length(network));
+    circular_graph = circularGraph(network,'Colormap',myColorMap);
+
+elseif iscell(network)
+
+    for i = 1:length(network)
+        figure;
+        myColorMap = lines(length(network{1,i}));
+        circular_graph = circularGraph(network{1,i},'Colormap',myColorMap);
+        title(['Circular Graph (Session ' num2str(i) ')']);
+    end
+end

subfunctions/circularGraph.m

@@ -0,0 +1,203 @@
+classdef circularGraph < handle
+% CIRCULARGRAPH Plot an interactive circular graph to illustrate connections in a network.
+%
+%% Syntax
+% circularGraph(X)
+% circularGraph(X,'PropertyName',propertyvalue,...)
+% h = circularGraph(...)
+%
+%% Description
+% A 'circular graph' is a visualization of a network of nodes and their
+% connections. The nodes are laid out along a circle, and the connections
+% are drawn within the circle. Click on a node to make the connections that
+% emanate from it more visible or less visible. Click on the 'Show All'
+% button to make all nodes and their connections visible. Click on the
+% 'Hide All' button to make all nodes and their connections less visible.
+%
+% Required input arguments.
+% X : A symmetric matrix of numeric or logical values.
+%
+% Optional properties.
+% Colormap : A N by 3 matrix of [r g b] triples, where N is the 
+%            length(adjacenyMatrix).
+% Label    : A cell array of N strings.
+%%
+% Copyright 2016 The MathWorks, Inc.
+  properties
+    Node = node(0,0); % Array of nodes
+    ColorMap;         % Colormap
+    Label;            % Cell array of strings
+    ShowButton;       % Turn all nodes on
+    HideButton;       % Turn all nodes off
+  end
+
+  methods
+    function this = circularGraph(adjacencyMatrix,varargin)
+      % Constructor
+      p = inputParser;
+
+      defaultColorMap = parula(length(adjacencyMatrix));
+      defaultLabel = cell(length(adjacencyMatrix));
+      for i = 1:length(defaultLabel)
+        defaultLabel{i} = num2str(i);
+      end
+
+      addRequired(p,'adjacencyMatrix',@(x)(isnumeric(x) || islogical(x)));
+      addParameter(p,'ColorMap',defaultColorMap,@(colormap)length(colormap) == length(adjacencyMatrix));
+      addParameter(p,'Label'   ,defaultLabel   ,@iscell);
+
+      parse(p,adjacencyMatrix,varargin{:});
+      this.ColorMap = p.Results.ColorMap;
+      this.Label    = p.Results.Label;
+
+      this.ShowButton = uicontrol(...
+        'Style','pushbutton',...
+        'Position',[0 40 80 40],...
+        'String','Show All',...
+        'Callback',@circularGraph.showNodes,...
+        'UserData',this);
+
+      this.HideButton = uicontrol(...
+        'Style','pushbutton',...
+        'Position',[0 0 80 40],...
+        'String','Hide All',...
+        'Callback',@circularGraph.hideNodes,...
+        'UserData',this);
+
+      fig = gcf;
+      set(fig,...
+        'UserData',this,...
+        'CloseRequestFcn',@circularGraph.CloseRequestFcn);
+
+      % Draw the nodes
+      delete(this.Node);
+      t = linspace(-pi,pi,length(adjacencyMatrix) + 1).'; % theta for each node
+      extent = zeros(length(adjacencyMatrix),1);
+      for i = 1:length(adjacencyMatrix)
+        this.Node(i) = node(cos(t(i)),sin(t(i)));
+        this.Node(i).Color = this.ColorMap(i,:);
+        this.Node(i).Label = this.Label{i};
+      end
+
+      % Find non-zero values of s and their indices
+      [row,col,v] = find(adjacencyMatrix);
+
+      % Calculate line widths based on values of s (stored in v).
+      minLineWidth  = 0.5;
+      lineWidthCoef = 5;
+      lineWidth = v./max(v);
+      if sum(lineWidth) == numel(lineWidth) % all lines are the same width.
+        lineWidth = repmat(minLineWidth,numel(lineWidth),1);
+      else % lines of variable width.
+        lineWidth = lineWidthCoef*lineWidth + minLineWidth;
+      end
+
+      % Draw connections on the Poincare hyperbolic disk.
+      %
+      % Equation of the circles on the disk:
+      % x^2 + y^2 
+      % + 2*(u(2)-v(2))/(u(1)*v(2)-u(2)*v(1))*x 
+      % - 2*(u(1)-v(1))/(u(1)*v(2)-u(2)*v(1))*y + 1 = 0,
+      % where u and v are points on the boundary.
+      %
+      % Standard form of equation of a circle
+      % (x - x0)^2 + (y - y0)^2 = r^2
+      %
+      % Therefore we can identify
+      % x0 = -(u(2)-v(2))/(u(1)*v(2)-u(2)*v(1));
+      % y0 = (u(1)-v(1))/(u(1)*v(2)-u(2)*v(1));
+      % r^2 = x0^2 + y0^2 - 1
+
+      for i = 1:length(v)
+        if row(i) ~= col(i)
+          if abs(row(i) - col(i)) - length(adjacencyMatrix)/2 == 0 
+            % points are diametric, so draw a straight line
+            u = [cos(t(row(i)));sin(t(row(i)))];
+            v = [cos(t(col(i)));sin(t(col(i)))];
+            this.Node(row(i)).Connection(end+1) = line(...
+              [u(1);v(1)],...
+              [u(2);v(2)],...
+              'LineWidth', lineWidth(i),...
+              'Color', this.ColorMap(row(i),:),...
+              'PickableParts','none');
+          else % points are not diametric, so draw an arc
+            u  = [cos(t(row(i)));sin(t(row(i)))];
+            v  = [cos(t(col(i)));sin(t(col(i)))];
+            x0 = -(u(2)-v(2))/(u(1)*v(2)-u(2)*v(1));
+            y0 =  (u(1)-v(1))/(u(1)*v(2)-u(2)*v(1));
+            r  = sqrt(x0^2 + y0^2 - 1);
+            thetaLim(1) = atan2(u(2)-y0,u(1)-x0);
+            thetaLim(2) = atan2(v(2)-y0,v(1)-x0);
+
+%             if u(1) >= 0 && v(1) >= 0 
+%               % ensure the arc is within the unit disk
+%               theta = [linspace(max(thetaLim),pi,50),...
+%                        linspace(-pi,min(thetaLim),50)].';
+%             else
+%               theta = linspace(thetaLim(1),thetaLim(2)).';
+%             end
+
+            sortedthetaLim=sort(thetaLim);
+            if sortedthetaLim(2)-sortedthetaLim(1)>pi
+            thetaLim(1)=sortedthetaLim(2)-(2*pi);
+            thetaLim(2)=sortedthetaLim(1);
+            end
+            theta = linspace(thetaLim(1),thetaLim(2)).';
+
+            this.Node(row(i)).Connection(end+1) = line(...
+              r*cos(theta)+x0,...
+              r*sin(theta)+y0,...
+              'LineWidth', lineWidth(i),...
+              'Color', this.ColorMap(row(i),:),...
+              'PickableParts','none');
+          end
+        end
+      end
+
+      axis image;
+      ax = gca;
+      for i = 1:length(adjacencyMatrix)
+        extent(i) = this.Node(i).Extent;
+      end
+      extent = max(extent(:));
+      ax.XLim = ax.XLim + extent*[-1 1];
+      fudgeFactor = 1.75; % Not sure why this is necessary. Eyeballed it.
+      ax.YLim = ax.YLim + fudgeFactor*extent*[-1 1];
+      ax.Visible = 'off';
+      ax.SortMethod = 'depth';
+
+      fig = gcf;
+      fig.Color = [1 1 1];
+    end
+
+  end
+
+  methods (Static = true)
+    function showNodes(this,~)
+      % Callback for 'Show All' button
+      n = this.UserData.Node;
+      for i = 1:length(n)
+        n(i).Visible = true;
+      end
+    end
+
+    function hideNodes(this,~)
+      % Callback for 'Hide All' button
+      n = this.UserData.Node;
+      for i = 1:length(n)
+        n(i).Visible = false;
+      end
+    end
+
+    function CloseRequestFcn(this,~)
+      % Callback for figure CloseRequestFcn
+      c = this.UserData;
+      for i = 1:length(c.Node)
+        delete(c.Node(i));
+      end
+      delete(gcf);
+    end
+
+  end
+
+end

subfunctions/connectome_inf.m

@@ -0,0 +1,34 @@
+function [net,edges] = connectome_inf(p)
+
+disp(['To set a threshold, please enter a value for '  (char(945)) ': '])
+alpha = input([ '(Note that, the threshold is defined as: Threshold = ' (char(956)) ' + ' (char(945)) ' * ' (char(963)) '): ' ]);
+
+if ~iscell(p)
+    net = p;
+    net = triu(net,1);
+    idx = net ~= 0;
+    th = mean(net(idx)) + alpha*std(net(idx));
+    net(net<th) = 0;
+    [row,col,weights] = find(net);
+    edges = [row,col,weights];
+    net = net+net';
+
+    figure;
+    imagesc(net); colormap('hot'); colorbar; xlabel('Neurons'); ylabel('Neurons'); title('Heatmap of Neuronal Connectome');
+
+elseif iscell(p)
+
+    for i = 1:length(p)
+        x = p{1,i};
+        x = triu(x,1);
+        idx = x ~= 0;
+        th = mean(x(idx)) + alpha*std(x(idx));
+        x(x<th) = 0;
+        [row,col,weights] = find(x);
+        edges = [row,col,weights];
+        net{1,i} = x+x';
+
+        figure;
+        imagesc(net{1,i}); colormap('hot'); colorbar; xlabel('Neurons'); ylabel('Neurons'); title(['Heatmap of Neuronal Connectome (Session ' num2str(i) ')']);
+    end
+end