mouse-imaging-centre/rminc

Statistics for MINC volumes: A library to integrate voxel-based statistics for MINC volumes into the R environment. Supports getting and writing of MINC volumes, running voxel-wise linear models, correlations, etc.; correcting for multiple comparisons using the False Discovery Rate, and more. With contributions from Jason Lerch, Chris Hammill, Jim Nikelski and Matthijs van Eede. Some additional information can be found here:

rminc's Issues

Convert to roxygen

1. Add explicit @export directives to roxygen blocks
2. Add roxygen blocks to all user visible functions
3. Remove old NAMESPACE to have it kept synchronised with current package state (we have 285 export statements but only 83 documented exported functions, 55 of which are generated by roxygen2).

provide better feedback on corrupt MINC files

When MINC files get corrupted, you can see error messages like this when loading them:

DF5-DIAG: Error detected in HDF5 (1.8.10-patch1) thread 0:
#000: H5Dio.c line 174 in H5Dread(): can't read data
major: Dataset
minor: Read failed
#1: H5Dio.c line 449 in H5D__read(): can't read data
major: Dataset
minor: Read failed
#2: H5Dchunk.c line 1735 in H5D__chunk_read(): unable to read raw data chunk
major: Low-level I/O
minor: Read failed
#3: H5Dchunk.c line 2766 in H5D__chunk_lock(): data pipeline read failed
major: Data filters
minor: Filter operation failed
#4: H5Z.c line 1120 in H5Z_pipeline(): filter returned failure during read
major: Data filters
minor: Read failed
#5: H5Zdeflate.c line 125 in H5Z_filter_deflate(): inflate() failed
major: Data filters
minor: Unable to initialize object
/projects/mice/share/arch/linux-3_2_0-36-generic-x86_64-eglibc-2_15/src/minc-toolkit/libminc/libsrc2/hyper.c:794 (from MINC2): HDF5 function H5Dread failed
Error in mincGetVolume(testfiles$filename[i]) : Could not get hyperslab.

Currently the RMINC functions (mincLm for instance) will die on this without indicating which file from your list of input files was the culprit.

make mincFDR work with output from mincLme

Add the ability for RMINC to install minc libraries if not present.

This was a suggestion by Vlad. If people receive some MINC files they want to analyze, but do not have MINC installed, this could be taken care off during the RMINC installation.

BUG: RMINC fails with latest lmerTest (2.0-20)

There was a code re-factoring in 2.0-20 which renamed some internal functions that RMINC is using from lmerTest. Version 2.0-11 works fine.

Particularly rhoInitthough there may be others, that was the first bug I ran into.

All .o and .so files should be removed prior to (re-)compiling/installing the RMINC library

anatGetAll call without "defs"

It should be possible to run anatGetAll (and other functions that use the same procdure) without a definitions file. This would make it a lot easier to provide a quickly made segmentation (or mask) to this function.

FEATURE REQUEST: vertexLmer like mincLmer

mincLmer is available for linear mixed effect models, it would be nice of there were a vertex version as well.

anatGetAll help displays without newlines

I get a solid blob where paragraph breaks should be:

 jacobians Each file contains log jacobians, and the volume for
 each atlas label is computed by multiplying the jacobian with the
 voxel volume at each voxel.  labels Each file contains integer
 labels (i.e. same as the atlas).  The volume is computed by
 counting the number of voxels with each label and multiplying by
 the voxel volume.  means Each file contains an arbitrary number
 and the mean of all voxels inside each label is computed.  sums
 Each file contains an aribtrary number and the sum of all voxels
 inside each label is computed.  text Each file is a comma
 separated values text file and is simply read in.

provide better feedback when mixing files with different dimension order

Currently if you run a mincLm call with input files have the zspace,yspace,xspace dimension order, and the mask you provide has a different order, you get a somewhat obscure error message:

Method: lm
Number of volumes: 20
Volume sizes: 172 340 225
N: 20 P: 2
In slice
HDF5-DIAG: Error detected in HDF5 (1.8.10-patch1) thread 0:
#000: H5Dio.c line 153 in H5Dread(): selection+offset not within extent
major: Dataspace
minor: Out of range
/projects/mice/share/arch/linux-3_2_0-36-generic-x86_64-eglibc-2_15/src/minc-toolkit/libminc/libsrc2/hyper.c:794 (from MINC2): HDF5 function H5Dread failed
Error in mincLm(gf$absolute_jacobians_0.5 ~ gf$Genotype, mask = "/some_mask.mnc") :
Error opening mask buffer.

We should do some file checking to overcome these kinds of error messages.

Option to ignore NA values in mincSummary functions

mincSummary should have an option to ignore NAs. The default is that this option will NOT be used. If it is used, then any NA values in the data will be ignored when mean, sum, sd, or var are calculated. This should mirror the default R behaviour.

If NAs are not ignored, and there are NAs in the data, mincSummary should give a message to the user, rather than simply ignoring the NA values.

Downdated VtV is not positive definite error causes mincLmer to crash without output

Running simple voxelwise mixed effect fitting:

model<-mincLmer(log_det~cohort+(1|subj),data=dbm,mask='model_nl_c/avg.016_mask_lr.mnc',parallel=c("snowfall", 4))

works for a while and then prints:

Error in checkForRemoteErrors(val) :
4 nodes produced errors; first error: Downdated VtV is not positive definite

and produces no output

FEATURE REQUEST: Provide residuals from vertexLm

Are the residuals from vertexLm stored in the object somewhere? I don't see anything immediately obvious.

make MINC compression work when writing out MINC files

Doesn't seem to work now in the presence of the MINC_COMPRESS environment variable.

the test for writeVertex does not actually test anything

or at least, it does not use any of the "test_that" tests

anatFDR (and possible others) fail when in the input contains (an) NA value(s)

If an anatLm call produces an NA in its output, and this is fed into anatFDR, the latter will return NA for all its thresholds, even though excluding the NA value might give you completely different answers. This is clearly a bug. The presence of NA values should be checked somewhere (either in anatLm or anatFDR) and dealt with by providing a proper error message

sfInit, should that really happen inside pMincApply?

Currently in the development branch of RMINC, sfInit() is called when the method for pMincApply is snowfall. However, it looks like this means that user doesn't have any control anymore over how many cpus are used. This could be remedied by feeding the number of workers (=cpus) to sfInit, however, it's still possible that you'll need to set things through sfExport and sfLibrary. I think all these snowfall related calls should be left to the user before calling pMincApply.

Print functions for output from anatGetAll (and others?)

The mincMultiDim classes have a print function. This should probably also be implemented for some of the outputs of anat* functions. Right now, printing the output of anatGetAll gives you an overwhelming amount of numbers.

Speed up mincFDR

Investigate slowness of mincFDR

problem using function miopen_volume

Sorry for asking this here but I don't have any other way to fix my problem.
After installing minc2 library when i try open a minc image using miopen_volume I always get an error opening the file.
I'm using a minc file from brainweb database from BIC. Here is my code (almost equal to tutorial 1):

include <minc2.h>

include <stdio.h>

int main(int argc, char *argv) {
mihandle_t minc_volume;
double voxel;
int result;
unsigned long location[3];
//char name[74]="/home/patrick/Desktop/06_20120908_114444/06_20120908_114444_201e1_mri.mnc";
/ open the volume - first and only command line argument */

result = miopen_volume("/home/patrick/t1.mnc", MI2_OPEN_READ, &minc_volume);
/* check for error on opening */
if (result != MI_NOERROR) {
fprintf(stderr, "Error opening input file: %d.\n", result);
}

/* get the first voxel in the image */
location[0] = location[1] = location[2] = 0;
miget_real_value(minc_volume, location, 3, &voxel);
printf("Hello MINC2 - first voxel was: %f\n", voxel);
//result=miclose_volume(minc_volume);
return(0);
}

Can the error opening the file happen because some mistake i did installing?
Again sorry for asking this here but can you help me with this?

Feature Request: Allow specification of blur level for civet.getAllFilenames

Right now 20mm is hardcoded. For backwards compatability between CIVET 1.1.12 and 1.1.10, we need to specify rms 28.28 because of the off by sqrt(2) bug in CIVET

See http://www.bic.mni.mcgill.ca/pipermail/minc-users/2015-January/003979.html

Defaults should be in /projects/mice/share if they aren't already

Example: in anatGetAll, the default location for the definitions is in /projects/mice/jlerch. This should be changed here, and anywhere else the default is not in /projects/mice/share.

Visualization tool: vertexHistogram

Just thought I'd dump this here, might be of use (or maybe addition).

Visualizes the distribution of statistics for a vertexLm:

vertexHistogram = function(linmodel, fdrmodel, column ){
  #Generates a nice histogram of vertexLm t-values, with FDR levels labelled
  legendline = c(c("FDR 5%", attributes(fdrmodel)$thresholds["0.05",column]), c("MAX", max(linmodel[,column])), c("MIN", min(linmodel[,column]) )   )
  hist(linmodel[,column],breaks = seq(floor(min(linmodel[,column])),ceiling(max(linmodel[,column])), by=0.25),
       xlab = paste( deparse(substitute(linmodel)), ",", column), main =  paste( deparse(substitute(linmodel)), ",", column))
  legend("topleft", legend = legendline)
  lim <- par("usr")
  rect(lim[1], lim[3], -attributes(fdrmodel)$thresholds["0.01",column], lim[4], border = '#00FF0055', col = '#00FF0077')
  rect(attributes(fdrmodel)$thresholds["0.01",column], lim[3], lim[2], lim[4], border = '#00FF0055', col = '#00FF0077')

  rect(-attributes(fdrmodel)$thresholds["0.01",column], lim[3], -attributes(fdrmodel)$thresholds["0.05",column], lim[4], border = '#00FF0055', col = '#00FF0055')
  rect(attributes(fdrmodel)$thresholds["0.05",column], lim[3], attributes(fdrmodel)$thresholds["0.01",column], lim[4], border = '#00FF0055', col = '#00FF0055')

  rect(-attributes(fdrmodel)$thresholds["0.05",column], lim[3], -attributes(fdrmodel)$thresholds["0.1",column], lim[4], border = '#FFFF0055', col = '#FFFF0055')
   rect(attributes(fdrmodel)$thresholds["0.1",column], lim[3], attributes(fdrmodel)$thresholds["0.05",column], lim[4], border = '#FFFF0055', col = '#FFFF0055')

  rect(-attributes(fdrmodel)$thresholds["0.1",column], lim[3], -attributes(fdrmodel)$thresholds["0.2",column], lim[4], border = '#FF000055', col = '#FF000055')
  rect(attributes(fdrmodel)$thresholds["0.2",column], lim[3], attributes(fdrmodel)$thresholds["0.1",column], lim[4], border = '#FF000055', col = '#FF000055')
}

Usage:

> vertexHistogram(vs.quad.left.gi, fdr.quad.left.gi, "tvalue-poly(Age, 2)1" )

Example:

masking (maskval) in minc* functions inconsistent

minc* functions seem to specify maskval = NULL in their default options list, however mincLm then checks for null and sets a minimum and maximum value of 1 and 99999. This is undocumented.

In the meantime, mincTtest also does maskval = NULL, and really does mean null in this case. Looking inside mincSummary however, NULL is then overridden to the values for mincLm.

This is a bit confusing looking at the code/docs.

check for and load lme4 library within mincLmer

anatGetAll should stop when shape of first input file and atlas are incompatible

remove inf, na, nan from volume before writing (mincWriteVolume)

Currently you get a very cryptic error message when you are trying to write a file out to disk where in R, the data contains either inf, na, or nan. This should be captures in mincWriteVolume. Either it should set these data points to 0, or tell the user what is happening and how it could be solved.

getRMINCTestData() can't access wiki

I think the site is long gone now, so the example data is no longer accessible. Can anyone provide an alternative link?

vertexAnova is broken

The F-statistics from a vertexLm and vertexAnova call in a simple two way test should be identical. This is currently not the case. Needs to be fixed.

Updates to mincApply and pMincApply

Currently, when creating a function to be used by mincApply or pMincApply, the function must have only one argument (x), which indicates that the function must be executed on every voxel within the specified mask.

If possible, and to increase the reuseability of functions created for mincApply and pMincApply, it would be nice if they could take additional arguments. This would reduce the amount of hard-coded parameters within each function.

This might not be feasible and/or may require modification of minc2_model, which is probably undesirable.

At the very least, we should think about if this is possible and how to accomplish it.

mincFDR - method "qvalue", "pFDR: and "fastqvalue" are broken-ish

The qvalue package has been removed from CRAN, and can not be installed anymore through install.packages("qvalue"). This currently breaks the above mentioned methods for mincFDR. Also, when mincFDR is called using a .mincSingleDim the df and statType arguments are not properly parsed/used. In the main mincFDR function they are overwritten right away by looking at the attributes "df" and "stat-type". If these are not present in the minc file object, mincFDR will fail.

Create verbose/quiet option for runRMINCTestbed

When we do our own testing, we should be able to see all output of the commands we run in order to catch strange warnings/errors. When an end user installs RMINC and runs the testbed, it's fine to sink all that output to /dev/null so they just see whether the tests return the correct values etc., but with a single switch we should be able to run the testbed and see all the output (i.e., not do any sink())

writeVertex Fails when header is set to TRUE

inconsistency in allowed masks for various functions

When using mincFDR or mincAnova, a set of labels (ie not a binary mask) can be specified as the mask.

When using mincLm, a binary mask is required or the calculation fails.

I think that this masking behaviour should be consistent across all functions. Either everything requires a binary mask or everything can use labels (or a binary mask). [Unless of course the current behaviour is deliberate, in which case this should be documented somewhere.]

mincGetVoxel and mincGetWorldVoxel fail on float MINC input files

Where fail means that values are retrieved, but they are incorrect.

anatGetAll(method="means") fails on some MINC input files

The problem seems to come from the program compute_values_across_segmentation:

atlas = volumeFromFile(args[1], dtype="ushort")

did not work for the Allen Brain atlas file.

rename runTestbed(), and remove temp files

The way to run the RMINC testbed is by calling the function runTestbed(). This function should be renamed to runRMINCTestbed() (or something along those lines). Also, all temp files and the data used for the test bed should be removed. We can use the unlink() function for that.

Additional mixed effects functions

Currently, RMINC has a mincLme function, which allows you to calculate a linear mixed effects model at every voxel, but only allows for one fixed effect and one random effect.

As we move forward with additional longitudinal studies and more complex modelling, it would be nice (if possible) to allow for multiple fixed and random effects in this (or another) function.

mincLme (using the nlme package) is broken, needs fixing

I don't have a full set of commands that reproduces this error, but a call to mincLme can currently return something like this:

mincLme(ss$filename, x ~ Group * Scan, random.effect=~1|Exp_ID, mask=mask)
Error in .C("matrixLog_pd", Factor = double(Ncol * Ncol), as.integer(Ncol), :
"matrixLog_pd" not available for .C() for package "nlme"

Feature Request: "mask" for vertexLm and related vertex operations

CIVET (for example) has a region of "known bad" measures which I would never interpret as being significant in a t-stat map.

Right now what I do is zero out the t-stat file RMINC creates with a mask (flat text file) file of "0" and "1" the length of the vertex file using bc.

It might be useful to allow for masking of vertexLm (and related functions) with such a mask file, and fill the masked region with either 0 (preferred) or NA (probably breaks some viewing software).

parallelize parametric bootstrap code

The parametric bootstrap code is quite slow - should be easy to parallelize using snowfall.

Add input filenames as an attribute to anatGetAll output

Finish historic information about version changes in the NEWS file

Purge R/minc_support, R/deprecated

As far as I can tell these files aren't used or documented, would anyone object to removing them?

pMincApply Enhancements

1. Add some sort of progress bar

2. Load mask file into memory rather than saving to file

Provide a more meaning full error when loading MINC1 files

Currently when a MINC1 file is loaded, the error message provided is:

Error opening volume: filename.mnc

Which is quite confusing, given that all other MINC tools can properly open this file. It would be better to have something along the lines:

Input file filename.mnc is MINC1. RMINC only reads MINC2, please use mincconvert to convert this file to MINC2

BUG: Can't access data from vertexFDR

vertexFDR will happily display it's results, but there is no clear way to extract the F-statistic threshold values from it's calculation

anatLm allows for only up to 2 terms

E.g, the following works:

anatLm(~ Genotype + Sex, data_frame, volume_data)

but this doesn't work:

anatLm(~ Genotype + Sex + Age, data_frame, volume_data)

The error message you'll get will be as follows:

Error in parse(text = rCommand) : :1:14: unexpected '+'
1: term <- data$+

Reading more than 1020 MNC files

An error occures when you try to read more than 1020 MNC files.

gf<-read.csv("/projects/muisjes/dfernandes/enrichment/scripts/scrap/file_for_R_analysis3.csv"); nr <- length(gf[,1]); for (i in 1:nr) {gene=mincGetVolume(gf[i,3]);print(i)}

The code will run fine for the first 1019 files and then refuse to open any subsequent files (MINC or otherwise).

If you try to open another non-MINC file after this error, you get a "cannot open the connection"/"Too many open files" error.

I tried the closeAllConnections in the R session, but it did not seem to help.

RMINC depends upon https://github.com/Mouse-Imaging-Centre/minc-stuffs

Particularly label_volumes_from_jacobians

Discovered following https://wiki.mouseimaging.ca/display/MICePub/Pydpiper+Virtual+Machine