-
Notifications
You must be signed in to change notification settings - Fork 3
/
getExpressionEnergy.m
53 lines (43 loc) · 1.8 KB
/
getExpressionEnergy.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
function expressionEnergy = getExpressionEnergy(geneEntrezID,structureIDs)
% Retrieves expression energy from the Allen SDK
%-------------------------------------------------------------------------------
%-------------------------------------------------------------------------------
% Input checking:
%-------------------------------------------------------------------------------
if nargin < 1
geneEntrezID = 20604;
end
if nargin < 2 || isempty(structureIDs)
C = load('Mouse_Connectivity_Data.mat','RegionStruct');
isCortex = ([C.RegionStruct.OhRegionIndex]==1);
% IDs of all cortical brain regions:
structureIDs = [C.RegionStruct(isCortex).id];
end
%-------------------------------------------------------------------------------
% Write the gene entrez ID to geneEntrezID.csv
csvwrite(fullfile('AllenSDK','geneEntrezID.csv'),geneEntrezID);
% Write the structure list to structureIDs.csv
csvwrite(fullfile('AllenSDK','structureIDs.csv'),structureIDs);
% Run the python script
% (seems like it's possible to do this through matlab without file i/o, but
% modules need to sit in the python search path?
% cf. https://au.mathworks.com/help/matlab/matlab_external/call-user-defined-custom-module.html)
%
% P = py.sys.path;
% if count(P,'AllenSDK') == 0
% insert(P,int32(0),'AllenSDK');
% end
cd('AllenSDK');
system('python RetrieveGene.py');
cd('../')
% Read in the csv output
csvFileName = sprintf('expressionEnergy_gene%u.csv',geneEntrezID);
expressionEnergy = csvread(csvFileName);
% Match output structures to input structures, just to be extra safe:
[~,ia,ib] = intersect(expressionEnergy(:,1),structureIDs);
if length(ia) < length(structureIDs)
error('Could not retrieve expression data for all structures??');
end
% Output the expression energy vector
expressionEnergy = expressionEnergy(ia,2);
end