Closed komalsrathi closed 6 years ago
I've reproduced the issue with your example and tracked down what code in chimeraviz causes the issue (the line edited in this commit). It looks like there's simply too many reads aligning to the fusion sequence so that Gviz in the end gets problems displaying them all when I've set min.height = 20
, leading to the error. When I change that line to min.height = 10
it seems to work fine. I've pushed the change, so let me know if that solves this problem for you.
I am unable to get the dev version using:
devtools::install_git("git@git.bioconductor.org:packages/chimeraviz")
My Rstudio just hangs and does not proceed.
Could you try again? I noticed chimeraviz was missing from the latest bioconductor build yesterday, but it seems to be back now. And I just tried installing with that command myself, and it worked.
I can't seem to understand what's going on -- when I try that command, it just hangs and doesn't do anything. Did I screw things up by upgrading Bioconductor to 3.5? Seems like everything was working fine before that :/
Hi Stian
The only way I am able to install chimeraviz
is:
source("https://bioconductor.org/biocLite.R")
biocLite("chimeraviz")
This is version chimeraviz_1.0.4
. Is this the correct one?
I'm sorry, I thought that would work. I'll get back to you on how you can install chimeraviz directly from the Bioconductor git repository. In the meantime you have two more common ways of installing the package:
I wanted you to install directly from the Bioconductor git repository so that you would get the very latest version, as there is a lag between me submitting changes and them appearing in the development version. But right now it doesn't matter since the fix for this Github issue has been synced to the devel version of chimeraviz. So go ahead and try to install that, or wait until the next Bioconductor release next week :)
This should work: devtools::install_git("https://git.bioconductor.org/packages/chimeraviz.git")
Thanks!! Closing this now. I am now testing chimeraviz_1.2.3
and will open new issues if any. Thanks again!!
Hi Stian,
I would like to reopen this again. This time trying with a different sample but getting the same error.
data is here: https://drive.google.com/open?id=0B-8gQV1WZcYdWFBVT0FlbjRWY2c
Script:
devtools::install_git("https://git.bioconductor.org/packages/chimeraviz.git")
library(chimeraviz)
# Get results file from fusioncatcher
fc <- importFusioncatcher(filename = 'sample_out/final-list_candidate-fusion-genes.GRCh37.txt', "hg19", 25)
# Load the results file into a list of fusion objects
fusion <- getFusionByGeneName(fc, geneName = 'RCC1')
# first look at 14 and then at 15
fusion <- getFusionById(fusion, id = '14')
# information about upstream and downstream partners
upstreamPartnerGene(fusion)
downstreamPartnerGene(fusion)
#############################################
# Fusion Reads Plot
#############################################
# read fastq files with reads for that fusion
fastq1 <- 'sample_R1.fq'
fastq2 <- 'sample_R2.fq'
# extract the fusion junction sequence
referenceFilename <- "reference.fa"
writeFusionReference(fusion = fusion, filename = referenceFilename)
# map fusion reads to the fusion junction sequence
# align the reads you have against the fusion junction sequence
source(system.file(
"scripts",
"rsubread.R",
package="chimeraviz"))
# create an index for our reference sequence, the fusion junction seq
rsubreadIndex(referenceFasta = referenceFilename)
# align the reads against the reference
# this will create fusionAlignment.bam file
rsubreadAlign(
referenceName = referenceFilename,
fastq1 = fastq1,
fastq2 = fastq2,
outputBamFilename = "fusionAlignment")
# this bam file has reads supporting the fusion event
bamfile <- 'fusionAlignment.bam'
# Add bamfile of fusion reads to the fusion oject
fusion <- addFusionReadsAlignment(fusion, bamfile)
plotFusionReads(fusion = fusion, showAllNucleotides = FALSE, nucleotideAmount = 15)
Error in data.frame(x1 = start(pairGaps) - 1, y1 = gy, x2 = end(pairGaps) + :
arguments imply differing number of rows: 0, 1
There are two fusions for the same genes, one is in-frame and other is out-of-frame. The code works with id: 15 but not with id: 14.
This is the session info:
> sessionInfo()
R version 3.4.2 (2017-09-28)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: OS X El Capitan 10.11.6
Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8
attached base packages:
[1] grid stats4 parallel stats graphics grDevices utils datasets methods base
other attached packages:
[1] biomaRt_2.32.1 chimeraviz_1.2.3 ensembldb_2.0.4 AnnotationFilter_1.0.0
[5] GenomicFeatures_1.28.5 AnnotationDbi_1.38.2 Biobase_2.36.2 Gviz_1.20.0
[9] GenomicRanges_1.28.6 GenomeInfoDb_1.12.3 Biostrings_2.44.2 XVector_0.16.0
[13] IRanges_2.10.5 S4Vectors_0.14.7 BiocGenerics_0.22.1 BiocInstaller_1.26.1
loaded via a namespace (and not attached):
[1] ProtGenerics_1.8.0 bitops_1.0-6 matrixStats_0.52.2
[4] devtools_1.13.3 bit64_0.9-7 RColorBrewer_1.1-2
[7] httr_1.3.1 rprojroot_1.2 Rgraphviz_2.20.0
[10] tools_3.4.2 backports_1.1.1 DT_0.2
[13] R6_2.2.2 rpart_4.1-11 Hmisc_4.0-3
[16] DBI_0.7 lazyeval_0.2.0 colorspace_1.3-2
[19] nnet_7.3-12 withr_2.0.0 gridExtra_2.3
[22] bit_1.1-12 curl_3.0 compiler_3.4.2
[25] git2r_0.19.0 graph_1.54.0 htmlTable_1.9
[28] DelayedArray_0.2.7 rtracklayer_1.36.6 scales_0.5.0.9000
[31] checkmate_1.8.5 readr_1.1.1 RCircos_1.2.0
[34] stringr_1.2.0 digest_0.6.12 Rsamtools_1.28.0
[37] foreign_0.8-69 rmarkdown_1.6 Rsubread_1.26.1
[40] base64enc_0.1-3 dichromat_2.0-0 pkgconfig_2.0.1
[43] htmltools_0.3.6 BSgenome_1.44.2 htmlwidgets_0.9
[46] rlang_0.1.2 RSQLite_2.0 shiny_1.0.5
[49] BiocParallel_1.10.1 acepack_1.4.1 VariantAnnotation_1.22.3
[52] RCurl_1.95-4.8 magrittr_1.5 GenomeInfoDbData_0.99.0
[55] Formula_1.2-2 Matrix_1.2-11 Rcpp_0.12.13
[58] munsell_0.4.3 stringi_1.1.5 yaml_2.1.14
[61] SummarizedExperiment_1.6.5 zlibbioc_1.22.0 org.Hs.eg.db_3.4.1
[64] plyr_1.8.4 AnnotationHub_2.8.3 blob_1.1.0
[67] lattice_0.20-35 splines_3.4.2 hms_0.3
[70] knitr_1.17 reshape2_1.4.2 XML_3.98-1.9
[73] evaluate_0.10.1 biovizBase_1.24.0 latticeExtra_0.6-28
[76] data.table_1.10.4-2 httpuv_1.3.5 gtable_0.2.0
[79] ggplot2_2.2.1.9000 mime_0.5 xtable_1.8-2
[82] ArgumentCheck_0.10.2 survival_2.41-3 tibble_1.3.4
[85] GenomicAlignments_1.12.2 memoise_1.1.0 cluster_2.0.6
[88] interactiveDisplayBase_1.14.0 BiocStyle_2.4.1
Thank you again @komalsrathi! I'm just writing to confirm that I get the same issue on my end. I'll try to fix this as soon as I can.
This should be taken care of in version 1.2.6. Let me know if it works :)
Hi Stian - I tried this with v1.2.6 but it is still throwing the same error.
R version 3.4.2 (2017-09-28)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: OS X El Capitan 10.11.6
Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8
attached base packages:
[1] grid stats4 parallel stats graphics grDevices utils
[8] datasets methods base
other attached packages:
[1] chimeraviz_1.2.6 ensembldb_2.0.4 AnnotationFilter_1.0.0
[4] GenomicFeatures_1.28.5 AnnotationDbi_1.38.2 Biobase_2.36.2
[7] Gviz_1.20.0 GenomicRanges_1.28.6 GenomeInfoDb_1.12.3
[10] Biostrings_2.44.2 XVector_0.16.0 IRanges_2.10.5
[13] S4Vectors_0.14.7 BiocGenerics_0.22.1 BiocInstaller_1.26.1
loaded via a namespace (and not attached):
[1] ProtGenerics_1.8.0 bitops_1.0-6
[3] matrixStats_0.52.2 devtools_1.13.3
[5] bit64_0.9-7 RColorBrewer_1.1-2
[7] httr_1.3.1 rprojroot_1.2
[9] tools_3.4.2 backports_1.1.1
[11] DT_0.2 R6_2.2.2
[13] rpart_4.1-11 Hmisc_4.0-3
[15] DBI_0.7 lazyeval_0.2.0
[17] colorspace_1.3-2 nnet_7.3-12
[19] withr_2.0.0 gridExtra_2.3
[21] bit_1.1-12 curl_3.0
[23] compiler_3.4.2 git2r_0.19.0
[25] htmlTable_1.9 DelayedArray_0.2.7
[27] rtracklayer_1.36.6 scales_0.5.0.9000
[29] checkmate_1.8.5 readr_1.1.1
[31] RCircos_1.2.0 stringr_1.2.0
[33] digest_0.6.12 Rsamtools_1.28.0
[35] foreign_0.8-69 rmarkdown_1.6
[37] Rsubread_1.26.1 base64enc_0.1-3
[39] dichromat_2.0-0 pkgconfig_2.0.1
[41] htmltools_0.3.6 BSgenome_1.44.2
[43] htmlwidgets_0.9 rlang_0.1.2
[45] RSQLite_2.0 shiny_1.0.5
[47] BiocParallel_1.10.1 acepack_1.4.1
[49] VariantAnnotation_1.22.3 RCurl_1.95-4.8
[51] magrittr_1.5 GenomeInfoDbData_0.99.0
[53] Formula_1.2-2 Matrix_1.2-11
[55] Rcpp_0.12.13 munsell_0.4.3
[57] stringi_1.1.5 yaml_2.1.14
[59] SummarizedExperiment_1.6.5 zlibbioc_1.22.0
[61] org.Hs.eg.db_3.4.1 plyr_1.8.4
[63] AnnotationHub_2.8.3 blob_1.1.0
[65] lattice_0.20-35 splines_3.4.2
[67] hms_0.3 knitr_1.17
[69] biomaRt_2.32.1 XML_3.98-1.9
[71] evaluate_0.10.1 biovizBase_1.24.0
[73] latticeExtra_0.6-28 data.table_1.10.4-2
[75] httpuv_1.3.5 gtable_0.2.0
[77] ggplot2_2.2.1.9000 mime_0.5
[79] xtable_1.8-2 ArgumentCheck_0.10.2
[81] survival_2.41-3 tibble_1.3.4
[83] GenomicAlignments_1.12.2 memoise_1.1.0
[85] cluster_2.0.6 interactiveDisplayBase_1.14.0
[87] BiocStyle_2.4.1
Hmm, that's weird. This code:
setwd("/home/stian/dev/chimeraviz_3")
devtools::install_git("https://git.bioconductor.org/packages/chimeraviz.git")
# Get results file from fusioncatcher
fc <- importFusioncatcher(filename = 'SRR1559136_out/final-list_candidate-fusion-genes.GRCh37.txt', "hg19", 25)
# Load the results file into a list of fusion objects
fusion <- getFusionByGeneName(fc, geneName = 'LMO1')
# first look at 14 and then at 15
fusion <- getFusionById(fusion, id = '14')
# information about upstream and downstream partners
upstreamPartnerGene(fusion)
downstreamPartnerGene(fusion)
#############################################
# Fusion Reads Plot
#############################################
# read fastq files with reads for that fusion
fastq1 <- 'SRR1559136_fusion_R1.fq'
fastq2 <- 'SRR1559136_fusion_R2.fq'
# extract the fusion junction sequence
referenceFilename <- "reference.fa"
writeFusionReference(fusion = fusion, filename = referenceFilename)
# map fusion reads to the fusion junction sequence
# align the reads you have against the fusion junction sequence
source(system.file(
"scripts",
"rsubread.R",
package="chimeraviz"))
# create an index for our reference sequence, the fusion junction seq
rsubreadIndex(referenceFasta = referenceFilename)
# align the reads against the reference
# this will create fusionAlignment.bam file
rsubreadAlign(
referenceName = referenceFilename,
fastq1 = fastq1,
fastq2 = fastq2,
outputBamFilename = "fusionAlignment")
# this bam file has reads supporting the fusion event
bamfile <- 'fusionAlignment.bam'
# Add bamfile of fusion reads to the fusion oject
fusion <- addFusionReadsAlignment(fusion, bamfile)
plotFusionReads(fusion = fusion, showAllNucleotides = FALSE, nucleotideAmount = 15)
Creates this plot for me:
My sessionInfo:
> sessionInfo()
R version 3.4.2 (2017-09-28)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.3 LTS
Matrix products: default
BLAS: /usr/lib/libblas/libblas.so.3.6.0
LAPACK: /usr/lib/lapack/liblapack.so.3.6.0
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=nb_NO.UTF-8 LC_TIME=nb_NO.UTF-8
[4] LC_COLLATE=en_US.UTF-8 LC_MONETARY=nb_NO.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=nb_NO.UTF-8 LC_NAME=C LC_ADDRESS=C
[10] LC_TELEPHONE=C LC_MEASUREMENT=nb_NO.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] grid stats4 parallel stats graphics grDevices utils datasets methods
[10] base
other attached packages:
[1] chimeraviz_1.2.6 ensembldb_2.0.4 AnnotationFilter_1.0.0 GenomicFeatures_1.28.5
[5] AnnotationDbi_1.38.2 Biobase_2.36.2 Gviz_1.20.0 GenomicRanges_1.28.6
[9] GenomeInfoDb_1.12.3 Biostrings_2.44.2 XVector_0.16.0 IRanges_2.10.5
[13] S4Vectors_0.14.7 BiocGenerics_0.22.1
loaded via a namespace (and not attached):
[1] colorspace_1.3-2 rprojroot_1.2 biovizBase_1.24.0
[4] htmlTable_1.9 ArgumentCheck_0.10.2 base64enc_0.1-3
[7] dichromat_2.0-0 rstudioapi_0.7 roxygen2_6.0.1
[10] DT_0.2 bit64_0.9-7 interactiveDisplayBase_1.14.0
[13] xml2_1.1.1 splines_3.4.2 knitr_1.17
[16] Formula_1.2-2 Rsamtools_1.28.0 cluster_2.0.6
[19] graph_1.54.0 shiny_1.0.5 readr_1.1.1
[22] compiler_3.4.2 httr_1.3.1 backports_1.1.1
[25] assertthat_0.2.0 Matrix_1.2-11 lazyeval_0.2.0
[28] acepack_1.4.1 htmltools_0.3.6 tools_3.4.2
[31] bindrcpp_0.2 gtable_0.2.0 glue_1.1.1
[34] GenomeInfoDbData_0.99.0 dplyr_0.7.4 Rcpp_0.12.13
[37] rtracklayer_1.36.5 stringr_1.2.0 testthat_1.0.2
[40] mime_0.5 devtools_1.13.3 XML_3.98-1.9
[43] AnnotationHub_2.8.2 org.Hs.eg.db_3.4.1 zlibbioc_1.22.0
[46] RCircos_1.2.0 scales_0.5.0 BiocStyle_2.4.1
[49] BSgenome_1.44.2 VariantAnnotation_1.22.3 BiocInstaller_1.26.1
[52] hms_0.3 ProtGenerics_1.8.0 SummarizedExperiment_1.6.5
[55] RColorBrewer_1.1-2 yaml_2.1.14 curl_3.0
[58] memoise_1.1.0 gridExtra_2.3 ggplot2_2.2.1
[61] biomaRt_2.32.1 rpart_4.1-11 latticeExtra_0.6-28
[64] stringi_1.1.5 RSQLite_2.0 checkmate_1.8.4
[67] BiocParallel_1.10.1 Rsubread_1.26.1 rlang_0.1.2
[70] pkgconfig_2.0.1 commonmark_1.4 matrixStats_0.52.2
[73] bitops_1.0-6 evaluate_0.10.1 lattice_0.20-35
[76] bindr_0.1 GenomicAlignments_1.12.2 htmlwidgets_0.9
[79] bit_1.1-12 plyr_1.8.4 magrittr_1.5
[82] R6_2.2.2 Hmisc_4.0-3 DelayedArray_0.2.7
[85] DBI_0.7 foreign_0.8-69 withr_2.0.0
[88] survival_2.41-3 RCurl_1.95-4.8 nnet_7.3-12
[91] tibble_1.3.4 crayon_1.3.4 rmarkdown_1.6
[94] data.table_1.10.4-2 blob_1.1.0 git2r_0.19.0
[97] Rgraphviz_2.20.0 digest_0.6.12 xtable_1.8-2
[100] httpuv_1.3.5 munsell_0.4.3
Could it be that something is cached on your machine? Just to exclude that as an option, could you delete your chimeraviz installation, exit R, start a new session and install chimeraviz again? (Running .libPaths()
will tell you where packages are installed on your system. Find the chimeraviz folder and delete it.)
Hi Stian,
I deleted chimeraviz, restarted R and installed chimeraviz again - but I am getting the same error. Is there anything else I can provide you so you can check what's going on?
This is the plotFusionReads
function for the version that I have installed (just in case it helps):
trace(plotFusionReads, edit = TRUE)
function (fusion, showAllNucleotides = TRUE, nucleotideAmount = 10)
{
.validatePlotFusionReadsParams(fusion, showAllNucleotides,
nucleotideAmount)
fusionSequence <- Biostrings::DNAStringSet(x = c(fusion@geneA@junctionSequence,
fusion@geneB@junctionSequence))
names(fusionSequence) <- "chrNA"
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, nchar(fusion@geneB@junctionSequence))))
names(fusionSequenceShort) <- "chrNA"
ref_chimera_short <- Gviz::SequenceTrack(fusionSequenceShort,
genome = fusion@genomeVersion, name = "Fusion Sequence")
ref_chimera <- Gviz::SequenceTrack(fusionSequence, genome = fusion@genomeVersion,
name = "Fusion Sequence")
displayParameters <- list(showTitle = FALSE, col = "blue")
Gviz::displayPars(ref_chimera) <- displayParameters
Gviz::displayPars(ref_chimera_short) <- displayParameters
Gviz::displayPars(fusion@fusionReadsAlignment) <- list(showTitle = FALSE,
showMismatches = FALSE)
if (showAllNucleotides) {
nf <- grid::grid.layout(nrow = 1, ncol = 1)
grid::grid.newpage()
grid::pushViewport(grid::viewport(layout = nf))
grid::pushViewport(grid::viewport(layout.pos.row = 1,
layout.pos.col = 1))
Gviz::plotTracks(c(fusion@fusionReadsAlignment, ref_chimera),
from = 1, to = nchar(fusionSequence), chromosome = fusion@fusionReadsAlignment@chromosome,
type = "pileup", sizes = c(6, 1), add = TRUE)
breakpoint <- nchar(fusion@geneA@junctionSequence)/nchar(fusionSequence)
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)
grid::popViewport(1)
}
else {
nf <- grid::grid.layout(nrow = 2, ncol = 1, heights = grid::unit(c(6,
1), "null"), widths = grid::unit(c(1, 1), "null"))
grid::grid.newpage()
grid::pushViewport(grid::viewport(layout = nf))
grid::pushViewport(grid::viewport(layout.pos.row = 1,
layout.pos.col = 1))
Gviz::plotTracks(fusion@fusionReadsAlignment, from = 1,
to = nchar(fusionSequence), chromosome = fusion@fusionReadsAlignment@chromosome,
type = "pileup", add = TRUE)
breakpoint <- nchar(fusion@geneA@junctionSequence)/nchar(fusionSequence)
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)
grid::popViewport(1)
grid::pushViewport(grid::viewport(layout.pos.row = 2,
layout.pos.col = 1))
Gviz::plotTracks(ref_chimera_short, from = 1, to = nchar(fusionSequenceShort),
chromosome = fusion@fusionReadsAlignment@chromosome,
add = TRUE)
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)
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)
grid::popViewport(1)
}
}
Thanks Komal
That looks like the most recent version of the function. Maybe this is another mac-specific issue, because I cannot reproduce it on my Ubuntu laptop. I'll check it on the mac I have at work as soon as I can.
Could you try doing this?
# Calculate image size based on supporting reads and lenght of junction
# sequence.
imageWidth <- (nchar(partnerGeneJunctionSequence(upstreamPartnerGene(fusion))) +
nchar(partnerGeneJunctionSequence(downstreamPartnerGene(fusion)))) * 15
imageHeight <- (fusionSplitReadsCount(fusion)+fusionSpanningReadsCount(fusion)) * 20
# Open device
png("plot.png", width = imageWidth, height = imageHeight)
# Now we can plot
plotFusionReads(fusion)
# Close device
dev.off()
That is, instead of using the standard graphics device, store the plot in a png file with a predefined height and width. The issue might have to do with the available space in the plot.
Hi Stian,
That did not work but I just took your idea of using a different plotting device and used pdf instead:
fusion <- addFusionReadsAlignment(fusion, bamfile)
pdf(file = pdfFileName, width = 6, height = 8)
plotFusionReads(fusion = fusion, showAllNucleotides = FALSE, nucleotideAmount = 15)
dev.off()
If you get a chance to test this on a mac then please do. But I am closing this now. Thanks a lot!
Hi Stian,
Tried pdf and your solutions but still facing this problem in a bunch of samples - will update this comment in a few hours with some examples.
Thanks Komal
Thank you! Examples along with the data that leads to the errors makes my job a lot easier :)
Hi Stian,
Here is the Rscript and associated data: https://drive.google.com/drive/folders/0B-8gQV1WZcYdWFBVT0FlbjRWY2c
Thanks!!
As far as I can understand there is something going wrong in the gviz package. I'll create an issue for gviz to see if Florian Hahne (the maintainer) has any tips. Is it alright if I use the files you provided in the issue I create for gviz?
Steps to reproduce the error:
library(chimeraviz)
fc <- importFusioncatcher(filename = 'SRR1559032_out/final-list_candidate-fusion-genes.GRCh37.txt', "hg19", 100)
# Load the results file into a list of fusion objects
fusion <- getFusionById(fc, id = '44')
fusion
# read fastq files with reads for that fusion
fastq1 <- 'SRR1559032_fusion_R1.fq'
fastq2 <- 'SRR1559032_fusion_R2.fq'
# extract the fusion junction sequence
referenceFilename <- "reference.fa"
writeFusionReference(fusion = fusion, filename = referenceFilename)
# map fusion reads to the fusion junction sequence
# align the reads you have against the fusion junction sequence
source(system.file(
"scripts",
"rsubread.R",
package="chimeraviz"))
# create an index for our reference sequence, the fusion junction seq
rsubreadIndex(referenceFasta = referenceFilename)
# align the reads against the reference
# this will create fusionAlignment.bam file
rsubreadAlign(
referenceName = referenceFilename,
fastq1 = fastq1,
fastq2 = fastq2,
outputBamFilename = "fusionAlignment")
bamfile <- 'fusionAlignment.bam'
# Create Biostrings::DNAStringSet of fusion sequence
fusionSequence <- Biostrings::DNAStringSet(
x = c(fusion@geneA@junctionSequence,
fusion@geneB@junctionSequence))
names(fusionSequence) <- "chrNA"
fusionReadsAlignment <- Gviz::AlignmentsTrack(
bamfile,
isPaired = TRUE,
# Setting chromosome to chrNA because this is a fusion sequence not found in
# any reference genome.
chromosome = "chrNA",
name="Fusion Reads",
genome = fusion@genomeVersion)
# This call leads to the error:
Gviz::plotTracks(
fusionReadsAlignment,
chromosome = "chrNA",
from = 1,
to = nchar(fusionSequence))
# Note that if I set isPaired to FALSE instead, then the plot works (but this is not what we want to do)
fusionReadsAlignment2 <- Gviz::AlignmentsTrack(
bamfile,
isPaired = FALSE,
# Setting chromosome to chrNA because this is a fusion sequence not found in
# any reference genome.
chromosome = "chrNA",
name="Fusion Reads",
genome = fusion@genomeVersion)
Gviz::plotTracks(
fusionReadsAlignment2,
chromosome = "chrNA",
from = 1,
to = nchar(fusionSequence))
Hi Stian,
Sure - but could you just share the data and Rscript in the same google drive folder instead of sending the files over? So that once this is resolved, I will delete the google drive folder as well. Let me know if this works for you.
Thanks Komal
Thanks! I've posted the issue here: https://support.bioconductor.org/p/103474/. Please upvote it if you're registered there.
Hi Stian, not to rush you in any way but any updates on this?
Thanks for the heads up, @komalsrathi! I hadn't noticed their answer. If the bug has indeed been fixed, you should only need to update your Bioconductor packages. Take a look here. I'll take a look myself as soon as I'm able :)
@komalsrathi, did you get past this error by updating your Bioconductor packages?
Hi Stian - sorry I just got around to testing this. It is working fine now. Thank you!!
Hi Stian,
I am facing a different issue with another sample:
The
plotFusionReads
function is giving me this error:I am attaching all required data and script here: https://drive.google.com/open?id=0B-8gQV1WZcYdNnJyWmROclI0UFE
Please look at this whenever you get some time.
Thanks for all your help!