umccr / RNAsum

Pipeline for generating RNAseq-based cancer patient reports
https://umccr.github.io/RNAsum/
Other
7 stars 4 forks source link

exprGroupsStats_geneWise.R: Subscript out of bounds #122

Closed pdiakumis closed 11 months ago

pdiakumis commented 1 year ago

During local debugging. Might have something to do with dependencies. Using SBJ04391_L2301279.

Quitting from lines 1177-1234 (rnasum.Rmd) 
Error in `dplyr::mutate()`:
ℹ In argument: `rnasum_rmd = RNAsum::rnasum_rmd(...)`.
ℹ In row 1.
Caused by error in `data.cum[, colnames(data)[targets$Target %in% group]]`:
! subscript out of bounds
Backtrace:
 1. RNAsum::exprTable(...)
 2. RNAsum::exprGroupsStats_geneWise(data, targets)
      at rnasum/R/exprTable.R:55:4
 3. base::cbind(...)
      at rnasum/R/exprGroupsStats_geneWise.R:97:6

The data.cum matrix generated in exprGroupsStats_geneWise only has rownames, so the colnames() doesn't find anything. data.cum is specified as:

##### Calculate cumulative sums and perform range standardization between 0 and 1
data.cum <- t(apply(data, 1, cumsum_ordered))

Will see if I can push a fix. Maybe the behaviour of apply() or t() has changed. Wonder how this worked before.


Debugging:

Browse[3]> class(data)
[1] "matrix" "array" 
Browse[3]> dim(data)
[1] 19497   371
Browse[3]> colnames(data) |> tail()
[1] "PAAD..UMCCR..Sample_36" "PAAD..UMCCR..Sample_37" "PAAD..UMCCR..Sample_38" "PAAD..UMCCR..Sample_39" "PAAD..UMCCR..Sample_40"
[6] "SBJ04391_L2301279"     
Browse[3]> data[1:5, 368:371]
         PAAD..UMCCR..Sample_38 PAAD..UMCCR..Sample_39 PAAD..UMCCR..Sample_40 SBJ04391_L2301279
TSPAN6                 3.918646               4.000431               4.897812          4.188811
DPM1                   4.499283               4.888805               4.314643          4.182681
SCYL3                  4.886067               4.882794               4.974602          3.989217
C1orf112               3.257823               3.230045               3.936439          2.699573
FGR                    2.727048               2.667036               2.352330          3.790203
Browse[3]> apply(data[1:10, 1:10], 1, cumsum_ordered) |> str()
 num [1:10, 1:10] 0.62459 0.00347 1 0.45969 0.05855 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:10] "PANCAN..TCGA..TCGA.OR.A5J1.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J2.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J3.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J5.01A.11R.A29S.07" ...
  ..$ : chr [1:10] "TSPAN6" "DPM1" "SCYL3" "C1orf112" ...
Browse[3]> apply(data, 1, cumsum_ordered) |> str()
 num [1:371, 1:19497] 0.437 0.184 0.769 0.429 0.239 ...
 - attr(*, "dimnames")=List of 2
  ..$ : NULL
  ..$ : chr [1:19497] "TSPAN6" "DPM1" "SCYL3" "C1orf112" ...

Comparing to a sample that works (SBJ04215), we see that there is a negative value that throws the cumsum_ordered off:

Browse[3]> rownames(data2) |> str()
 chr [1:19497] "TSPAN6" "DPM1" "SCYL3" "C1orf112" "FGR" "CFH" "FUCA2" "GCLC" "NFYA" "STPG1" "NIPAL3" "LAS1L" "ENPP4" "SEMA3F" "CFTR" ...
Browse[3]> rownames(data) |> str()
 chr [1:21251] "TSPAN6" "TNMD" "DPM1" "SCYL3" "C1orf112" "FGR" "CFH" "FUCA2" "GCLC" "NFYA" "STPG1" "NIPAL3" "LAS1L" "ENPP4" "SEMA3F" "CFTR" ...
Browse[3]> class(data)
[1] "matrix" "array" 
Browse[3]> class(data2)
[1] "matrix" "array" 
Browse[3]> attributes(data) |> str()
List of 2
 $ dim     : int [1:2] 21251 371
 $ dimnames:List of 2
  ..$ : chr [1:21251] "TSPAN6" "TNMD" "DPM1" "SCYL3" ...
  ..$ : chr [1:371] "PANCAN..TCGA..TCGA.OR.A5J1.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J2.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J3.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J5.01A.11R.A29S.07" ...
Browse[3]> attributes(data2) |> str()
List of 2
 $ dim     : int [1:2] 19497 371
 $ dimnames:List of 2
  ..$ : chr [1:19497] "TSPAN6" "DPM1" "SCYL3" "C1orf112" ...
  ..$ : chr [1:371] "PANCAN..TCGA..TCGA.OR.A5J1.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J2.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J3.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J5.01A.11R.A29S.07" ...
Browse[3]> apply(data2[1:19497, 1:371], 1, cumsum_ordered) |> str()
 num [1:371, 1:19497] 0.437 0.184 0.769 0.429 0.239 ...
 - attr(*, "dimnames")=List of 2
  ..$ : NULL
  ..$ : chr [1:19497] "TSPAN6" "DPM1" "SCYL3" "C1orf112" ...
Browse[3]> apply(data2[1:10, 1:20], 1, cumsum_ordered) |> str()
 num [1:20, 1:10] 0.185941 0.000727 0.326461 0.117036 0.020571 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:20] "PANCAN..TCGA..TCGA.OR.A5J1.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J2.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J3.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J5.01A.11R.A29S.07" ...
  ..$ : chr [1:10] "TSPAN6" "DPM1" "SCYL3" "C1orf112" ...
Browse[3]> apply(data2[1:1000, 1:20], 1, cumsum_ordered) |> str()
 num [1:20, 1:1000] 0.185941 0.000727 0.326461 0.117036 0.020571 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:20] "PANCAN..TCGA..TCGA.OR.A5J1.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J2.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J3.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J5.01A.11R.A29S.07" ...
  ..$ : chr [1:1000] "TSPAN6" "DPM1" "SCYL3" "C1orf112" ...
Browse[3]> apply(data2[1:2000, 1:20], 1, cumsum_ordered) |> str()
 num [1:20, 1:2000] 0.185941 0.000727 0.326461 0.117036 0.020571 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:20] "PANCAN..TCGA..TCGA.OR.A5J1.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J2.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J3.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J5.01A.11R.A29S.07" ...
  ..$ : chr [1:2000] "TSPAN6" "DPM1" "SCYL3" "C1orf112" ...
Browse[3]> apply(data2[1:18000, 1:20], 1, cumsum_ordered) |> str()
 num [1:20, 1:18000] 0.185941 0.000727 0.326461 0.117036 0.020571 ...
 - attr(*, "dimnames")=List of 2
  ..$ : NULL
  ..$ : chr [1:18000] "TSPAN6" "DPM1" "SCYL3" "C1orf112" ...
Browse[3]> apply(data2[1:16000, 1:20], 1, cumsum_ordered) |> str()
 num [1:20, 1:16000] 0.185941 0.000727 0.326461 0.117036 0.020571 ...
 - attr(*, "dimnames")=List of 2
  ..$ : NULL
  ..$ : chr [1:16000] "TSPAN6" "DPM1" "SCYL3" "C1orf112" ...
Browse[3]> apply(data2[1:13000, 1:20], 1, cumsum_ordered) |> str()
 num [1:20, 1:13000] 0.185941 0.000727 0.326461 0.117036 0.020571 ...
 - attr(*, "dimnames")=List of 2
  ..$ : NULL
  ..$ : chr [1:13000] "TSPAN6" "DPM1" "SCYL3" "C1orf112" ...
Browse[3]> apply(data2[1:11000, 1:20], 1, cumsum_ordered) |> str()
 num [1:20, 1:11000] 0.185941 0.000727 0.326461 0.117036 0.020571 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:20] "PANCAN..TCGA..TCGA.OR.A5J1.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J2.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J3.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J5.01A.11R.A29S.07" ...
  ..$ : chr [1:11000] "TSPAN6" "DPM1" "SCYL3" "C1orf112" ...
Browse[3]> apply(data2[1:12000, 1:20], 1, cumsum_ordered) |> str()
 num [1:20, 1:12000] 0.185941 0.000727 0.326461 0.117036 0.020571 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:20] "PANCAN..TCGA..TCGA.OR.A5J1.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J2.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J3.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J5.01A.11R.A29S.07" ...
  ..$ : chr [1:12000] "TSPAN6" "DPM1" "SCYL3" "C1orf112" ...
Browse[3]> apply(data2[1:12500, 1:20], 1, cumsum_ordered) |> str()
 num [1:20, 1:12500] 0.185941 0.000727 0.326461 0.117036 0.020571 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:20] "PANCAN..TCGA..TCGA.OR.A5J1.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J2.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J3.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J5.01A.11R.A29S.07" ...
  ..$ : chr [1:12500] "TSPAN6" "DPM1" "SCYL3" "C1orf112" ...
Browse[3]> apply(data2[1:127500, 1:20], 1, cumsum_ordered) |> str()
Error in data2[1:127500, 1:20] : subscript out of bounds
Browse[3]> apply(data2[1:12750, 1:20], 1, cumsum_ordered) |> str()
 num [1:20, 1:12750] 0.185941 0.000727 0.326461 0.117036 0.020571 ...
 - attr(*, "dimnames")=List of 2
  ..$ : NULL
  ..$ : chr [1:12750] "TSPAN6" "DPM1" "SCYL3" "C1orf112" ...
Browse[3]> apply(data2[1:12700, 1:20], 1, cumsum_ordered) |> str()
 num [1:20, 1:12700] 0.185941 0.000727 0.326461 0.117036 0.020571 ...
 - attr(*, "dimnames")=List of 2
  ..$ : NULL
  ..$ : chr [1:12700] "TSPAN6" "DPM1" "SCYL3" "C1orf112" ...
Browse[3]> apply(data2[1:12600, 1:20], 1, cumsum_ordered) |> str()
 num [1:20, 1:12600] 0.185941 0.000727 0.326461 0.117036 0.020571 ...
 - attr(*, "dimnames")=List of 2
  ..$ : NULL
  ..$ : chr [1:12600] "TSPAN6" "DPM1" "SCYL3" "C1orf112" ...
Browse[3]> apply(data2[1:12550, 1:20], 1, cumsum_ordered) |> str()
 num [1:20, 1:12550] 0.185941 0.000727 0.326461 0.117036 0.020571 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:20] "PANCAN..TCGA..TCGA.OR.A5J1.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J2.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J3.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J5.01A.11R.A29S.07" ...
  ..$ : chr [1:12550] "TSPAN6" "DPM1" "SCYL3" "C1orf112" ...
Browse[3]> apply(data2[1:12570, 1:20], 1, cumsum_ordered) |> str()
 num [1:20, 1:12570] 0.185941 0.000727 0.326461 0.117036 0.020571 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:20] "PANCAN..TCGA..TCGA.OR.A5J1.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J2.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J3.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J5.01A.11R.A29S.07" ...
  ..$ : chr [1:12570] "TSPAN6" "DPM1" "SCYL3" "C1orf112" ...
Browse[3]> apply(data2[1:12580, 1:20], 1, cumsum_ordered) |> str()
 num [1:20, 1:12580] 0.185941 0.000727 0.326461 0.117036 0.020571 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:20] "PANCAN..TCGA..TCGA.OR.A5J1.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J2.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J3.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J5.01A.11R.A29S.07" ...
  ..$ : chr [1:12580] "TSPAN6" "DPM1" "SCYL3" "C1orf112" ...
Browse[3]> apply(data2[1:12590, 1:20], 1, cumsum_ordered) |> str()
 num [1:20, 1:12590] 0.185941 0.000727 0.326461 0.117036 0.020571 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:20] "PANCAN..TCGA..TCGA.OR.A5J1.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J2.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J3.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J5.01A.11R.A29S.07" ...
  ..$ : chr [1:12590] "TSPAN6" "DPM1" "SCYL3" "C1orf112" ...
Browse[3]> apply(data2[1:12595, 1:20], 1, cumsum_ordered) |> str()
 num [1:20, 1:12595] 0.185941 0.000727 0.326461 0.117036 0.020571 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:20] "PANCAN..TCGA..TCGA.OR.A5J1.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J2.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J3.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J5.01A.11R.A29S.07" ...
  ..$ : chr [1:12595] "TSPAN6" "DPM1" "SCYL3" "C1orf112" ...
Browse[3]> apply(data2[1:12596, 1:20], 1, cumsum_ordered) |> str()
 num [1:20, 1:12596] 0.185941 0.000727 0.326461 0.117036 0.020571 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:20] "PANCAN..TCGA..TCGA.OR.A5J1.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J2.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J3.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J5.01A.11R.A29S.07" ...
  ..$ : chr [1:12596] "TSPAN6" "DPM1" "SCYL3" "C1orf112" ...
Browse[3]> apply(data2[1:12597, 1:20], 1, cumsum_ordered) |> str()
 num [1:20, 1:12597] 0.185941 0.000727 0.326461 0.117036 0.020571 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:20] "PANCAN..TCGA..TCGA.OR.A5J1.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J2.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J3.01A.11R.A29S.07" "PANCAN..TCGA..TCGA.OR.A5J5.01A.11R.A29S.07" ...
  ..$ : chr [1:12597] "TSPAN6" "DPM1" "SCYL3" "C1orf112" ...
Browse[3]> apply(data2[1:12598, 1:20], 1, cumsum_ordered) |> str()
 num [1:20, 1:12598] 0.185941 0.000727 0.326461 0.117036 0.020571 ...
 - attr(*, "dimnames")=List of 2
  ..$ : NULL
  ..$ : chr [1:12598] "TSPAN6" "DPM1" "SCYL3" "C1orf112" ...
Browse[3]> data2[12596:12600, 1:3]
          PANCAN..TCGA..TCGA.OR.A5J1.01A.11R.A29S.07 PANCAN..TCGA..TCGA.OR.A5J2.01A.11R.A29S.07 PANCAN..TCGA..TCGA.OR.A5J3.01A.11R.A29S.07
LINC00174                                   3.831145                                   1.632413                                   2.816775
GEMIN4                                      4.424853                                   3.973434                                   3.909745
HNRNPCL4                                   -7.617169                                  -7.617169                                  -7.617169
FJX1                                        1.247121                                   3.283797                                   4.053823
KLHL28                                      4.253055                                   4.216638                                   5.032185
pdiakumis commented 11 months ago

Closing. Issue potentially related to local setup of BLAS/LAPACK accelerated libraries which do not generate NaNs for genes with 0 sd/variance.