The step for DGE and differential isoform usage is interrupted.

[Choeun-fran]

Thank you for your effort for this impressive tool. I'm running this tool with multiple 10x genomics single cell data. However, the process is suddenly halted before the step for DGE and differential isoform usage.

This is my full script to run SCALPEL:

#!/bin/bash 

# Set directory to save temporal files 
# (if these are not set, then nextflow make those files in home directory - memory issue)
export NXF_HOME=/node200data/dlchdms0109/.nextflow
export NXF_JAVA_TMPDIR=/node200data/dlchdms0109/tmp
export NXF_CONDA_CACHEDIR=/node200data/dlchdms0109/conda

# clean cache first
rm -rf /node200data/dlchdms0109/.nextflow/assets/

/home/dlchdms0109/.local/bin/nextflow run /home/dlchdms0109/Isoform_imputation/sr_tx_analysis/scalpel/SCALPEL/main.nf \
  -work-dir /node200data/dlchdms0109/run_SCALPEL/work \
  --sequencing chromium \
  --samplesheet /node200data/dlchdms0109/run_SCALPEL/samplesheet_for_SCALPEL.csv \
  --transcriptome /home/dlchdms0109/Isoform_imputation/Reference_genome/gencode.v49.transcripts.fa.gz \
  --gtf /home/dlchdms0109/Isoform_imputation/Reference_genome/gencode.v49.primary_assembly.annotation.gtf \
  --ipdb /home/dlchdms0109/Isoform_imputation/Reference_genome/ipdb/hg38_ipriming_sites.bed.tar.gz \
  --cpus 30 \
  -with-conda /home/dlchdms0109/Isoform_imputation/sr_tx_analysis/scalpel/SCALPEL/requirements.yml \
  --outputDir /node200data/dlchdms0109/run_SCALPEL/scalpel_output

I describe my sample sheet file below:

PO13,/node200data/dhkim1674/SC_SR_Samsung_Alzheimer/fastq/PO_/Samsung_hBrain_PO13_NC_GEX/Samsung_hBrain_PO13_NC_GEX_S1_R1_001.fastq.gz,/node200data/dhkim1674/SC_SR_Samsung_Alzheimer/fastq/PO_/Samsung_hBrain_PO13_NC_GEX/Samsung_hBrain_PO13_NC_GEX_S1_R2_001.fastq.gz,/node200data/dhkim1674/SC_SR_Samsung_Alzheimer/cellranger_output/PO13/

When I checked the output directory, iDGE_seurat.RDS and DIU_table.csv are not made.

This is the error state when I run the tool:

Error in `vec_rep_each()`:
  ! `times` can't be missing. Location 1 is missing.
  Backtrace:
      ▆
   1. ├─... %>% ...
   2. ├─tidyr::pivot_wider(...)
   3. ├─tidyr:::pivot_wider.data.frame(...)
   4. │ └─tidyr::pivot_wider_spec(...)
   5. │   └─vctrs::vec_rep_each(fill, nrow * ncol)
   6. └─rlang::abort(message = message, call = call)
  Warning message:
  In nrow * ncol : NAs produced by integer overflow
  Execution halted

I didn't understand why those processes are suddenly interrupted

[mschilli87]

How big is your dataset? I.e. how many gene isoforms and cells does your iDGE cover? My best guess is that the R script generating the files you are missing uses some inefficient approach of converting the wide matrix into a long table that will have # isoforms x # cells rows and that your dataset exceed the integer limit of R (thus the In nrow * ncol : NAs produced by integer overflow warning that might create the NA missing in vec_rep_each's times) argument. If that's the case, we might have to do that conversion in chunks for big dataset.

@Franzx7: Could you have a look at the code and check if my therory is plausible?

@Choeun-fran: Could you try sub-sampling your dataset to, say a few thousand cells, and repeating running the same script on the same machine to confirm it works with smaller inputs?

[Choeun-fran]

Thank you for your replying.
I need to run SCALPEL with 8 samples. There are about 60,000 cells and 30,000 genes. I will try to subsampling of my data. Additionally, I have one more question related to this problem. Is it possible that the problem is related to the R version? Because when I checked the packages that are needed to run this tool, some packages need later than 4.0. But my R version is 3.6. So if it's right, then if I update my R, then is it working correctly?

[mschilli87]

I am not sure if the integer limit was raised between R 3.x and 4.0. Either way, SCALPEL should be written in a way that scales with larger dataset. So knowing if it works for you using 6,000 cells (10 %) would be helpful.

---

PS: I just did a quick check. On my system (64-bit Linux, R 4.3.1), .Machine$integer.max is 2147483647 (*i.e. 2^23 / 2 - 1, so signed 32-bit integers). Thus 6e5 * 3e5 is larger than this limit, and even 6e4 * 3e5 would still not fit. So for your test, please try to also subsample the genes to 3,000 (10 %), as 6e4 * 3e4 is 1.8e9 and thus < .Machine$integer.max.

[mschilli87]

@Franzx7: I suspect this operation to be the culprit:

SCALPEL/src/merge_DGE_PRED.R

Lines 30 to 31 in 3c31ffc

	distinct(gene_transcript, bc, tr_counts) %>%
	pivot_wider(names_from=bc, values_from=tr_counts, values_fill = 0)

If n_distinct(gene_transcript, bc, tr_counts) > 2e31 - 1, pivot_wider(..., values_fill = 0) will eventually try to construct a vector longer than 2e31 - 1 full of 0s (the fill value) via vctrs::vec_rep_each(fill, nrow * ncol).
This in turn just calls ffi_vec_rep_each, passing the whole operation over from R to C, but eventually casting the times argument to int. This integer overflow results in an NA (and a warning), but as repeating a vector NA times is illdefined, vec_rep_each throws the error that @Choeun-fran reported.

So either we find a way around this pivot_wider call, or we should consider 'chunking' the long table into pieces < 1e31 rows, calling pivot_wider on each chunk separately and bind_rows the resulting wide tables back together.

But I'd need your input regarding what exactly is done here to come up with a chunking strategy that ensures the sub-tables are independent. I assume as long as we don't split in-between isoforms of the same gene (or even ensure that all instances of a single bc entry are within the same 'chunk'), this is probably quite indendent. But I'd love to hear your thoughts before digging further.

[Choeun-fran]

Thank you for your advice.
I will try again running SCALPEL with subsampling. I have few more questions about this process. I wil try use --subsample parameter. Is it correct way? Additional question is, if SCALPEL will be running with subsampled data, then should I run it repeatedly? Because doing only one time is quiet difficult to see that the subsampled data represent the original one.

[mschilli87]

@Choeun-fran: I meant to just run the failing script and modifying the input to contain a random subset of 10% of cells and 10 % genes. This will abviously not give you the results you want but know if this makes the error go away would help us to make a new SCALPEL version that 'just works' on your full dataset, using the command you tried initially.
If SCALPEL's own --subsample parameter can be used for this is probably better answered by @Franzx7.

[Choeun-fran]

Hi, thank you for your help to solve this problem.I tried to run SCALPEL with subsampled data (around 5,000 cells) as you mentioned. I only modified input data using the same script that previously failed.

Command error:
 
  Attaching package: ‘dplyr’
 
  The following objects are masked from ‘package:stats’:
 
      filter, lag
 
  The following objects are masked from ‘package:base’:
 
      intersect, setdiff, setequal, union
 
 
  Attaching package: ‘data.table’
 
  The following objects are masked from ‘package:dplyr’:
 
      between, first, last
 
  Warning message:
  In melt.data.table(dge) :
    id.vars and measure.vars are internally guessed when both are 'NULL'. All non-numeric/integer/logical type columns are considered id.vars, which in this case are columns [GENE, ...]. Consider providing at least one of 'id' or 'measure' vars in future.
  Error in `vec_rep_each()`:
  ! `times` can't be missing. Location 1 is missing.
  Backtrace:
      ▆
   1. ├─... %>% ...
   2. ├─tidyr::pivot_wider(...)
   3. ├─tidyr:::pivot_wider.data.frame(...)
   4. │ └─tidyr::pivot_wider_spec(...)
   5. │   └─vctrs::vec_rep_each(fill, nrow * ncol)
   6. └─rlang::abort(message = message, call = call)
  Warning message:
  In nrow * ncol : NAs produced by integer overflow
  Execution halted

Even I provided smaller samples, I got exactly same error with previous error. I hope this trial can help with the development of this impressive tool.

[mschilli87]

@Choeun-fran: How many features (isoforms) did you use for the subsampled datasets? If you used 60,000 features for 5,000 cells, your matrix should have 300,000,000 cells, which is way below the maximim signed 32 bit integer. So this is certainly helpful.

Could you try to modify

SCALPEL/src/merge_DGE_PRED.R

Line 30 in 3c31ffc

distinct(gene_transcript, bc, tr_counts) %>%

to

  distinct(gene_transcript, bc, tr_counts) %T>% print %>%

and add

  library(magrittr)

after

SCALPEL/src/merge_DGE_PRED.R

Line 5 in 3c31ffc

library(tidyr)

and re-run your smallest subsample that reproduces the error and share the output with us?

[Choeun-fran]

I'm using 5,000 cells and 35,000 features and I ran again SCALPEL with the modified merge_DGE_PRED.R. But I got an error again.
Below is the error statement:

 Attaching package: ‘magrittr’
 
  The following object is masked from ‘package:tidyr’:
 
      extract
 
  Warning message:
  In melt.data.table(dge) :
    id.vars and measure.vars are internally guessed when both are 'NULL'. All non-numeric/integer/logical type columns are considered id.vars, which in this case are columns [GENE, ...]. Consider providing at least one of 'id' or 'measure' vars in future.
              gene_transcript                 bc tr_counts
                       <char>             <char>     <num>
         1: ACADM***ACADM-203 AAACCAAAGCGATCTA-1     2.904
         2: ACADM***ACADM-205 AAACCAAAGCGATCTA-1     0.003
         3: ACADM***ACADM-208 AAACCAAAGCGATCTA-1     0.000
         4: ACADM***ACADM-209 AAACCAAAGCGATCTA-1     0.000
         5: ACADM***ACADM-212 AAACCAAAGCGATCTA-1     0.003
        ---
  47626625:   UBA1***UBA1-241 TGTGTTGAGCTGGCAT-1     0.390
  47626626:   UBA1***UBA1-251 TGTGTTGAGCTGGCAT-1     0.125
  47626627:   UBA1***UBA1-254 TGTGTTGAGCTGGCAT-1     0.000
  47626628:   UBA1***UBA1-256 TGTGTTGAGCTGGCAT-1     0.002
  47626629: ZFP92***ZFP92-201 TGTGTTGAGCTGGCAT-1     1.000
  Error in `vec_rep_each()`:
  ! `times` can't be missing. Location 1 is missing.
  Backtrace:
      ▆
   1. ├─... %>% ...
   2. ├─tidyr::pivot_wider(...)
   3. ├─tidyr:::pivot_wider.data.frame(...)
   4. │ └─tidyr::pivot_wider_spec(...)
   5. │   └─vctrs::vec_rep_each(fill, nrow * ncol)
   6. └─rlang::abort(message = message, call = call)
  Warning message:
  In nrow * ncol : NAs produced by integer overflow
  Execution halted

[mschilli87]

@Choeun-fran: Thank you!
Yes - I was expecting to still see the error. But now we have confirmation that the data.table has indeed <= 175,000,000 rows for your 5,000 times 35,000 matrix. I will have to look a bit deeper into the code to figure out where you might be hitting the integer limit then. I will report back.

@Franzx7: If you have any ideas, let me know.

[franzake]

Hi @Choeun-fran , Sorry for the late reply, and thank you for your valuable report.

Like @mschilli87 , my first suspicion was that this might be related to pivot_wider. However, since you’re seeing the same error even when using a subset of the cells, it seems likely the root cause lies elsewhere.

Could you please share a very small subset of the input files (just the minimum needed to reproduce the issue)? That would allow us to replicate the error on our side and investigate further...

Best,
F