About Error in "Example analysis - pbmc8k"

[Renown-TAL]

I try the "Using polyApipe" pipeline and get a error in "Further processing in R" as following:
################ code ################
polyApiper:::do_ensembl_organism(
out_path="human_ens100",
species="Homo sapiens",
version="100")

sne <- read.csv("/home/mjy/pertub_seq_APA/polyApipe_result/test_data/tsne/2_components/projection.csv")
cells_to_use <- tsne$Barcode

cell_name_func <- function(batch,cell) paste0(batch,cell,"-1")

polyApiper:::do_pipeline(
out_path ="/home/mjy/pertub_seq_APA/polyApipe_result",
counts_files =c("/home/mjy/pertub_seq_APA/polyApipe_result/pbmc8k_counts.tab.gz"),
batch_names =c(""),
peak_info_file="/home/mjy/pertub_seq_APA/polyApipe_result/pbmc8k_polyA_peaks.gff",
organism ="human_ens100",
cells_to_use =cells_to_use,
cell_name_func=cell_name_func
)
############## output info ####################
-- 1/4 load --
Loading /home/mjy/pertub_seq_APA/polyApipe_result/pbmc8k_counts.tab.gz
Loaded 118500 x 8381 matrix of counts
-- 2/4 assign --
'getOption("repos")' replaces Bioconductor standard repositories, see 'help("repositories", package = "BiocManager")' for details.
Replacement repositories:
CRAN: https://mirror-hk.koddos.net/CRAN/
-- 3/4 weitrices --
119278 peaks in 8381 cells
Kept 21594 peaks in 8381 cells
10442 genes
4519 genes with two or more peaks
Compute normalized, log transformed counts
Preclustering
15 preclusters
Compute APA shifts and peak proportions
Error in counts_shift(assay(peaks_se_final, "counts"), grouping, typecast = realize_RleArray) :
unused argument (typecast = realize_RleArray)
In addition: Warning message:
In .Seqinfo.mergexy(x, y) :
Each of the 2 combined objects has sequence levels not in the other:

• in 'x': GL000008.2, GL000221.1, GL000224.1, KI270706.1, KI270708.1, KI270719.1
• in 'y': CHR_HG107_PATCH, CHR_HG109_PATCH, CHR_HG126_PATCH, CHR_HG1277_PATCH, CHR_HG1298_PATCH, CHR_HG1309_PATCH, CHR_HG1311_PATCH, CHR_HG1320_PATCH, CHR_HG1342_HG2282_PATCH, CHR_HG1362_PATCH, CHR_HG1384_PATCH, CHR_HG1395_PATCH, CHR_HG1398_PATCH, CHR_HG142_HG150_NOVEL_TEST, CHR_HG1445_PATCH, CHR_HG1485_PATCH, CHR_HG151_NOVEL_TEST, CHR_HG1524_PATCH, CHR_HG1531_PATCH, CHR_HG1535_PATCH, CHR_HG1651_PATCH, CHR_HG1708_PATCH, CHR_HG1815_PATCH, CHR_HG1832_PATCH, CHR_HG1_PATCH, CHR_HG2002_PATCH, CHR_HG2021_PATCH, CHR_HG2022_PATCH, CHR_HG2023_PATCH, CHR_HG2030_PATCH, CHR_HG2046_PATCH, CHR_HG2047_PATCH, CHR_HG2057_PATCH, CHR_HG2058_PATCH, CHR_HG2060_PATCH, CHR_HG2062_PATCH, CHR_HG2063_PATCH, CHR_HG2066_PATCH, CHR_HG2067_PATCH, CHR_HG2072_PATCH, CHR_HG2087_PATCH, CHR_HG2088_PATCH, CHR_HG2095_PATCH, CHR_HG2104_PATCH, CHR_HG2111_PATCH, CHR_HG2114_ [... truncated]
  ##############################

Could you tell me how to solve this problam? Thanks in advanced.

[Renown-TAL]

Another example data having other error:
############## code ################
do_ensembl_organism(
out_path="mouse_ens100",
species="Mus musculus",
version="100")

tsne <- read.csv("/home/mjy/pertub_seq_APA/polyApipe_result/SRR5259422_demo_counts.tab.gz", sep = "\t")
cells_to_use <- tsne[!duplicated(tsne[2]),]$cell
cells_to_use <- paste0(cells_to_use,"-1")

cell_name_func <- function(batch,cell) paste0(batch,cell,"-1")

polyApiper:::do_pipeline(
out_path ="/home/mjy/pertub_seq_APA/polyApipe_result",
counts_files =c("/home/mjy/pertub_seq_APA/polyApipe_result/SRR5259422_demo_counts.tab.gz"),
batch_names =c(""),
peak_info_file="/home/mjy/pertub_seq_APA/polyApipe_result/SRR5259422_demo_polyA_peaks.gff",
organism ="mouse_ens100",
cells_to_use =cells_to_use,
cell_name_func=cell_name_func
)

############ output info#########
-- 1/4 load --
Loading /home/mjy/pertub_seq_APA/polyApipe_result/SRR5259422_demo_counts.tab.gz
Loaded 116 x 2839 matrix of counts
-- 2/4 assign --
'getOption("repos")' replaces Bioconductor standard repositories, see 'help("repositories", package = "BiocManager")' for details.
Replacement repositories:
CRAN: https://mirror-hk.koddos.net/CRAN/
-- 3/4 weitrices --
117 peaks in 2839 cells
Kept 5 peaks in 2839 cells
4 genes
1 genes with two or more peaks
Compute normalized, log transformed counts
Preclustering
Error in .local(x, ...) : size factors should be positive

[PoroProcessor]

@Renown-TAL, I also ran into the same issue at this step.

Preprocessing

I have been using the demo data for the processing, here are the codes before I ran into the error:

# In Linux
demo_dir="${main_dir}/demo_group"
polyApipe/polyApipe.py \
-i /~/scRNA_polyA/polyApipe/data/demo/ \
-o ${demo_dir}

# In R
> counts_file_dir
[1] "/~/scRNA_polyA/demo_group_counts"
list.files(counts_file_dir)
[1] "SRR5259354_demo.tab.gz" "SRR5259422_demo.tab.gz"

First Error and temporary solution:

I then ran the source code of do_pipeline() line by line since I encountered this error. I found the error to be coming from
result_gene <- se_logcounts(result_gene, do_logNormCounts=do_logNormCounts, do_computeSumFactors=do_computeSumFactors)

Error message:

Preclustering
Error: BiocParallel errors
2 remote errors, element index: 1, 2
0 unevaluated and other errors
first remote error:
Error in .local(x, ...): size factors should be positive

Cause of Error?

Which is coming from Line 319 in pipeline.R, which is using a function defined in the same file se_logcounts, and I believe Line 178 in pipeline.R is the cause of such error, which is using scran::quickCluster() (source code: https://rdrr.io/bioc/scran/src/R/quickCluster.R#sym-.quick_cluster). And the "size factors should be positive" appears to be coming from:

sf <- librarySizeFactors(x, subset_row=subset.row)
y <- normalizeCounts(x, size_factors=sf, subset_row=subset.row)

Temporary Solution

After making sure that there were no negative values in sf; there were just some zeros, I tried to surpass this error by adding a pseudocount to make all size factors positive by adding this between the two commands above:

sf <- librarySizeFactors(x, subset_row=subset.row)

epsilon <- 0.000001
cat("Adjusting size factors by adding", epsilon, "\n")
sf <- sf + epsilon

y <- normalizeCounts(x, size_factors=sf, subset_row=subset.row)
cat("Done normalizing counts with adjusted size factors.", "\n")

Second Error and no solution yet...

Although I can surpass this error, I encountered another error from the same se_logcounts line:

> result_gene <- se_logcounts(result_gene, do_logNormCounts=do_logNormCounts, do_computeSumFactors=do_computeSumFactors)
Preclustering
Adjusting size factors by adding 1e-06 
[1] "Done normalizing counts with adjusted size factors."

Adjusting size factors by adding 1e-06 
[1] "Done normalizing counts with adjusted size factors."

Error: BiocParallel errors
  2 remote errors, element index: 1, 2
  0 unevaluated and other errors
  first remote error:
Error in assign(".effective_grid", effective_grid, envir = envir): invalid 'envir' argument

Error traceback

I am attaching the error traceback here (error.txt), generated immediately after I encounter this error:

sink("error3.txt", append = TRUE)
print("Traceback:")
traceback()
sink()

Final notes

@pfh I deeply appreciate the tool and the Python portion worked seamlessly for me, which was really nice! I noticed that the same error was reported in July by @cihaterdogan #3 (Error in .local(x, ...) : size factors should be positive), but no solution was provided specifically for this error. I am not sure if my second error (...invalid 'envir' argument) was due to the temporary solution I used for my first error. Perhaps there could be a more systematic fix to this issue?

Again, thank you so much for making this tool. I also appreciate anyone who may have a solution for this.

[pfh]

The second error, with the included demo data, is probably a bit small for the R part of the pipeline to work properly. Not great as a demo, sorry!

Negative size factors from computeSumFactors is a thing that happens. When I wrote this, I thought using computeSumFactors was the right thing to do. I am no longer so convinced. There's an option to disable it.

Zero size factors is a new one. Hmm.

There really should be a working end-to-end example.

Really, the whole code needs more testing. I'm not sure if that is going to happen.

If the python part written by Sarah Williams is giving you useful results at least, I'm pleased to hear that.

[shimi-ly]

Hi，unused argument (typecast = realize_RleArray)，I'm also experiencing this problem, it shows that this does not exist for this parameter, have you solved this problem yet?

[pfh]

@shimi-ly, I have to apologize, this package is not really being actively maintained at the moment.

I think this might be a different error from the other ones. I'm guessing that is from a call to a weitrix function such as counts_shift or counts_proportion, but recent versions should have this argument. What version of weitrix do you have?