Skip to content

Commit b3be667

Browse files
lingminhaolingminhao
authored
Merge pull request #565 from GoekeLab/tidy_readClassFile_quantData
Tidy readClassList and quantData
2 parents 6dcb1d5 + d112d3a commit b3be667

10 files changed

Lines changed: 313 additions & 170 deletions

‎R/bambu-assignDist.R‎

Lines changed: 31 additions & 76 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,13 @@ assignReadClasstoTranscripts <- function(readClassList, annotations, rcAssignmen
1010
readClassList <- splitReadClassFiles(readClassList)
1111
readClassDt <- genEquiRCs(metadata(readClassList)$readClassDist, annotations, verbose)
1212
readClassDt$eqClass.match = match(readClassDt$eqClassById,metadata(readClassList)$eqClassById)
13+
14+
# Add columnIds and columnCounts columns
15+
readClassDt[, `:=`(columnIds = vector("list", .N), columnCounts = vector("list", .N))]
16+
observedIdx <- which(!is.na(readClassDt$eqClass.match))
17+
readClassDt$columnIds[observedIdx] <- metadata(readClassList)$columnIds[readClassDt$eqClass.match[observedIdx]]
18+
readClassDt$columnCounts[observedIdx] <- metadata(readClassList)$columnCounts[readClassDt$eqClass.match[observedIdx]]
19+
1320
readClassDt <- simplifyNames(readClassDt)
1421
readClassDt <- readClassDt %>% group_by(eqClassId, gene_sid) %>%
1522
mutate(multi_align = length(unique(txid))>1) %>%
@@ -18,93 +25,41 @@ assignReadClasstoTranscripts <- function(readClassList, annotations, rcAssignmen
1825
data.table()
1926
#return non-em counts
2027
ColData <- generateColData(readClassList, sampleMetadata, extractBarcodeUMI)
21-
quantData <- SummarizedExperiment(assays = SimpleList(
22-
counts = generateUniqueCounts(readClassDt, metadata(readClassList)$countMatrix, annotations)),
23-
rowRanges = annotations,
24-
colData = ColData)
25-
colnames(quantData) <- ColData$id
26-
if(sum(metadata(readClassList)$incompatibleCountMatrix)==0){
27-
metadata(quantData)$incompatibleCounts <- NULL
28-
}else{
29-
metadata(quantData)$incompatibleCounts <- generateIncompatibleCounts(metadata(readClassList)$incompatibleCountMatrix, annotations)
30-
}
31-
metadata(quantData)$nonuniqueCounts <- generateNonUniqueCounts(readClassDt, metadata(readClassList)$countMatrix, annotations)
32-
metadata(quantData)$readClassDt <- readClassDt
33-
metadata(quantData)$countMatrix <- metadata(readClassList)$countMatrix
34-
metadata(quantData)$incompatibleCountMatrix <- metadata(readClassList)$incompatibleCountMatrix
35-
metadata(quantData)$sampleName <- metadata(readClassList)$sampleData$sampleName
36-
if(returnDistTable)
37-
metadata(quantData)$distTable <- metadata(metadata(readClassList)$readClassDist)$distTableOld
38-
39-
if(trackReads)
40-
metadata(quantData)$readToTranscriptMap <-
41-
generateReadToTranscriptMap(readClassList,
42-
metadata(readClassList)$readClassDist,
43-
annotations)
28+
29+
incompatibleCountMatrix <- metadata(readClassList)$incompatibleCountMatrix
30+
incompatibleCounts <- generateIncompatibleCounts(incompatibleCountMatrix, annotations)
4431

45-
return(quantData)
32+
distTable <- if(returnDistTable) {
33+
metadata(metadata(readClassList)$readClassDist)$distTableOld
34+
} else{
35+
NULL
36+
}
4637

47-
}
38+
readToTranscriptMap <- if(trackReads) {
39+
generateReadToTranscriptMap(readClassList, metadata(readClassList)$readClassDist,annotations)
40+
} else{
41+
NULL
42+
}
4843

49-
#' Generate unique counts
50-
#' @noRd
51-
generateUniqueCounts <- function(readClassDt, countMatrix, annotations){
52-
x <- readClassDt %>% filter(!multi_align & !is.na(eqClass.match))
53-
uniqueCounts <- countMatrix[x$eqClass.match,]
54-
uniqueCounts.tx <- sparse.model.matrix(~ factor(x$txid) - 1)
55-
uniqueCounts <- t(uniqueCounts.tx) %*% uniqueCounts
56-
rownames(uniqueCounts) <- names(annotations)[match(as.numeric(levels(factor(x$txid))),mcols(annotations)$txid)]
57-
counts <- sparseMatrix(length(annotations), ncol(uniqueCounts), x = 0)
58-
rownames(counts) <- names(annotations)
59-
counts[rownames(uniqueCounts),] <- uniqueCounts
60-
return(counts)
61-
62-
# these three lines appear after return, so it's not used, is this used for debug only?
63-
# counts.total = colSums(countMatrix) + colSums(incompatibleCountMatrix)
64-
# counts.total[counts.total==0] = 1
65-
# counts.CPM = counts/counts.total * 10^6
44+
quantData <- constructQuantData(
45+
sampleData = data.frame(ColData),
46+
readClassDt = readClassDt,
47+
incompatibleCounts = incompatibleCounts,
48+
distTable = distTable,
49+
readToTranscriptMap = readToTranscriptMap
50+
)
6651

52+
return(quantData)
6753
}
6854

69-
7055
#' Generate incompatible counts
7156
#' @noRd
7257
generateIncompatibleCounts <- function(incompatibleCountMatrix, annotations){
7358
genes <- levels(factor(unique(mcols(annotations)$GENEID)))
7459
rownames(incompatibleCountMatrix) <- genes[as.numeric(rownames(incompatibleCountMatrix))]
7560
geneMat <- sparseMatrix(length(genes), ncol(incompatibleCountMatrix), x = 0)
7661
rownames(geneMat) <- genes
77-
geneMat[rownames(incompatibleCountMatrix),] <- incompatibleCountMatrix
78-
return(geneMat)
79-
}
80-
81-
82-
#' Generate non-unique counts
83-
#' @noRd
84-
generateNonUniqueCounts <- function(readClassDt, countMatrix, annotations){
85-
#fuse multi align RCs by gene
86-
x <- readClassDt %>% filter(multi_align & !is.na(eqClass.match))
87-
x <- x %>% distinct(eqClassId, .keep_all = TRUE)
88-
nonuniqueCounts <- countMatrix[x$eqClass.match,, drop = FALSE]
89-
if(nrow(x)>1 & length(unique(x$gene_sid))>1){
90-
nonuniqueCounts.gene <- sparse.model.matrix(~ factor(x$gene_sid) - 1)
91-
nonuniqueCounts <- t(nonuniqueCounts.gene) %*% nonuniqueCounts
92-
} else{
93-
warning("The factor variable 'gene_sid' has only one level. Adjusting output.")
94-
nonuniqueCounts.gene <- Matrix(1, nrow = nrow(x), ncol = 1, sparse = TRUE)
95-
nonuniqueCounts <- t(nonuniqueCounts.gene) %*% nonuniqueCounts
96-
}
97-
#covert ids into gene ids
98-
geneids <- as.numeric(levels(factor(x$gene_sid)))
99-
geneids <- x$txid[match(geneids, x$gene_sid)]
100-
geneids <- mcols(annotations)$GENEID[as.numeric(geneids)]
101-
rownames(nonuniqueCounts) <- geneids
102-
#create matrix for all annotated genes
103-
genes <- levels(factor(unique(mcols(annotations)$GENEID)))
104-
geneMat <- sparseMatrix(length(genes), ncol(nonuniqueCounts), x = 0)
105-
rownames(geneMat) <- genes
106-
if(!is.null(rownames(nonuniqueCounts))){
107-
geneMat[rownames(nonuniqueCounts),] <- nonuniqueCounts
108-
}
109-
return(geneMat)
62+
colnames(geneMat) <- colnames(incompatibleCountMatrix)
63+
geneMat[match(rownames(incompatibleCountMatrix), rownames(geneMat)), ] <- incompatibleCountMatrix
64+
return(geneMat[unique(mcols(annotations)$GENEID), , drop = FALSE])
11065
}

‎R/bambu-classes.R‎

Lines changed: 78 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,78 @@
1+
#' @import methods
2+
#' @importFrom Matrix Matrix
3+
#' @importFrom data.table data.table
4+
#' @importClassesFrom Matrix Matrix
5+
#' @importClassesFrom data.table data.table
6+
7+
# Valid values for metadata(se)$seType, which identifies the SE variant:
8+
# EMCounts — transcript-level SE with EM-estimated counts (main bambu output)
9+
# geneCounts — gene-level SE derived by collapsing transcript counts
10+
# uniqueCounts — transcript-level SE with uniquely-assigned counts only, no EM
11+
SE_TYPES <- c(
12+
EMCounts = "EMCounts",
13+
geneCounts = "geneCounts",
14+
uniqueCounts = "uniqueCounts"
15+
)
16+
17+
# quantData holds per-sample intermediate results from the assignDist step.
18+
# distTable and readToTranscriptMap are optional: populated only when
19+
# returnDistTable=TRUE or trackReads=TRUE respectively, and lifted into
20+
# metadata(countsSe) at the end of bambu().
21+
setClass("quantData",
22+
slots = c(
23+
sampleData = "data.frame", # per-sample metadata (id, sampleName)
24+
readClassDt = "data.table", # read-class-level count and assignment data
25+
incompatibleCounts = "sparseMatrix", # counts of reads incompatible with any annotation
26+
distTable = "ANY", # read-class-to-transcript compatibility table (DataFrame or NULL)
27+
readToTranscriptMap = "ANY" # per-read assignment to transcripts (tibble or NULL)
28+
),
29+
prototype = list(
30+
sampleData = data.frame(id = integer(), sampleName = character()),
31+
readClassDt = data.table::data.table()
32+
),
33+
validity = function(object) {
34+
errs <- character()
35+
if (!all(c("id", "sampleName") %in% names(object@sampleData)))
36+
errs <- c(errs, "sampleData must have columns 'id' and 'sampleName'")
37+
if (!is.null(object@distTable) && !is(object@distTable, "DataFrame"))
38+
errs <- c(errs, "distTable must be a DataFrame or NULL")
39+
if (!is.null(object@readToTranscriptMap) && !is(object@readToTranscriptMap, "tbl_df"))
40+
errs <- c(errs, "readToTranscriptMap must be a tibble or NULL")
41+
if (length(errs)) errs else TRUE
42+
})
43+
44+
#' Construct a quantData object
45+
#' @noRd
46+
constructQuantData <- function(sampleData, readClassDt,
47+
incompatibleCounts,
48+
distTable = NULL,
49+
readToTranscriptMap = NULL) {
50+
new("quantData",
51+
sampleData = sampleData,
52+
readClassDt = readClassDt,
53+
incompatibleCounts = incompatibleCounts,
54+
distTable = distTable,
55+
readToTranscriptMap = readToTranscriptMap)
56+
}
57+
58+
#' @noRd
59+
setGeneric("getSampleData", function(x) standardGeneric("getSampleData"))
60+
#' @noRd
61+
setGeneric("getReadClassDt", function(x) standardGeneric("getReadClassDt"))
62+
#' @noRd
63+
setGeneric("getIncompatibleCounts", function(x) standardGeneric("getIncompatibleCounts"))
64+
#' @noRd
65+
setGeneric("getDistTable", function(x) standardGeneric("getDistTable"))
66+
#' @noRd
67+
setGeneric("getReadToTranscriptMap", function(x) standardGeneric("getReadToTranscriptMap"))
68+
69+
#' @noRd
70+
setMethod("getSampleData", "quantData", function(x) x@sampleData)
71+
#' @noRd
72+
setMethod("getReadClassDt", "quantData", function(x) x@readClassDt)
73+
#' @noRd
74+
setMethod("getIncompatibleCounts", "quantData", function(x) x@incompatibleCounts)
75+
#' @noRd
76+
setMethod("getDistTable", "quantData", function(x) x@distTable)
77+
#' @noRd
78+
setMethod("getReadToTranscriptMap", "quantData", function(x) x@readToTranscriptMap)

‎R/bambu-processReads.R‎

Lines changed: 25 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -56,6 +56,7 @@ bambu.processReads <- function(reads, annotations, genomeSequence,
5656
processByChromosome = processByChromosome, trackReads = trackReads, fusionMode = fusionMode,
5757
extractBarcodeUMI = extractBarcodeUMI, dedupUMI = dedupUMI, index = 1)},
5858
BPPARAM = bpParameters)
59+
names(readClassList) <- names(reads)
5960
} else {
6061
readGrgList <- bplapply(seq_along(reads), function(i) {
6162
bambu.readsByFile(bam.file = reads[i],
@@ -162,22 +163,23 @@ bambu.processReadsByFile <- function(bam.file, genomeSequence, annotations,
162163

163164
mcols(readGrgList)$id <- seq_along(readGrgList)
164165

165-
if(extractBarcodeUMI){
166-
mcols(readGrgList)$sampleID <- as.numeric(mcols(readGrgList)$CB)
166+
if(extractBarcodeUMI){
167+
mcols(readGrgList)$columnID <- as.numeric(mcols(readGrgList)$CB)
167168
} else {
168-
mcols(readGrgList)$sampleID <- index
169+
mcols(readGrgList)$columnID <- index
169170
}
170171

172+
runName <- names(bam.file)[1]
171173
# construct read classes for each chromosome seperately
172174
if(processByChromosome){
173175
se <- lowMemoryConstructReadClasses(readGrgList, genomeSequence,
174-
annotations, stranded, verbose,bam.file)
176+
annotations, stranded, verbose, bam.file)
175177
} else{
176178
unlisted_junctions <- unlistIntrons(readGrgList, use.ids = TRUE)
177179
uniqueJunctions <- isore.constructJunctionTables(unlisted_junctions,
178180
annotations,genomeSequence, stranded = stranded, verbose = verbose)
179181
se <- isore.constructReadClasses(readGrgList,
180-
unlisted_junctions, uniqueJunctions, runName = "TODO",
182+
unlisted_junctions, uniqueJunctions, runName = runName,
181183
annotations, stranded, verbose)
182184

183185
}
@@ -297,18 +299,18 @@ bambu.readsByFile <- function(bam.file, genomeSequence, annotations,
297299
constructReadClasses <- function(readGrgList, genomeSequence, annotations,
298300
stranded = FALSE, min.readCount = 2,
299301
fitReadClassModel = TRUE, min.exonOverlap = 10, defaultModels = NULL, returnModel = FALSE,
300-
verbose = FALSE, processByChromosome = FALSE, trackReads = FALSE, fusionMode = FALSE){
302+
verbose = FALSE, processByChromosome = FALSE, trackReads = FALSE, fusionMode = FALSE, runName = "sample"){
301303

302304
if(processByChromosome){
303305
# construct read classes for each chromosome seperately
304306
se <- lowMemoryConstructReadClasses(readGrgList, genomeSequence,
305-
annotations, stranded, verbose,"TODO", fusionMode)
307+
annotations, stranded, verbose, runName, fusionMode)
306308
} else{
307309
unlisted_junctions <- unlistIntrons(readGrgList, use.ids = TRUE)
308310
uniqueJunctions <- isore.constructJunctionTables(unlisted_junctions,
309311
annotations,genomeSequence, stranded = stranded, verbose = verbose)
310312
se <- isore.constructReadClasses(readGrgList,
311-
unlisted_junctions, uniqueJunctions, runName = "TODO",
313+
unlisted_junctions, uniqueJunctions, runName = runName,
312314
annotations, stranded, verbose)
313315

314316
}
@@ -336,21 +338,22 @@ constructReadClasses <- function(readGrgList, genomeSequence, annotations,
336338
#' Low memory mode for construct read classes (processByChromosome)
337339
#' @noRd
338340
lowMemoryConstructReadClasses <- function(readGrgList, genomeSequence,
339-
annotations, stranded, verbose,bam.file, fusionMode = FALSE){
341+
annotations, stranded, verbose, bam.file, fusionMode = FALSE){
340342
if(fusionMode){
341343
readGrgList <- list(readGrgList)
342344
names(readGrgList) <- c("fusion")
343345
} else{
344346
readGrgList <- split(readGrgList, getChrFromGrList(readGrgList))
345347
}
348+
runName <- names(bam.file)[1]
346349
se <- lapply(names(readGrgList),FUN = function(i){
347350
if(length(readGrgList[[i]]) == 0) return(NULL)
348351
# create error and strand corrected junction tables
349352
unlisted_junctions <- unlistIntrons(readGrgList[[i]], use.ids = TRUE)
350353
uniqueJunctions <- isore.constructJunctionTables(unlisted_junctions,
351354
annotations,genomeSequence, stranded = stranded, verbose = verbose)
352355
se.temp <- isore.constructReadClasses(readGrgList[[i]],
353-
unlisted_junctions, uniqueJunctions, runName = "TODO",
356+
unlisted_junctions, uniqueJunctions, runName = runName,
354357
annotations, stranded, verbose)
355358
return(se.temp)
356359
})
@@ -385,10 +388,12 @@ splitReadClassFiles = function(readClassFile){
385388
distTable <- metadata(metadata(readClassFile)$readClassDist)$distTable
386389
eqClasses <- distTable %>% group_by(eqClassById) %>%
387390
distinct(eqClassById, readCount,GENEID, totalWidth, firstExonWidth, .keep_all = TRUE)
388-
eqClasses$sampleIDs <- rowData(readClassFile)$sampleIDs[match(eqClasses$readClassId, rownames(readClassFile))]
391+
eqClasses$columnIds <- rowData(readClassFile)$columnIds[match(eqClasses$readClassId, rownames(readClassFile))]
389392
eqClasses <- eqClasses %>% summarise(nobs = sum(readCount),
390-
sampleIDs = list(unlist(sampleIDs)))
391-
counts.table <- tableFunction(eqClasses$sampleIDs)
393+
columnIds = list(unlist(columnIds)))
394+
counts.table <- tableFunction(eqClasses$columnIds)
395+
metadata(readClassFile)$columnIds <- lapply(counts.table, function(x) as.numeric(names(x)))
396+
metadata(readClassFile)$columnCounts <- lapply(counts.table, function(x) as.numeric(x))
392397
counts <- sparseMatrix(
393398
i = rep(seq_along(counts.table), lengths(counts.table)),
394399
j = as.numeric(names(unlist(counts.table))),
@@ -397,14 +402,14 @@ splitReadClassFiles = function(readClassFile){
397402
#incompatible counts
398403
distTable <- metadata(metadata(readClassFile)$readClassDist)$distTable.incompatible
399404
if(nrow(distTable)==0) {
400-
counts.incompatible <- sparseMatrix(i= 1, j = 1, x = 0,
401-
dims = c(1, length(metadata(readClassFile)$sampleData$id)))
402-
rownames(counts.incompatible) <- "TODO"
405+
counts.incompatible <- sparseMatrix(i= integer(0), j = integer(0), x = numeric(0),
406+
dims = c(0, length(metadata(readClassFile)$sampleData$id)))
407+
rownames(counts.incompatible) <- character(0)
403408
} else{
404-
distTable$sampleIDs <- rowData(readClassFile)$sampleIDs[match(distTable$readClassId, rownames(readClassFile))]
409+
distTable$columnIds <- rowData(readClassFile)$columnIds[match(distTable$readClassId, rownames(readClassFile))]
405410
distTable <- distTable %>% group_by(GENEID.i) %>% summarise(counts = sum(readCount),
406-
sampleIDs = list(unlist(sampleIDs)))
407-
counts.table <- lapply(distTable$sampleIDs, FUN = function(x){table(x)})
411+
columnIds = list(unlist(columnIds)))
412+
counts.table <- lapply(distTable$columnIds, FUN = function(x){table(x)})
408413
counts.incompatible <- sparseMatrix(
409414
i = rep(seq_along(counts.table), lengths(counts.table)),
410415
j = as.numeric(names(unlist(counts.table))),
@@ -416,7 +421,6 @@ splitReadClassFiles = function(readClassFile){
416421
colnames(counts) <- metadata(readClassFile)$sampleData$id
417422
metadata(readClassFile)$eqClassById <- eqClasses$eqClassById
418423
#rownames(counts) = eqClasses$eqClassById
419-
metadata(readClassFile)$countMatrix <- counts
420424
metadata(readClassFile)$incompatibleCountMatrix <- counts.incompatible
421425
return(readClassFile)
422426
}
@@ -426,7 +430,7 @@ splitReadClassFiles = function(readClassFile){
426430
#' @importFrom Matrix
427431
#' @noRd
428432
splitReadClassFilesByRC <- function(readClassFile){
429-
counts.table <- tableFunction(rowData(readClassFile)$sampleIDs)
433+
counts.table <- tableFunction(rowData(readClassFile)$columnIds)
430434
counts <- sparseMatrix(
431435
i = rep(seq_along(counts.table), lengths(counts.table)),
432436
j = as.numeric(names(unlist(counts.table))),

0 commit comments

Comments
 (0)