Correlation of thalamic functional organization disturbances and genetic architecture in motor subtypes of Parkinson's disease.

1. %%step1
2. %individual matrix
3. clc;clear;
4. Timepoints=200; % change based on your data
5. inputall=[uigetdir(pwd,'please choose analysis file') filesep];
6. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
8. pathdata=['F:\FunImgTARWDCF_PPMI\20240807\PPMI\TD' filesep]; % input data
9. outputpath_affinity=[inputall 'affinity_matrix' filesep];%path of output for affinity_matrix
10. outputpath_others=[inputall 'intermediate_files' filesep];% intermediate_files
11. outputpath_emb=[inputall 'emb' filesep];%path of output for emb
12. outputpath_res=[inputall 'res' filesep];%path of output for res
13. outputpath_res_exp=[inputall 'exp' filesep];%path of output for exp
14. outputpath_nii=[inputall 'result_nii' filesep];%path of output for nii
15. outputpath_aligned=[inputall 'alignedgradient' filesep];%path of output for aligned
16. outputpath_groupresults=[inputall 'groupresults' filesep];%path of output for aligned
17. mask_A=[inputall 'mask' filesep 'Tha.nii'];% A mask
18. mask_B=[inputall 'mask' filesep 'Tha_without.nii'];% B mask
19. mkdir(outputpath_affinity,'A_B');% to store affinity matrix of A_B
21. filenames=dir(pathdata);
22. filenames(1:2,:)=[];
24. %connectivity
25. V_A=spm_vol(mask_A);
26. Y_A=spm_read_vols(V_A);
27. V_B=spm_vol(mask_B);
28. Y_B=spm_read_vols(V_B);
29. all_subjects_A_timeseries=zeros(Timepoints*size(filenames,1),size(find(Y_A>0),1));
30. all_subjects_B_timeseries=zeros(Timepoints*size(filenames,1),size(find(Y_B>0),1));
31. % affinity matrix for single subject
32. for m=1:size(filenames,1)
33. filenames11=[pathdata filenames(m).name filesep filenames(m).name 'F.nii'];
34. V1=spm_vol(filenames11);
35. Y1=spm_read_vols(V1);
36. data_R=reshape(Y1,size(Y1,1)*size(Y1,2)*size(Y1,3),size(V1,1));
37. clear filenames11 V1 Y1;
38. data_A=data_R(Y_A>0,:)';
39. data_B=data_R(Y_B>0,:)';
40. clear data data_R;
41. data_N_A=zscore(data_A);
42. data_N_B=zscore(data_B);
43. all_subjects_A_timeseries(((m-m)*size(data_A,1)+1):1*size(data_A,1),:)=data_N_A;
44. all_subjects_B_timeseries(((m-m)*size(data_B,1)+1):1*size(data_B,1),:)=data_N_B;
45. ['finished_affnity' num2str((m)) '/' num2str(size(filenames,1))]
46. %% find zeros
47. number_of_nan0_A=zeros(1,size(all_subjects_A_timeseries,2));
48. for i=1:size((all_subjects_A_timeseries),2);
49. number_of_nan0A=length(find(all_subjects_A_timeseries(:,i)~=0));
50. number_of_nan0_A(1,i)=number_of_nan0A;
51. end
52. number_of_nan0_B=zeros(1,size(all_subjects_B_timeseries,2));
53. for k=1:size((all_subjects_B_timeseries),2);
54. number_of_nan0B=length(find(all_subjects_B_timeseries(:,k)~=0));
55. number_of_nan0_B(1,k)=number_of_nan0B;
56. end
57. if length(find(number_of_nan0_A==0))==0 && length(find(number_of_nan0_B==0))==0;
58. ['there is no all zero data']
59. else
60. ['there has all zero data']
61. end
62. clear k i number_of_nan0_A number_of_nan0_B number_of_nan0A number_of_nan0B
63. con_w3=corr(all_subjects_A_timeseries, all_subjects_B_timeseries);
64. con_w3= atanh(con_w3);
65. %con_w3= fisherz(con_w3);
66. con_w3=reshape(con_w3,size(find(Y_A>0),1),size(find(Y_B>0),1));
67. bb=prctile(con_w3',90);% top 10%
68. for kk=1:size(con_w3,1)
69. con_w3(kk,con_w3(kk,:)<bb(1,kk))=0;
70. end
71. con_w3(con_w3<0)=0;
72. clear kk bb;
73. Cosine_w=1-pdist(con_w3,'cosine');
74. index = 1:size(con_w3,1)*size(con_w3,1);
75. index = reshape(index,size(con_w3,1),size(con_w3,1))';
76. index = triu(index,1);
77. index=index';
78. index=index(:)';
79. index(index==0)=[];
80. Cosine_w_F = zeros(size(con_w3,1),size(con_w3,1));
81. Cosine_w_F(index) = Cosine_w;
82. Cosine_w_F=Cosine_w_F+Cosine_w_F';
83. for kk=1:size(con_w3,1)
84. Cosine_w_F(kk,kk)=1;
85. end
86. clear kk;
87. save([outputpath_affinity 'A_B' filesep filenames(m).name 'Ave_A_B'],'Cosine_w_F');
88. clear index con_w3 Cosine_w Cosine_w_F;
89. clear all_subjects_A_timeseries all_subjects_B_timeseries;
90. end
92. %%Éú³Éÿ¸ö±»ÊÔµÄÌݶȾØÕó
93. root_dir = 'E:\tha_gradient\Gradient_thalamus\TD1\affinity_matrix\A_B';
94. out_path = 'E:\tha_gradient\Gradient_thalamus\TD1\gradients';
95. %%
96. if ~exist(out_path,'dir')
97. mkdir(out_path);
98. end
99. data_all = dir([root_dir filesep '*.mat']);
100. for sub = 1:numel(data_all)
101. sub_func = load([root_dir filesep data_all(sub).name]);
102. %
103. %GrandientMaps
104. gm = GradientMaps('kernel','normalizedAngle','approach','diffusionEmbedding');
105. gm = gm.fit(sub_func.Cosine_w_F);
106. out_name = [out_path filesep data_all(sub).name];
107. mkdir(out_name);
108. save([out_name '/gredient.mat'],'gm');
110. %if~exist(out_name,'dir')
111. %mkdir(out_name);
112. %end
113. sub %command windows
114. end
116. %%step2
117. %labeling = load('D:\4BNTG_day4\scheafer_400_7network_order\label.mat');
118. root_dir = 'E:\tha_gradient\Gradient_thalamus\top5\local\gradients';
119. out_path = 'E:\tha_gradient\Gradient_thalamus\top5\local\gradients_allign';
120. %%
121. if ~exist(out_path,'dir')
122. mkdir(out_path);
123. end
124. data_all = dir([root_dir filesep '*.mat']);%sub-*
125. %
126. n_sub= numel(data_all);
127. for sub = 1:n_sub
128. sub_func = load([root_dir filesep data_all(sub).name '\gredient.mat']);
129. all_sub_gradients{sub} = sub_func.gm.gradients{1};
130. sub
131. end
132. %
133. [aligned, xfms] = procrustes_alignment(all_sub_gradients,'nIterations', 100);
134. for sub = 1:n_sub
135. aligned_gradient = aligned{sub};
136. aligned_gradient_all(:,:,sub)=aligned_gradient;
137. out_name = [out_path filesep data_all(sub).name];
138. if ~exist(out_name,'dir')
139. mkdir(out_name)
140. end
141. save([out_name,'/aligned_gradient.mat'],'aligned_gradient');
142. sub
143. end
144. %%
145. group_dir = 'E:\tha_gradient\Gradient_thalamus\top5\local\group_results';
146. if ~exist(group_dir,'dir')
147. mkdir(group_dir);
148. end
149. group1_gradient= mean(aligned_gradient_all(:,:,1:46),3); %HC
150. group2_gradient= mean(aligned_gradient_all(:,:,47:118),3); %PIGD
151. group3_gradient= mean(aligned_gradient_all(:,:,119:189),3);%TD
152. save([group_dir '\group1_mean_gradient.mat'],'group1_gradient');
153. save([group_dir '\group2_mean_gradient.mat'],'group2_gradient');
154. save([group_dir '\group3_mean_gradient.mat'],'group3_gradient');
155. save([group_dir '\all_subjects_aligned_gradient.mat'],'aligned_gradient_all');
156. %
157. figure;h1=histfit(group1_gradient(:,1),[],'kernel');hold on; h2=histfit(group2_gradient(:,1),[],'kernel');h3=histfit(group3_gradient(:,1),[],'kernel');
158. h1(1).FaceAlpha = 0.3;h2(1).FaceAlpha = 0.3;h3(1).FaceAlpha = 0.3;h1(2).Color = [.9 .0 .0];h2(2).Color = [0.2 0.2 0.5];h3(2).Color = [.4 .2 .3];
160. figure;h1=histfit(group1_gradient(:,2),[],'kernel');hold on; h2=histfit(group2_gradient(:,2),[],'kernel');
161. h1(1).FaceAlpha = 0.2;h2(1).FaceAlpha = 0.2;h1(2).Color = [.2 .2 .2];h2(2).Color = [.2 .2 .2];
163. figure;h1=histfit(group1_gradient(:,3),[],'kernel');hold on; h2=histfit(group2_gradient(:,3),[],'kernel');
164. h1(1).FaceAlpha = 0.2;h2(1).FaceAlpha = 0.2;h1(2).Color = [.2 .2 .2];h2(2).Color = [.2 .2 .2];
165. %% range and std
166. emb_range = zeros(n_sub,3);%emb_range
167. for i = 1:3
168. for j = 1:n_sub
169. emb_range(j,i) = max(aligned_gradient_all(:,i,j)) - min(aligned_gradient_all(:,i,j));
170. end
171. end
172. save([group_dir,'/emb_range.mat'],'emb_range');
174. % calculate variation
175. emb_std = zeros(n_sub,3);%emb_std
176. for i = 1:3
177. for j = 1:n_sub
178. emb_std(j,i) = std(aligned_gradient_all(:,i,j));
179. end
180. end
181. save([group_dir,'/emb_std.mat'],'emb_std');

Untitled.m, no license · at the source