-
Notifications
You must be signed in to change notification settings - Fork 2
/
Replicating_LB_2D_imco_started_20200911.Rmd
4646 lines (3322 loc) · 241 KB
/
Replicating_LB_2D_imco_started_20200911.Rmd
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
---
title: "R Notebook for IMCO"
output:
html_notebook
---
```{r graphic_setup, echo=TRUE}
library(knitr) # For knitting document and include_graphics function
library(ggplot2) # For plotting
library(png) # For grabbing the dimensions of png files
```
This project was started on 20200911
IMCO, or intermodal coupling (or "coupling"), was a project that begun in 2018 and has had numerous data anaylists/managers on the project. It is based on work by Simon 's work in 2016 (Paper: Vandekar, S. N., Shinohara, R. T., Raznahan, A., Hopson, R. D., Roalf, D. R., Ruparel, K., … Satterthwaite, T. D. (2016). Subject-level Measurement of Local Cortical Coupling. NeuroImage, 133, 88–97.), which looked at the relationship between cortical thickness and sulcal depth, which demonstrated regional differences, non-linear trajectories, as well as sex differences. His work used surface gray matter maps. The general process was the brain was divided into vertices and a neighborhood around each vertex was defined. For each neighborhood, the relationship between ct and sulcal depth was calculated with a variety of covariates.
In 2018, Lauren Beard replicated this work, studying cbf (ASL) with ALFF (rs-fmri). This was all done in 2D at the surface. Ali and Kristin then developed a 3D method (using volumes rather than surface), but ran into significant partial volume effects given white matter flanking edges of GM. As such, we have decided to move forward with Lauren Beard's main analyses in 2D space.
#Relevant Wikis:
1) how to generate 2D coupling maps: https://github.com/PennBBL/tutorials/wiki/Surface-Coupling
2) CBF-ALFF specifically (single subject): https://github.com/PennBBL/imcoScripts/wiki/Rest-CBF-Coupling
*Links to* https://github.com/PennBBL/tutorials/wiki/3D-Volume-to-Surface-Projection-(FS)
3) group analysis: https://github.com/PennBBL/imcoScripts/wiki/AlffCbf-RehoCbf-GLM-Summary
4) Vertex wise group analysis: https://github.com/PennBBL/tutorials/wiki/FreeSurfer-Vertex-Wise-Group-Level-Analyses
5) results:https://github.com/PennBBL/imcoScripts/blob/master/alffCbf_rehoCbf_effects_20180607%20(1).pdf
Replicated by Graham Baum
---
Lab notebook by date:
#### **20200911:**
#####1) Moving through Surface-Coupling Wiki:
*Source code*: **https://bitbucket.org/simonvandekar/coupling**
Input for this analysis consists of a csv that lists bblid/scanid: **/data/joy/BBL/tutorials/exampleData/surfCoupling/subjList.csv**
It also requires that the subjects have been processed using freesurfer. Specifically, **these files must be present:** lh.sphere.reg lh.sulc lh.thickness rh.sphere.reg rh.sulc rh.thickness.The **fsaverage5 directory must be present** as well.
*Code location*: **/data/joy/BBL/tutorials/code/surfCoupling**
--- *\*This is all on sample data, not on actual data --- skipped replicating this step and moved onto second Wiki*
---
#####2) CBF-ALFF Wiki:
Surface projections are carried out using **freesurfer version 6**. Surface coupling carried out with **freesurfer version 5.3**. Please refer to this wiki for more information on surface projections:
*Code Location*: **/data/jux/BBL/projects/coupling/imcoScripts/surfProjection/\***
**Steps: **
* 1. copy LB's stuff to PMACS from chead - *completed 20200911*
+ Remote connect via F5
+ ssh -Y eballer@chead
+ cd /data/jux/BBL/projects/
+ in a second terminal, ssh -Y [email protected]
+ cd /project/imco
+ mkdir from_lb #denoting data copied from lb, on chead
+ scp -r coupling [email protected]:/project/imco/from_lb
* 2. Making changes so that paths are PMACS and not chead
+ \*9/15: Ran into issues here - the scripts not only have chead paths, but they also require different freesurfer loads. These loads will need to be copied to PMACS in addition to having all scripts changed. This will require more effort to replicate than originally expected.
#### **20200915:**
- Decided to change course from original plan to move everything from chead to PMACS and rerun
- We will reprocess images with ASLPrep and FMRIPREP+XCP on cubic
- Figuring out whether we
+ Rerun but use existing container, project to surface later
+ Rerun and wait for new container that projects to surface
#### **20200916:**
- Will proceed in two steps
+ Will rerun in ASLPrep/fmriPrep + xcp ( a little later)
+ But first, will start with datafreeze data
- move data from datafreeze alff and cbf onto PMACS
- redo surface projections for a few people
- run coupling maps for a few people
- this will let me confirm that this stuff is good/her data is okay
#### **20200917:**
- Trying to find the input files for surface projection based on wikis
- issue is that code is in example data, but not clear where original files are
- need to start with freesurfer files
- to meet with Azeez 9/18 to discuss more about freesurfer
#### **20200918:**
- Seems to be something possibly useful in: /data/joy/BBL/studies/pnc/processedData/structural/freesurfer53/\*/fsaverage/
- WAIT, I think I found the files (these have mgh files in them)!!!!!:
/data/joy/BBL/studies/pnc/processedData/structural/freesurfer53/\*/\*/surf
- pcasl and restbold have asl and alff subcategories, I think this is where things came from
- made a directory:
```{bash}
#get into directory
cd /data/jux/BBL/projects/coupling/coupling_test_eb_20200918/
#make new subdirectories
mkdir code subjectLists surfaceProjection surfaceProjection/cbf surfaceProjection/alff
#go into code dir
cd code
#copy LB's surf projection code
cp /data/jux/BBL/projects/coupling/imcoScripts/surfProjection/* .
#go into subjectLists
cd /data/jux/BBL/projects/coupling/coupling_test_eb_20200918/subjectLists/
#copy LB's alff surface projection list
cp /data/jux/BBL/projects/coupling/subjectsLists/n2166_alff_surfProj.csv .
#extract first 5
more n2166_alff_surfProj.csv | head -5 >! test_n5_alff_surfProg.csv
```
- Testing whether I can replicate vol2surf
```{bash}
cd /data/jux/BBL/projects/coupling/coupling_test_eb_20200918
cd /data/jux/BBL/projects/coupling/coupling_test_eb_20200918/code
./vol2surf_wrapper.sh -i /data/jux/BBL/projects/coupling/coupling_test_eb_20200918/subjectLists/test_n5_alff_surfProg.csv -s /data/joy/BBL/studies/pnc/processedData/structural/freesurfer53/ -o /data/jux/BBL/projects/coupling/coupling_test_eb_20200918/surfaceProjection/alff
qstat #to see progress
```
**Check output**
```{bash}
cd /data/jux/BBL/projects/coupling/coupling_test_eb_20200918/surfaceProjection
#double check that 5 directories of BBLIDs are created
cd 100031/20100918x3818
more 100031_20100918x3818_projection.log
```
**Error**
- line 1: /opt/sge/default/spool/compute-0-20/job_scripts/921984: line 26: lta_convert: command not found
regio_read_register(): No such file or directory -> **make sure I have right freesurfer libraries**
- more error.log
+ mghRead(/data/jux/BBL/projects/coupling/coupling_test_eb_20200918/surfaceProjection//100031/20100918x3818/100031_20100918x3818_rh_surf.mgh, -1): could not open file **likely related to failing step 1***
**looked at ~/.bash_profile and realized freesurfer 5.3.0 was version being used**
- per Surface Coupling wiki, coupling requires you change the bash_profile to
"export FREESURFER_HOME=$CFNAPPS/freesurfer/5.3.0"
```{bash}
gedit ~/.bash_profile
# added export FREESURFER_HOME=$CFNAPPS/freesurfer/6.0.0 to profile and commented out 5.3.0 version and saved
source ~/.bash_profile
```
**then compared - and I can replicate surface projection!**
```{bash}
diff 100088_20120514x6853_matrixFinalTransform.dat /data/jux/BBL/projects/coupling/surfaceMaps/alff/100088/20120514x6853/100088_20120514x6853_matrixFinalTransform.dat
```
Output: Same... i.e. diff returns **nothing**
- tried this for a bunch of the files, for a bunch of the people, and it worked
- also compared files of different types to verify, and would get errors:
```{bash}
diff 100031_20100918x3818_lh_surf.mgh ../../../../surfaceMaps/alff/100031/20100918x3818/100031_20100918x3818_lh_fs5_surf.mgh
```
Result: \"Binary files 100031_20100918x3818_lh_surf.mgh and ../../../../surfaceMaps/alff/100031/20100918x3818/100031_20100918x3818_lh_fs5_surf.mgh differ\"
---
**From meeting with Azeez @ 1:30pm**
- First step of preprocessing is to run freesurfer on T1 (i.e. in T1 space)
+ scp to get cbf, reho, alff, all in native space (\*std means standard or pnc space)
+ seq2struct -> changes reho, cbf, alff into T1 space
- after transform, they are in T1 space
+ then, you have to do vol2surf
+ lta_convert changes to fsl format
+ R script transformMatrix.R is important because fsl matrix is in a different format than freesurfer, so you need to do some manipulation,multiply some column by -1. It is a whole to-do
+ mri_vol2surf actually goes from fsl to freesurfer
- reg -> takes matrixFinalTransform.dat, which is output of transformMatrix.R
- now you have actually moved from native space to freesurfer space
+ fsaverage5 (fs5) - freesurfer space that is most similar to 2x2x2mm **most often used**
+ fsaverage6 (fs6) - more similar to 1x1x1mm, takes a lot more processing speed so people often use fs5
- surface to surface - moves from one surface to another
---
**Visualization to view mgh images**
```{bash}
freeview -f /share/apps/freesurfer/6.0.0/subjects/fsaverage5/surf/lh.pial:overlay=100031_20100918x3818_lh_fs5_surf.mgh
#you can add lh.inflated if you'd like to see an inflated map by clicking the green + next to the brain
```
**Now make sure I can project cbf into fs5 space**
```{bash}
#LB does not have a subject list for asl surface projection, so I will adapt my most recent one
#move into directory
cd /data/jux/BBL/projects/coupling/coupling_test_eb_20200918/subjectLists
#copy alff to cbf list
cp test_n5_alff_surfProg.csv test_n5_cbf_surfProg.csv
#change directory name from restbold to pacsl, greedy
more test_n5_cbf_surfProg.csv | perl -pe 's/restbold\/restbold_201607151621/pcasl\/pcasl_201607291423/g' > temp1
#change alff to cbf, greedy
more temp1 | perl -pe 's/alff/asl/g'> temp2
#change _asl to _aslQuantSST1 which refers to native space asl
more temp2 | perl -pe 's/_asl/_aslQuantSST1/' > temp3
#Change temp3 name back to original file
mv temp3 test_n5_cbf_surfProg.csv
#delete temp files
rm temp1 temp2
#go back to code
cd /data/jux/BBL/projects/coupling/coupling_test_eb_20200918/code
./vol2surf_wrapper.sh -i /data/jux/BBL/projects/coupling/coupling_test_eb_20200918/subjectLists/test_n5_cbf_surfProg.csv -s /data/joy/BBL/studies/pnc/processedData/structural/freesurfer53/ -o /data/jux/BBL/projects/coupling/coupling_test_eb_20200918/surfaceProjection/cbf
```
**great news! Was able to make all of the appropriate files**
Checking to see if they match LB's cbf surface projections
```{bash}
#check if surface projections match
diff /data/jux/BBL/projects/coupling/surfaceMaps/cbf/100088/20120514x6853/100088_20120514x6853_matrixFinalTransform.dat /data/jux/BBL/projects/coupling/coupling_test_eb_20200918/surfaceProjection/cbf/100088/20120514x6853/100088_20120514x6853_matrixFinalTransform.dat #match
diff /data/jux/BBL/projects/coupling/surfaceMaps/cbf/100088/20120514x6853/100088_20120514x6853_lh_fs5_surf.mgh /data/jux/BBL/projects/coupling/coupling_test_eb_20200918/surfaceProjection/cbf/100088/20120514x6853/100088_20120514x6853_lh_fs5_surf.mgh #match
#trying different person to be sure
diff /data/jux/BBL/projects/coupling/surfaceMaps/cbf/100050/20130810x8309/100050_20130810x8309_rh_fs5_surf.mgh /data/jux/BBL/projects/coupling/coupling_test_eb_20200918/surfaceProjection/cbf/100050/20130810x8309/100050_20130810x8309_rh_fs5_surf.mgh #match
```
**We can now be confident that surface projections in both alff and cbf (from native space through to freesurfer space) are good **
#### 20200922
Will proceed with attempting coupling today on surface maps generated last week
Step 1: Make subject list to be used for coupling.
- For this, I edited the surface projection list
```{bash}
#go into directory
cd /data/jux/BBL/projects/coupling/coupling_test_eb_20200918/subjectLists
#take alff subject list (though could use cbf as well), grab BBLID and SCANID, and put them in new file format <BBLID>/<SCANID> by getting rid of , between, stores in n5_surfCoupling_list
more test_n5_alff_surfProg.csv | perl -pe 's/(.*),(\d{8}x\d{4}),.*/$1\/$2/' > n5_surfCoupling_list
```
**Comparing /data/joy/BBL/tutorials/code/surfCoupling/grid_submit.sh with /data/jux/BBL/projects/coupling/imcoScripts/surfCoupling/R/grid_submit.sh **
- Essentially identical. Big difference
+ /data/jux/BBL/projects/coupling/imcoScripts/surfCoupling/R/grid_submit.sh calls cbfAlffgrid.sh in qsub command, whereas /data/joy/BBL/tutorials/code/surfCoupling/grid_submit.sh calls surfCoup_wrapper.sh
- Minor differences
+ paths (obviously)
+ commented out line #sublist=$scripts/bblid_scanid.txt in /data/jux/BBL/projects/coupling/imcoScripts/surfCoupling/R/grid_submit.sh
+ logdir=$scripts/logdir2 in /data/jux/BBL/projects/coupling/imcoScripts/surfCoupling/R/grid_submit.sh (versus logdir1 in tutorials example)
**Comparing /data/jux/BBL/projects/coupling/imcoScripts/surfCoupling/R/cbfAlffgrid.sh with /data/joy/BBL/tutorials/code/surfCoupling/surfCoupling_wrapper.sh**
- Major Difference
- flags at the end of the call
- surfCoupling_wrapper: -m thickness,sulc -f 15 -t fsaverage5
- cbfAlffgrid.sh: -m alff,cbf
- Minor Difference
+ both call coupling_v2.R
+ both call different paths/subject lists (of course)
- Assessment:
+ my read on this is that alff,cbf should be our calls (this indicates what coupling measures we actually want to use)
+ flags for fwhm at 15 and fsaverage5 are correct, and seem to be defaults for cbfAlffgrid.sh
**Will use surfCoupling_wrapper.sh and ammend thickness,sulc to be alff,cbf **
---
**Comparing /data/joy/BBL/tutorials/code/surfCoupling/coupling_v2.R and /data/jux/BBL/projects/coupling/imcoScripts/surfCoupling/R/coupling_v2.R **
1) In first section
- /data/joy/BBL/tutorials/code/surfCoupling/coupling_v2.R has a lot of commented out stuff that is not there in the cleaner example code script
2) Next section is identical
3) In section of \#### GET ENVIRONMENT VARIABLES ####:
- Only difference is the signal to change scripts path
4) In section of \#### NOW FIND THE NEAREST NEIGHBORS FOR THIS TEMPLATE ####
- chunk after this is only different in the paths and instructions for what paths to replace
5) Both are the same below the line #### CREATE LABELS IF THEY DON'T EXIST YET ####
- line 90 in /data/joy/BBL/tutorials/code/surfCoupling/coupling_v2.R
- line 78 /data/joy/BBL/tutorials/code/surfCoupling/coupling_v2.R
**Will proceed with updating/changing coupling_v2.R and subsequent code copied from tutorial**
---
Step 2: Copy appropriate code into my code directory
```{bash}
#go into my directory
/data/jux/BBL/projects/coupling/coupling_test_eb_20200918/code
#copy tutorial code
cp /data/joy/BBL/tutorials/code/surfCoupling/* .
```
Step 3: Change directories within grid_submit.sh
- SUBJECTS_DIR=/data/joy/BBL/studies/pnc/processedData/structural/freesurfer53/
- scripts=/data/jux/BBL/projects/coupling/coupling_test_eb_20200918/code
- sublist=/data/jux/BBL/projects/coupling/coupling_test_eb_20200918/subjectLists/n5_surfCoupling_list
- add line: outdir=/data/jux/BBL/projects/coupling/coupling_test_eb_20200918/surfCoupling
- add output directories (added **\$outdir** to the end of both lines):
+ qsub -V -N coupling_\$time -t 1 -cwd -o \$logdir -e \$logdir -q \$queue \$scripts/surfCoup_wrapper.sh -s \$sublist -f \$fwhm -o **\$outdir**
+ qsub -V -t 2-\$nsub -cwd -hold_jid coupling_\$time -o \$logdir -e \$logdir -q \$queue \$scripts/surfCoup_wrapper.sh -s \$sublist -f \$fwhm -o **\$outdir**
- kept fwhm at 15, and wrapper as surfCoup_wrapper
- saved
- *info on what to replace is in the script
- *have kept grid_submit_eb.sh as reference
---
Step 4: Change directories within surfCoup_wrapper.sh
- change /data/joy/BBL/tutorials/code/surfCoupling/coupling_v2.R to /data/jux/BBL/projects/coupling/coupling_test_eb_20200918/code/coupling_v2.R
- change /data/joy/BBL/tutorials/exampleData/surfCoupling/subjList.csv to /data/jux/BBL/projects/coupling/coupling_test_eb_20200918/subjectLists/n5_surfCoupling_list
- change thickness,sulc to alff,cbf
- verified that the changes actually point to real directories that are correctly formatted
- saved
---
Step 5: Change directories within /data/jux/BBL/projects/coupling/coupling_test_eb_20200918/code/coupling_v2.R
- line 41, changed script directory to: /data/jux/BBL/projects/coupling/coupling_test_eb_20200918/code
- line 75, changed kth neighbor path to: /data/jux/BBL/projects/coupling/coupling_test_eb_20200918/code/kth_neighbors_v3.R
- verified these scripts/locations exist
- saved
---
Step 6: Verified that no changes need to be made to kth_neighbors_v3.R
- no changes needed to be made
---
Step 7: Change environment settings per LB (surfCoupling will only work on FS 5.3)
- open bash profile (gedit ~/.bash_profile)
- verify FS version
+ line should read: FREESURFER_HOME=$CFNAPPS/freesurfer/5.3.0
- if commented out, uncomment
- if not present, add it
- comment out line with 6.0.0
- source ~/.bash_profile
- set SUBJECTS_DIR (this is confusing because it looks like it is set in grid_submit.sh), but did it anyway
```{bash}
#export subjects dir, have gone back and forth with this
export SUBJECTS_DIR=/data/joy/BBL/studies/pnc/processedData/structural/freesurfer53/
#export SUBJECTS_DIR=/data/jux/BBL/projects/coupling/coupling_test_eb_20200918_surfCoupling
```
- after completed, **WILL NEED TO** comment out the 5.3.0 line and bring back the 6.0.0 line
---
Step 8: Ran: qsub -q himem.q -V -N baller_imco_replication ./grid_submit.sh
**ERROR**
In /data/jux/BBL/projects/coupling/coupling_test_eb_20200918/code/logdir :
coupling_161336.o925032.1, surfCoup_wrapper.sh.o925033.[1-5], all the same:
Loading required package: optparse
Warning message:
In library(package, lib.loc = lib.loc, character.only = TRUE, logical.return = TRUE, :
there is no package called ‘optparse’
Error in make_option(c("-s", "--subjectlist"), action = "store", default = NULL, :
could not find function "make_option"
Execution halted
---
#### 20200923
- seems like an issue with python optparse - is this accurate? How to fix?
+ looks like /data/jux/BBL/projects/coupling/imcoScripts/surfCoupling/setup_test.R does this
```{bash}
# copied it into my directory
cp /data/jux/BBL/projects/coupling/imcoScripts/surfCoupling/setup_test.R /data/jux/BBL/projects/coupling/coupling_test_eb_20200918/code/.
```
- evaluated, looks like it sets up the scripts, main pieces are to add a new SHEBANG and install package.
- Added these directly at the top of my grid_submit.sh
+ \#!/share/apps/R/R-3.2.3/bin/Rscript --vanilla --no-init-file
+ install.packages('optparse', repo='http://cran.us.r-project.org')
- Planned to run grid_submit.sh again, but not R script
- decided to add it to top of coupling_v2.R
+ \#!/share/apps/R/R-3.2.3/bin/Rscript --vanilla --no-init-file
+ install.packages('optparse', repo='http://cran.us.r-project.org')
---
Found a new script that looks more promising, and uses the directories that were previously made
- /data/jux/BBL/projects/coupling/imcoScripts/surfCoupling/old/alffCbf_surfCoupling.R
**Will evaluate and amend this one**
---
Step 1: Copied this into my directory
```{bash}
#copy into my directory
cp /data/jux/BBL/projects/coupling/imcoScripts/surfCoupling/old/alffCbf_surfCoupling.R /data/jux/BBL/projects/coupling/coupling_test_eb_20200918/code
```
Step 2: Changed directory paths **have not done this yet**
- alffCbf_subjects <- read.csv("/data/jux/BBL/projects/coupling/coupling_test_eb_20200918/subjectLists/n5_surfCoupling_list")
- demos <- read.csv("/data/joy/BBL/studies/pnc/n1601_dataFreeze/demographics/n1601_demographics_go1_20161212.csv")
- setwd("/data/jux/BBL/projects/coupling/coupling_test_eb_20200918/couplingSurfaceMaps/alffCbf/lh/stat")
**This requires the surface coupling was already done, this just looks at covariates **
**Don't know why this is considered \"old\" **
---
Per meeting with Azeez:
**Big problems**
- Chaining of multiple scripts whose inputs don't feed directly into each other
- R version used by script is outdated (may be a minor problem)
- Input and output directory for LB's stuff writes to the same directory
- Seems to be some manual moving and renaming of files, which makes it very hard to replicate
- Tutorial files/info was also moved manually into a strange structure that is not present in actual data.
**Minor issues**
- I am guessing the actual algorithm works, but figuring out input/outputs very difficult
- Lots of multi levels of naming stuff weirdly in coupling_v2.R
+ plan to avoid environmental variables and just name things locally- making it all super confusiong
**Actionable ideas**
1) don't do grid_submit, could use surfCoupling_wrapper.sh or just start with coupling_v2.R
2) She did actually take files from /data/joy/BBL/studies/pnc/processedData/structural/freesurfer53/ and dumbed it back into the same directory. Looks like she then manually copied these files
3) There are many iterations of files in /data/joy/BBL/studies/pnc/processedData/structural/freesurfer53/, some with doubles
4) Good news is that
- diff ../surfaceProjection/alff/100031/20100918x3818/100031_20100918x3818_lh_fs5_surf.mgh /data/joy/BBL/studies/pnc/processedData/structural/freesurfer53/100031/20100918x3818/surf/lh.alff.fsaverage5.mgh
- these are the same files
5) Will need to take the coupling_v2.R script (added SHEBANG and optparse load) and rewrite so it grabs data from my directories and writes to mine
- INPUT: /data/jux/BBL/projects/coupling/coupling_test_eb_20200918/surfaceProjection/[alff,cbf]/*
- OUTPUT: /data/jux/BBL/projects/coupling/coupling_test_eb_20200918/couplingSurfaceMaps/alffCbf/lh/
6) Once this is done/verified, will need to copy the actual files to PMACS (or preferably CUBIC) and do correlations with age, sex, cognition there
7) Will also need to figure out which files that are in the datafreeze actually match what I generate
8) containerize?
**surface files just has numbers, asc is same info, just different format**
Afternoon:
1) Started rewriting coupling script. Will be called **coupling_v2_eb.R**
- For now, I am going to get rid of the wrappers and just assign things directly
#### 20200924
- Continuing to rewrite the code in coupling_v2_eb.R
+ kth_neighbor_v3.R makes the neighborhoods based on whatever size you want to use
- will need to give this a better output file
- Found whole new directory with seemingly important stuff:
+ /data/jux/BBL/projects/coupling/imcoScripts/restCbf/
+ unclear how this relates, but within dependencies, there is an imco script
+ looks like this is done at volume and not surface level, but should theoretically be a similar approach
+ **will pass for now, coupling seems to be contained in coupling_v2.R script**
- also starting to work on kth_neighbor_v3_eb.R
+ Am going to send it an additional argument so that I can specify output directory
- sending neigh_outdir (/data/jux/BBL/projects/coupling/coupling_test_eb_20200918/surfCoupling/neighborhood_info) to the kth_neighbor_v3_eb.R script
#### 20200925
- Start at line 15 for kth_neighbor_v3_eb.R
+ make sure you are coping the asc files (that is the command entered yesterday)
To do:
- [x] go through the script and make sure indirectory/out directory are ok
- [x] suggest copying files from structural into my local directory and doing everything from within that
- [x] compare with neighborhoods from structural dir
+ x <-readRDS("/data/joy/BBL/studies/pnc/processedData/structural/freesurfer53/fsaverage5/surf/rh.10neighbors_degree.rds")
+ y <- readRDS("/data/jux/BBL/projects/coupling/coupling_test_eb_20200918/surfCoupling/neighborhood_info/surf/rh.10neighbors_degree.rds")
+ plot(x[,5], y[,5]) -> identical #randomly selected
+ cor.test(x[,5], y[,5]) : **t = inf, df = 10248, p-val <2.2 x 10^-16, cor = 1**
+ identical(x,y) = TRUE
- **kth_neighbor_v3_eb.R has been amended and replicated, now writing neighborhood tables to my local directory**
- [x] we are using structural kth neighbor - seems like those vertices should be okay as it is not actually running coupling but getting neighborhoods, which should be derived structurally
+ [x] double check with azeez
- azeez says that this is some sort of partial volume correction but I am not sure
#### 20200929
To do
- [x] read simon's paper again to understand the neighborhood business, as well as the imco manuscript. My sense is that the defining of the neighborhoods on the structural inflated brain is okay, because alff and cbf surface projections were based on those structural brains. - still think this is true
- [o] surface projections were in native space, so that group map might be useless
#### 20200930
- [x] will keep working through script. Important part is that it replicates. We are just trying to make sure the output files can be used
+ having a lot of issues getting the MRI preprocessing command to work
+ it is set up to use structural data and structural paths
+ please see beginning of script where I have to set environmental variable to run with freesurfer
```{bash}
#must set subject directory or freesurfer will not know where to look
#currently set to /data/joy/BBL/studies/pnc/processedData/structural/freesurfer53
Sys.setenv
```
- [x] what is .asc (looks like suffix to indicate ascii) vs mgh? Do i need to take my single subject surface projections (in mgh) and convert to asc and then do kth neighbor on those files?
+ https://cran.r-project.org/web/packages/freesurferformats/freesurferformats.pdf (page 40)
+ http://brainvis.wustl.edu/pipermail/caret-users/2014-June/006150.html
+ **it is ascii abbreviation and no, don't need to take indiv surface for template b/c the surface projection was based on structural template, so it should similarly define neighborhoods**
-**Spoke with Azeez, there are so many dependencies with respect to using freesurfer stuff in joy that this is essentially super challenging to replicate bottom up.**
-**Will instead just run wrapper /data/jux/BBL/projects/coupling/imcoScripts/surfCoupling/RcbfAlffgrid.sh and change outasc in coupling_v2.R to see if I can make the files and if I can, just go on from there**
-Went back to original coupling_v2.R in my personal coupling directory, and saved a copy
coupling_v2_eb_justoutasc_change.R
+ In this, just changing where files are written to, to see if I can replicate. Not going to rewrite the whole thing. Feels like a waste.
```{r}
# Added lines 140 to 147
localdir = "/data/jux/BBL/projects/coupling/coupling_test_eb_20200918/surfCoupling/couplingMaps/surf"
subdir = file.path(localdir, id)
#added by EB, if direc doesn't exist, make it
if (!dir.exists(subdir)) {
system(paste0("mkdir ", subdir))
}
```
- Then realized that this was not the problem, the script writes everything to a single directory. So I deleted that piece and instead added two lines (163 and 169)
```{r}
#dded by EB
localdir = "/data/jux/BBL/projects/coupling/coupling_test_eb_20200918/surfCoupling/couplingMaps/"
localsubj_dir = file.path(localdir, id, surf)
if (!dir.exists(localsubj_dir)) {
system(paste0("mkdir ", file.path(localsubj_dir)))
}
outasc = file.path(localsubj_dir, paste(hemi, paste('coupling_', retcoef, '_', paste(measures, collapse='_'), '.fwhm', fwhm, sep=''), template, 'asc', sep='.') )
```
- **spoke with Azeez, will take a break so he can work on some stuff and we can look together next week**
#### 20201007
- spoke with Azeez. he stated:
"I look at the code, to make it work for coupling computation, the projected surfaces ike rh.alff.fsaverage5.mgh etc has to be copied to surf/rh.alff . A lot of steps were skipped that are not included in the wiki for all and cbf coupling. the original code (gitlab) by Simeons will work on structural only ( curv,thickness)"
He plans to move this all to python soon.
##### To do:
[x] can replicate sample construction
- using /data/jux/BBL/projects/coupling/imcoScripts/restCbf/n831_alff_cbf_makeSample.R
- **question - why is ltnExclude used rather than healthRatingExclude? **
+ *answer: ltn is low threshold normal (or dirty normal). healthExclude is major medical issues, squeaky clean is TDs. fewest # meet criteria for squeaky clean, then LTN, then health Exclude
- **question - there is no original file output from this, but it does seem to make logical sense and gets to 831 **
- ltnExclude takes out 340 people, healthRating exclude only 154
[x] figure out which high level script she used and replicate backwards (i.e. from highest level 2nd level analysis to lower)
+ starting in /data/jux/BBL/projects/coupling/imcoScripts/restCbf
- /data/jux/BBL/projects/coupling/imcoScripts/restCbf/alffCbf_surfCoupling.Rmd (this is most up-to-date script)
+good news is it runs
+not quite sure what the actual models show
+looks like she loops through (each voxel? 10242) and runs a linear model, relatng the coupling to age, sex and motion.
+ every time she does the loop though, it should rewrite the model.
+ it is showing huge age effect, not sex effect
+ the big missing chunk right now is actually making the coupling maps, but may be able to replicate group analysis and sample construction
- /data/jux/BBL/projects/coupling/imcoScripts/surfCoupling/* scripts will be the last to replicate
#### 20201008
- Stepping through /data/jux/BBL/projects/coupling/imcoScripts/restCbf/alffCbf_surfCoupling.Rmd to see if I can make sense of what is happening. Code below
```{r}
---
title: "Alff-Cbf Surface Coupling"
author: "Erica Baller"
date: '20201007'
output:
word_document: default
pdf_document: default
html_document: default
---
### ALFF-CBF SURF COUPLING
# read in packages
library(plyr)
# read in demographics, this is directly from LB's stuff
alffCbf_subjDemos <- read.csv("/data/jux/BBL/projects/coupling/subjectsLists/n831_rest_cbf_finalSample_imageOrder.csv")
#factorize sex
alffCbf_subjDemos$sex <- as.factor(alffCbf_subjDemos$sex)
# read in all files and get in the right format
#change working directory within R - this only works if you knit, not if you run line by line
setwd("/data/jux/BBL/projects/coupling/couplingSurfaceMaps/alffCbf/lh/stat")
#pth <- paste0("/data/jux/BBL/projects/coupling/couplingSurfaceMaps/alffCbf/lh/stat/")
#store all files with .asc in lh_alffCbf_files - this is the part that only works if you knit - otherwise, you stay in your current working directory and these local files without full paths don't exist. Tried to add full paths but could not do this. If I wanted to try again, could do a string replace but thought better of it.
lh_alffCbf_files = list.files(pattern="*.asc")
#goes through each file, reads in the table, and then makes each table a row in the lh_alff_Cbf_data data frame. This makes one massive data frame.
#do.call actually makes every element of the list into a row and stacks them together. So each person should be represented by 10242 rows - this is # of voxels?
lh_alffCbf_data = do.call(rbind, lapply(lh_alffCbf_files, function(x) read.table(x, stringsAsFactors = FALSE)))
#extracts the 5th column, which looks like T values
lh_alffCbf_coupling <- as.data.frame(lh_alffCbf_data$V5)
#takes the data frame and transposes it. Specifically, it will divide the # of rows total by 831 and group (this basically re-separates the lines by person), and then transposes it
#the end result is a data frame with 831 rows, with one column per voxel, and the value of each column is the T value associated with the T-test at that voxel
lh_alffCbf_coupling_n <- t(as.data.frame(split(lh_alffCbf_coupling,1:831)))
#do same thing for right side
setwd("/data/jux/BBL/projects/coupling/couplingSurfaceMaps/alffCbf/rh/stat")
rh_alffCbf_files = list.files(pattern="*.asc")
rh_alffCbf_data = do.call(rbind, lapply(rh_alffCbf_files, function(x) read.table(x, stringsAsFactors = FALSE)))
rh_alffCbf_coupling <- as.data.frame(rh_alffCbf_data$V5)
rh_alffCbf_coupling_n <- t(as.data.frame(split(rh_alffCbf_coupling,1:831)))
# run model. For each voxel, run linear model regressing age and sex on T. This confuses me, because by the ned, each of the models only contain the last i (10242), and not all of the interim models
for (i in 1:10242) {
lh_alffCbf_agemodel <- lm(lh_alffCbf_coupling_n[,i] ~ ageAtScan1, data=alffCbf_subjDemos)
lh_alffCbf_ageSexmodel <- lm(lh_alffCbf_coupling_n[,i] ~ ageAtScan1 + sex, data=alffCbf_subjDemos)
lh_alffCbf_ageSex_qaModel <- lm(lh_alffCbf_coupling_n[,i] ~ ageAtScan1 + sex + pcaslRelMeanRMSMotion + restRelMeanRMSMotion, data=alffCbf_subjDemos)
}
for (i in 1:10242) {
rh_alffCbf_agemodel <- lm(rh_alffCbf_coupling_n[,i] ~ ageAtScan1, data=alffCbf_subjDemos)
rh_alffCbf_ageSexmodel <- lm(rh_alffCbf_coupling_n[,i] ~ ageAtScan1 + sex, data=alffCbf_subjDemos)
rh_alffCbf_ageSex_qaModel <- lm(rh_alffCbf_coupling_n[,i] ~ ageAtScan1 + sex, pcaslRelMeanRMSMotion + restRelMeanRMSMotion, data=alffCbf_subjDemos)
}
# print out summary results
#summary(lh_alffCbf_agemodel)
#summary(lh_alffCbf_ageSexmodel)
summary(lh_alffCbf_ageSex_qaModel)
#summary(rh_alffCbf_agemodel)
#summary(rh_alffCbf_ageSexmodel)
summary(rh_alffCbf_ageSex_qaModel)
```
- started testing models with psychopathology and cognition
- mood and overall psychopathology - no diff
- accuracy, speed, efficiency - no diff
- age continues to be. by far, the most significant factor
- Azeez will figure out if model that LB ran makes sense, or just rewrites the same piece over and over again
#### 20201027
Doing LRP work now. Need to figure out number of subjects with good ASL and rest data.
- Final N = 969 (after excluding pcaslExclude, restExclude, t1Exclude, smry_dep = NA, healthExcludev2)
```{r}
aslQA <- read.csv("/Users/eballer/BBL/from_chead/ballerDepHeterogen/data/neuroimaging/neuroimaging/n1601_PcaslQaData_20170403.csv", sep = ",", header = TRUE)
aslQA <- aslQA[which(aslQA$pcaslExclude == 0),] #n = 1508
restQA <-read.csv("/Users/eballer/BBL/from_chead/ballerDepHeterogen/data/neuroimaging/rest/n1601_RestQAData_20170714.csv", sep = ",", header = TRUE)
restQA <- restQA[which(restQA$restExclude == 0),] #n = 1123
structQA <- read.csv("/Users/eballer/BBL/from_chead/ballerDepHeterogen/data/neuroimaging/t1struct/n1601_t1QaData_20170306.csv", sep = ",", header = TRUE)
structQA <- restQA[which(structQA$t1Exclude == 0),] #n = 1540
psych <- read.csv("/Users/eballer/BBL/from_chead/ballerDepHeterogen/data/n9498_goassess_psych_summary_vars_20131014.csv", sep = ",", header = TRUE)
psych <- psych[!is.na(psych$smry_dep),] #n = 9411
med <- read.csv("/Users/eballer/BBL/from_chead/ballerDepHeterogen/data/n1601_health_with_meds_20170421.csv", sep = ",", header = TRUE)
med <- med[which(med$healthExcludev2 == 0), ] #n = 1447
aslQA_subset <- subset(aslQA, by=c("bblid","scanid","pcaslExclude"))
restQA_subset <- subset(restQA, by=c("bblid","scanid","restExclude"))
structQA_subset <- subset(structQA, by=c("bblid","scanid","t1Exclude"))
merged <- merge(aslQA_subset, restQA_subset, by = c("bblid","scanid"))
merged <- merge(merged, structQA_subset, by = c("bblid","scanid"))
merged <- merge(merged, psych, by = "bblid")
merged <- merge(merged, med, by = "bblid") #final n = 969
print(paste0("Summary #s from smry_dep: Dep = 4: ",
length(which(merged$smry_dep == 4)), #n = 146
", Dep = 3: ", length(which(merged$smry_dep == 3)), #n = 7
", Dep = 2: ", length(which(merged$smry_dep == 2)), #n = 49
", Dep = 1: ", length(which(merged$smry_dep == 1)),# n = 277
", Dep = 0: ", length(which(merged$smry_dep == 0)))) #n = 490
```
#### 20201029
Will use anxious-misery factor score rather than depression to define groups. See below:
```{r}
aslQA <- read.csv("/Users/eballer/BBL/from_chead/ballerDepHeterogen/data/neuroimaging/neuroimaging/n1601_PcaslQaData_20170403.csv", sep = ",", header = TRUE)
aslQA <- aslQA[which(aslQA$pcaslExclude == 0),] #n = 1508
restQA <-read.csv("/Users/eballer/BBL/from_chead/ballerDepHeterogen/data/neuroimaging/rest/n1601_RestQAData_20170714.csv", sep = ",", header = TRUE)
restQA <- restQA[which(restQA$restExclude == 0),] #n = 1123
structQA <- read.csv("/Users/eballer/BBL/from_chead/ballerDepHeterogen/data/neuroimaging/t1struct/n1601_t1QaData_20170306.csv", sep = ",", header = TRUE)
structQA <- restQA[which(structQA$t1Exclude == 0),] #n = 1540
psych <- read.csv("/Users/eballer/BBL/from_chead/ballerDepHeterogen/data/n9498_goassess_psych_summary_vars_20131014.csv", sep = ",", header = TRUE)
#psych <- psych[!is.na(psych$smry_dep),] #n = 9411
#factor scores
goassess <- read.csv("/Users/eballer/BBL/from_chead/ballerDepHeterogen/data/n9498_goassess_itemwise_bifactor_scores_20161219.csv", sep = ",", header = TRUE)
med <- read.csv("/Users/eballer/BBL/from_chead/ballerDepHeterogen/data/n1601_health_with_meds_20170421.csv", sep = ",", header = TRUE)
med <- med[which(med$healthExcludev2 == 0), ] #n = 1447
aslQA_subset <- subset(aslQA, by=c("bblid","scanid","pcaslExclude"))
restQA_subset <- subset(restQA, by=c("bblid","scanid","restExclude"))
structQA_subset <- subset(structQA, by=c("bblid","scanid","t1Exclude"))
merged <- merge(aslQA_subset, restQA_subset, by = c("bblid","scanid"))
merged <- merge(merged, structQA_subset, by = c("bblid","scanid"))
merged <- merge(merged, psych, by = "bblid")
merged <- merge(merged, med, by = "bblid") #final n = 969
merged <- merge(merged, goassess, by = "bblid") #final n = 969
print(paste0("Summary #s from anxious-misery: Patients (moodv2 > 1): ",
length(which(merged$mood_4factorv2 > 1)), #n = 157
", TD: ", length(which(merged$mood_4factorv2 <= 1)))) #TD 812
number_of_depressive_pts <- length(which((merged$smry_dep == 4) | (merged$smry_gad == 4) | (merged$smry_ptd == 4) | merged$smry_man == 4))
number_of_tds <- dim(merged)[1] - number_of_depressive_pts
print(paste0("Summary #s for meeting full criteria (score =4: MDD + BAD + GAD + PTSD = ", number_of_depressive_pts, #n=231
", TD = ", number_of_tds, #n = 738
", mood_cat: ", length(which(merged$smry_mood_cat == 4)))) #n=154
#liberalize to people with 3/4
number_of_depressive_pts <- length(which((merged$smry_dep >= 3) | (merged$smry_gad >= 3) | (merged$smry_ptd >= 3) | merged$smry_man >= 3))
number_of_tds <- dim(merged)[1] - number_of_depressive_pts
print(paste0("Summary #s for liberalizing to 3/4: MDD + BAD + GAD + PTSD = ", number_of_depressive_pts, #n=259
", TD = ", number_of_tds, #n = 710
", mood_cat: ", length(which(merged$smry_mood_cat >= 3)))) #n=162
```
Based on this review, using moodv2 >2 as a cutoff, n = 157 of affected.
1) If you include people who have mdd, bad, gad, ptsd, you get n = 231 affected, 738 unaffected. Mood disorders specifically, n = 154
2) If you liberalize and look at people who have 3 or 4 for those measures, things do not change much
mdd+bad+gad+ptsd = 259, td = 710, mood_cat = 162
3) As such, **will do mdd+bad+gad+ptsd == 4**
#### 20201117
Rechanged criteria
- will do dimensional analysis across all mood states and HYDRA, use the criteria from my depression paper
Making figures from:
/data/jux/BBL/projects/coupling/fsGLM/lhAlffCbf_oneSampletTest/lh.n831.coupling_coef_alff_cbf.fwhm15.fsaverage5.mean.fwhm15.fsaverage5/contrast1/z.mgh
Contrast 1: Age
Contrast 2: Sex
Contrast 3: asl motion
Contrast 4: rest motion
To view:
```{bash}
freeview -f /share/apps/freesurfer/6.0.0/subjects/fsaverage5/surf/lh.pial:overlay=/data/jux/BBL/projects/coupling/fsGLM/lhAlffCbf/lh.n831.coupling_coef_alff_cbf.fwhm15.fsaverage5.ageAtScanmean_sex_pcaslRelMeanRMSMotion_restRelMeanRMSMotion.fwhm15.fsaverage5/contrast1/z.mgh
```
To change threshold, can add:
```{bash}
example visualization command: freeview -f /share/apps/freesurfer/6.0.0/subjects/fsaverage5/surf/lh.pial:overlay=/data/jux/BBL/projects/coupling/fsGLM/lhAlffCbf/lh.n831.coupling_coef_alff_cbf.fwhm15.fsaverage5.ageAtScanmean_sex_pcaslRelMeanRMSMotion_restRelMeanRMSMotion.fwhm15.fsaverage5/contrast1/z.mgh:overlay_threshold=3.09,30 #or whatever the range you would like
```
I am thresholding at 3.09 for my grant.
#### 20201118
Spoke with Azeez. Big hole in data is that the directory/naming format that the surface projection maps are saved under is not what the coupling script requires. The surface projection maps were actually manually renamed, which is why things have been so complicated to replicate.
A few notes from our meetingL
Coupling flow -> native space -> T1 -> surface (using fmri vox_to_surf), using fsaverage5
To make the coupling run,
1) Have to move surfaces into the freesurfer directory (which requires writing to frozen)
2) In coupling_v2
- mri_preproc is doing regression from one surface to another
- line 183 - mris_preproc
- saves file as mgh
- then converts to ascii which only saves vertex number, vertex location (x,y,z) and beta weigtht
- line 205 - if any regions w/NA, they will be replaced by 0s. This would account for areas that were masked out
- lr: local regression
- line 11, uncomment
- line 18, must change to cbf, alff. Files must be saved in different format
- this is the part of the script I'll need to write
- format needs to be subj/surf/rh.alff.fsaverage5.mgh and needs to be copied to freesurfer directory
----
For the grant,
Going through LB's scripts
- the nice lhAlffCBF script has a fatal mistake - does not store LMs
-take this script, get average of betas and correlate with age
- removed vertices with NA
- did not work, spread is all within 3-4 for betas. Got an essentially flat line.
- fsGLM_lhAlffCbf_eb.sh - useful?
- azeez will look for where the code is that generated the glm - that seems useful
- this directory that she points to in the script doesn't exist-
homedir=/data/jux/BBL/projects/coupling/imcoScripts/surfCoupling/groupAnalysis/lhAlffCbf
-----
#### 20201124
Working now on fsGLM, have identified script to work on:
I have changed the
- home directory: /data/jux/BBL/projects/coupling/coupling_test_eb_20200918/fsGLM
- data directory:/data/jux/BBL/projects/coupling/coupling_test_eb_20200918/subjectLists/n831_cbf_alff_w_mood_and_cognition.csv
- scripts directory: /data/jux/BBL/projects/coupling/coupling_test_eb_20200918/code/fsGLM/groupAnalysis
I have copied the design matrix locally
```{bash}
cp /data/jux/BBL/projects/coupling/imcoScripts/fsGLM/designMatrix.R /data/jux/BBL/projects/coupling/coupling_test_eb_20200918/code/fsGLM/groupAnalysis/.
```
I have constructed the new data frame using
/data/jux/BBL/projects/coupling/imcoScripts/fsGLM/make_new_subj_list.R
```{r}
#read in LB's n831 script and concatenate cognition and mood data for GLM
#read in LB's df
df <- read.csv("/data/jux/BBL/projects/coupling/subjectsLists/n831_rest_cbf_finalSample_imageOrder.csv", sep = ",", header = TRUE)
#read in clinical and cog dfs
clinical <- read.csv("/data/joy/BBL/studies/pnc/n1601_dataFreeze/clinical/n1601_goassess_itemwise_corrtraits_scores_20161219.csv", sep = ",", header = TRUE)
cog <- read.csv("/data/joy/BBL/studies/pnc/n1601_dataFreeze/cnb/n1601_cnb_factorscores.csv", sep = ",", header = TRUE)
#subset
#combine dfs
final_df <- merge(df, clinical, by = c("bblid", "scanid"))
final_df <- merge(final_df, cog, by = c("bblid", "scanid"))
write.csv(x = final_df, file = "/data/jux/BBL/projects/coupling/coupling_test_eb_20200918/subjectLists/n831_cbf_alff_w_mood_and_cognition.csv", quote = FALSE, row.names=FALSE)
```
Experimented with the glm models by adding mood_corrtraitsv2, mood_corrtraitsv2 + psychosis_corrtraitsv2 to see if each contrast was independent or if the F depended on how many terms in the model
F DEPENDS ON HOW MANY TERMS IN THE MODEL
- went back to just mood, right and left, this can be okay for sobp I believe
Output can be found in:
/data/jux/BBL/projects/coupling/coupling_test_eb_20200918/fsGLM/lh.n831.coupling_coef_alff_cbf.fwhm15.fsaverage5.ageAtScanmean_sex_pcaslRelMeanRMSMotion_restRelMeanRMSMotion_mood_corrtraitsv2.fwhm15.fsaverage5
/data/jux/BBL/projects/coupling/coupling_test_eb_20200918/fsGLM/rh.n831.coupling_coef_alff_cbf.fwhm15.fsaverage5.ageAtScanmean_sex_pcaslRelMeanRMSMotion_restRelMeanRMSMotion_mood_corrtraitsv2.fwhm15.fsaverage5
We now can replicate
-age effects (stronger coupling posteriorly) - contrast 1
- sex effects (stronger coupling in females > males) - contrast 2
And new mood effects (anterior, posterio cingulate, temporal) - contrast 3 but much weaker
Probs, nothing corrected
#### 20201201
- make pngs for ted
- abstract for sobp - age effects, sex effects, mood effects
- what is the F cutoff for "sig"
#### 20201202
- something is up
- I run fsglm with models including mood and psychosis together, i get different results than if I run a model with either mood or psychosis
- I also seem to be getting something funny for the right side of the brain
- right and left do NOT matcj
- Design files are also not hemisphere specific
- also confusing because l&r look ok in age and sex, just not mood and psychosis.
#### 20201203
plan:
[x] review script to run R and L sides of the brain
- mood one identical to each other with exception of lh/rh
[x] rerun mood independent analysis?
[x] talk to ted, am I using the right corr traits folder? -yes, correct to use non-age regressed because they are regressed on full sample, not my sample
[x] do I want bifactor? - YES
[x] which abstract to use? I think just development and sex differences one, can add other stuff once I get it working.- development and sex differences
[x] talk to adam on how to analyze surfaces
- he is going to get me some surface code for visualization.
[] tinashe working on general analysis tool to do models. This would be a good testing dataset
#### 20201204
[] Go into R data, store the results of each regression in a vector. FDR correct. Then he will help with visualization
- first, I made anew directory that will hold group level analyses in R, then copied the following scripts into it:
```{bash}
mkdir /data/jux/BBL/projects/coupling/coupling_test_eb_20200918/code/restCbf_R
cp /data/jux/BBL/projects/coupling/imcoScripts/restCbf/alffCbf_surfCoupling_eb.* /data/jux/BBL/projects/coupling/coupling_test_eb_20200918/code/restCbf_R/.
```
- Now, will go into alffCbf_surfCoupling_eb.Rmd and make some code changes:
- The plan is to take all subjects' list of vertices (10242 per person, 831 people) and make a matrix that is 831 x 10242 (each row is a person, each column is the beta at that vertex from cbf ~ alff)
- Next, I will run a gam model at every vertex across subjects: gam(coupling ~ s(age) + sex + motionalff + motioncbf + mood)
- Each model will be stored in a vector of 10242 (one per vertex)
- I will then make 5 other vectors that correspond to each p value
- Next, I will run an FDR correction for each p
- lastly, I will use adam's data to visualize
#### 20201208
-Crashed the rstudio on chead when I tried to make a matrix that stores the results of the linear models at each of the vertices. Brought it to my local machine to see if that would help
- Will write matrix to one data frame and send it locally to work on further