diff --git a/DESCRIPTION b/DESCRIPTION index 0695c58..727272a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -2,6 +2,8 @@ Package: chimeraviz Type: Package Title: Visualization tools for gene fusions Version: 1.5.4 +Author: Stian Lågstad +Maintainer: Stian Lågstad Authors@R: c( person( given = "Stian", diff --git a/NAMESPACE b/NAMESPACE index 62d5594..b30ab43 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -19,6 +19,7 @@ export(getFusionById) export(getTranscriptsEnsembldb) export(importDefuse) export(importEricscript) +export(importFusionInspector) export(importFusioncatcher) export(importFusionmap) export(importInfusion) @@ -31,6 +32,7 @@ export(partnerGeneJunctionSequence) export(plotCircle) export(plotFusion) export(plotFusionReads) +export(plotFusionReadsSimple) export(plotFusionSeparate) export(plotFusionTogether) export(plotFusionTranscript) diff --git a/R/importStarfusion.R b/R/importStarfusion.R index 10ff00e..f69d6a8 100644 --- a/R/importStarfusion.R +++ b/R/importStarfusion.R @@ -1,3 +1,53 @@ +#' read the results from runnning FusionInspector +#' +#' Note: there may not be perfect match between the fusions found by STAR-Fusion +#' and FusionInspector (since the latter reruns STAR on the mini contig; also, +#' selfie fusions are ignored). +#' If no match can be found, the whole fusion is skipped +#' +#' @export +importFusionInspector <- function (filename='FusionInspector-inspect/finspector.igv.FusionJuncSpan', + limit=Inf) { + ## Try to read the FusionInspector report. + ## Within one 'scaffold', only keep the leftmost and rightmost coordinate + report <- withCallingHandlers({ + col_types_fusioninspector = readr::cols_only( + "#scaffold" = col_character(), + "fusion_break_name" = col_skip(), + "break_left" = col_integer(), + "break_right" = col_integer(), + "num_junction_reads" = col_skip(), + "num_spanning_frags" = col_skip(), + "spanning_frag_coords" = col_skip() + ) + if (missing(limit)) { + readr::read_tsv( + file = filename, + col_types = col_types_fusioninspector) + } else { + readr::read_tsv( + file = filename, + col_types = col_types_fusioninspector, + n_max = limit) + } + }, + error = function(cond) { + message(paste0("Reading ", filename, " caused an error: ", cond[[1]])) + stop(cond) + }, + warning = function(cond) { + message(paste0("Reading ", filename, " caused a warning: ", cond[[1]])) + warning(cond) + }) + colnames(report)[1] <- 'id' + d <- + data.frame(id=with(report, id[ !duplicated(id) ]), + start=with(report, tapply(break_left, id, min)), + end=with(report, tapply(break_right, id, max))) + rownames(d) <- d$id + d +} #importFusionInspector + #' Import results from a STAR-Fusion run into a list of Fusion objects. #' #' A function that imports the results from a STAR-Fusion run, typically from @@ -20,8 +70,8 @@ #' # This should import a list of 3 fusions described in Fusion objects. #' #' @export -importStarfusion <- function (filename, genomeVersion, limit) { - +importStarfusion <- function (filename, genomeVersion, limit=Inf, + useFusionInspector=FALSE) { # Is the genome version valid? validGenomes <- c("hg19", "hg38", "mm10") if (is.na(match(tolower(genomeVersion), tolower(validGenomes)))) { @@ -39,7 +89,8 @@ importStarfusion <- function (filename, genomeVersion, limit) { report <- withCallingHandlers( { col_types_starfusion = readr::cols_only( - "#FusionName" = col_skip(), +## "#FusionName" = col_skip(), + "#FusionName" = col_character(), "JunctionReadCount" = col_integer(), "SpanningFragCount" = col_integer(), "SpliceType" = col_skip(), @@ -52,7 +103,8 @@ importStarfusion <- function (filename, genomeVersion, limit) { "LeftBreakEntropy" = col_number(), "RightBreakDinuc" = col_character(), "RightBreakEntropy" = col_number(), - "FFPM" = col_number() + "FFPM" = col_number(), + "PROT_FUSION_TYPE" = col_character() ) if (missing(limit)) { # Read all lines @@ -79,10 +131,15 @@ importStarfusion <- function (filename, genomeVersion, limit) { } ) + if(useFusionInspector) { + l <- limit+10 + fi.table <- importFusionInspector(limit=l) + } + # Set variables id <- NA inframe <- NA - fusionTool <- "starfusion" + fusionTool <- ifelse(useFusionInspector, "starfusion+fusioninspector", "starfusion") spanningReadsCount <- NA splitReadsCount <- NA junctionSequence <- NA @@ -96,6 +153,7 @@ importStarfusion <- function (filename, genomeVersion, limit) { # Import starfusion-specific fields fusionToolSpecificData <- list() + fusionToolSpecificData[["FusionName"]] = report[[i, "#FusionName"]] fusionToolSpecificData[["LargeAnchorSupport"]] = report[[i, "LargeAnchorSupport"]] fusionToolSpecificData[["LeftBreakDinuc"]] = report[[i, "LeftBreakDinuc"]] fusionToolSpecificData[["LeftBreakEntropy"]] = report[[i, "LeftBreakEntropy"]] @@ -103,6 +161,17 @@ importStarfusion <- function (filename, genomeVersion, limit) { fusionToolSpecificData[["RightBreakEntropy"]] = report[[i, "RightBreakEntropy"]] fusionToolSpecificData[["FFPM"]] = report[[i, "FFPM"]] + if(!is.null(report$PROT_FUSION_TYPE)) { + ## only available when loading from coding_effect file (i.e. + ## star-fusion.fusion_predictions.abridged.annotated.coding_effect.tsv) + ft <- report[[i, "PROT_FUSION_TYPE"]] + fusionToolSpecificData[["inframe"]] <- ft + if (ft == 'INFRAME') + inframe <- TRUE + if (ft == 'FRAMESHIFT') + inframe <- FALSE + } + # id for this fusion id <- as.character(i) @@ -139,6 +208,31 @@ importStarfusion <- function (filename, genomeVersion, limit) { ensemblIdA <- geneNames1[2] ensemblIdB <- geneNames2[2] + if (useFusionInspector) { + ## override things to go with the FI's virtual 'mini genome', but + ## first save them in fusionToolSpecificData: + fusionToolSpecificData[['orig_chromosomeA']] <- chromosomeA + fusionToolSpecificData[['orig_breakpointA']] <- breakpointA + fusionToolSpecificData[['orig_strandA']] <- strandA + fusionToolSpecificData[['orig_chromosomeB']] <- chromosomeB + fusionToolSpecificData[['orig_breakpointB']] <- breakpointB + fusionToolSpecificData[['orig_strandB']] <- strandB + + fusion.name <- report[[i, "#FusionName"]] + if(fusion.name %in% fi.table$id) { + chromosomeA <- fusion.name + chromosomeB <- fusion.name + strandA <- '+' + strandB <- '+' + breakpointA <- fi.table[fusion.name, 'start'] + breakpointB <- fi.table[fusion.name, 'end'] + } else { + warning(sprintf("Could not find location for fusion %s +on virtual contig, ignoring it (line %d)", fusion.name, i+1)) + next + } + } + # PartnerGene objects geneA <- new(Class = "PartnerGene", name = nameA, @@ -168,8 +262,8 @@ importStarfusion <- function (filename, genomeVersion, limit) { geneB = geneB, inframe = inframe, fusionToolSpecificData = fusionToolSpecificData) - } + } #for i - # Return the list of Fusion objects - fusionList -} + ## when useFusionInspector was used, some list elements can be NULL + Filter(function(elt)!is.null(elt), fusionList) +} #importStarFusion diff --git a/R/plotCircle.R b/R/plotCircle.R index fc72d5f..9504cda 100644 --- a/R/plotCircle.R +++ b/R/plotCircle.R @@ -180,18 +180,18 @@ plotCircle <- function(fusionList) { stop("Invalid input. genomeVersion must be either \"hg19\", \"hg38\" or \"mm10\".") } - if (any(rapply(fusionList, function(fusion) fusion@geneA@chromosome == "chrM" || fusion@geneB@chromosome == "chrM"))) { - indexesMitochondrialGenes <- - rapply(fusionList, function(fusion) fusion@geneA@chromosome == "chrM" || fusion@geneB@chromosome == "chrM") - fusionList <- fusionList[!indexesMitochondrialGenes] - message( - paste0( - "Removing ", - length(fusionList[indexesMitochondrialGenes]), - " fusions involving mitochondrial genes as they cannot be plotted in the circle plot." - ) - ) - } +## if (any(rapply(fusionList, function(fusion) fusion@geneA@chromosome == "chrM" || fusion@geneB@chromosome == "chrM"))) { +## indexesMitochondrialGenes <- +## rapply(fusionList, function(fusion) fusion@geneA@chromosome == "chrM" || fusion@geneB@chromosome == "chrM") +## fusionList <- fusionList[!indexesMitochondrialGenes] +## message( +## paste0( +## "Removing ", +## length(fusionList[indexesMitochondrialGenes]), +## " fusions involving mitochondrial genes as they cannot be plotted in the circle plot." +## ) +## ) +## } # Read cytoband information depending on genome version if (fusionList[[1]]@genomeVersion == "hg19") { diff --git a/R/plotFusionReads.R b/R/plotFusionReads.R index dae8983..a870975 100644 --- a/R/plotFusionReads.R +++ b/R/plotFusionReads.R @@ -229,7 +229,7 @@ plotFusionReads <- function( grid::popViewport(1) } -} +} #plotFusionReads .validatePlotFusionReadsParams <- function( fusion, diff --git a/R/plotFusionReadsSimple.R b/R/plotFusionReadsSimple.R new file mode 100644 index 0000000..eac19c5 --- /dev/null +++ b/R/plotFusionReadsSimple.R @@ -0,0 +1,26 @@ +#' Simple versin of plotFusionReads +#' +#' +#' +#' @export +plotFusionReadsSimple <- function(fusion, leftFlank=10000, rightFlank=leftFlank) { + Gviz::displayPars(fusion@fusionReadsAlignment) <- list( + showTitle = FALSE, + showMismatches = FALSE # Show mismatched reads? + ) + + nf <- grid::grid.layout( + nrow = 1, + ncol = 1, + heights = grid::unit(c(6), "null"), + widths = grid::unit(c(1), "null")) + + grid::grid.newpage() + Gviz::plotTracks( + fusion@fusionReadsAlignment, + from = fusion@geneA@breakpoint-leftFlank, + to = fusion@geneB@breakpoint+rightFlank, + chromosome = fusion@fusionReadsAlignment@chromosome, + type = "pileup", + add = TRUE) +} #plotFusionReadsSimple diff --git a/R/utilities.R b/R/utilities.R index 60cd956..5c95e91 100644 --- a/R/utilities.R +++ b/R/utilities.R @@ -701,7 +701,7 @@ getTranscriptsEnsembldb <- function(fusion, edb) { #' fusion <- addFusionReadsAlignment(fusion, bamfile5267) #' #' @export -addFusionReadsAlignment <- function(fusion, bamfile) { +addFusionReadsAlignment <- function(fusion, bamfile, chromosome="chrNA") { # Check if we got a fusion object if (class(fusion) != "Fusion") { @@ -713,7 +713,7 @@ addFusionReadsAlignment <- function(fusion, bamfile) { isPaired = TRUE, # Setting chromosome to chrNA because this is a fusion sequence not found in # any reference genome. - chromosome = "chrNA", + chromosome = chromosome, name="Fusion Reads", genome = fusion@genomeVersion) @@ -847,15 +847,17 @@ selectTranscript <- function( # If the user has chosen one of the four transcript categories, then we want to check whether or not such # transcripts exist. If they exist, simply return them. If they don't exist, go on to try the other categories. Try - # the wanted category first: - if (length(genePartner@transcripts[mcols(genePartner@transcripts)$transcriptCategory == whichTranscripts[[1]] ]) > 0) { - message(paste0("..found transcripts of type ", whichTranscripts[[1]])) + # the wanted category first: + l <- length(genePartner@transcripts[mcols(genePartner@transcripts)$transcriptCategory == whichTranscripts[[1]] ]) + if (l > 0) { + message(paste0("... found ", l, " transcripts of type ", whichTranscripts[[1]])) return(genePartner@transcripts[mcols(genePartner@transcripts)$transcriptCategory == whichTranscripts[[1]] ]) } # Check the remaining categories for(transcriptCategory in transcriptCategories[transcriptCategories != whichTranscripts[[1]]]) { - if (length(genePartner@transcripts[mcols(genePartner@transcripts)$transcriptCategory == transcriptCategory ]) > 0) { - message(paste0("..found transcripts of type ", transcriptCategory)) + l <- length(genePartner@transcripts[mcols(genePartner@transcripts)$transcriptCategory == transcriptCategory ]) + if (l > 0) { + message(paste0(" ... found ",l, " transcripts of type ", transcriptCategory)) return(genePartner@transcripts[mcols(genePartner@transcripts)$transcriptCategory == transcriptCategory ]) } } @@ -1154,13 +1156,19 @@ is.nucleotideAmount.valid <- function(argument_checker, nucleotideAmount, fusion length(fusion@geneA@junctionSequence) + length(fusion@geneB@junctionSequence) + if (fusionJunctionSequenceLength == 0 ) { + ArgumentCheck::addWarning( + msg = "length of junction sequence is 0, have they been determined?", + argcheck = argument_checker + ) + } + if (class(nucleotideAmount) != "numeric" || nucleotideAmount <= 0 || nucleotideAmount > fusionJunctionSequenceLength) { - ArgumentCheck::addError( - msg = paste0("'nucleotideAmount' must be a numeric bigger than or equal ", - "to 0 and less than or equal to the fusion junction ", - "sequence length."), + ArgumentCheck::addWarning( + msg = paste0("'nucleotideAmount' must be a numeric larger > 0 ", + "and <= fusion junction sequence length."), argcheck = argument_checker ) } diff --git a/inst/extdata/UCSC.HG38.Human.CytoBandIdeogram.txt b/inst/extdata/UCSC.HG38.Human.CytoBandIdeogram.txt index 82abd9a..8da32e0 100644 --- a/inst/extdata/UCSC.HG38.Human.CytoBandIdeogram.txt +++ b/inst/extdata/UCSC.HG38.Human.CytoBandIdeogram.txt @@ -859,4 +859,5 @@ chrY 12400000 17100000 q11.221 gpos50 chrY 17100000 19600000 q11.222 gneg chrY 19600000 23800000 q11.223 gpos50 chrY 23800000 26600000 q11.23 gneg -chrY 26600000 57227415 q12 gvar \ No newline at end of file +chrY 26600000 57227415 q12 gvar +chrM 0 16569 m0 gneg diff --git a/man/addFusionReadsAlignment.Rd b/man/addFusionReadsAlignment.Rd index 2855213..2b2b02d 100644 --- a/man/addFusionReadsAlignment.Rd +++ b/man/addFusionReadsAlignment.Rd @@ -4,7 +4,7 @@ \alias{addFusionReadsAlignment} \title{Add fusion reads alignment to fusion object} \usage{ -addFusionReadsAlignment(fusion, bamfile) +addFusionReadsAlignment(fusion, bamfile, chromosome = "chrNA") } \arguments{ \item{fusion}{The fusion object to add a genomic alignment to.} diff --git a/man/importFusionInspector.Rd b/man/importFusionInspector.Rd new file mode 100644 index 0000000..bdbfb99 --- /dev/null +++ b/man/importFusionInspector.Rd @@ -0,0 +1,16 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/importStarfusion.R +\name{importFusionInspector} +\alias{importFusionInspector} +\title{read the results from runnning FusionInspector} +\usage{ + + importFusionInspector(filename = "FusionInspector-inspect/finspector.igv.FusionJuncSpan", + limit = Inf) +} +\description{ +Note: there may not be perfect match between the fusions found by STAR-Fusion +and FusionInspector (since the latter reruns STAR on the mini contig; also, +selfie fusions are ignored). +If no match can be found, the whole fusion is skipped +} diff --git a/man/importStarfusion.Rd b/man/importStarfusion.Rd index 3f225ae..ebaa697 100644 --- a/man/importStarfusion.Rd +++ b/man/importStarfusion.Rd @@ -4,7 +4,7 @@ \alias{importStarfusion} \title{Import results from a STAR-Fusion run into a list of Fusion objects.} \usage{ -importStarfusion(filename, genomeVersion, limit) +importStarfusion(filename, genomeVersion, limit = Inf) } \arguments{ \item{filename}{Filename for the STAR-Fusion diff --git a/man/plotFusionReadsSimple.Rd b/man/plotFusionReadsSimple.Rd new file mode 100644 index 0000000..1fd7df1 --- /dev/null +++ b/man/plotFusionReadsSimple.Rd @@ -0,0 +1,11 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plotFusionReadsSimple.R +\name{plotFusionReadsSimple} +\alias{plotFusionReadsSimple} +\title{Simple versin of plotFusionReads} +\usage{ +plotFusionReadsSimple(fusion) +} +\description{ +Simple versin of plotFusionReads +}