FRASER Error in check_returned_array

[lemdcock]

Dear C. Mertes,

I'm running FRASER on 61 PBMC samples for the detection of aberrant splicing events in these samples. Counting of the spliced and none spliced reads and parameter optimisation worked perfectly. However, during the model fitting step I run into the following error:

Wed Feb  1 12:08:32 2023: Fit step for: 'psi5'.
Wed Feb  1 12:08:34 2023: Running fit with correction method: PCA
Wed Feb  1 12:08:36 2023: Computing PCA ...
Wed Feb  1 12:10:39 2023: Fitting rho ...
Thu Feb  2 05:31:11 2023: rho fit:
             Min.   1st Qu.    Median      Mean   3rd Qu.      Max.
        0.0000000 0.0000001 0.0000001 0.0047353 0.0000001 0.5999270
Thu Feb  2 05:36:43 2023: Compute p values for: 'psi5'.
Thu Feb  2 22:18:56 2023: Adjust p values for: 'psi5'.
Thu Feb  2 22:19:42 2023: Compute Z scores for: 'psi5'.
Error in check_returned_array(ans, expected_dim, "extract_array", class(x)) :
  The "extract_array" method for HDF5ArraySeed objects returned an array
  with incorrect dimensions. Please contact the author of the
  HDF5ArraySeed class (defined in the HDF5Array package) about this and
  point him/her to the man page for extract_array() in the DelayedArray
  package (?extract_array).
Calls: FRASER ... extract_array -> extract_array -> check_returned_array
Execution halted

As the error occurs in the middle of the fitting step, I believe that it is not directly related to my input files and thus I was wondering whether you had an idea what might be causing the error?

This is the script I use for the model fit:

############ NVQ_Run093-NVQ_Run642 Untreated Samples FRASER Analysis: Fit Model and Call Outliers ############

#### Load Required Libraries ####
library(FRASER)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
library(org.Hs.eg.db)

#### Define Global Parameters ####
WORK_DIR <- '/data/gent/shared/001/gvo00101/RNA-Seq_in_Diagnostics/Aberrant_Splicing/FRASER/work/'
OUT_DIR <- '/data/gent/shared/001/gvo00101/RNA-Seq_in_Diagnostics/Aberrant_Splicing/FRASER/output/'
register(MulticoreParam(workers = 10))

#### Read in FRASER Object ####
fds <- readRDS(paste(OUT_DIR,'NVQ_Run093-NVQ_Run642_FRASER_Untreated_Filtered_Optim_psi5_psi3_and_theta.rds', sep = ''))

#### Detection of Aberrant Splicing Events ####
# Fit the Splicing Model 
fds <- FRASER(fds, q=c(psi5=5, psi3=5, theta=5))

# Annotate Introns
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
orgDb <- org.Hs.eg.db

fds <- annotateRangesWithTxDb(fds, txdb=txdb, orgDb=orgDb)

# Call Splicing Outliers
res <- results(fds, zScoreCutoff=NA, padjCutoff = 0.999999, deltaPsiCutoff=0.0000001)

#### Save Objects ####
# Processed FRASER Object
FDS_File <- paste(OUT_DIR,'NVQ_Run093-NVQ_Run642_FRASER_Untreated_Processed.rds', sep = '')

saveRDS(fds, FDS_File)

# Results Object
RES_File <- paste(OUT_DIR,'NVQ_Run093-NVQ_Run642_FRASER_Untreated_Results.rds', sep = '')

saveRDS(res, RES_File)

#### Generate Results File for Each Sample ####
Samples <- unique(res$sampleID)

for (s in Samples){
  Results_Samples <- res[res$sampleID == s]
  
  write.xlsx(Results_Samples, paste(OUT_DIR, 'Untreated/','NVQ_Run093-NVQ_Run642_FRASER_Untreated_Results_deltaPsi_0_padj_1_', s, '.xlsx', sep = ''))
} 

Kind regards,
Laurenz De Cock