STAGER error: subscript out of bounds

[jma1991]

Hello,

We recently integrated the stageR package into the nf-core/rnasplice workflow. However, a user has reported a technical issue that we're struggling to resolve. When the user sets the alpha parameter in the stageWiseAdjustment function to the standard value of 0.05, they receive a 'subscript out of bounds' error.

Error in `[<-`(`*tmp*`, idCon, 1, value = unlist(txLevelAdjustments)) : 
  subscript out of bounds

Interestingly, adjusting alpha to 1 triggers a different error.

Error in FUN(X[[i]], ...) : 
  Some genes have only one transcript; this is incompatible with DTU correction. Remove these transcripts.

We're seeking guidance on how to address this issue. The user has shared their input files for reference, and the relevant code is located in the run_stager.R script, which you can find in the bin folder of our GitHub repository.

Thanks,
James

[koenvandenberge]

Hi @jma1991

Thank you for reporting. I will take some time to look into this in the coming days.

[SethMagnusJarvis]

I am also having this error. It's when using DRIMSeq though. Someone previously said that it was due to NA values but it would be fixed in the next patch. Is there a quick fix e.g. an easy way to remove NA values so it doesn't throw errors?

[RPdavies]

This seems to occur when there are NA values in the pScreen vector. Filtering these genes out (and their corresponding transcripts) seems to avoid this error. I can't imagine doing so having any unintended consequences but it would be great if one of the devs could confirm that. Here's the augmented workflow integrating DRIMseq and stageR:

# extract results, test_x is the dmFit() output
gene_res <- results(test_x, level = "gene")      # gene‑level results
tx_res   <- results(test_x, level = "feature")   # transcript‑level results

# filter out genes with NA p‑values
gene_res_filt <- gene_res[!is.na(gene_res$pvalue), ]

pScreen <- gene_res_filt$pvalue
names(pScreen) <- gene_res_filt$gene_id

# grab the gene IDs you want to keep
keep_genes <- gene_res_filt$gene_id

# subset the transcript table to just those genes
tx_res_filt <- tx_res[ tx_res$gene_id %in% keep_genes, ]

# now gene_res_filt and tx_res_filt have no NAs at the gene level,
# and tx_res_filt only contains transcripts from genes in gene_res_filt

## Assign transcript-level pvalues to the confirmation stage
pConfirmation <- matrix(tx_res_filt$pvalue, ncol = 1)
rownames(pConfirmation) <- tx_res_filt$feature_id

## Create the gene-transcript mapping
tx2gene <- tx_res_filt[, c("feature_id", "gene_id")]

## Create the stageRTx object and perform the stage-wise analysis
stageRObj <- stageRTx(pScreen = pScreen, pConfirmation = pConfirmation,
                      pScreenAdjusted = FALSE, tx2gene = tx2gene)

stageRObj <- stageWiseAdjustment(object = stageRObj, method = "dtu",
                                 alpha = 0.05, allowNA=TRUE)

padj <- getAdjustedPValues(stageRObj, order = TRUE,
                           onlySignificantGenes = FALSE)