What can cause `dim(Red_grid) != dim(detP_sink_grid)` in detectionP

[rajwanir]

Hi,

I have a large dataset to analyze (27,000+ scans). To speed things up and also constraining memory, I split the comple list of files in batches and use read.metharray2 to create hd5 for each batch. Eventually, I tried simply loading the DelayedArray for each file and run cbind to create a combined DelaydMatrix per assay to create an RgChannelSet. The RgChannelSet creation works fine and some of the other functions. However, detectionP raises an error at the following:

minfi/R/detectionP.R

Lines 87 to 88 in 6e9859a

	stopifnot(dim(Red_grid) == dim(Green_grid),
	dim(Red_grid) == dim(detP_sink_grid))

Any idea what could cause this dim(Red_grid) != dim(detP_sink_grid) , the block arrangement during writes and chunking? Does the idea for creating merged DelayedMatrix backed by multiple datasets is sound? Or it is logically incompatible for how minfi expects a DelayedMatrix for each assay.

Thank you for your help.

[PeteHaitch]

Are you able to create a minimal reproducible example (reprex)?
It's many years since I've used minfi, but with a reprex I might be able to help.

[rajwanir]

Thanks a lot @PeteHaitch. I could resolve this issue. The dim(Red_grid) != dim(detP_sink_grid) was simply my bad. I didn't align the order of the intermediate batches with read.metharray2 before doing a do.call(cbind,allDelayedArrayBatches).

I have a different general question though if I may, do you know if many of the functions in minfi that uses DelayedArray are compatible with a MultiCoreParam or SnowParam instead of the default SerialParam? The reason I asked is that I see many of the functions written very nicely to use DelayedArray and looks it should Parallelize but they don't. Probably because the HDF5Array backend itself does not support parallel writes unless the much of the related packages (rhd5,DelayedArray) are compiled with MPI. I wonder if you recall testing the minfi with Parallel support and if it was with MPI?

With the dataset size (27,000+), I probably would need both memory efficiency and parallel support. So I am sort of hacking my way through different functions in minfi but seem like a wrong approach if by default those functions are supposed to work with Parallel processing.

For instance, in detectionP BPPARAM is not exposed. It calls on an internal .detectionP based on the signature where BPPARAM is fixed to SerialParam. Even if try to export internal .dectectionP and try multiCoreParam, it doesn't work. Probably because it tries to call write_block within blockMapplyWithRealization which results in output matrix being silently NA when trying to process in Parallel.

Here is the internal .detectionP exported:

internal_detectionP <- getFromNamespace(".detectionP", "minfi")
.isRGOrStop <- getFromNamespace(".isRGOrStop", "minfi")
detectionP_Parallel <- function(rgSet,
                              BPPARAM = MulticoreParam(workers=4)) {
  # Input validation
  .isRGOrStop(rgSet)
  

  # Extract required components
  locusNames  <- getManifestInfo(rgSet, "locusNames")
  controlIdx  <- getControlAddress(rgSet, controlType = "NEGATIVE")
  Red         <- getRed(rgSet)
  Green       <- getGreen(rgSet)
  TypeI.Red   <- getProbeInfo(rgSet, type = "I-Red")
  TypeI.Green <- getProbeInfo(rgSet, type = "I-Green")
  TypeII      <- getProbeInfo(rgSet, type = "II")
  
  # Coerce to DelayedMatrix
  Red   <- as(Red, "DelayedMatrix")
  Green <- as(Green, "DelayedMatrix")
  
  setAutoRealizationBackend(NULL)
  # Call the internal detectionP function with custom BPPARAM
  detP <- internal_detectionP(
    Red         = Red,
    Green       = Green,
    locusNames  = locusNames,
    controlIdx  = controlIdx,
    TypeI.Red   = TypeI.Red,
    TypeI.Green = TypeI.Green,
    TypeII      = TypeII,
    BPPARAM     = BPPARAM
  )
  
  return(detP)
}

Here is my calls:

detectionP_Parallel(rgsetAll,MulticoreParam(workers=1,progressbar = TRUE))
  |=========================================================================| 100%

<866238 x 40> DelayedMatrix object of type "double":
           206693320168_R01C01 ... 206458350079_R08C01
cg18478105       4.582797e-175   .        8.962750e-77
cg09835024        3.823120e-03   .        1.894521e-39
cg14361672        5.549582e-15   .        1.208187e-07
cg01763666        5.846476e-23   .        4.314336e-03
cg12950382        9.119806e-03   .        1.231708e-01
       ...                   .   .                   .
cg23079522       2.627541e-306   .        6.011207e-59
cg16818145        0.000000e+00   .        0.000000e+00
cg14585103        0.000000e+00   .        0.000000e+00
cg10633746        0.000000e+00   .        0.000000e+00
cg12623625        0.000000e+00   .        0.000000e+00
> detectionP_Parallel(rgsetAll,MulticoreParam(workers=2,progressbar = TRUE))
  |=========================================================================| 100%

<866238 x 40> DelayedMatrix object of type "double":
           206693320168_R01C01 ... 206458350079_R08C01
cg18478105                  NA   .                  NA
cg09835024                  NA   .                  NA
cg14361672                  NA   .                  NA
cg01763666                  NA   .                  NA
cg12950382                  NA   .                  NA
       ...                   .   .                   .
cg23079522                  NA   .                  NA
cg16818145                  NA   .                  NA
cg14585103                  NA   .                  NA
cg10633746                  NA   .                  NA
cg12623625                  NA   .                  NA

Do have any suggestions or general advice on working with large datasets (27,000+) with minfi? Do you have recommendation on any configuration (e.g. the rhd5+MPI build) and everything might work out of the box? Or hacking modifying different functions like I am trying is the best bet or perhaps throwing a large memory and just being patient with SerialParam might be the most effective solution?

Thank you very much.

[PeteHaitch]

Hi @rajwanir,

No, I'm afraid I don't have any experience running minfi on large datasets.
You are right, though, that HDF5 files don't play nice with parallel access.

We explored implementing DelayedArray-backed minfi objects at the time I was adding these to bsseq, but never got as far as we did with bsseq.
Unfortunately, other priorities meant that the work in minfi basically stopped as far as I know.

@kasperdanielhansen might have more recent experience using minfi on large datasets and advice on strategies.
But this might be a better question for https://support.bioconductor.org/

Cheers,
Pete