[Question]: how can I get the effect size in tracts?

[yutinghe99]

Hi, @smeisler ,
Thanks for your sharing, this code really helps a lot!
After running fixelcfestats and mrthreshold, I got the significant fixels among groups, now I want to assign these fixels to tracts so that I can extract the mean FBA metrics within tracts to run the post-hoc test. I saw your scripts 10_itersect_with_tracts.sh can help me achieve this goal. But I noticed you ran GAMs with ModelArray instead of fixelcfestats. The outputs from fixelcfestats didn't include the effect size. How can I get the max effect size achieved in tract?

Background: I run the one-way ANOVA with fixelcfestats, here are my degisn and contrast matrix.
For one-way ANOVA:
Design matrix: HCs subgroup1 subgroup2 sex age TIV (sex, age, and TIV demeaned)
Contrast matrix: C1: 1 -1 0 0 0 0
C2: 0 1 -1 0 0 0
F matrix: [1 1]
Code: fixelcfestats log_fc_smooth/ files.txt design_matrix_F.txt contrast_matrix_F.txt matrix/ stats_log_fc_F/ -ftests ftest.txt -fonly

Below are my outputs

One-way ANOVA didn't get the abs_effect.mif, not like the T-test. So how can I calculate the max effect size achieved in tract? Could you please guide me?

Best,
Yuting

[yutinghe99]

New questions:
Follow your bash code 8b and 10, I ran the following lines

#!/bin/bash
echo 'generate fixels from tract tcks'
ts_out=template/tractseg
#mkdir -p $ts_out/tck_fixels
for tck in $ts_out/TOM_trackings/*.tck
do tract=$(basename "$tck")
$mrtrix tck2fixel $tck template/fixel_mask $ts_out/tck_fixels $tract.mif
done

echo 'binarizing tract fixel masks'
for tract in $ts_out/tck_fixels/*.mif
do $mrtrix mrthreshold $tract ${tract//.mif}_bin.mif -abs 1 -force
done

These run successfully, I got the binarizing tract fixel masks. Then I ran the bash code 10, but delete the max effect part.

#!/bin/bash
# creat output file
out_file=template/stats_fd_anova/tract_fixels.txt
touch $out_file
# loop over tract masks
for track in template/tractseg/tck_fixels/*bin*;
	# calculate proportion of tract occupied by significant fixels
	do proportion=$($mrtrix mrstats template/stats_fd_anova/sig_fixels_0.95.mif -mask $track -output mean);
	# get number of significant fixels per tract
	count_tract=$($mrtrix mrstats template/stats_fd_anova/sig_fixels_0.95.mif -mask $track -output count);
	# get number of total fixels per tract
	num_fixels=$(echo "$proportion * $count_tract" | bc);
	# add information to text file
	echo $(basename $track) >> $out_file;
	printf "%.0f\n" $num_fixles >> $out_file;
	echo $proportion >> $out_file;
done

However errors occurred:
mrstats: [ERROR] dimension mismatch between "template/tractseg/tck_fixels/directions_bin.mif" and "template/stats_fd_anova/sig_fixels_0.95.mif" between axes 0 and 2 (804320,3,1 vs. 804320,1,1) mrstats: [ERROR] dimension mismatch between "template/tractseg/tck_fixels/directions_bin.mif" and "template/stats_fd_anova/sig_fixels_0.95.mif" between axes 0 and 2 (804320,3,1 vs. 804320,1,1) (standard_in) 1: syntax error mrstats: [ERROR] dimension mismatch between "template/tractseg/tck_fixels/index_bin.mif" and "template/stats_fd_anova/sig_fixels_0.95.mif" between axes 0 and 2 (152,181,147,2 vs. 804320,1,1) mrstats: [ERROR] dimension mismatch between "template/tractseg/tck_fixels/index_bin.mif" and "template/stats_fd_anova/sig_fixels_0.95.mif" between axes 0 and 2 (152,181,147,2 vs. 804320,1,1) (standard_in) 1: syntax error

Do you have any idea of the error?

[smeisler]

Hi @yutinghe99 and thank you for opening this issue. I'm happy the code has worked well for you up to this point!

For the effect size, you would use the betaX.mif file, where X corresponds to the relevant column in your design matrix. So if your term of interest was the first column of the design matrix, you would use beta0.mif, second column would be beta1.mif, and so on.

I'll try to get a solution to the syntax error soon. I've never tried this code with fixelcfestats, maybe there is something different about output formats?

Best,
Steven

[yutinghe99]

Thanks for your quick reply.

For the effect size, you would use the betaX.mif file, where X corresponds to the relevant column in your design matrix.

Sorry, I'm bad at statistics. Could you please explain this more? The first three columns of my design matrix are my independent variable (i.e. three different groups: healthy controls, patients subgroup 1, and patients subgroup 2). I performed one-way ANOVA to test the FBA metrics differences among these three groups. If I calculate the max effect size achieved in tract with beta0.mif (i.e. the healthy controls column), what does the output mean?

[yutinghe99]

Hi, @smeisler

I found a post related to the syntax error, he said that he misused the mask output from tckmap. But with my understanding of the bash scripts 10, we have already used the mask generated by tck2fixel, right?

If you need more background information on my analysis, please let me know. Looking forward to your reply.

Cheers,
Yuting

[smeisler]

Oh, about the errors, you can ignore the ones for index and directions mif files, you only need the outputs for the actual metric files e.g., log_fc.

[smeisler]

For the ANOVA effect size, the F value might be a good proxy. As the test statistic, it should be effect size normalized by sample variance.

[yutinghe99]

Thanks for your sincere help! :)

To make sure that I understanding correctly, I will do some confirmations:

1. I can temporarily remove the index.mif and directions.mif from the template/tractseg/tck_fixels file to solve the syntax errors, right?

2. In my case, I should choose the Fvalue.mif to calculate the max effect size, the code should be:

#!/bin/bash
# creat output file
out_file=template/stats_fd_anova/tract_fixels.txt
touch $out_file
# loop over tract masks
for track in template/tractseg/tck_fixels/*bin*;
	# calculate proportion of tract occupied by significant fixels
	do proportion=$($mrtrix mrstats template/stats_fd_anova/sig_fixels_0.95.mif -mask $track -output mean);
	# get number of significant fixels per tract
	count_tract=$($mrtrix mrstats template/stats_fd_anova/sig_fixels_0.95.mif -mask $track -output count);
	# get number of total fixels per tract
	num_fixels=$(echo "$proportion * $count_tract" | bc);
        # Get max effect size achieved in tract
	max_eff_tract=$($mrtrix mrstats template/stats_fd_anova/Fvalue.mif -mask $track -output max);
	# Add information to text file
	echo $(basename $track) >> $out_file;
	printf "%.0f\n" $num_fixles >> $out_file;
	echo $proportion >> $out_file;
        echo $max_eff_tract >> $out_file;
done

If I am incorrect, please point out. Thank you!

[yutinghe99]

Hi @smeisler ,

After I ran the bash code above, it went without error. But the tract_files.txt seems weird.

Based on my understanding of this script, the tract_files.txt contains 4 rows for each tract.:
The first row is the name of the tract, the second is the number of significant fixels per tract, the thrid is the proportion of tract occupied by significant fixels, and the last one is the max effect size achieved in the tract.

But in my tract_files.txt the second row (i.e. num_fixles) is always 0 for all tracts, and some of them don't have information on the fourth row (i.e. max_eff_tract). Is this reasonable?

I have attached my tract_files.txt here. Could you help me check it? tract_fixels.txt

I noticed that in your elife paper, you said that

To ascribe significant fixels to tracts, we intersected significant fixels (qFDR < 0.05) and the binarized tract fixel masks.

How do you achieve this? Based on the max effect size to decide the tract?

Cheers,
Yuting

[smeisler]

Hi,

No, it's not expected to have 0 for num_fixels. Have you tried visualizing the fixel masks in mrview?

And to answer the last question, the script we are talking about now is what I meant by that intersection. Looking at the results table you can see I report the number of fixels significant at pFDR < 0.05 and max effect size achieved within a fixel mask.

For the code you provide above, have you tried running the individual parts line by line to see if the outputs are what you expect?

[yutinghe99]

No, it's not expected to have 0 for num_fixels. Have you tried visualizing the fixel masks in mrview?

Yes, I have tried visualizing the fixel masks in mrview, but I am not sure if I visualized it correctly...

Additionally, I just found that my sig_fixels_0.95.mif was also a fixel mask.

Is this the reason for the 0 num_fixel?
I created the sig_fixels_0.95.mif by performing:
mrthreshold stats_fd_anova/fwe_1mpvalue.mif stats_fd_anova/sig_fixels_0.95.mif -abs 0.95 -comparison gt

[smeisler]

What are values you get for proportion and count_tract?

[yutinghe99]

Hi @smeisler ,

I'm sorry for not replying to you sooner. Following your suggestion, I output the proportion and count_tract respectively. I found both of them looked sensible. Then I scrutinized my code...there is a typo between num_fixels and num_fixles....
Now I can get the correct num_fixels. Sorry for my negligence.

But now I have another question. I want to intersect the sig_fixels.mif with tract_bin.mif to get the tract mask containing only the significant part. Do you have any idea about which command I should use?

Cheers,
Yuting

[smeisler]

From https://community.mrtrix.org/t/fba-creating-intersection-fixel-mask/1917, mrcalc mask1.mif mask2.mif -min intersectionmask.mif

[smeisler]

Also, now that the typo was fixed and the original question addressed, mind if I close this issue?

[yutinghe99]

Of course not. Thank you for your help!