Skip to content

NF: AFNI 3dNetCorr as afni.NetCorr #3263

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 5 commits into from
Mar 12, 2021
Merged
Show file tree
Hide file tree
Changes from 2 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions nipype/interfaces/afni/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@
LFCD,
Maskave,
Means,
NetCorr,
OutlierCount,
QualityIndex,
ROIStats,
Expand Down
175 changes: 175 additions & 0 deletions nipype/interfaces/afni/preprocess.py
Original file line number Diff line number Diff line change
Expand Up @@ -2556,6 +2556,181 @@ def _format_arg(self, name, trait_spec, value):
return super(TCorrMap, self)._format_arg(name, trait_spec, value)


class NetCorrInputSpec(AFNICommandInputSpec):
in_file = File(
desc="input time series file (4D data set)",
exists=True,
argstr="-inset %s",
mandatory=True)
in_rois = File(
desc="input set of ROIs, each labelled with distinct integers",
exists=True,
argstr="-in_rois %s",
mandatory=True)
mask = File(
desc="can include a whole brain mask within which to "
"calculate correlation. Otherwise, data should be "
"masked already",
exists=True,
argstr="-mask %s")
weight_ts = File(
desc="input a 1D file WTS of weights that will be applied "
"multiplicatively to each ROI's average time series. "
"WTS can be a column- or row-file of values, but it "
"must have the same length as the input time series "
"volume. "
"If the initial average time series was A[n] for "
"n=0,..,(N-1) time points, then applying a set of "
"weights W[n] of the same length from WTS would "
"produce a new time series: B[n] = A[n] * W[n]",
exists=True,
argstr="-weight_ts %s")
fish_z = traits.Bool(
desc="switch to also output a matrix of Fisher Z-transform "
"values for the corr coefs (r): "
"Z = atanh(r) , "
"(with Z=4 being output along matrix diagonals where "
"r=1, as the r-to-Z conversion is ceilinged at "
"Z = atanh(r=0.999329) = 4, which is still *quite* a "
"high Pearson-r value",
argstr="-fish_z")
part_corr = traits.Bool(
desc="output the partial correlation matrix",
argstr="-part_corr")
ts_out = traits.Bool(
desc="switch to output the mean time series of the ROIs that "
"have been used to generate the correlation matrices. "
"Output filenames mirror those of the correlation "
"matrix files, with a '.netts' postfix",
argstr="-ts_out")
ts_label = traits.Bool(
desc="additional switch when using '-ts_out'. Using this "
"option will insert the integer ROI label at the start "
"of each line of the *.netts file created. Thus, for "
"a time series of length N, each line will have N+1 "
"numbers, where the first is the integer ROI label "
"and the subsequent N are scientific notation values",
argstr="-ts_label")
ts_indiv = traits.Bool(
desc="switch to create a directory for each network that "
"contains the average time series for each ROI in "
"individual files (each file has one line). "
"The directories are labelled PREFIX_000_INDIV/, "
"PREFIX_001_INDIV/, etc. (one per network). Within each "
"directory, the files are labelled ROI_001.netts, "
"ROI_002.netts, etc., with the numbers given by the "
"actual ROI integer labels",
argstr="-ts_indiv")
ts_wb_corr = traits.Bool(
desc="switch to create a set of whole brain correlation maps. "
"Performs whole brain correlation for each "
"ROI's average time series; this will automatically "
"create a directory for each network that contains the "
"set of whole brain correlation maps (Pearson 'r's). "
"The directories are labelled as above for '-ts_indiv' "
"Within each directory, the files are labelled "
"WB_CORR_ROI_001+orig, WB_CORR_ROI_002+orig, etc., with "
"the numbers given by the actual ROI integer labels",
argstr="-ts_wb_corr")
ts_wb_Z = traits.Bool(
desc="same as above in '-ts_wb_corr', except that the maps "
"have been Fisher transformed to Z-scores the relation: "
"Z=atanh(r). "
"To avoid infinities in the transform, Pearson values "
"are effectively capped at |r| = 0.999329 (where |Z| = 4.0). "
"Files are labelled WB_Z_ROI_001+orig, etc",
argstr="-ts_wb_Z")
ts_wb_strlabel = traits.Bool(
desc="by default, '-ts_wb_{corr,Z}' output files are named "
"using the int number of a given ROI, such as: "
"WB_Z_ROI_001+orig. "
"With this option, one can replace the int (such as '001') "
"with the string label (such as 'L-thalamus') "
"*if* one has a labeltable attached to the file",
argstr="-ts_wb_strlabel")
nifti = traits.Bool(
desc="output any correlation map files as NIFTI files "
"(default is BRIK/HEAD). Only useful if using "
"'-ts_wb_corr' and/or '-ts_wb_Z'",
argstr="-nifti")
output_mask_nonnull = traits.Bool(
desc="internally, this program checks for where there are "
"nonnull time series, because we don't like those, in "
"general. With this flag, the user can output the "
"determined mask of non-null time series.",
argstr="-output_mask_nonnull")
push_thru_many_zeros = traits.Bool(
desc="by default, this program will grind to a halt and "
"refuse to calculate if any ROI contains >10 percent "
"of voxels with null times series (i.e., each point is "
"0), as of April, 2017. This is because it seems most "
"likely that hidden badness is responsible. However, "
"if the user still wants to carry on the calculation "
"anyways, then this option will allow one to push on "
"through. However, if any ROI *only* has null time "
"series, then the program will not calculate and the "
"user will really, really, really need to address their masking",
argstr="-push_thru_many_zeros")
ignore_LT = traits.Bool(
desc="switch to ignore any label table labels in the "
"'-in_rois' file, if there are any labels attached",
argstr="-ignore_LT")
out_file = File(
desc="output file name part",
name_template="%s_netcorr",
argstr="-prefix %s",
position=1,
name_source="in_file",
)

class NetCorrOutputSpec(TraitedSpec):
out_matrix = File(desc="output text file for correlation stats")

class NetCorr(AFNICommand):
"""Calculate correlation matrix of a set of ROIs (using mean time series of
each). Several networks may be analyzed simultaneously, one per brick.

For complete details, see the `3dNetCorr Documentation
<https://afni.nimh.nih.gov/pub/dist/doc/program_help/3dNetCorr.html>`_.

Examples
--------
>>> from nipype.interfaces import afni
>>> ncorr = afni.NetCorr()
>>> ncorr.inputs.in_file = 'functional.nii'
>>> ncorr.inputs.mask = 'mask.nii'
>>> ncorr.inputs.in_rois = 'maps.nii'
>>> ncorr.inputs.ts_wb_corr = True
>>> ncorr.inputs.ts_wb_Z = True
>>> ncorr.inputs.fish_z = True
>>> ncorr.inputs.out_file = 'sub0.tp1.ncorr'
>>> ncorr.cmdline
'3dNetCorr -prefix sub0.tp1.ncorr -inset functional.nii -mask mask.nii -in_rois maps.nii -ts_wb_corr -ts_wb_Z -fish_z'
>>> res = ncorr.run() # doctest: +SKIP

"""

_cmd = "3dNetCorr"
input_spec = NetCorrInputSpec
output_spec = NetCorrOutputSpec

def _list_outputs(self):
outputs = self.output_spec().get()

if not isdefined(self.inputs.out_file):
prefix = self._gen_fname(self.inputs.in_file, suffix="_netcorr")
else:
prefix = self.inputs.out_file

# All outputs should be in the same directory as the prefix
out_dir = os.path.dirname(os.path.abspath(prefix))

outputs["out_matrix"] = (
fname_presuffix(prefix, suffix="_000", use_ext=False, newpath=out_dir) + ".netcc"
)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is the output file called? I would expect, if I set netcorr.inputs.out_file, the output file is called exactly that.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@effigies : the netcorr.inputs.out_file defines the prefix for the output.
Afni 3dNetCorr doesn't have actual output options, but it uses the prefix to write out the correlation matrix and correlation maps to the output directory.
This added function finds and returns the convention filename for the output correlation matrix.

I am also interested to return the correlation maps return to a subdirectory in the output directory. The number of these correlation maps depends on the number of regions in the input ROI label map. Do you have any idea how I can those to the output as well?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@effigies : Thanks for the approval.
also, do you have any idea on this?
Please let me know if my question needs more clarification

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah, sorry, missed this. Assuming these files exist, you can use glob to get them. See an example here:

class FEAT(FSLCommand):
"""Uses FSL feat to calculate first level stats
"""
_cmd = "feat"
input_spec = FEATInputSpec
output_spec = FEATOutputSpec
def _list_outputs(self):
outputs = self._outputs().get()
is_ica = False
outputs["feat_dir"] = None
with open(self.inputs.fsf_file, "rt") as fp:
text = fp.read()
if "set fmri(inmelodic) 1" in text:
is_ica = True
for line in text.split("\n"):
if line.find("set fmri(outputdir)") > -1:
try:
outputdir_spec = line.split('"')[-2]
if os.path.exists(outputdir_spec):
outputs["feat_dir"] = outputdir_spec
except:
pass
if not outputs["feat_dir"]:
if is_ica:
outputs["feat_dir"] = glob(os.path.join(os.getcwd(), "*ica"))[0]
else:
outputs["feat_dir"] = glob(os.path.join(os.getcwd(), "*feat"))[0]
print("Outputs from FEATmodel:", outputs)
return outputs

Also, can you run make specs and commit the new file?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks! That was very helpful.

return outputs


class TCorrelateInputSpec(AFNICommandInputSpec):
xset = File(
desc="input xset",
Expand Down