Pipeline crash in compute_prob.R when transcriptomic_distance is positive (10x 3' data)

[Chris-Kato]

Hi,

I encountered a crash in the SCALPEL pipeline at the step:

isoform_quantification:probability_distribution

The error is:

Error in seq.default(0, read_tab$transcriptomic_distance[1], -BINS) :
wrong sign in 'by' argument

This occurs in compute_prob.R at:

part_neg = c(seq(0,read_tab$transcriptomic_distance[1],-BINS), ...)

From debugging, the issue arises because the pipeline assumes that
transcriptomic_distance values are negative (i.e. upstream of the 3' end).

However, in my dataset (10x Chromium 3' snRNA-seq), many distances are positive.

Example from all_unique_reads.txt:

chr21 44335310 44335372 + 478 ...
chr21 44335399 44335489 + 361 ...
chr21 44335490 44335553 + 297 ...

Here the transcript coordinates are:

transcript: 44335251–44335851 (+ strand)

So the read lies upstream of the 3' end, but the computed dist_END is positive.

Because of this, the following call fails:

seq(0, positive_value, -BINS)

which leads to the crash.

As a temporary workaround, I inverted the sign of dist_END in compute_prob.R:

reads = reads %>%
dplyr::filter(gene_name %in% gene.tokeep) %>%
mutate(dist_END = -as.numeric(dist_END))

After this change, the pipeline proceeds normally.

My questions are:

1. Should transcriptomic_distance be expected to be negative upstream of the 3' end?
2. Is the current sign convention in mapping_filtering.R intended?
3. Should compute_prob.R handle both positive and negative distances more robustly (e.g. using min() instead of [1])?

Thank you for developing SCALPEL.

Best regards

[mschilli87]

I think that this is not a bug in SCALPEL. SCALPEL is annotation based. As such, the pipeline assumes that the transcriptome annotation you provide contains all transcripts that can occur in the sample. Thus, I'd suggest you to extend your annotation based on your data first and then to re-run SCALPEL with your extended reference. I successfully used such a strategy some years back in the following paper: Diag & Schilling et al., Dev. Cell (2018).

The only biologically valid reason for such 'downstream' reads to occur that I can think of, specifically in the context of nuclei-seequencing, would be primary transcripts pre-cleave and poly-adenylation that also contain long enough stretches of 'A's to capture them with oligo-(dT)s. But those should definitely be very rare and I'd argue it would be acceptable for SCALPEL to neglect them.

That said, SCALPEL should certainly handle this situation more gracefully and give an error message that explains what is going wrong and how to potentially fix it.

@Franzx7: Do you think it would be feasible to extent SCALPEL to handle read exceeding beyon the annotated 3' end? If not, what do you think about providing an (optional) extra step to extend the annotation (pre-salmon)?

[Chris-Kato]

Hi, thank you very much for the quick response and for the helpful explanation.

I understand your point that reads extending beyond the annotated 3′ end should generally be rare and could reflect primary transcripts prior to cleavage and polyadenylation, especially in the context of nuclei sequencing.

However, in my dataset the situation appears somewhat different. In the all_unique_reads.txt produced by SCALPEL, the distribution of dist_END values is not dominated by a small number of downstream reads. Instead, most distances are positive:

positive distances: 99.1%
negative distances: 0.8%
zero: 0.1%

So in this case the positive values do not seem to represent a small set of reads extending beyond the annotated 3′ end. Rather, the majority of reads fall into this category.

For example, one of the entries is:

transcript: chr21:44335251–44335851 (+ strand)
read: chr21:44335310–44335372
dist_END: 478

Here the read is upstream of the annotated 3′ end, yet the computed dist_END is positive.

Because of this, the call in compute_prob.R

seq(0, read_tab$transcriptomic_distance[1], -BINS)

fails when the first value is positive.

From this behaviour it seems possible that either:

the sign convention of dist_END differs from what compute_prob.R expects, or

the transcript end used for distance calculation is systematically offset from the effective capture positions in my reference/data.

Thank you again for taking the time to look into this. I appreciate the suggestion regarding annotation extension and would be interested to hear whether this interpretation of the distance distribution also seems plausible from your perspective.

[mschilli87]

Did you map those genes to the genome? An interesting quick check would be to pool all reads (across nuclei and samples) for all single-isoform(!) genes and plot the genomic coverage centered around the annotated 3' end ('metagene plot').
Also, is your dataset publicly available? How well annotated is the organism your are dealing with? Is it possible that your samples are dominated by unannotated isoforms using more distal PAS sites? Or do you suspect an error in SCALPEL's processing? Can you reproduce this behaviour with any other (publicly available!) dataset?

[Chris-Kato]

Hi, thank you for the suggestions — these are very helpful points.

Genome mapping:
Yes, the reads were mapped to the genome prior to running SCALPEL.
Looking at the genomic coordinates of the reads in IGV, most reads appear to fall upstream of the annotated 3′ end, as expected for 3′ capture data.

Annotation:
The analysis uses a custom transcriptome annotation derived from our own analysis combined with a reference genome annotation.

Unannotated isoforms vs processing issue:
It is certainly possible that some reads correspond to unannotated distal PAS usage, especially in snRNA-seq where pre-mRNA may still be present.

After investigating the intermediate files and the code more closely, I think the issue may be related to an implementation assumption in compute_prob.R.

In the raw all_unique_reads.txt, both positive and negative values of dist_END are present. However, after the filtering steps inside compute_prob.R (gene filtering, distinct(), and grouping), the resulting distribution of read_tab$transcriptomic_distance can sometimes become one-sided.

The current code implicitly assumes that both negative and positive distances exist:

part_neg = c(seq(0, read_tab$transcriptomic_distance[1], -BINS),
read_tab$transcriptomic_distance[1]) %>% unique() %>% rev()

part_pos = c(seq(0, max(read_tab$transcriptomic_distance), BINS),
max(read_tab$transcriptomic_distance)) %>% unique()

If the filtered distribution happens to be only positive or only negative, this results in a sign mismatch in seq(), for example:

seq(0, positive_value, -BINS)

or

seq(0, negative_value, +BINS)

which produces the error:

Error in seq.default(...): wrong sign in 'by' argument

Importantly, this does not necessarily mean that the raw reads are one-sided, but rather that the distribution becomes one-sided after the filtering and collapsing steps within compute_prob.R.

To address this, I modified only the interval construction step so that it checks whether negative or positive distances exist before generating the corresponding bins. This avoids the crash while keeping the definition and sign of dist_END unchanged.

It might therefore be useful for SCALPEL to make this part of the code slightly more robust to one-sided transcriptomic_distance distributions.

[mschilli87]

@Chris-Kato: Thank you for your detailed follow-up. This does indeed sound like an oversight on our side. Could you by any chance provide a PR and maybe even a small test case that reproduces this error in the current code?