Error in xj[i] : invalid subscript type 'list' when running jDUreport

[bertapozzi]

Hi all,

I'm using ASpli (mosly copied the example pipeline and replaced with my files) and encountered the following error:

jdur <- jDUreport(asd, contrast=c(-1,1)) Running junctionsPJU test Running junctionsPIR test Error in xj[i] : invalid subscript type 'list'

This doesn't happen when I run it with the example data, and cannot find a way around it

[estepi]

Hi Berta, thanks for your post! Can you paste your complete script? I need more info to decipher from where the error comes...thanks!!!!

[bertapozzi]

Hi estepi, thanks fro the reply.

I tried with different samples and did not encounter that error. However, exportIntegratedSignals did not work. This is what I get:

BAMFiles [1] "g1.bam" "g2.bam" "n1.bam" "n2.bam" targets <- data.frame(row.names = paste0('Sample',c(1:4)),

• bam = BAMFiles[1:4],
• f1 = c( 'gfp','gfp','ns5','ns5'),
• stringsAsFactors = FALSE) mBAMs <- data.frame( bam = sub("[12]","",targets$bam[c(1,3)]),
• condition = c("gfp","ns5")) getConditions( targets ) [1] "gfp" "ns5" targets bam f1 Sample1 g1.bam gfp Sample2 g2.bam gfp Sample3 n1.bam ns5 Sample4 n2.bam ns5 mBAMs bam condition 1 g.bam gfp 2 n.bam ns5 gbcounts <- gbCounts(features=features, targets=targets,
• minReadLength = 100, maxISize = 50000) Summarizing Sample1 ETA: 4 min Summarizing Sample2 ETA: 3 min Summarizing Sample3 ETA: 1 min Summarizing Sample4 ETA: 0 min gbcounts Object of class ASpliCounts Gene counts: 61544 genes analysed. Access using countsg(object) Gene RD: 61544 genes analysed. Access using rdsg(object) Bin counts: 1044348 bins analysed. Access using countsb(object) Bin RD: 1044348 bins analysed. Access using rdsb(object) Junction counts: 82803 junctions analysed. Access using countsj(object) asd <- jCounts(counts=gbcounts, features=features, minReadLength=100) Junctions PJU completed Junctions PIR completed Junctions IR PIR completed Junctions AltSS PSI completed Junctions ES PSI completed asd Object of class ASpliAS IR PIR: 398429 intron bins analysed. 8206 intron bins passed the filters. Access using irPIR(object) ES PSI: 233035 exon bins analysed. 1687 exon bins passed the filters. Access using esPSI(object) AltSS PSI: 35499 exon bins analysed. 1016 exon bins passed the filters. Access using altPSI(object) Junctions PIR: 10873 junctions analysed. Access using junctionsPIR(object) Junctions PJU: 58926 junctions analysed. Access using junctionsPJU(object) gb <- gbDUreport(gbcounts, contrast = c(-1,1)) Filtering genes: Filtering by reads. Filtering using gfp,ns5 conditions Filtering any condition with mean minimum value 10 Filtering by read density. Filtering using gfp,ns5 conditions Filtering any condition with mean minimum value 0.05 Filtering genes done Genes differential expression: Contrast:-1gfp 1ns5
  Running GLM LRT Genes differential expression done Genes DE completed Bins DE completed gb Object of class ASpliDU Gene DE: 4119 genes analysed. Access using genesDE(object) Bins DU: 5236 bins analysed. Access using binsDU(object) jdur <- jDUreport(asd, contrast=c(-1,1)) Running junctionsPJU test Running junctionsPIR test Running irPIR test Running esPSI test Running altPSI test jdur Object of class ASpliJDU Locale Junctions: 1934 junctions. Access using localej(object) Locale Clusters: 793 clusters Access using localec(object) Anchor Junctions: 315 junctions. Access using anchorj(object) Anchor Clusters: 315 clusters Access using anchorc(object) Intron Retention Junctions: 193 junctions. Access using jir(object) Exon Skipping Junctions: 81 junctions. Access using jes(object) Alternative Splicing Junctions: 212 junctions. Access using jalt(object) sr <- splicingReport(gb, jdur, counts=gbcounts) is <- integrateSignals(sr,asd) exportSplicingReports( sr, output.dir="resultsns5") exportIntegratedSignals(is,sr=sr,
• output.dir = "aspliNS5",
• counts=gbcounts,features=features,asd=asd,
• mergedBams = mBAMs) Generating graphs... | | 0%[E::hts_open_format] Failed to open file "g.bam" : No such file or directory samtools depth: Cannot open input file "g.bam": No such file or directory |============================ | 25%[E::hts_open_format] Failed to open file "g.bam" : No such file or directory samtools depth: Cannot open input file "g.bam": No such file or directory |======================================================== | 50%[E::hts_open_format] Failed to open file "g.bam" : No such file or directory samtools depth: Cannot open input file "g.bam": No such file or directory |===================================================================================== | 75%[E::hts_open_format] Failed to open file "g.bam" : No such file or directory samtools depth: Cannot open input file "g.bam": No such file or directory |=================================================================================================================| 100%[E::hts_open_format] Failed to open file "g.bam" : No such file or directory samtools depth: Cannot open input file "g.bam": No such file or directory

[estepi]

Hi Berta, Cool! Happy your error has gone

Can you check this line: mBAMs bam condition 1 g.bam gfp 2 n.bam ns5

Did you merge bams files into a bigger one (g.bam and n.bam)? You can try just using "single bam files"... g1.bam and n1.bam

thanks!

[bertapozzi]

Hi estepi, thanks for your reply. I now tried without merging BAMs and this is what I get, maybe .bai are missing?

mBAMs bam condition 1 g1.bam gfp 2 n1.bam ns5 exportIntegratedSignals(is,sr=sr,

• output.dir = "aspliNS5",
• counts=gbcounts,features=features,asd=asd,
• mergedBams = mBAMs) Generating graphs... | | 0%[E::idx_find_and_load] Could not retrieve index file for 'g1.bam' samtools depth: cannot load index for "g1.bam" |============================ | 25%[E::idx_find_and_load] Could not retrieve index file for 'g1.bam' samtools depth: cannot load index for "g1.bam" |======================================================== | 50%[E::idx_find_and_load] Could not retrieve index file for 'g1.bam' samtools depth: cannot load index for "g1.bam" |===================================================================================== | 75%[E::idx_find_and_load] Could not retrieve index file for 'g1.bam' samtools depth: cannot load index for "g1.bam" |=================================================================================================================| 100%[E::idx_find_and_load] Could not retrieve index file for 'g1.bam' samtools depth: cannot load index for "g1.bam"

[estepi]

Thanks! Yes please, you should index bam files before and save them in the same folder

[reck999]

I also got the same initial error. Here is my code, session info, and traceback. I've been able to get Aspli to work with other bams with this same code. Any advice? Thank you for a great package!

RNA-seq BAM files

flsall <- dir(getwd(),".bam") flsall<-paste0(getwd(),'/',flsall) names(flsall)<-gsub('.bam','',dir(getwd(),".bam")) BAMFiles <- flsall

RNA-seq GTF file

flsgtf <- dir(getwd(),".gtf") genomeTxDb <- makeTxDbFromGFF(file = flsgtf) features <- binGenome(genomeTxDb)

Combine BAM and target files

targets <- data.frame(row.names = paste0('Sample',c(1:9)), bam = BAMFiles[1:9], f1 = c( 'treatment','control','treatment','treatment','treatment','treatment','control','control','control'), stringsAsFactors = FALSE) mBAMs <- data.frame( bam = sub("_[012]","",targets$bam[c(1,4)]), condition = c("control","treatment"))

Read counting

gbcounts <- gbCounts(features=features, targets=targets, minReadLength = 295, maxISize = 50000, libType="PE")

Junction-based de-novo counting

asd <- jCounts(counts=gbcounts, features=features, minReadLength=295, libType="PE")

Differential gene expression

gb <- gbDUreport(gbcounts, contrast = c(-1,1))

Differential junction usage

jdur <- jDUreport(asd, contrast=c(-1,1), strongFilter=TRUE, useSubset = FALSE)

library(ASpli) library(GenomicFeatures)

Differential junction usage

jdur <- jDUreport(asd, contrast=c(-1,1), strongFilter=TRUE, useSubset = FALSE) Running junctionsPJU test Running junctionsPIR test Error in xj[i] : invalid subscript type 'list' sessionInfo() R version 4.3.3 (2024-02-29 ucrt) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows 10 x64 (build 19045)

Matrix products: default

locale: [1] LC_COLLATE=English_United States.utf8 LC_CTYPE=English_United States.utf8 LC_MONETARY=English_United States.utf8 [4] LC_NUMERIC=C LC_TIME=English_United States.utf8

time zone: America/Los_Angeles tzcode source: internal

attached base packages: [1] stats4 parallel stats graphics grDevices utils datasets methods base

other attached packages: [1] BiocManager_1.30.22 GenomicFeatures_1.54.4 GenomicRanges_1.54.1 GenomeInfoDb_1.38.8 ASpli_2.12.0
[6] AnnotationDbi_1.64.1 IRanges_2.36.0 S4Vectors_0.40.2 Biobase_2.60.0 BiocGenerics_0.48.1
[11] edgeR_4.0.16 limma_3.58.1

loaded via a namespace (and not attached): [1] RColorBrewer_1.1-3 rstudioapi_0.16.0 magrittr_2.0.3 rmarkdown_2.26
[5] BiocIO_1.12.0 zlibbioc_1.48.2 vctrs_0.6.5 memoise_2.0.1
[9] Rsamtools_2.18.0 RCurl_1.98-1.14 base64enc_0.1-3 htmltools_0.5.8
[13] S4Arrays_1.2.1 progress_1.2.3 curl_5.2.1 SparseArray_1.2.4
[17] Formula_1.2-5 htmlwidgets_1.6.4 plyr_1.8.9 Gviz_1.46.1
[21] cachem_1.0.8 GenomicAlignments_1.38.2 igraph_2.0.3 lifecycle_1.0.4
[25] pkgconfig_2.0.3 Matrix_1.6-5 R6_2.5.1 fastmap_1.1.1
[29] GenomeInfoDbData_1.2.11 MatrixGenerics_1.14.0 digest_0.6.35 colorspace_2.1-0
[33] Hmisc_5.1-2 RSQLite_2.3.5 filelock_1.0.3 fansi_1.0.6
[37] httr_1.4.7 abind_1.4-5 compiler_4.3.3 bit64_4.0.5
[41] htmlTable_2.4.2 backports_1.4.1 BiocParallel_1.34.2 DBI_1.2.2
[45] UpSetR_1.4.0 biomaRt_2.58.2 MASS_7.3-60.0.1 rappdirs_0.3.3
[49] DelayedArray_0.28.0 rjson_0.2.21 tools_4.3.3 foreign_0.8-86
[53] nnet_7.3-19 glue_1.7.0 restfulr_0.0.15 grid_4.3.3
[57] checkmate_2.3.1 cluster_2.1.6 generics_0.1.3 gtable_0.3.5
[61] BSgenome_1.70.2 tidyr_1.3.1 ensembldb_2.26.0 data.table_1.15.2
[65] hms_1.1.3 xml2_1.3.6 utf8_1.2.4 XVector_0.42.0
[69] pillar_1.9.0 stringr_1.5.1 splines_4.3.3 dplyr_1.1.4
[73] BiocFileCache_2.10.2 lattice_0.22-6 deldir_2.0-4 rtracklayer_1.62.0
[77] bit_4.0.5 biovizBase_1.50.0 tidyselect_1.2.1 locfit_1.5-9.9
[81] Biostrings_2.70.3 knitr_1.46 gridExtra_2.3 ProtGenerics_1.34.0
[85] SummarizedExperiment_1.32.0 xfun_0.43 statmod_1.5.0 matrixStats_1.2.0
[89] DT_0.33 stringi_1.8.3 lazyeval_0.2.2 yaml_2.3.8
[93] evaluate_0.23 codetools_0.2-19 interp_1.1-6 tibble_3.2.1
[97] cli_3.6.2 rpart_4.1.23 pbmcapply_1.5.1 munsell_0.5.1
[101] dichromat_2.0-0.1 Rcpp_1.0.12 dbplyr_2.5.0 png_0.1-8
[105] XML_3.99-0.16.1 ggplot2_3.5.1 blob_1.2.4 prettyunits_1.2.0
[109] jpeg_0.1-10 latticeExtra_0.6-30 AnnotationFilter_1.26.0 bitops_1.0-7
[113] VariantAnnotation_1.48.1 scales_1.3.0 purrr_1.0.2 crayon_1.5.2
[117] BiocStyle_2.30.0 rlang_1.1.3 KEGGREST_1.42.0

traceback() 6: [.data.frame(J3, reliables, ) 5: J3[reliables, ] 4: .makeJunctions(data, targets, start_J1, start_J2, start_J3, minAvgCounts, filterWithContrasted, contrast, strongFilter) 3: .junctionDUreportExt(asd, minAvgCounts, contrast, filterWithContrasted, runUniformityTest, mergedBams, maxPValForUniformityCheck, strongFilter, maxConditionsForDispersionEstimate, formula, coef, maxFDRForParticipation, useSubset) 2: jDUreport(asd, contrast = c(-1, 1), strongFilter = TRUE, useSubset = FALSE) 1: jDUreport(asd, contrast = c(-1, 1), strongFilter = TRUE, useSubset = FALSE)