|
| 1 | +#' Create equivilence classes and assign to transcripts |
| 2 | +#' @inheritParams bambu |
| 3 | +#' @import data.table |
| 4 | +#' @noRd |
| 5 | +assignReadClasstoTranscripts <- function(readClassList, annotations, isoreParameters, |
| 6 | + verbose, demultiplexed, spatial, |
| 7 | + returnDistTable = FALSE, trackReads = TRUE) { |
| 8 | + if (is.character(readClassList)) readClassList <- readRDS(file = readClassList) |
| 9 | + metadata(readClassList)$readClassDist <- calculateDistTable(readClassList, annotations, isoreParameters, verbose, returnDistTable) |
| 10 | + readClassList <- splitReadClassFiles(readClassList) |
| 11 | + readClassDt <- genEquiRCs(metadata(readClassList)$readClassDist, annotations, verbose) |
| 12 | + readClassDt$eqClass.match = match(readClassDt$eqClassById,metadata(readClassList)$eqClassById) |
| 13 | + readClassDt <- simplifyNames(readClassDt) |
| 14 | + readClassDt <- readClassDt %>% group_by(eqClassId, gene_sid) %>% |
| 15 | + mutate(multi_align = length(unique(txid))>1) %>% |
| 16 | + ungroup() %>% |
| 17 | + mutate(aval = 1) %>% |
| 18 | + data.table() |
| 19 | + #return non-em counts |
| 20 | + ColData <- generateColData(colnames(metadata(readClassList)$countMatrix), clusters = NULL, demultiplexed, spatial) |
| 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)$sampleNames <- metadata(readClassList)$sampleNames |
| 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) |
| 44 | + |
| 45 | + return(quantData) |
| 46 | + |
| 47 | +} |
| 48 | + |
| 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 |
| 66 | + |
| 67 | +} |
| 68 | + |
| 69 | + |
| 70 | +#' Generate incompatible counts |
| 71 | +#' @noRd |
| 72 | +generateIncompatibleCounts <- function(incompatibleCountMatrix, annotations){ |
| 73 | + genes <- levels(factor(unique(mcols(annotations)$GENEID))) |
| 74 | + rownames(incompatibleCountMatrix) <- genes[as.numeric(rownames(incompatibleCountMatrix))] |
| 75 | + geneMat <- sparseMatrix(length(genes), ncol(incompatibleCountMatrix), x = 0) |
| 76 | + 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) |
| 110 | +} |
0 commit comments