Error in quality control barplots, if sample number is divisible by sample.batch.size

[FranziskaKimmig]

Dear RnBeads developers,

Thank you so much for all the work you put into this package!

I have encountered the following error during a preprocessing run of my newest dataset. (The full log file is linked below.)

> ##### Run Analysis
> rnb.run.analysis(dir.reports=report.dir,
+                  sample.sheet=sample.annotation,
+                  data.dir=idat.dir,
+                  data.type=data.type)
2025-01-23 13:26:24     1.7  STATUS STARTED RnBeads Pipeline
2025-01-23 13:26:24     1.7    INFO     Analysis Title: xxx
2025-01-23 13:26:24     1.7    INFO     Initialized report index and saved to index.html
2025-01-23 13:26:24     1.7  STATUS     STARTED Loading Data
2025-01-23 13:26:24     1.7    INFO         Number of cores: 1
2025-01-23 13:26:24     1.7    INFO         Loading data of type "idat.dir"
2025-01-23 13:26:24     1.7  STATUS         STARTED Loading Data from IDAT Files
2025-01-23 13:26:24     1.7    INFO             Added column barcode to the provided sample annotation table
2025-01-23 13:26:27     1.7    INFO             Detected platform: MethylationEPIC
2025-01-23 13:47:12     2.3  STATUS         COMPLETED Loading Data from IDAT Files
2025-01-23 14:46:12     3.2  STATUS         Loaded data from /xxx/
2025-01-23 14:46:28     9.7  STATUS         Added data loading section to the report
2025-01-23 14:46:28     9.7  STATUS         Loaded 500 samples and 866895 sites
2025-01-23 14:46:28     9.7    INFO         Output object is of type RnBeadRawSet
2025-01-23 14:46:28     9.7  STATUS     COMPLETED Loading Data
2025-01-23 15:21:37     3.2    INFO     Initialized report index and saved to index.html
2025-01-23 15:21:54     3.2  STATUS     STARTED Quality Control
2025-01-23 15:21:54     3.2    INFO         Number of cores: 1
2025-01-23 15:21:54     3.2  STATUS         STARTED Quality Control Section
2025-01-23 15:22:21     3.2  STATUS             Added quality control box plots
Error in qc(rnb.set)$Cy3[, sample.subset, drop = FALSE] : 
  subscript out of bounds
Calls: rnb.run.analysis ... FUN -> lapply -> lapply -> FUN -> lapply -> lapply -> FUN
Execution halted

I believe that this error occurs when add.qc.barplots() calls the function rnb.plot.control.barplot() in qualityControl.R, L412-L456. My dataset consists of exactly 500 samples and in this case the sample subsets, that are computed by add.qc.barplots(), are out of bounds considering my sample numbers with the default sample.batch.size==50

RnBeads/R/qualityControl.R

Lines 412 to 422 in 3ab7825

	nsamp<-length(samples(object))
	plot.names<-NULL
	if(nsamp %% sample.batch.size==1){
	sample.batch.size<-sample.batch.size-5
	}
	portion.starts<-0:(nsamp %/% sample.batch.size)*sample.batch.size+1
	portion.ends<-portion.starts+sample.batch.size-1
	portion.ends[length(portion.ends)]<-nsamp
	portions<-paste(portion.starts, portion.ends, sep="-")

# nsamp==500,
# sample.batch.size==50
> portions
[1] "1-50"    "51-100"  "101-150" "151-200" "201-250" "251-300" "301-350" "351-400" "401-450" "451-500" "501-500"

In the following code (lines 424 to 456) rnb.plot.control.barplot() is called for all portions in turn, which gives the subscript out of bounds error above for sample.subset == portions[11]:

RnBeads/R/controlPlots.R

Line 410 in 3ab7825

green<-qc(rnb.set)$Cy3[,sample.subset, drop=FALSE]

Looking at the code this error would occur for all sample sizes divisible by sample.batch.size, in this case 50. My suggestion would be to change that part for example to

nsamp<-length(samples(object))
plot.names<-NULL

if(nsamp %% sample.batch.size == 0){
  portion.starts<-0:(nsamp %/% sample.batch.size - 1)*sample.batch.size+1
}else if (nsamp %% sample.batch.size != 0){
  portion.starts<-0:(nsamp %/% sample.batch.size)*sample.batch.size+1
}
portion.ends<-portion.starts+sample.batch.size-1
portion.ends[length(portion.ends)]<-nsamp
portions<-paste(portion.starts, portion.ends, sep="-")

I have checked that all the iDat files are present and are loaded properly (expected dimensions and no NAs in qc(rnb.set)$Cy3 in rnbSet_unnormalized). The error is reproduced when I split my samples into two batches of 250 samples each for preprocessing, but it runs without error when I use 499 out of my 500 samples. So I am confident that there is no faulty file among the raw data files. As you can see in sessionInfo.txt I used the older version RnBeads 2.20.0, but the code snippets are taken from the most recent version.

Thank you in advance for your help!

sessionInfo.txt
Log_RnBeads_preprocessing.txt

[lutsik]

Hi Franziska, Thanks for spotting the bug, your fix makes sense. I created a dedicated branch: https://github.com/epigen/RnBeads/tree/fix/qc_plot_batches Would you be so kind to create a pull request? Kind regards, Pavlo

…

On Thu, 30 Jan 2025 at 11:07, FranziskaKimmig ***@***.***> wrote: Dear RnBeads developers, Thank you so much for all the work you put into this package! I have encountered the following error during a preprocessing run of my newest dataset. (The full log file is linked below.) > ##### Run Analysis > rnb.run.analysis(dir.reports=report.dir, + sample.sheet=sample.annotation, + data.dir=idat.dir, + data.type=data.type) 2025-01-23 13:26:24 1.7 STATUS STARTED RnBeads Pipeline 2025-01-23 13:26:24 1.7 INFO Analysis Title: xxx 2025-01-23 13:26:24 1.7 INFO Initialized report index and saved to index.html 2025-01-23 13:26:24 1.7 STATUS STARTED Loading Data 2025-01-23 13:26:24 1.7 INFO Number of cores: 1 2025-01-23 13:26:24 1.7 INFO Loading data of type "idat.dir" 2025-01-23 13:26:24 1.7 STATUS STARTED Loading Data from IDAT Files 2025-01-23 13:26:24 1.7 INFO Added column barcode to the provided sample annotation table 2025-01-23 13:26:27 1.7 INFO Detected platform: MethylationEPIC 2025-01-23 13:47:12 2.3 STATUS COMPLETED Loading Data from IDAT Files 2025-01-23 14:46:12 3.2 STATUS Loaded data from /xxx/ 2025-01-23 14:46:28 9.7 STATUS Added data loading section to the report 2025-01-23 14:46:28 9.7 STATUS Loaded 500 samples and 866895 sites 2025-01-23 14:46:28 9.7 INFO Output object is of type RnBeadRawSet 2025-01-23 14:46:28 9.7 STATUS COMPLETED Loading Data 2025-01-23 15:21:37 3.2 INFO Initialized report index and saved to index.html 2025-01-23 15:21:54 3.2 STATUS STARTED Quality Control 2025-01-23 15:21:54 3.2 INFO Number of cores: 1 2025-01-23 15:21:54 3.2 STATUS STARTED Quality Control Section 2025-01-23 15:22:21 3.2 STATUS Added quality control box plots Error in qc(rnb.set)$Cy3[, sample.subset, drop = FALSE] : subscript out of bounds Calls: rnb.run.analysis ... FUN -> lapply -> lapply -> FUN -> lapply -> lapply -> FUN Execution halted I believe that this error occurs when *add.qc.barplots()* calls the function *rnb.plot.control.barplot()* in qualityControl.R, L412-L456 <https://github.com/epigen/RnBeads/blob/3ab7825d89b142191cf43ee75bb0071ff8645717/R/qualityControl.R#L412C3-L456>. My dataset consists of exactly 500 samples and in this case the sample subsets, that are computed by *add.qc.barplots()*, are out of bounds considering my sample numbers with the default *sample.batch.size==50* https://github.com/epigen/RnBeads/blob/3ab7825d89b142191cf43ee75bb0071ff8645717/R/qualityControl.R?plain=1#L412-L422 # nsamp==500, # sample.batch.size==50 > portions [1] "1-50" "51-100" "101-150" "151-200" "201-250" "251-300" "301-350" "351-400" "401-450" "451-500" "501-500" In the following code (lines 424 to 456) <https://github.com/epigen/RnBeads/blob/3ab7825d89b142191cf43ee75bb0071ff8645717/R/qualityControl.R#L424C3-L456> *rnb.plot.control.barplot()* is called for all portions in turn, which gives the subscript out of bounds error above for *sample.subset == portions[11]*: https://github.com/epigen/RnBeads/blob/3ab7825d89b142191cf43ee75bb0071ff8645717/R/controlPlots.R?plain=1#L410 Looking at the code this error would occur for all sample sizes divisible by *sample.batch.size*, in this case 50. My suggestion would be to change that part for example to nsamp<-length(samples(object)) plot.names<-NULL if(nsamp %% sample.batch.size == 0){ portion.starts<-0:(nsamp %/% sample.batch.size - 1)*sample.batch.size+1 }else if (nsamp %% sample.batch.size != 0){ portion.starts<-0:(nsamp %/% sample.batch.size)*sample.batch.size+1 } portion.ends<-portion.starts+sample.batch.size-1 portion.ends[length(portion.ends)]<-nsamp portions<-paste(portion.starts, portion.ends, sep="-") I have checked that all the iDat files are present and are loaded properly (expected dimensions and no NAs in *qc(rnb.set)$Cy3* in rnbSet_unnormalized). The error is reproduced when I split my samples into two batches of 250 samples each for preprocessing, but it runs without error when I use 499 out of my 500 samples. So I am confident that there is no faulty file among the raw data files. As you can see in *sessionInfo.txt* I used the older version RnBeads 2.20.0, but the code snippets are taken from the most recent version. Thank you in advance for your help! sessionInfo.txt <https://github.com/user-attachments/files/18560443/sessionInfo.txt> Log_RnBeads_preprocessing.txt <https://github.com/user-attachments/files/18600209/All_RnBeads_preprocessing.txt> — Reply to this email directly, view it on GitHub <#61>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACQTCLQ64VANI62VJAO6TJT2NH2UTAVCNFSM6AAAAABWE3N3IGVHI2DSMVQWIX3LMV43ASLTON2WKOZSHAZDAMZZHAZTOMY> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[lutsik]

Dear Franziska, please test whether your issue has been fixed in this branch. If so, we can merge it into devel.

[FranziskaKimmig]

Hi Pavlo,

Please excuse the late reply.
I installed RnBeads from fix/qc_plot_batches and ran the preprocessing step for my samples again with the same setup as before (500 samples and rnb.options(qc.sample.batch.size = 50)) to test if everything works fine now. However, I got the following error during loading even before the qc step (Log_RnBeads_preprocessing_test.txt, sessionInfo.txt):

> ##### Run Analysis
> rnb.run.analysis(dir.reports=report.dir,
+                  sample.sheet=sample.annotation,
+                  data.dir=idat.dir,
+                  data.type=data.type)
2025-03-06 17:04:05     1.7  STATUS STARTED RnBeads Pipeline
2025-03-06 17:04:05     1.7    INFO     Analysis Title: xxx
2025-03-06 17:04:05     1.7    INFO     Initialized report index and saved to index.html
2025-03-06 17:04:05     1.7  STATUS     STARTED Loading Data
2025-03-06 17:04:05     1.7    INFO         Number of cores: 1
2025-03-06 17:04:05     1.7    INFO         Loading data of type "idat.dir"
2025-03-06 17:04:05     1.7  STATUS         STARTED Loading Data from IDAT Files
2025-03-06 17:04:05     1.7    INFO             Added column barcode to the provided sample annotation table
2025-03-06 17:04:08     1.7    INFO             Detected platform: MethylationEPIC
2025-03-06 17:04:08     1.7    INFO             Annotation package genome assembly version: hg19
2025-03-06 17:23:21     2.2  STATUS         COMPLETED Loading Data from IDAT Files
2025-03-06 18:05:23     3.3  STATUS         Loaded data from xxx
2025-03-06 18:05:40     9.8   ERROR         Undefined platform; unexpected or missing IDAT files
Error in logger.error(txt) : 
  Undefined platform; unexpected or missing IDAT files
Calls: rnb.run.analysis ... rnb.detect.infinium.platform -> rnb.error -> logger.error
Execution halted

As far as I could see rnb.detect.infinium.platform() is called twice, but the second time without the argument recursive=TRUE
Compare the first call

RnBeads/R/dataImport.R

Line 679 in 460c356

fn.available<-list.files(base.dir, pattern="idat", recursive=TRUE, full.names=TRUE)

with the second one

RnBeads/R/loading.R

Line 421 in 460c356

platform <- rnb.detect.infinium.platform(list.files(path = data.source[[1]], full.names = TRUE, pattern = 'idat'))

After I added recursive=TRUE in the second call, everything worked fine, so I think this was the issue. I also saw that the same error during loading was reported previously in issue #32, so maybe this fixes it as well.

Otherwise, the fix for indexing the sample batches seems to work without problems.

Kind regards,
Franziska