Merge pull request #401 from GoekeLab/bugFixIncompatibleCounts

fix overcounting issue for incompatible counts

Author: cying111

.github/workflows/check-bioc.yml

@@ -54,9 +54,9 @@ jobs:
       fail-fast: false
       matrix:
         config:
-          - { os: ubuntu-latest, r: '4.2', bioc: '3.16', cont: "bioconductor/bioconductor_docker:RELEASE_3_16", rspm: "https://packagemanager.rstudio.com/cran/__linux__/focal/latest" }
-          - { os: macOS-latest, r: '4.2', bioc: '3.16'}
-          - { os: windows-latest, r: '4.2', bioc: '3.16'}
+          - { os: ubuntu-latest, r: '4.3', bioc: '3.18', cont: "bioconductor/bioconductor_docker:RELEASE_3_18", rspm: "https://packagemanager.rstudio.com/cran/__linux__/focal/latest" }
+          - { os: macOS-latest, r: '4.3', bioc: '3.18'}
+            ##- { os: windows-latest, r: '4.3', bioc: '3.18'}
           ## Check https://github.com/r-lib/actions/tree/master/examples
           ## for examples using the http-user-agent
     env:

R/bambu-processReads_utilityConstructReadClasses.R

@@ -246,11 +246,12 @@ constructUnsplicedReadClasses <- function(reads.singleExon, annotations,
     counts = as.data.frame(reads.singleExon) %>% 
         mutate(id = mcols(reads.singleExon)$id) %>% 
         group_by(seqnames,start,end,strand) %>% 
-        summarise(n=n(), id = list(id)) %>% 
+        mutate(n=n(), id = list(id)) %>%  # change summarise to mutate as summarise will reorder the table
+        ungroup() %>%
         as.data.frame()
-    reads.singleExon = unique(reads.singleExon)
     mcols(reads.singleExon)$counts <- counts$n
     mcols(reads.singleExon)$id <- counts$id
+     reads.singleExon = unique(reads.singleExon)
 
     rcUnsplicedAnnotation <- getUnsplicedReadClassByReference(
         granges = reads.singleExon, grangesReference = referenceExons,

R/bambu-quantify_utilityFunctions.R

@@ -20,8 +20,8 @@ modifyIncompatibleAssignment <- function(distTable){
 #' Process incompatible counts
 #' @noRd
 processIncompatibleCounts <- function(readClassDist){
-  distTable <- data.table(as.data.frame(metadata(readClassDist)$distTable))[, 
-                                                                            .(readClassId, annotationTxId, readCount, GENEID, dist,equal)]
+  distTable <- unique(data.table(as.data.frame(metadata(readClassDist)$distTable))[, 
+               .(readClassId, annotationTxId, readCount, GENEID, equal)], by = NULL)
   distTableIncompatible <- distTable[grep("unidentified", annotationTxId)]
   # filter out multiple geneIDs mapped to the same readClass using rowData(se)
   geneRCMap <- as.data.table(as.data.frame(rowData(readClassDist)),