-
Notifications
You must be signed in to change notification settings - Fork 2
/
dbs_analysis_subject_Q.m
170 lines (127 loc) · 5.64 KB
/
dbs_analysis_subject_Q.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
%DBS Individual Subject Connectomics: rapid Creates a DBS connectome using
%MSN's (no figures / imags)
%
% Does Morphometrics Similarity Networks (MNS's) based on Freesurfer
% reconstructions of clinical DBS data
%
% Based on paper and primary research by:
% Siedlitz et al, Bioarchiv 2017, doi: 10.1101/135855
%
% As per dbs_analysis_subject_F.m but without individual metrics or images
%
% Dependencies: DBS_connectome toolbox (includes: importfile, make_zscores, make_networks, DBS_link_viewer, DBS_hub_viewer)
% BCT, L1precision, pwling, fdr_bh
%
% Inputs: Freesurfer .annot file
% Uses standard BrainNet design and Freesurfer xyz
% co-ordinates
%
% NB: need to set subject at initiation
%
% Outputs: BrainNet views (.tiff)
% Custom connectome views of hubs and links (.tiff)
%
% NB: saves all output files in working directory under standard names
%
% Version: 1.0
%
% Michael Hart, University of Cambridge, May 2017
% Overview
% 1. Take 7 metrics from .annot file
%.2. Add metrics to matlab variable
% 3. Transform to z-scores
%.4. Correlate measures across regions
% 5. Threshold
% Pipeline
% 1. Define subject directory
% 2. Compute & load data
% 3. Make matric table
% 4. Transform to z-scores
% 5. Make networks
% 6. Compare networks
% 7. Threshold
% Michael Hart, University of Cambridge, June 2018
%% 0. Choose subject
%clear variables
%clear global
clc
%This should be the only bit that needs changed for each subject
%mysubject = 'stn_20171127';
%% 1. Define subject directory
mysubjectpath = sprintf('/home/mgh40/scratch/functional/all/%s', mysubject);
cd(mysubjectpath); %stores all files here
mkdir dbs_connectome
%% 2 . Create & load data from .annot file
%first print outputs to stats files
setenv('SUBJECTS_DIR', mysubjectpath);
system('mris_anatomical_stats -a ./FS/label/lh.aparc.annot -b /FS lh > dbs_connectome/lh_msn_stats.txt');
system('mris_anatomical_stats -a ./FS/label/rh.aparc.annot -b /FS rh > dbs_connectome/rh_msn_stats.txt');
cd dbs_connectome
%path to file
lh_annot_file = strcat(mysubjectpath, '/dbs_connectome/lh_msn_stats.txt');
rh_annot_file = strcat(mysubjectpath, '/dbs_connectome/rh_msn_stats.txt');
%load data using custom loader
lh_msn_stats = importfile(lh_annot_file);
rh_msn_stats = importfile(rh_annot_file);
%% 3. Add metrics identifiers & make to single table
%scores
scores = {'number_vertices'; 'total_surface_area'; 'total_grey_matter'; 'average_cortical_thickness'; 'cortical_thickness_SD'; 'integrated_rectified_mean_curvature'; 'integrated_rectified_Gaussian_curvature'; 'folding_index'; 'intrinsic_curvature'};
%parcels
parcels = {'bankssts'; 'caudalanteriorcingulate'; 'caudalmiddlefrontal'; 'cuneus'; 'entorhinal'; 'fusiform'; 'inferiorparietal'; 'inferiortemporal'; 'isthmuscingulate'; 'lateraloccipital'; 'lateralorbitofrontal'; 'lingual'; 'medialorbitofrontal'; 'middletemporal'; 'parahippocampal'; 'paracentral'; 'parsopercularis'; 'parsorbitalis'; 'parstriangularis'; 'pericalcarine'; 'postcentral'; 'posteriorcingulate'; 'precentral'; 'precuneus'; 'rostralanteriorcingulate'; 'rostralmiddlefrontal'; 'superiorfrontal'; 'superiorparietal'; 'superiortemporal'; 'supramarginal'; 'frontalpole'; 'temporalpole'; 'transversetemporal'; 'insula'};
lh_parcels = cell(length(parcels), 1);
for i = 1:length(parcels)
grot = 'lh_';
lh_parcels{i} = sprintf('%s %s', grot, parcels{i});
end
rh_parcels = cell(length(parcels), 1);
for i = 1:length(parcels)
grot = 'rh_';
rh_parcels{i} = sprintf('%s %s', grot, parcels{i});
end
parcels = [lh_parcels; rh_parcels];
%combine (first 34 are left, next 34 are right)
msn_scores = [lh_msn_stats; rh_msn_stats];
msn_scores_table = array2table(msn_scores, 'VariableNames', scores, 'RowNames', parcels);
writetable(msn_scores_table, 'table_msn_scores.txt', 'delimiter', 'tab');
%% 4. Transform to z-scores
%subtract mean and divide by standard deviation
%creates zero mean and unit SD
[msn_zscores] = dbs_make_zscores(msn_scores);
%% 5. Make networks
msn_zscores = msn_zscores'; %flip to conventional format (scores x nodes)
msn_zscores(5, :) = []; %remove certain scores
msn_zscores(1, :) = []; %remove certain scores
%Pearson (1), Partial(2), L1(3), L2(4)
msn_networks = dbs_make_networks(msn_zscores); %approximately 1 minute
%% 6. Compare networks
%some quality control
nNetworks = size(msn_networks, 3);
R = zeros(nNetworks, 1); pos = zeros(nNetworks, 1);
maxmin = zeros(nNetworks, 1);
for iLevel = 1:nNetworks
%correlation
R(iLevel, 1) = mean(mean(triu(msn_networks(:, :, iLevel))));
%negative correlations
pos(iLevel, 1) = nnz(triu(msn_networks(:, :, iLevel))>0);
%range
maxmin(iLevel, 1) = range(range(triu(msn_networks(:, :, iLevel))));
end
network_stats = [R pos maxmin];
network_check_codes = {'Mean_correlation'; 'Negative_correlations';
'Range_of_correlations'};
network_types = {'Pearson'; 'Partial'; 'L1'; 'L2'};
%write table
network_stats_table = array2table(network_stats, 'VariableNames', network_check_codes, 'RowNames', network_types); %Matlab R2015 onwards
%% 7. Threshold
nNodes = size(msn_networks, 2);
%FWER with Bonferroni
fwer_thresh = 0.05 ./ (nNodes * nNodes);
[~, grot] = corr(msn_zscores); %get pvalues for Pearson net
fwer_thresh_mask = grot < fwer_thresh; %binarise if pvalue better than threshold
pearson_net_fwer_thresh = msn_networks(:, :, 1) .* fwer_thresh_mask; %1 per cent cost
%FDR
[fdr_thresh_mask, crit_p, adj_p] = fdr_bh(grot, 0.05, 'pdep', 'yes');
pearson_net_fdr_thresh = msn_networks(:, :, 1) .* fdr_thresh_mask; %~20 per cent cost
%% Close up
filename = sprintf('DBS_analysisQ_%s%s', mysubject, '.mat');
save(filename);