Import fusion inspector

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 <stianlagstad@gmail.com>
 Authors@R: c(
   person(
     given = "Stian",

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)

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,13 +153,25 @@ 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"]]
     fusionToolSpecificData[["RightBreakDinuc"]] = report[[i, "RightBreakDinuc"]]
     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

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") {

R/plotFusionReads.R

@@ -229,7 +229,7 @@ plotFusionReads <- function(
     grid::popViewport(1)
   }
 
-}
+}                                       #plotFusionReads
 
 .validatePlotFusionReadsParams <- function(
   fusion,

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

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
     )
   }

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
+chrY  26600000  57227415  q12 gvar
+chrM  0  16569 m0  gneg