-
Notifications
You must be signed in to change notification settings - Fork 0
/
MarmosetReadDataV1.m
1791 lines (1588 loc) · 73.6 KB
/
MarmosetReadDataV1.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
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
% B Pailthorpe, U. Sydney. V1 (27/9/19-24)
% Marmoset brain: Load data: connectivity; Coords;
% Dependencies: data files: AdjMarmoset.csv, AdjFullMarmoset.csv, AdjMarmoset_Rescale2b.csv,
% MarmosetAcrn116.csv, MarmosetAcrn139to116.csv, CoordsMarmoset.csv, MarmosetLobes116.csv,
% VolsMarmoset.csv, MarmosetPairListRescale2b.clu % InfoMap (c code) module ID list (.clu file).
% Marmoset_NodeFlow, Marmoset_FlowWt % Prob'y flow over Nodes & Links (.csv):
% for original or extra details & 3D vis.:
% Maybe: MarmosetGrayImageSlice.csv % saved image volume data, midslice,
% from Atlas volume (atlas.nii); nb. large files (86MB) - cf Appx.
% atlas_segmentation.nii (26MB); % needs NIFTI tools (Matlab), load_nii, etc.
% Do these things:
% #1, import link data & other info: Acrn; ??Coords (Atlas,
% & midslice of vol image
% Flows, nodes & links
% Infomap module ID, in array idx
% 1.8 d(i,j) dist matrix, from coords
% Link List ; top 40% (line 322)
% 1.9 Investigate node vol effects on FLNe wts (line 406)
% & 1.3a #Links: un-wt Deg (line 525 - out of order)
% 1.8.2 Calc LinkList of pos Wts & plot hist(log10)
% 2.0 LNe,1 (normalised to 1 mm^3 target vol) data (line 617 )
% 2.1 calc LN { * denom : tot#LN - LN,i } for ea target
% 2.2 calc LNe : new adj & wts (Line (778
% 3.0 "Symmetrise" the 116 x 55 Adj (line 692 +)
% 4.0 examine spread of Injection Volumes [line 741]
% 4.1 re-examine LN data (3 exp measurements): Denom ~ 0.2% variation: OK
% Apendices: % A.0 get RGB colours for areas: from Atlas files (line 958)
% & test & devmt codes
% A.1 get Atlas files { line 983
% calc Node volumes, CM, from bvol 3D data
% A.1.1 3D vol of areas ID# (in atlas_labels.txt) {line
% A.1.2a: mid slice grayscale (mm scale) from image data (line ; And 1193 [clean code])
% >>>>
%% 1.0 Read data files & setup
% may use matlab toolboxes for NIFTI files etc - in sub-paths
% nb. simpler path on this comp - vs laptop
% nb. read from InfoMap flow calc in code: PlotInfoMapV2.m
% & calc flows over i-j links in MB_flows.m (at #2.)
clear all; %close all
addpath(genpath('/bap_working/MatLabfiles/MatlabFiles/MarmosetBrain')); % Data files; include sub-directories
%addpath('/bap_working/MatLabfiles/MarmosetBrain/CodesOther'); % for other codes & dependencies
%
fprintf('\n >> Vis Marmoset Brain: : Load (raw) linklist data \n')
Actx=csvread('AdjMarmoset.csv'); % read square A(55 x 55); selection of raw A (both S&T); no self-loops
nsquare = size(Actx,1) % number of nodes: Sources AND Targets
Afull= csvread('AdjFullMarmoset.csv'); % calc prev. (116 x 116) % read raw A (all S, some T); no self-loops
[n, ~] = size(Afull) % number of Sources,
% set one of these Adj to A, for calc below:
figure; spyc(Actx) % 1854 links:
title('Marmoset Adj: Ctx-Ctx: 55 x 55 ')
figure; spyc(Afull) % 3474 links
title('Marmoset Adj: Ctx-Ctx: 116 x 55 ')
%plotArray(A) % coarser grid, larger squares, but only 1379 squares?
% set up larger figs for later
figwidth = 1024; figheight = 896; figposition = [100, 100, figwidth, figheight];
% A = Afull; % easier for other codes
% clear Afull
%% 1.01 fully connected 55x55 subnetwork (FLNe)
% count links (<~ 50-150)
n55Links = length(find(Actx(:)))
frn55 = n55Links/(55*54)
[sii sjj v]=find(Actx);
Alogical=full(sparse(sii, sjj, 1)); % "mask" of selected entries, to re-form sq array of 1's.
% nb. need to use sparse to form the array
A1 = Alogical.*1; % pick out the selected entries; yields a "full" array
figure; spy(A1) % debug
clear Alogical sii sjj v
length(find(A1(:))) % check 3474 links:
clear sii sjj v n55Links frn55
% length(find(A1(:)) )
% count links (<~ 50-150)
DegIn1=sum(A1,1)'; % i<-j; col sum (1st index): produces a row % show as Col
DegOut1=sum(A1,2); % i->j; row sum
figure; hist([DegIn1,DegOut1], 50); legend('k-In', 'k-Out')
xlabel('Degree (link count)'); ylabel('counts')
title('Marmoset, Degree of 55x55 binarised links')
% wt. degree (node strength)
DegIn55=sum(Actx,1)'; % i<-j; col sum (1st index): produces a row % show as Col
DegOut55=sum(Actx,2); % i->j; row sum
figure; hist([DegIn55,DegOut55], 50); legend('wtDeg-In', 'wtDeg-Out')
xlabel('weight (FLNe)'); ylabel('counts')
title('Marmoset Deg of 55x55 frn. wt. links')
[max(DegIn55) max(DegOut55) ]
% nb. scale doesnt make sense, since FLN is fraction amongst 116 sources (not 55).
%% 1.1d wt Degree (FLNe)
% unwt deg at 1.3a - far below
DegIn=sum(Afull,1)'; % i<-j; col sum (1st index): produces a row % show as Col
DegOut=sum(Afull,2); % i->j; row sum
% figure; hist(DegIn); title('Marmoset wt DegIn ')
%figure; hist(DegOut); hold on; title('Marmoset wt DegOut '); hold off
hFig=figure; hist([DegIn, DegOut]);
legend('DegIn', 'DegOut'); title(' Marmoset brain, FLN,e: Deg-In,-Out ')
% cf un-wt Deg (#Links) at #1.3a, line 480
kTot = length(find(Afull(:)))
clear n55Links frn55
%% 1.1a read V2b ADj, based on LNe
Anew= csvread('AdjMarmoset_Rescale2b.csv'); % read re-scaled Adj (116 x 116)
[n, ~] = size(Anew) % number of Sources
%% 1.1 read Node Acrn, Labels & Coord (CM)
% read the 116 Acrn for the Adj (Inj Source)
% reader skips blanks, so load 'na' to preserve numbering
filename = '/bap_working/MatLabfiles/MatlabFiles/MarmosetBrain/MarmosetAcrn116.csv';
delimiter = ','; fileID = fopen(filename,'r'); formatSpec = '%s%[^\n\r]';
dataArray = textscan(fileID, formatSpec, 'Delimiter', delimiter, 'ReturnOnError', false); fclose(fileID); % these are Cells, containing strings
Acrn = dataArray{1, 1}(1:end); % these are cells;
clearvars filename delimiter formatSpec fileID dataArray ans
Acrn{107} % check this is 'V1': ok
% 1.1a read the 139 Atlas labels
% reader skips blanks, so load 'na' to preserve numbering
filename = '/bap_working/MatLabfiles/MatlabFiles/MarmosetBrain/MarmosetLabels139.csv';
delimiter = ','; fileID = fopen(filename,'r'); formatSpec = '%s%[^\n\r]';
dataArray = textscan(fileID, formatSpec, 'Delimiter', delimiter, 'ReturnOnError', false); fclose(fileID); % these are Cells, containing strings
Labels = dataArray{1, 1}(1:end); % these are cells;
clearvars filename delimiter formatSpec fileID dataArray ans
Labels{126} % check this is 'V1': ok
NodeID=csvread('MarmosetAcrn139to116.csv'); % for the 116 Nodes (cf Links) vs 136 Labels
% translates the full 139 List (Atlas) to the 116 nodes (Source Injn)
% 1.1b coords - prev calc from Atlas volumes; cf Appx. A.1 below
% from 3D vol: 'atlas_segmentation.nii'
NodeCoord=csvread('CoordsMarmoset.csv'); % nb. some absent: NaN
% adjust z-origin, to match Paxinos atlas??
% NodeCoord(:,3)=NodeCoord(:,3)-20; % origin at V (top) & points down
max(NodeCoord(:, 3)) % z
NodeVols=csvread('VolsMarmoset.csv'); % (mm^3 - calc from bv.img[.nii]))
figure; hist(NodeVols,50); title('Marmoset: Node Vol Distrn (mm^3) ')
% NodeVoxels=csvread('VoxelsMarmoset.csv'); % get file to laptop
%% 1.1c read the 116 Lobes assigned (in Paxinos Atlas)
% reader skips blanks, so load 'na' to preserve numbering
filename = '/bap_working/MatLabfiles/MatlabFiles/MarmosetBrain/MarmosetLobes116.csv';
delimiter = ','; fileID = fopen(filename,'r'); formatSpec = '%s%[^\n\r]';
dataArray = textscan(fileID, formatSpec, 'Delimiter', delimiter, 'ReturnOnError', false); fclose(fileID); % these are Cells, containing strings
Lobes = dataArray{1, 1}(1:end); % these are cells;
clearvars filename delimiter formatSpec fileID dataArray ans
Lobes{108} % check this is 'Vis': ok
%LobesTxt=char(Lobes); % use char array to search the 116 Lobe names - faster
%% 1.1d wt Degree (FLNe)
% unwt deg at 1.3a - far below
DegIn=sum(Afull,1)'; % i<-j; col sum (1st index): produces a row % show as Col
DegOut=sum(Afull,2); % i->j; row sum
% figure; hist(DegIn); title('Marmoset wt DegIn ')
%figure; hist(DegOut); hold on; title('Marmoset wt DegOut '); hold off
hFig=figure; hist([DegIn, DegOut]);
legend('DegIn', 'DegOut'); title(' Marmoset brain, FLN,e: Deg-In,-Out ')
% cf un-wt Deg (#Links) at #1.3a, line 480
kIn = length(find(Afull(:)))
%% 1.1d.1 wt Degree (LNe)
DegIn=sum(Anew,1)'; % i<-j; col sum (1st index): produces a row % show as Col
DegOut=sum(Anew,2); % i->j; row sum
% figure; hist(DegIn); title('Marmoset wt DegIn ')
%figure; hist(DegOut); hold on; title('Marmoset wt DegOut '); hold off
hFig=figure; hist([DegIn, DegOut]);
legend('DegIn', 'DegOut'); title(' Marmoset brain, (LN,e): Deg-In,-Out ')
% #Links: un-wt Deg - binarise 116x116 A: UN-wt links
[sii sjj v]=find(Anew);
Alogical=full(sparse(sii, sjj, 1)); % "mask" of selected entries, to re-form sq array of 1's.
% nb. need to use sparse to form the array
A1 = Alogical.*1; % pick out the selected entries; yields a "full" array
figure; spy(A1)
clear Alogical sii sjj v
length(find(A1(:))) % check 3474 links:
% count links (<~ 50-150)
DegIn1=sum(A1,1)'; % i<-j; col sum (1st index): produces a row % show as Col
DegIn1=[DegIn1; 0; 0]; % append 2 missing 0's
DegOut1=sum(A1,2); % i->j; row sum
figure; hist(DegIn1, 50)
title('Marmoset DegIn of binarised links')
axis([0 110 0 10]) % nb. big # at 0: not sampled
%
figure; hist(DegOut1, 50); hold on
title('Marmoset DegOut of binarised links')
%axis([0 80 0 130]);
figure; stem (DegIn1, '.'); hold on
stem (DegOut1, 'r.');
legend('DegIn1', 'DegOut1')
title('Marmoset Deg of binarised links')
% joint plot
figure; hist([DegIn1, DegOut1], 50); xlim([2 Inf]); hold on
title('Marmoset DegIn, Out of binarised links')
legend('DegIn1', 'DegOut1')
%% node vol dependence of Deg, k
figure; subplot(2,2,1); plot(NodeVols, DegIn, '.', 'MarkerSize', 12); title('Marmoset: wt-DegIn')
subplot(2,2,2); plot(NodeVols, DegIn1, '.', 'MarkerSize', 12); title(' #Links, k-In')
subplot(2,2,3); plot(NodeVols, DegOut, '.', 'MarkerSize', 12); title(' wt-DegOut')
subplot(2,2,4); plot(NodeVols, DegOut1, '.', 'MarkerSize', 12); title(' #Links, k-Out')
% & expand scale:
figure; subplot(2,2,1); plot(NodeVols, DegIn, '.', 'MarkerSize', 12); title('Marmoset: wt-DegIn')
axis([0 40 -Inf Inf])
subplot(2,2,2); plot(NodeVols, DegIn1, '.', 'MarkerSize', 12); title(' #Links, k-In')
axis([0 40 -Inf Inf]);
subplot(2,2,3); plot(NodeVols, DegOut, '.', 'MarkerSize', 12); title(' wt-DegOut')
axis([0 40 -Inf Inf]); xlabel('Node Vol (mm^3) ')
subplot(2,2,4); plot(NodeVols, DegOut1, '.', 'MarkerSize', 12); title(' #Links, k-Out')
axis([0 40 -Inf Inf]); xlabel('Node Vol (mm^3) ')
% fits
figure; plot(NodeVols, DegIn1, '.', 'MarkerSize', 12); title(' #Links, k-In')
axis([0 40 -Inf Inf]); % use GUI for linear fit; avoid bias by many 0's
nodes1=find(DegIn1 >0);
DegIn11=DegIn1(nodes1); % need non-zero,
vols1=NodeVols(nodes1);
figure; plot(vols1, DegIn11, '.', 'MarkerSize', 12); title(' #Links, k-In')
axis([0 40 -Inf Inf]); %
[mean(DegIn11) std(DegIn11)] % stsats of non-zero k
clear nodes1 DegIn11 vols1
% Out links vs 2D perimeter
NodeRad=(3*NodeVols/(4*pi)); NodePerim =2*pi*NodeRad.^2; % assume cyl regions
figure; subplot(1,2,1); plot(NodePerim, DegIn1, '.', 'MarkerSize', 14); title(' #Marmoset: Links, k-In')
xlabel('Node perimeter (mm) '); axis([0 70 -Inf Inf]) % omit 2 big outliers
subplot(1,2,2); plot(NodePerim, DegOut1, '.', 'MarkerSize', 14); title(' #Links, k-Out')
xlabel('Node perimeter (mm) '); axis([0 70 -Inf Inf])
clear NodeRad NodePerim
%% 1.1e Strongest links
length(find(Afull(:)>= 0.1)) % strongest Out links
length(find(Afull(:)>= 0.4)) % V strongest Out links
%% 1.1g Deg, k - vol sampling (LNe)
figure; subplot(2,2,1); plot(NodeVols,
%% ++++>>> 1.4.1 NEEDED: set mid-transverse(x,y)-plane slice : as reference plane in 3D plots (mm space)
% needs the Atlas volume (.nii) - cf Appx , below
% read saved image volume data, midslice
gray_imageflip=csvread('MarmosetGrayImageSlice.csv'); % prev calc RGB; x-y flipped
gray_imageflip=uint8(gray_imageflip); % Need integer; nb image was x-y flipped
gray_imageflipR=fliplr(gray_imageflip); % reflect LHS to get RHS image
% set x-y box & desired z position of the image plane.
zpix=350; pixsize=[0.0400 0.5000 0.0400]; % set manually, avoid reloading large image
min_x = 1; max_x = 825; min_y = 1; max_y = 63;
min_z = 1; max_z = 550; imgzposition = zpix; % "half" height, pixel space
% pix-mm conversion, % apply vox -> mm tranform
origin = [411,24, 57]; % from bv.hdr.hist.originator
pixsize=[0.0400 0.5000 0.0400]; % set manually, avoid reloading large image
min_xmm=(min_x-origin(1))*pixsize(1); max_xmm=(max_x-origin(1))*pixsize(1);
min_ymm=(min_y-origin(2))*pixsize(2); max_ymm=(max_y-origin(2))*pixsize(2);
min_zmm=(min_z-origin(3))*pixsize(3); max_zmm=(max_z-origin(3))*pixsize(3);
imgzpositionmm=(imgzposition-origin(3))*pixsize(3); % centered, in mm-coords
clear origin pixsize min_x max_x min_y max_y min_z max_z imgzposition zpix
% test plot in mm space the image plane using surf.
figure; hold on; grid on % nb. need to eliminate "box"
surf([min_xmm max_xmm],[min_ymm max_ymm],repmat(imgzpositionmm, [2 2]),...
gray_imageflip,'facecolor','texture', 'FaceAlpha', 0.2,'edgecolor', 'none') %
colormap(gray); view(-30, 30)
plot3(0, 0, 20, '+k','MarkerSize', 30);
title('Marmoset Brain: Mid Slice z= -5.72 mm ')
ylabel('y (mm)'); xlabel('x (mm)'); % nb. swap
axis equal; axis([-10 10 -10 20 0 20])
% add RHS - symmetric about mid-line
surf([min_xmm max_xmm],[min_ymm max_ymm],repmat(imgzpositionmm, [2 2]),...
gray_imageflipR,'facecolor','texture', 'FaceAlpha', 0.2, 'EdgeAlpha', 0.0)
hold off
%% 1.7 ++NEEDED: Read Module IDs - InfoMap calc; Assign colours (ipsi-lateral modules)
% use A; etc, needs LinkList - calc from A116x52, A52x52, Aunwt52x52, etc
% needs idx: InfoMap module ID list = output of InfoMap calc, in .clu file
%idx=dlmread('MarmosetPairListRescale.clu'); % V2 of data: rescale by LNe/mm^3
%idx=dlmread('MarmosetPairList.clu'); % wt >=1; V1, based on FLNe
idx=dlmread('MarmosetPairListRescale2b.clu');
% idx=dlmread('MarmosetPairList55.clu'); % 55 x 55 Ctx only
% aggregate few orphans, based on linkages, tbd; 7 modules #8-14 orphans
% idx(find(idx>=7))=7; % orphans
idx(find(idx>=8))=8; % orphans: "last" module[s] {might be to #1 - majority of in/out links
%idx(51)=8; % reassign Apir as orphan [has no detected links in or out]
% & set Node, Module colors
%nm=max(idx) % # modules to use (from InfoMap)
fprintf(' >>>> use idx: Infomap of links, Assign Colors \n')
% %fprintf(' >>>> use idx: Infomap of links, bottom 10pc, Assign Colors \n')
nodecolors=zeros(n,3); modulecolors=zeros(8,3);
% ** update these (9/1/18) - to match modules members to Allen Atlas colours
% cf. http://atlas.brain-map.org/ for assigned colors
% >>> need to revise colours to match marmoset atlas
for i = 1:n
if idx(i) == 1
%nodecolors(i,:)=[1, 102/255, 216/255]; % for mouse: bright pink: MidBrain;
nodecolors(i,:)=[0.9412 1.0000 0.2314]; %use average of Atlas Node Colors
%modulecolors(1,:)=[1, 102/255, 216/255];% for mouse
modulecolors(1,:)= [0.9412 1.0000 0.2314]; % adjust for more contrast
%[0.9294 0.818 0.2839]; %use average of Atlas Node Colors
elseif idx(i) == 2
%nodecolors(i,:) =[2/255, 162/255, 88/255]; % mid-green: Vis
nodecolors(i,:) =[0.7037, 0.9495, 0.2913]; %use average of Atlas Node Colors
%modulecolors(2,:)=[2/255, 162/255, 88/255]; % for mouse
modulecolors(2,:)=[0.7037, 0.9495, 0.2913]; %use average of Atlas Node Colors
elseif idx(i) == 3
%nodecolors(i,:) = [250/255, 20/255, 20/255]; % red-ish/pink % DlpFC
nodecolors(i,:) = [1.0000 0.3255 0.0000]; % adjust for more contrast
%[0.9001 0.505 0.2539]; %use av
%modulecolors(3,:)= [250/255, 20/255, 20/255];
modulecolors(3,:)= [1.0000 0.3255 0.0000]; % use average of Atlas Node Colors
elseif idx(i) == 4
%nodecolors(i,:) =[121/255, 213/255, 190/255]; % : Tan
nodecolors(i,:) =[0.9592 0.4400 0.2694]; %use av
%modulecolors(4,:)=[0, 213/255, 190/255];
modulecolors(4,:)=[0.9592 0.4400 0.2694]; %use av
elseif idx(i) == 5
%nodecolors(i,:) = [0 100/255 137/255]; % : brownish : VL & DL pFC
nodecolors(i,:) = [0.9131 0.7196 0.4680]; % use av
%modulecolors(5,:)=[0 100/255 137/255];
modulecolors(5,:)=[0.9131 0.7196 0.4680]; % use av
elseif idx(i) == 6
%nodecolors(i,:) =[30/255, 240/255, 60/255]; % Yellow: Cing
nodecolors(i,:) =[0.9000 0.8961 0.2319]; % use av
% modulecolors(6,:)=[30/255, 240/255, 60/255];
modulecolors(6,:)=[0.9000 0.8961 0.2319]; % use av
elseif idx(i) == 7
%nodecolors(i,:) =[0.9,112/255, 160/255]; % yellow/orange: LTemp, Insul
nodecolors(i,:) =[0.9508 0.7109 0.2050]; % use av
%modulecolors(7,:)=[0.9,112/255, 160/255];
modulecolors(7,:)=[0.9508 0.7109 0.2050]; % use av
else % Gray - (grp 8+ are orpans)
nodecolors(i,:) =[0.5,0.5,0.5]; % others, minor modules : assign colours later?
modulecolors(8,:)=[0.5,0.5,0.5];
end
end
%% 1.7.1 Flow: Prob'y FLOWS over links - InfoMap
% use code PlotInfoMapV2.m to parse the InfoMap .map output file
% read flow over Nodes & Links:
NodeFlow= csvread('Marmoset_NodeFlow');
FlowWt=csvread('Marmoset_FlowWt');
%% 1.7.2 Flow: sort lists
list=FlowWt(find(FlowWt>0) ); % 1854 links
figure; hist(list, 1090)
listb=find(FlowWt(:)>1e-3); % 261 of
% make & sort list
[ii jj]=ind2sub(n,listb);
flowList=[ii jj zeros(length(ii),1)];
%
for i=1:length(ii)
flowList(i,3) =FlowWt(ii(i), jj(i) );
end
%
flowList3=flowList(:,3);
[degsort, degindex] = sort(flowList3); % sorted list & get new index
flowListSort=flowList(degindex,:);
flowListSort=flipud(flowListSort); % list are largest... smaller
sum(flowListSort(:,3)) % is 0.7669, ie. 77% of all flows
clear degsort degindex flowList3
%% 1.7.3 Flow: list top 10, 20, 30 of Deg, Flow
%NodeFlow= csvread('MBnodeflows.csv');
[degsort, degindex] = sort(NodeFlow); % sorted list & get new index
hifi = degindex(n-19:n); % list top 20, 30
hifi=flipud(hifi);
clear degsort degindex
fprintf('\n >> Marmoset Brain: hi (In-) InfoFlow nodes \n')
for ii=1:20 %length(hifi)
i=hifi(ii);
fprintf(' %4.0f %s %s %7.4e \n', i, Acrn{i}, Lobes{i}, NodeFlow(i) )
end
clear degsort degindex % use hifi
%% 1.7.3 Calc flows over wt links [or read prev calc - at 1.6.1 above]
fprintf('\n > Marmoset Brain, flows In: thru Nodes & over wt-links (V2b) \n')
FlowWt=zeros(n,n); % allocate array
A = Anew; % v2b wts
% length(find(A(:) ) ) % #links = 3474 ok
for ii=1:n
%OUT: list of 1st NN
% jj = find(A(ii,:)); % find only strong Out links; > this cut off (>0
%jj = find(A(ii,:)>=1); % find only Out links
jj = find(A(:, ii)>=1); % find only In links
wts = A(ii,jj); % check the i-{nn} Out weights
TotOutWt=sum(wts); % sum =1 ok
%wts = round(A(ii,jj)) % check the Out weights
for j=1:length(jj) % loop over NN of ii
jjj= jj(j);
FlowWt(ii,jjj)= NodeFlow(ii)*A(ii,jjj)/TotOutWt; % partiton out flow to this link
end % loop over linked nodes
end % loop over all nodes
sum(FlowWt(:) ) % = 1.0 ok
% csvwrite('Marmoset_FlowWt', FlowWt);
clear jj wts TotOutWt
% checks
[sum(FlowWt(:, 3)) sum(FlowWt(3, :))] % : 0.0172 0.0174 % sum in & out
NodeFlow(3) % 0.0174 ok
find( FlowWt == max(FlowWt(:)) )
[ii jj] = ind2sub(116, 4794) % is (38, 42)
% sort flows over links:
flowLinkList=adj2edgeL(FlowWt); % 16864 links (no self loops)
% csvwrite('MBrawLinkList.csv', rawLinkList); % save file
% sort by wts: weak to strong
wts=flowLinkList(:,3);
[wtsort, wtindex] = sort(wts);
flowLinkListSort=flowLinkList(wtindex,:); % now sorted by raw wt (incrs)
flowListSort=flipud(flowLinkListSort); % Decrs. order
% clear wts wtsort wtindex flowLinkList flowLinkListSort
%
fprintf('top 20 flows, NodeFlow(i:out)')
for ii=1:20
i=flowListSort(ii,1) ; j=flowListSort(ii,2) ;
fprintf(' %4.0f %s %s %7.4e %4.0f %s %s \n', i, Acrn{i}, Lobes{i}, NodeFlow(i), j, Acrn{j}, Lobes{j} )
end
%
fprintf('\n top 20 flows, i -> j: OUT links \n')
for ii=1:20
i=flowListSort(ii,1) ; j=flowListSort(ii,2) ;
fprintf(' %4.0f %s %s %7.6f %4.0f %s %s \n', i, Acrn{i}, Lobes{i}, flowListSort(ii,3), j, Acrn{j}, Lobes{j} )
end
% form dist matrix
%% 1.8 Calc i-j Distance Matrix - for links that exist
% needs Coords array - read up at #1.1, or here:
NodeCoord=csvread('CoordsMarmoset.csv'); % nb. some absent: NaN
Adist=csvread( 'MarmosetDistPairs.csv'); % read prev calc
% prev calc from bv.img [.nii] below at A.1.2 [line ~400]
fprintf('\n Calc dist array, from calc Atlas Coords (for links >=0) \n')
Adist=zeros(n); % initial;ise array: d(i,j) in mm - from Atlas area coords
for jjj=1:n
for kk=1:n
if Afull(jjj,kk)>0 % for non-zero wts only
xj=NodeCoord(jjj,1); yj=NodeCoord(jjj,2); zj=NodeCoord(jjj,3);
dxjk=abs(NodeCoord(kk,1)-xj); %col vec % nb. need to use actual x-coord, not "idx"
dyjk=abs(NodeCoord(kk,2)-yj); %
dzjk=abs(NodeCoord(kk,3)-zj);
drjk=sqrt(dxjk*dxjk+ dyjk*dyjk+ dzjk*dzjk); %clear dxs dys dzs
Adist(jjj,kk)=drjk;
end
end
end
clear xj yj zj dxjk dyjk dzjk drjk
length(find(Adist(:)>0)) % check 3223
%
% csvwrite( 'MarmosetDistPairs.csv',Adist); % save file
%
figure; hist(Adist(:), 100); hold on % ~ skewedGaussian, peaks at 7, 11, 13 mm
title('MarmosetBrain: Link distances','Fontname','Times New Roman','FontSize',12,'FontWeight','Bold');
axis([0.2 22.5 0 100])
xlabel('link length (mm)','Fontname','Times','FontSize',12,'FontWeight','Bold')
ylabel('counts','Fontname','Times New Roman','FontSize',12,'FontWeight','Bold'); hold off
% looks log Normal
%% 1.8.1 Calc LinkList of Wts (FLNe) & plot hist(log10)
fprintf('\n Calc LinkList from rescaled Adj-new (LNe) \n')
% LinkList=csvread('MarmosetLinkList.csv'); % i-j-wt; 6325 entries, incl many 0's
% Or, calc from Adj:
LinkList=adj2edgeL(Anew); % Rescaled; finds 3474 links for 116 x 166 Adj
% eliminate wt:0
posWts=find(LinkList(:,3)>0); % find the non-zero entries % 3474 of
LinkList=LinkList(posWts,:);
clear posWts
[nl ~] =size(LinkList) % now 3474 non-0 wts
% nb. Acrn-i, Acrn-j, wt] from .txt file & in LinkListACrn.csv
wt=LinkList(:,3);
figure; hist(wt,50); title('Marmoset rescaled LNe distribution ')
wtlog=log10(wt);
figure; hist(wtlog,50)
%axis([-5.3 0 0 150]) % omit large peak at "0")
title('Marmoset original log-10 wt (LNe) distribution ')
xlabel('log10 (weight, LNe)'); ylabel('Counts')
% hh =gcf; print(hh, 'FigS3.tif', '-dtiff' , '-r300');
% sort list of wt
wts=LinkList(:,3);
[wtsort, wtindex] = sort(wts);
LinkListSort=LinkList(wtindex,:); % now sorted by raw wt (incrs)
clear wts wtsort wtindex
% Top 40% wt>= 0.0033 ~10^-2.48 (#2093-3474: 1382 Links)
%LinkListT40=LinkList(2093:end, :);
%A40=edgeL2adj(LinkListT40); % no diag, Ok / missing 3 nodes?
%S=LinkList(:, 1); S=unique(S); % 116 ok
%T=LinkList(:, 2); T=unique(T); % 55, as before -
% clear LinkList
%% 1.8.2 Alt. Calc LinkList of Wts (FLNe) & plot hist(log10)
fprintf('\n Calc LinkList from orig Adj-full \n')
% LinkList=csvread('MarmosetLinkList.csv'); % i-j-wt; 6325 entries, incl many 0's
% Or, calc from Adj:
LinkList=adj2edgeL(Afull); % finds 3474 links for 116 x 166 Adj
% eliminate wt:0
posWts=find(LinkList(:,3)>0); % find the non-zero entries % 3474 of
LinkList=LinkList(posWts,:);
clear posWts
[nl ~] =size(LinkList) % now 3474 non-0 wts
% nb. Acrn-i, Acrn-j, wt] from .txt file & in LinkListACrn.csv
wt=LinkList(:,3);
figure; hist(wt,50); title('Marmoset original FLNe distribution ')
wtlog=log10(wt);
figure; hist(wtlog,50)
axis([-5.3 0 0 150]) % omit large peak at "0")
title('Marmoset original log-10 (FLNe) distribution ')
xlabel('log10 (weight, FLNe)'); ylabel('Counts')
% hh =gcf; print(hh, 'FigS3.tif', '-dtiff' , '-r300');
clear LinkList
% sort list of wt
wts=LinkList(:,3);
[wtsort, wtindex] = sort(wts);
LinkListSort=LinkList(wtindex,:); % now sorted by raw wt (incrs)
clear wts wtsort wtindex
% Top 40% wt>= 0.0033 ~10^-2.48 (#2093-3474: 1382 Links)
%LinkListT40=LinkList(2093:end, :);
%A40=edgeL2adj(LinkListT40); % no diag, Ok / missing 3 nodes?
%S=LinkList(:, 1); S=unique(S); % 116 ok
%T=LinkList(:, 2); T=unique(T); % 55, as before -
%% need to do long-hand
A40=zeros(116); % set up array
for i=1:length(LinkListT40)
A40(LinkListT40(i, 1), LinkListT40(i, 2)) = LinkListT40(i, 3);
end
%% 1.8.3 Calc Dist; & wt vs Dist
% LinkList calc at #
DistCol=zeros(length(LinkList),1);
for i=1:nl
dist=( (NodeCoord(LinkList(i,1),1)- NodeCoord(LinkList(i,2),1) )^2 ...
+ (NodeCoord(LinkList(i,1),2)- NodeCoord(LinkList(i,2),2) )^2 ...
+ (NodeCoord(LinkList(i,1),3)- NodeCoord(LinkList(i,2), 3) )^2 );
DistCol(i)=sqrt(dist);
end
clear i dist
LinkListDist=[LinkList DistCol];
clear DistCol LinkList
% sort list by wt
wt=LinkListDist(:,3);
[wtsort, wtindex] = sort(wt);
LinkListSortWt=LinkListDist(wtindex,:); % now sorted by raw wt (incrs);
% nb 2851 wt entries are 0 ?? ; min ~ 4.59e-6
clear wtsort wtindex wt
% LinkList92=LinkListSortWt(284:end,:); % top 92% : wt> 5e-5
% other
DistCol=LinkListSortWt(:,4);
wts=LinkListSortWt(:,3);
wtslog=log10(wts);
mean(wtslog) % now -2.80
% Next: sort list by dist
DistCol=LinkListDist(:,4);
[wtsort, wtindex] = sort(DistCol);
LinkListSortDist=LinkListDist(wtindex,:); % now sorted by d(i-j) (incrs);
% max Dist is 22.65mm; nb #5941:6325 are NaN :no Injn ??
clear wtsort wtindex % DistCol {is needed for plot}
% now in dist-sorted order
DistCol=LinkListSortDist(:,4);
wts=LinkListSortDist(:,3);
wtslog=log10(wts);
mean(wtslog) % now -2.80
% length(find(wtslog < -100) ) % man wt: 0 > -Inf for log
% log10(mean(wts)) % -2.79 : too many 0!!
% wtsnz=wts(wts>0); log10(mean(wtsnz)) % 3474 non-zero elements; mean -1.80
% eliminate wt:0 % Dist: nan
posWts=find(LinkListDist(:,3)>0); % find the non-zero entries % 3474 of
LinkListDist=LinkListDist(posWts,:);
clear posWts
goodDist=find(LinkListDist(:,4) < 30); % 3223 of
% max(LinkListDist(:,4)) % 21.04 mm
LinkListDist=LinkListDist(goodDist,:);
clear goodDist
% save for InfoMap calc
%dlmwrite( 'MarmosetPairList.txt',LinkList,'Delimiter',' '); % InfoMap code wants .txt , space delim
%LinkListSmall=LinkListDist(:,1:3);
%dlmwrite( 'MarmosetSmallPairList.txt',LinkListSmall,'Delimiter',' ');
%
%% 1.9 Investigate node vol effects on FLNe
LinkListVol=zeros(length(LinkList),5);
LinkListVol(:,1) = LinkList(:,1); LinkListVol(:,3) = LinkList(:,2); % Load Source, Target IDs
LinkListVol(:,5) = LinkList(:,3); % and wts
for i=1:length(LinkList)
LinkListVol(:,2) = NodeVols(LinkList(:,1)); % Source vol
LinkListVol(:,4) = NodeVols(LinkList(:,2)); % Target vol
end
%
wtslog=log10(LinkListVol(:,5));
%volslog=log10(LinkListVol(:,2)) %+log10(LinkListVol(:,4)); % log(Vs*Vt)
volslog=log10(LinkListVol(:,4));
figure; hold on; %plot(LinkListVol(:,2), LinkListVol(:,5), '.')
%plot(volslog, wtslog, '.') % log-lin plot
%plot(LinkListVol(:,4), wtslog, '.') % log-lin plot
plot(volslog, wtslog, '.') % log-log plot
% title('Marmoset, compare FLNe vs Vol-Source')
%title('Marmoset, compare FLNe vs Vol-Target')
title('Marmoset, compare FLNe vs Vols-Source*Target')
xlabel('log(vol-T) (mm^3)') %xlabel('vol-Source (mm^3)')
ylabel('log(FLNe)')
%% 1.2.1 Plot LinkWts, LNe or FLNe vs Dist {for exp
% calc DistCol at # 1.8.3, from LinkList [calc at # 1.8.1 (LNe) or 1.8.2 (FLNe)
fprintf('\n Calc LinkList from rescaled Adj-new (LNe) \n')
figure; %figure('position',figposition, 'units','pixels'); hold on; % big
plot(DistCol, wtslog, '.', 'MarkerSize', 16); hold on
title('Marmoset Brain log10(LNe) vs ListDist(mm) ')
%title('Marmoset Brain, link weight vs dist ','Fontname','Times New Roman','FontSize',14,'FontWeight','Bold')
xlabel('dist(mm)','Fontname','Times','FontSize',14,'FontWeight','Bold')
ylabel('log10(LNe wt) ','Fontname','Times New Roman','FontSize',14,'FontWeight','Bold');
% highlight weakest wt (< 9e-6)
for i=1:17 % fitst 17 pts in wt-ordered list
plot(DistCol(i), wtslog(i), '.', 'Color', 'red','MarkerSize', 16)
end
tmp=LinkListSortWt(1:17,:);
% save S-T
tmpS=cell(17,1);
for i=1:17
tmpS=Acrn{tmp(i,1)};
tmpT=Acrn{tmp(i,2)};
[tmpS tmpT]
end
clear tmp
% highlight next group of weaker, wt (< 5e-5)
for i=18:283 % fitst 17 pts in wt-ordered list
plot(DistCol(i), wtslog(i), '.', 'Color', [1, 0.65, 0],'MarkerSize', 16) % orange
end
% top 40%, by wt is at cutoff 10^0.0033 : 1.0076
% ie. rows 2085:3474 : ie. 1390 links
% drop 10% - ie. top 90% is rows 348:3474
%% & Curve Fitting tool
% use Matlab Curve Fitting Toolbox
cftool % invokes GUI, import data
% ** refit - so edit below **
% linear fit, to logwt: -2.019 - 0.0863*dist ; mean 0 about that line; std : xx
plot([0 22], [-2.019 -3.9176], 'k-', 'LineWidth', 2, 'LineSmoothing', 'on') % along the mean of the data
% & 3*sigma CI: nb. log scale
plot([0 22], [-2.019+3*0.967 -3.9176+3*0.967], 'r--', 'LineWidth', 1.5, 'LineSmoothing', 'on') % 3*sigma
plot([0 22], [-2.019-3*0.967 -3.9176-3*0.967], 'r--', 'LineWidth', 1.5, 'LineSmoothing', 'on')
% get residuals - in log-space
[nll ~]=size(LinkListSortDist);
resLogWt=zeros(nll,1);
resLogWt(:)=wtslog(:) - (-2.019 - 0.0863*DistCol(:) );
% [mean(resLogWt) std(resLogWt)] % -4.6944e-04, 0.9665
% short < 10mm
% [mean(resLogWt(1:2116)) std(resLogWt(1:2116)) ]
% short < 10mm
% [mean(resLogWt(2117:end)) std(resLogWt(2117:end)) ]
% explore wt cut-off
% log-wt < 1e-5
%mean(LinkListSortWt(1:17, 4)) % link dist: 13.5
% log-wt < 2e-5
%mean(LinkListSortWt(1:65, 4)) % link dist: 11.23
% log-wt < 3e-5
% mean(LinkListSortWt(1:129, 4)) % link dist: 10.7
% residuals - in original-space (FLNe)
% reswt=10.^resLogWt;
%% 1.2.1 Plot LinkWts, LNe or FLNe vs Dist {for exp fit
% calc DistCol at # 1.8.3, from LinkList [calc at # 1.8.1 (LNe) or 1.8.2 (FLNe)
fprintf('\n Plot Wt-Dist from rescaled Adj-new (LNe) \n')
figure; %figure('position',figposition, 'units','pixels'); hold on; % big
plot(DistCol, wtslog, '.', 'MarkerSize', 16); hold on
title('Marmoset Brain log10(LNe) vs ListDist(mm) ')
%title('Marmoset Brain, link weight vs dist ','Fontname','Times New Roman','FontSize',14,'FontWeight','Bold')
xlabel('dist(mm)','Fontname','Times','FontSize',14,'FontWeight','Bold')
ylabel('log10(LNe wt) ','Fontname','Times New Roman','FontSize',14,'FontWeight','Bold');
%% highlight weakest wt (< 9e-6) : for FLNe
for i=1:17 % fitst 17 pts in wt-ordered list
plot(DistCol(i), wtslog(i), '.', 'Color', 'red','MarkerSize', 16)
end
tmp=LinkListSortWt(1:17,:);
% save S-T
tmpS=cell(17,1);
for i=1:17
tmpS=Acrn{tmp(i,1)};
tmpT=Acrn{tmp(i,2)};
[tmpS tmpT]
end
clear tmp
% highlight next group of weaker, wt (< 5e-5)
for i=18:283 % fitst 17 pts in wt-ordered list
plot(DistCol(i), wtslog(i), '.', 'Color', [1, 0.65, 0],'MarkerSize', 16) % orange
end
% top 40%, by wt is at cutoff 10^0.0033 : 1.0076
% ie. rows 2085:3474 : ie. 1390 links
% drop 10% - ie. top 90% is rows 348:3474
%% 1.2.1a Log-Log Plot LinkWts, LNe or FLNe vs Dist {for power law fit
% nb. LinkListDist(:,3) is wt; LinkListDist(:,4) is dist
% calc DistCol at # 1.8.2or3, from LinkList [calc at # 1.8.1 (LNe) or 1.8.2 (FLNe)
% wtslog=log10(wts); mean(wtslog) % now -2.80
fprintf('\n Plot Wt-Dist from rescaled Adj-new (LNe) \n')
% figure; plot(DistCol, wts) % debug
DistLn=log(DistCol); WtsLn =log(wts);
figure; %figure('position',figposition, 'units','pixels'); hold on; % big
% plot(DistCol, wtslog, '.', 'MarkerSize', 16); hold on % was log 10
%loglog(DistCol, wts); hold on
plot(DistLn, WtsLn, '.', 'MarkerSize', 16); hold on % now ln = log_e
%xlim([0.5 25])
title('Marmoset Cortex log_e(LNe) vs log_e(ListDist (mm)) ')
%title('Marmoset Cortex, link weight vs dist ','Fontname','Times New Roman','FontSize',14,'FontWeight','Bold')
xlabel('log_e(dist) (mm)','Fontname','Times','FontSize',14,'FontWeight','Bold')
ylabel('log_e(LNe wt) ','Fontname','Times New Roman','FontSize',14,'FontWeight','Bold');
xlim([-0.2 3.25])
% cftool: GUI for curve fitting
line ([0 3], [6.9 (6.9 - 2.0*3) ], 'Color', 'k', ...
'LineWidth', 1.5); % linear fit 6.9 -2.0*d(mm)
%% 1.2.1b Plot LinkWts, LNe vs Dist {+ exp fit)
% calc DistCol at # 1.8.3, from LinkList [calc at # 1.8.1 (LNe) or 1.8.2 (FLNe)
fprintf('\n Plot Wt-Dist from rescaled Adj-new (LNe) \n')
figure; %figure('position',figposition, 'units','pixels'); hold on; % big
plot(DistCol, WtsLn, '.', 'MarkerSize', 16); hold on
title('Marmoset Cortex: log_e(LNe) vs ListDist(mm) ')
%title('Marmoset Cortex, link weight vs dist ','Fontname','Times New Roman','FontSize',14,'FontWeight','Bold')
xlabel('dist(mm)','Fontname','Times','FontSize',14,'FontWeight','Bold')
ylabel('log_e(LNe, wt) ','Fontname','Times New Roman','FontSize',14,'FontWeight','Bold');
line ([1 24], [4.77- 0.219 5.77 - 0.219*24], 'Color', 'k', ...
'LineWidth', 1.5); % linear fit 4.77 - 0.219d
%% 1.3 node wt-Degree / Strength :: nb. numbering out of order
DegIn=sum(Afull,1)'; % i<-j; col sum (1st index): produces a row % show as Col
DegOut=sum(Afull,2); % i->j; row sum
%sum(Afull(:,1)) % check : these are "1" ??
figure; hist(DegIn,50)
title('Marmoset Brain, wt Degree-In, histogram','FontSize',14)
% get stats on Deg.
binrange=linspace(0,3.4E4,100); % get range from hist plot already
bincounts=histc(DegIn,binrange);
figure; bar(binrange, bincounts,'histc')
% appears to be exp fall off up to ~ 4xE7; only greater; & are 1 off
figure; semilogy(binrange, bincounts, '.','MarkerSize',14)
hold on; title('Marmoset Brain, wt Degree-In, log histogram','FontSize',14)
xlabel('Degree (bins)','FontSize',14)
ylabel('log counts','FontSize',14); hold off
% clear binrange bincounts
%
[degsort, degindex] = sort(DegIn); % sorted list & get new index
hdin = degindex(end-9:end); % list top 10, Hi Deg nodes (Deg = 113 .. 16 )
hdin=flipud(hdin); % list are largest... smaller
%mdi = degindex(n-43:n-20)'; % next 24, Medium deg (Deg = 15 .. 10 )
% list top 10, 20
md=mean(DegIn)
sd=std(DegIn)
for ii=1:length(hdin)
i=hdin(ii);
fprintf(' %4.0f %s %7.3f %5.2f \n', i, Acrn{i}, DegIn(i), (DegIn(i)-md)/sd )
end
clear md sd
%% 1.3a #Links: un-wt Deg - binarise full 116x116 A: UN-wt links
[sii sjj v]=find(Afull);
Alogical=full(sparse(sii, sjj, 1)); % "mask" of selected entries, to re-form sq array of 1's.
% nb. need to use sparse to form the array
A1 = Alogical.*1; % pick out the selected entries; yields a "full" array
%figure; spy(A1)
clear Alogical sii sjj v
length(find(A1(:))) % check 3474 links:
% alt calc
%Afull= csvread('AdjFullMarmoset.csv'); % was calc from expt data
LinkList=adj2edgeL(Afull); %LinkList=adj2edgeL(Afull);
% use "filtered" LinklList (of 3462 links ~0)
% LinkList=csvread('MarmosetLinkList.csv'); % i-j-wt
[nl ~]=size(LinkList)
wt1=ones(nl);
LinkList1=[LinkList(:,1:2) wt1]; % append wt=1
A1=edgeL2adj(LinkList1);
% nb. need to use sparse to form the array
%figure; spy(A1)
%
clear wt1 LinkList LinkList1
% figure; spy(A1) % 1948 entries: banded & v sparse
length(find(A1(:)) )
% count links (<~ 50-150)
DegIn1=sum(A1,1)'; % i<-j; col sum (1st index): produces a row % show as Col
DegIn1=[DegIn1; 0; 0]; % append 2 missing 0's
DegOut1=sum(A1,2); % i->j; row sum
figure; hist(DegIn1, 50)
title('Marmoset DegIn of binarised links')
axis([0 110 0 10]) % nb. big # at 0: not sampled
%
figure; hist(DegOut1, 50); hold on
title('Marmoset DegOut of binarised links')
%axis([0 80 0 130]);
figure; stem (DegIn1, '.'); hold on
stem (DegOut1, 'r.');
legend('DegIn1', 'DegOut1')
title('Marmoset Deg of binarised links')
% all
figure; hist(DegIn1, 50); xlim([2 Inf]); hold on;
hist(DegOut1, 50, 'Color', 'r') % magenta
ylim([0 20]);
% all - joint p[lot
figure; hist([DegIn1, DegOut1], 50); hold on % Fails, diff # in & out
title('Marmoset Distribution of in- and out-links'); legend('k-In', 'k-Out')
%axis([1 110 0 20]) % omit large peaks at 0
%% plot & fit distrn
DegIn11=DegIn1(find(DegIn1 >0)); % need non-zero, eg for logNorm
DegOut11=DegOut1(find(DegOut1 >0));
figure; hist([DegIn11, DegOut11], 50); hold on % Fails, diff # in & out
title('Marmoset Distribution of in- and out-links'); legend('k-In', 'k-Out')
axis([1 110 0 20]) % omit large peaks at 0
clear DegIn11 DegOut11
% csvwrite('Marmoset_Adj1unwt.csv', A1); % save binarises Adj (UN-wt Afull)
%% 1.3a.1 probe missing In-link info
linkRatio=DegIn1./DegOut1;
% figure; stem(linkRatio, '.'); xlabel('Node # (all nodes')
title('Marmoser (LNe) k-in/ k-out');
figure; plot(linkRatio, '.','MarkerSize',18) % lin fit dominated by 0's
xlabel('Node # (all nodes'); title('Marmoser (LNe) k-in/ k-out');
linkRatNon=linkRatio(find(linkRatio)); % get the 55 non-zero ratios
figure; plot(linkRatNon, '.','MarkerSize',18); xlabel('Node # (injection sites only'); hold on
title('Marmoset: #linksIn/#Out (non-0) ');
text(11+1, 4.25, Acrn{11}); text(32+1, 4.68, Acrn{32}); text(33+1, 6.33, Acrn{33}); % label outliers
figure; hist(linkRatNon, 50); % skewed normal
clear linkRatio linkRatNon
%% 1.3.2 wt/link distrn
% wt/Link-In, per anatomical area (node)
WtLinkIn=DegIn./DegIn1;
length(find(isnan(WtLinkIn))) % check on 0 #Links:
WtLinkIn(find(isnan(WtLinkIn)))=0; % omit the 0/0: NaN
figure; hist(WtLinkIn, 50); title('Marmoset, Distribution of Wt-Deg-In/ k-In ')
figure; %stem(WtLinkIn)
plot(WtLinkIn , '.','MarkerSize',12); title('Marmoset, Wt-Deg-In/ k-In vs Node-ID '); hold on
% this plot shows av ~ xx wt/Link;
[mean(WtLinkIn) mode(WtLinkIn) std(WtLinkIn)] % finds mode=low, since exp distrn
length(find(WtLinkIn==0)) % 61 of; only 55 targets measured
WtLinkIn1=WtLinkIn(find(WtLinkIn >0)); % need non-zero, eg for logNorm
% typical range is xx wt/Link
% wt/Link-Out, per area (node)
WtLinkOut=DegOut./DegOut1;
length(find(isnan(WtLinkOut))) % check on 0 #Links:
WtLinkIn(find(isnan(WtLinkOut)))=0; % omit the 0/0: NaN
figure; hist(WtLinkOut, 50); title('Marmoset, Distribution of Wt-Deg-Out/ k-Out ')
% use dfittool % title('Marmoset, Distribution of Weight-In/link-In ')
figure; %stem(WtLinkOut)
plot(WtLinkOut , '.','MarkerSize',12); title('Marmoset, Wt-Deg-Out/ k-Out vs Node-ID '); hold on
% this plot shows av ~ xx wt/Link;
[mean(WtLinkOut) mode(WtLinkOut) std(WtLinkOut)] % finds mode=xxx, since exp distrn
% typical range is xx wt/Link
%% 1.4 form other LinkLists, with CutOffs, from Adj - for InfoMap calc
% cf noTarget & noTarget lists at A.2.1, line 785.
% for Ctx: 55 x 55 Source AND Target sites
% cf. A.2.1, below, to find many missing Targets (116-52 of)
% ie. omit Cols 4-9, etc
% nb. this discards a few Sources also: ie. Rows 4-9, ...
LinkList55=adj2edgeL(A);
length(find(A(:))) % 1854 links
dlmwrite( 'MarmosetPairList52.txt',LinkList55,'Delimiter',' '); % needs space deliminter
% un-wt
tmp=ones(length(LinkList55),1);
LinkList55(:,3)=tmp;
clear tmp
dlmwrite( 'MarmosetPairList55u.txt',LinkList55,'Delimiter',' '); % needs space deliminter
% check A-unwt
Au=adjLu2adj(LinkList52); % un-wt
figure; spy(Au)
title('Marmoset: Adj 55x55 un-wt')
length(find(Au(:)) ) % 1854, ok
% count links (<~ 50)
DegIn1=sum(Au,1)'; % i<-j; col sum (1st index): produces a row % show as Col
DegOut1=sum(Au,2); % i->j; row sum
figure; hist(DegIn1)
title('Marmoset 55x55 DegIn of binarised links')
figure; hist(DegOut1); hold on
title('Marmoset 55x55 DegOut of binarised links')
%axis([0 80 0 130]);
figure; stem (DegIn1, '.'); hold on
stem (DegOut1, 'r.');
legend('DegIn1', 'DegOut1')
title('Marmoset 55x55 Deg of binarised links')
%% 2.0 LNe,1: raw data, Test: normalised to 1 mm^3 target vol
% 55 injn sites (targets)
% cf XL workings %& calc: Marmoset_raw1_2_sort.xlx
LNe1=csvread('LNE1.csv');
wt=LNe1(:,3);
figure; hist(wt,50); title('Marmoset original LNe1 distribution ')
wtlog=log10(wt);
figure; hist(wtlog,50)
%axis([-5.3 0 0 150]) % omit large peak at "0")
title('Marmoset original log-10 (LNe1) distribution ')
%% 2.1 calc LN from LNe link list & LNt,e,i data; cf #4.1 below
fprintf('\n Calc #LN(i-j) from raw data \n')
% data is [ Node Id, TotLN, LNe, LMi], with ~2-9 repeats
% calc using the average LN & LNi for the 55 targets - in cols [5, 6]
% as calc by hand in XL. see Marmoset_raw1_2_sort.xlsx : combines 2 data sets
rawLN1=csvread('rawLN1.csv');
LNtot1=sum(rawLN1(:,5)); % tot of all mean(LN-repeats); now 905k
FLNdenom = zeros(116,2);
FLNav = zeros(116,3); % to store av values [FLNtot, FLNe FLNi] for ea target
count=1; % nb. need to skip over repeats - just get the av.
for i =1:143
if rawLN1(i,5) > 0 % then av for these Targets recorded
ii=rawLN1(i,1); % index recorded on col 1 (w repeats)
[ii rawLN1(i,6)] % debug
FLNdenom(ii,1)=ii; % target index [only 55 of the 116]
FLNdenom(ii,2)=(LNtot1-rawLN1(i,6)); % av LN, after repeats
FLNav = [rawLN1(i,5), rawLN1(i,6), (rawLN1(i,5)-rawLN1(i,6))]; % calc FLNe also
count=count+1;
end
end
count % debug - check
clear count
sum(FLNdenom(:,2))/55 % to get mean
[min(FLNdenom(:,2)) sum(FLNdenom(:,2))/55 max(FLNdenom(:,2))]
% variation about mean: (9.05-9.01)/9.01 : 0.0044; ie. 0.4% : Trivial !
%
figure; hold on; plot(FLNdenom(:,2), '.', 'MarkerSize', 11);
% axis([0 150 2.13e6 2.18e6])
title('Marmoset Brain: Tot#LN - LNi (average at 55 targets) ')
ylabel('All(LN) - LNi '); xlabel('Node ID # ');
%% 2.2 scan LinkList & undo normalisation by denom: get LN, Anew
% LinkList of [i j Wt] calc at #1.8.2 [line 322], above, from raw wts
[nl ~] = size(LinkList)
extraCol=zeros(nl,1); % append to the 3-col LinkList
%LinkList= [LinkList, extraCol]; % 4th col wil be LNe(i-j); 5th LNe/Tot#LN
for i=1:nl
LinkList(i,4)=LinkList(i,3)*FLNdenom(LinkList(i,2),2); % * denom(target)
LinkList(i,5)=LinkList(i,3)*FLNdenom(LinkList(i,2),2)/LNtot1; % and as frn of tot#LN
end
clear extraCol
figure; plot(LinkList(:,3), LinkList(:,5), '.', 'MarkerSize', 11); % good st line!!
% calc new Adj - need 3 cols in link list now:
LinkListNew=[LinkList(:,1:2) LinkList(:,4)]; % is [i j LN(i-j)/TotLN ]
Anew=zeros(n);
for i=1:nl
Anew(LinkListNew(i,1), LinkListNew(i,2))=LinkListNew(i,3);
end
%% 2.3 new (LNe) wt Deg-In,Out (Strength)
DegIn=sum(Anew,1)'; % i<-j; col sum (1st index): produces a row % show as Col
DegOut=sum(Anew,2); % i->j; row sum
figure; hist(DegIn); title('Marmoset rescaled wt-DegIn ')
figure; hist(DegOut); hold on; title('Marmoset rescaled wt-DegOut '); hold off
figure; plot(DegOut, DegIn, '.', 'MarkerSize', 9 ); title('marmoset: rescaled DegOut vs DegIn ')
figure; stem (DegIn, '.'); hold on % ~ "wt-DegIn"