Import fusion inspector support

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",

R/importStarfusion.R

@@ -20,8 +20,7 @@
 #' # This should import a list of 3 fusions described in Fusion objects.
 #'
 #' @export
-importStarfusion <- function (filename, genomeVersion, limit) {
-
+importStarfusion <- function (filename, genomeVersion, limit=Inf) {
   # Is the genome version valid?
   validGenomes <- c("hg19", "hg38", "mm10")
   if (is.na(match(tolower(genomeVersion), tolower(validGenomes)))) {
@@ -39,7 +38,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 +52,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
@@ -96,13 +97,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)
 

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,179 @@ plotFusionReads <- function(
     grid::popViewport(1)
   }
 
-}
+}                                       #plotFusionReads
+
+plotFusionReadsSimple <- function(fusion) {
+###  .validatePlotFusionReadsParams(fusion, showAllNucleotides, nucleotideAmount)
+
+  # Create Biostrings::DNAStringSet of fusion sequence
+  fusionSequence <- Biostrings::DNAStringSet(
+    x = c(fusion@geneA@junctionSequence,
+          fusion@geneB@junctionSequence))
+  names(fusionSequence) <- "chrNA"
+  # And also of a "zoomed-in" version only displaying some nucleotides on each
+  # end of the fusion junction
+  fusionSequenceShort <- Biostrings::DNAStringSet(
+    x = c(
+      substr(
+        fusion@geneA@junctionSequence,
+        nchar(fusion@geneA@junctionSequence)-nucleotideAmount+1,
+        nchar(fusion@geneA@junctionSequence)),
+      substr(
+        fusion@geneB@junctionSequence,
+        nchar(fusion@geneB@junctionSequence)-nucleotideAmount+1 -1,
+        # Added -1 again on the line above because Gviz does not show the last
+        # nucleotide. Ref: http://permalink.gmane.org/gmane.science.biology.informatics.conductor.devel/5726
+        # Or if that link is down, google
+        # "Gviz doesn't provide enough horizontial space"
+        nchar(fusion@geneB@junctionSequence))))
+  names(fusionSequenceShort) <- "chrNA"
+
+  # Load fusion sequences to tracks
+  ref_chimera_short <- Gviz::SequenceTrack(
+    fusionSequenceShort,
+    genome = fusion@genomeVersion,
+    name = "Fusion Sequence")
+  ref_chimera <- Gviz::SequenceTrack(
+    fusionSequence,
+    genome = fusion@genomeVersion,
+    name = "Fusion Sequence")
+
+  # DISPLAY PARAMETERS
+
+  # Set display parameters for the fusion sequence
+  displayParameters <- list(
+    showTitle = FALSE,
+    col = "blue" # The color of the line when no indiviual letters can be plotted due to size limitations.
+  )
+  Gviz::displayPars(ref_chimera) <- displayParameters
+  Gviz::displayPars(ref_chimera_short) <- displayParameters
+
+  # Set display parameters for the fusion reads
+  Gviz::displayPars(fusion@fusionReadsAlignment) <- list(
+    showTitle = FALSE,
+    showMismatches = FALSE # Show mismatched reads?
+  )
+
+  # Do some plotting
+
+  if (showAllNucleotides) {
+    # Create the grid we're gonna use for the plot
+    nf <- grid::grid.layout(
+      nrow = 1,
+      ncol = 1)
+
+    # Open a new page to plot it
+    grid::grid.newpage()
+
+    # Divide the page according to the grid we created
+    grid::pushViewport(grid::viewport(layout = nf))
+
+    # Open row 1, column 1
+    grid::pushViewport(grid::viewport(layout.pos.row = 1, layout.pos.col = 1))
+    # Plot the fusion read alignment
+    Gviz::plotTracks(
+      c(fusion@fusionReadsAlignment, ref_chimera),
+      from = 1,
+      to = nchar(fusionSequence),
+      chromosome = fusion@fusionReadsAlignment@chromosome,
+      type = "pileup",
+      sizes = c(6, 1),
+      add = TRUE
+    )
+    # Calculate where in the alignment the breakpoint is
+    breakpoint <- nchar(fusion@geneA@junctionSequence)/nchar(fusionSequence)
+    # Plot line indicating fusion breakpoint
+    grid::grid.lines(x = grid::unit(breakpoint, "npc"),
+               y = grid::unit(c(0, 1), "npc"),
+               default.units = "npc",
+               arrow = NULL,
+               name = NULL,
+               gp = grid::gpar(
+                 col = "black",
+                 lty = "dotted"),
+               draw = TRUE,
+               vp = NULL)
+    # Close row 1, column 1
+    grid::popViewport(1)
+  } else {
+    # Create the grid we're gonna use for the plot
+    nf <- grid::grid.layout(
+      nrow = 2,
+      ncol = 1,
+      heights = grid::unit(c(6, 1), "null"),
+      widths = grid::unit(c(1, 1), "null"))
+
+    # Open a new page to plot it
+    grid::grid.newpage()
+
+    # Divide the page according to the grid we created
+    grid::pushViewport(grid::viewport(layout = nf))
+
+    # Open row 1, column 1
+    grid::pushViewport(grid::viewport(layout.pos.row = 1, layout.pos.col = 1))
+    # Plot the fusion read alignment
+    Gviz::plotTracks(
+      fusion@fusionReadsAlignment,
+      from = 1,
+      to = nchar(fusionSequence),
+      chromosome = fusion@fusionReadsAlignment@chromosome,
+      type = "pileup",
+      add = TRUE
+    )
+    # Calculate where in the alignment the breakpoint is
+    breakpoint <- nchar(fusion@geneA@junctionSequence)/nchar(fusionSequence)
+    # Plot line indicating fusion breakpoint
+    grid::grid.lines(x = grid::unit(breakpoint, "npc"),
+               y = grid::unit(c(0, 1), "npc"),
+               default.units = "npc",
+               arrow = NULL,
+               name = NULL,
+               gp = grid::gpar(
+                 col = "black",
+                 lty = "dotted"),
+               draw = TRUE,
+               vp = NULL)
+    # Close row 1, column 1
+    grid::popViewport(1)
+
+    # Open row 2, column 1
+    grid::pushViewport(grid::viewport(layout.pos.row = 2, layout.pos.col = 1))
+    # Plot a part of the fusion junction sequence
+    Gviz::plotTracks(
+      ref_chimera_short,
+      from = 1,
+      to = nchar(fusionSequenceShort),
+      chromosome = fusion@fusionReadsAlignment@chromosome,
+      add = TRUE
+    )
+    # Plot line indicating fusion breakpoint
+    grid::grid.lines(x = grid::unit(0.5, "npc"),
+               y = grid::unit(c(0, 1), "npc"),
+               default.units = "npc",
+               arrow = NULL,
+               name = NULL,
+               gp = grid::gpar(
+                 col = "black",
+                 lty = "dotted"),
+               draw = TRUE,
+               vp = NULL)
+    # Plot line between lines
+    grid::grid.lines(x = grid::unit(c(0.5, breakpoint), "npc"),
+               y = grid::unit(1, "npc"),
+               default.units = "npc",
+               arrow = NULL,
+               name = NULL,
+               gp = grid::gpar(
+                 col = "black",
+                 lty = "dotted"),
+               draw = TRUE,
+               vp = NULL)
+    # Close row 1, column 2
+    grid::popViewport(1)
+  }
+
+}                                       #plotFusionReadsSimple
 
 .validatePlotFusionReadsParams <- function(
   fusion,

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