phanstiellab / sushi Goto Github PK

this is my code

bdata <-read.table("gene.bed",header=T)
head(bdata)
chrom            = "SM_V7_1"
chromstart       = 33135620
chromend         = 33143037
plotGenes(bdata, chrom=chrom,chromstart = chromstart,chromend = chromend, types=bdata$type, maxrows=50,bheight=0.08,plotgenetype="box",bentline=FALSE,col="brown", labeloffset=.2,fontsize=0.9,arrowlength = 0.025,labeltext=TRUE)
labelgenome( chrom, chromstart,chromend,side=1,scipen=20,n=3,scale="Mb",line=.18,chromline=.5,scaleline=0.5)

Run this code in jutyter notebook. Only labelgenome() function work

chrom | start | stop | gene | score | strand | type
-- | -- | -- | -- | -- | -- | --
SM_V7_1 | 33135620 | 33135696 | Smp_012860.1 | . | 1 | CDS
SM_V7_1 | 33135730 | 33135847 | Smp_012860.1 | . | 1 | CDS
SM_V7_1 | 33135886 | 33136029 | Smp_012860.1 | . | 1 | CDS
SM_V7_1 | 33136062 | 33136237 | Smp_012860.1 | . | 1 | CDS
SM_V7_1 | 33137617 | 33137839 | Smp_012860.1 | . | 1 | CDS
SM_V7_1 | 33139243 | 33139362 | Smp_012860.1 | . | 1 | CDS
SM_V7_1 | 33140579 | 33140630 | Smp_012860.1 | . | 1 | CDS
SM_V7_1 | 33141423 | 33141616 | Smp_012860.1 | . | 1 | CDS
SM_V7_1 | 33142706 | 33142927 | Smp_012860.1 | . | 1 | CDS
SM_V7_1 | 33134853 | 33135619 | Smp_012860.1 | . | 1 | utr
SM_V7_1 | 33142926 | 33143037 | Smp_012860.1 | . | 1 | utr
SM_V7_1 | 33134851 | 33135030 | Smp_012860.2 | . | 1 | CDS
SM_V7_1 | 33136062 | 33136237 | Smp_012860.2 | . | 1 | CDS
SM_V7_1 | 33137617 | 33137839 | Smp_012860.2 | . | 1 | CDS
SM_V7_1 | 33139243 | 33139362 | Smp_012860.2 | . | 1 | CDS
SM_V7_1 | 33140579 | 33140630 | Smp_012860.2 | . | 1 | CDS
SM_V7_1 | 33141423 | 33141616 | Smp_012860.2 | . | 1 | CDS
SM_V7_1 | 33142706 | 33142927 | Smp_012860.2 | . | 1 | CDS
SM_V7_1 | 33134803 | 33134850 | Smp_012860.2 | . | 1 | utr
SM_V7_1 | 33142926 | 33143037 | Smp_012860.2 | . | 1 | utr
SM_V7_1 | 33135620 | 33135696 | g789.t1 | . | 1 | CDS
SM_V7_1 | 33135730 | 33135847 | g789.t1 | . | 1 | CDS
SM_V7_1 | 33135886 | 33136029 | g789.t1 | . | 1 | CDS
SM_V7_1 | 33136062 | 33136237 | g789.t1 | . | 1 | CDS
SM_V7_1 | 33137617 | 33137839 | g789.t1 | . | 1 | CDS
SM_V7_1 | 33139243 | 33139362 | g789.t1 | . | 1 | CDS
SM_V7_1 | 33140579 | 33140630 | g789.t1 | . | 1 | CDS
SM_V7_1 | 33141423 | 33141616 | g789.t1 | . | 1 | CDS
SM_V7_1 | 33142706 | 33142927 | g789.t1 | . | 1 | CDS
SM_V7_1 | 33134853 | 33135619 | g789.t1 | . | 1 | utr
SM_V7_1 | 33142926 | 33143037 | g789.t1 | . | 1 | utr

Hi,

Is there a way to select the colors when using plotBedpe function for the loops?

I am able to generate a plot with different colors for the different samples.. However, I don't have any control over the colors that are chosen. I tried different palettes (topo.colors, terrain.colors, heat.colors, cm.colors and rainbow) but am not happy with the choices.. When I tried to specify a vector of colors for the colorbycol, it says it expects a matrix and not a vector for col. Please let me know if the feature already exists to choose our own colors! If not, I think it will be a nice addition!

pbde = plotBedpe(inputFile, chrom, chromstart, chromend, heights=inputFile$cc, plottype="loops", colorby=inputFile$Sample, colorbycol=heat.colors, lwd=1, bty='n', border="black")

Thanks!

figure out how to write vignettes

Dear Author,

Thank you for developing this package.
It is incredible to visualise some data.
However, I was wondering if it is possible to plot refseq genes below the x coordinates of the
plotBedpe plot. I was thinking if calling layout() function would be only possible to achieve that. Do you know if there is a simpler way to do the same?

Thank you for your time,

Vidal

Dear authors,

Thanks for developing this package.
Recently I asked about plotGenes with below plotBedpe to help visualising the data,
however, I keep receving the following error:

[1] "yes"
Error in FUN(X[[i]], ...) : 
  only defined on a data frame with all numeric variables

These are my specifications:

> chrom <- "chr9"
> chromstart <- 18377770
> chromend <- 20274437
> pg=plotGenes(RNAseq_in_position_chr9.bed,chrom,chromstart,chromend,
+              types=RNAseq_in_position_chr9.bed$type,
+              colorby=RNAseq_in_position_chr9.bed$score+2,
+              colorbycol= SushiColors2(7),colorbyrange=c(0,1),
+              labeltext=TRUE,maxrows=50,height=0.4,plotgenetype="box")

head(RNAseq_in_position_chr9.bed)

 chrom     start      stop  gene      score strand type
1    chr9  99403532  99417599 AAED1  0.4323268    -1 exon
164  chr9 135600964 135754198   AK8 -1.8297926    -1 exon
165  chr9 135600964 135754198   AK8 -1.8297926    -1 exon
166  chr9 135600964 135754198   AK8 -1.8297926    -1 exon
193  chr9 101978706 101984246  ALG2  0.3562529    -1 exon
194  chr9 101978706 101984246  ALG2  0.3562529    1 exon

Would you know why I'm having this error?

Thank you in advance

Hello @dphansti

I have questions with respect plotBed function.
I have been trying to do a graphic's set, like the facet function, but i can't. is possible to do it with plotBed?
Moreover, I have to Bed archive with only 2 datas to graph, and i can't to do a graphic with it, which is the minimum number of datas i need to use to graph ?

I will grateful for your response
Danitza

Hi @dphansti

Thanks again for developing this fantastic tool. It is so useful to visualize my data. However, it seems the plotbedpe function can only plot the loops with both ends located within the defined genome range. I'm wondering how to draw part of the loop when there's only one end within the range and the other one outside the range.

Many thanks in advance!
Fan

I'm trying to plot a bedfile with chromstart = 4077 and chromend = 2153224, but it never completes. It seems to hang somewhere chromend = 440,000 and 450,000. values lower than this take around 1 min on my pc. I left it overnight running the full amount but it never finished.

Is there any way around this?

Thanks!

We you want to zoom into a gene at exon level you don't see the exon anymore. As shown by the following figure:

I believe this has to do with a transcriptinfo filtering step in plotGenes:

transcriptinfo = transcriptinfo[which(
    (transcriptinfo[, 2] > chromstart & transcriptinfo[, 2] < chromend) | 
    (transcriptinfo[, 3] > chromstart & transcriptinfo[, 3] < chromend)
), ]

I believe you should add the situation where the zoom window falls completely inside a gene:

transcriptinfo = transcriptinfo[which(
    (transcriptinfo[, 2] > chromstart & transcriptinfo[, 2] < chromend) | 
    (transcriptinfo[, 3] > chromstart & transcriptinfo[, 3] < chromend) |
    (transcriptinfo[, 2] <= chromstart & transcriptinfo[, 3] >= chromend)
), ]

Which gives the figure I expected:

Here is the code to reproduce the plot

library(Sushi)
## Get data
Sushi_data = data(package = 'Sushi')
data(list = Sushi_data$results[,3])

## Setup position
chrom = "chr11"
chromstart = 1945000
chromend = 1960000

## Get genes
mart = useMart(biomart = "ensembl", dataset = "hsapiens_gene_ensembl")
geneinfo = getBM(attributes = c("chromosome_name", "exon_chrom_start","exon_chrom_end", "external_gene_name", "strand"), filters = c("chromosome_name", "start", "end"), values = list(gsub("chr","",chrom), chromstart, chromend), mart = mart)
names(geneinfo) = c("chrom", "start", "stop", "gene","strand")
geneinfo$score = "."
geneinfo = geneinfo[, c(1, 2, 3, 4, 6, 5)]
geneinfo$chrom = paste0("chr", sub("chr","",geneinfo$chrom))

## Set graphic parameters
par(mar=c(3,4,1,1))
layout(matrix(c(1,2,3,4), 4, 1, byrow = TRUE))

## Plot coverage
plotBedgraph(Sushi_ChIPSeq_CTCF.bedgraph,
   chrom,chromstart,chromend, 
   transparency=.50,color=SushiColors(2)(2)[1]
)                                                                                                                                                                                       
axis(side=2,las=2,tcl=.2)
labelgenome(chrom,chromstart,chromend, n=5,scale="bp")

## Plot gene
plotGenes(geneinfo, chrom,chromstart,chromend)

## Create zoom box
chrom = "chr11"
chromstart = 1956000
chromend = 1958000
zoomsregion(c(chromstart,chromend), extend=c(0.01,0.13), wideextend=0.05, offsets=c(0,0))

## Plot coverage
plotBedgraph(Sushi_ChIPSeq_CTCF.bedgraph,
   chrom,chromstart,chromend, 
   transparency=.50,color=SushiColors(2)(2)[1]
)                                                                                                                                                                                       
axis(side=2,las=2,tcl=.2)
labelgenome(chrom,chromstart,chromend, n=5,scale="bp")

## Plot gene
plotGenes(geneinfo, chrom,chromstart,chromend)

issue was resolved

Hi Doug, I've sent an email to report a issue on plotGenes, but I think maybe it's easier for you to track issues here.

I’m trying to use the Sushi package for some plotting, and I really like the tools in this package!

Just one little bug, probably, in plotGenes(). I’m trying to plot genes without giving a bed file, and I had this error:

Error in getBM(attributes = c("chromosome_name", "exon_chrom_start", "exon_chrom_end", : Invalid attribute(s): external_gene_id 
Please use the function 'listAttributes' to get valid attribute names

It looks like "external_gene_id” is not a valid attributes in biomaRt. At least when I did listAttributes(), it’s not in there.
I noticed it does have “external_gene_name” in the attribute list, though.

So I manually change “ external_gene_id” to “ external_gene_name” in the plotGenes() function, and another error occurred:

Error: could not find function "convertstrandinfo"

It looks like the function is in the package by looking at library(help = Sushi), but maybe it’s not being appropriately added in the package?

Here’s the version of package I have:
Sushi_1.10.0 biomaRt_2.28.0

Could you please let me know how can I fix these issues?

Thanks a lot!

Best,
Hurley

Dear authors,

Is there any way to remove the arrows that indicate the gene strand leaving at the same time labeltext=TRUE in the function plotGenes?
I had a look on the source code, but I couldn't figure out how to sort this out - it is only because it gets quite noise if you have a lot of genes.

Many thanks in advance;

Hello! I have been trying to use plotGenes with geneinfo = NULL and I think there might be a problem with the biomart attributes called in the function.

The function calls "external_gene_id" as attribute to retrieve from biomart, but in the version of biomart that is being called with my R session, the correct attribute name is "external_gene_name". Is this problem happening only to me? I can workaround it by calling biomart myself, then running plotGenes with my own geneinfo, but I was wondering whether the source code for plotGenes had to be corrected.

Kind regards,

Silvia

---- sessionInfo() output ------------------

R version 3.6.1 (2019-07-05)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Mojave 10.14.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.6/Resources/lib/libRlapack.dylib

Random number generation:
RNG: Mersenne-Twister
Normal: Inversion
Sample: Rounding

locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8

attached base packages:
[1] parallel stats4 stats graphics grDevices utils datasets
[8] methods base

other attached packages:
[1] org.Hs.eg.db_3.8.2
[2] TxDb.Hsapiens.UCSC.hg38.knownGene_3.4.6
[3] GenomicFeatures_1.36.4
[4] AnnotationDbi_1.46.1
[5] Biobase_2.44.0
[6] annotatr_1.10.0
[7] Sushi_1.22.0
[8] biomaRt_2.40.5
[9] zoo_1.8-8
[10] rtracklayer_1.44.4
[11] forcats_0.5.0
[12] stringr_1.4.0
[13] dplyr_1.0.2
[14] purrr_0.3.4
[15] readr_1.4.0
[16] tidyr_1.1.2
[17] tibble_3.0.4
[18] ggplot2_3.3.2
[19] tidyverse_1.3.0
[20] GenomicRanges_1.36.1
[21] GenomeInfoDb_1.20.0
[22] IRanges_2.18.3
[23] S4Vectors_0.22.1
[24] BiocGenerics_0.30.0

loaded via a namespace (and not attached):
[1] bitops_1.0-6 matrixStats_0.57.0
[3] fs_1.5.0 lubridate_1.7.9
[5] bit64_4.0.5 progress_1.2.2
[7] httr_1.4.2 tools_3.6.1
[9] backports_1.2.0 utf8_1.1.4
[11] R6_2.5.0 DBI_1.1.0
[13] colorspace_1.4-1 withr_2.3.0
[15] tidyselect_1.1.0 prettyunits_1.1.1
[17] curl_4.3 bit_4.0.4
[19] compiler_3.6.1 cli_2.1.0
[21] rvest_0.3.6 xml2_1.3.2
[23] DelayedArray_0.10.0 labeling_0.4.2
[25] scales_1.1.1 rappdirs_0.3.1
[27] digest_0.6.27 Rsamtools_2.0.3
[29] XVector_0.24.0 htmltools_0.5.0
[31] pkgconfig_2.0.3 BSgenome_1.52.0
[33] regioneR_1.16.5 fastmap_1.0.1
[35] dbplyr_2.0.0 rlang_0.4.8
[37] readxl_1.3.1 rstudioapi_0.11
[39] RSQLite_2.2.1 shiny_1.5.0
[41] generics_0.1.0 farver_2.0.3
[43] jsonlite_1.7.1 BiocParallel_1.18.1
[45] RCurl_1.98-1.2 magrittr_1.5
[47] GenomeInfoDbData_1.2.1 Matrix_1.2-18
[49] Rcpp_1.0.5 munsell_0.5.0
[51] fansi_0.4.1 lifecycle_0.2.0
[53] yaml_2.2.1 stringi_1.5.3
[55] SummarizedExperiment_1.14.1 zlibbioc_1.30.0
[57] plyr_1.8.6 BiocFileCache_1.8.0
[59] AnnotationHub_2.16.1 grid_3.6.1
[61] blob_1.2.1 promises_1.1.1
[63] crayon_1.3.4 lattice_0.20-41
[65] Biostrings_2.52.0 haven_2.3.1
[67] hms_0.5.3 pillar_1.4.6
[69] reshape2_1.4.4 reprex_0.3.0
[71] XML_3.99-0.3 glue_1.4.2
[73] BiocManager_1.30.10 modelr_0.1.8
[75] httpuv_1.5.4 vctrs_0.3.4
[77] cellranger_1.1.0 gtable_0.3.0
[79] assertthat_0.2.1 xfun_0.19
[81] mime_0.9 xtable_1.8-4
[83] broom_0.7.2 later_1.1.0.1
[85] GenomicAlignments_1.20.1 tinytex_0.27
[87] memoise_1.1.0 interactiveDisplayBase_1.22.0
[89] ellipsis_0.3.1

Hi authors, thanks for this exciting tool!

I'm trying to plot raw reads from different small RNA libraries using the plotBed function. My problem is I want to maintain fidelity for the y-axis between different libraries. I've tried common arguments (like ymax) and looked thoroughly through the recommended arguments to no avail. Perhaps I just missed it.

This is a feature available with common genome browsers (i.e. igv) - is it possible to lock the Y-axis using the plotBed function in sushi?

I linked an IGV screenshot of what I'm hoping to achieve with sushi:

Thanks for your attention and your great software.
Nate

Dear author;
can the Sushi software plot the gene density on two or more chromosome?I have many chromosomes to plot and want to plot it all into a pictures.Is there a method to solve it ?any help will be appreciate.

Thanks
Alex

This is an amazing tool! But I'm wondering is it possible to plot personal bed graph files?

chrom = "chr11" chromstart = 1650000 chromend = 2350000 plotBedgraph(my.bedgraph, chrom, chromstart, chromend, colorbycol = SushiColors(5))

object 'my.bedgraph' not found

I'm using the R "Sushi package to plot CTCF ChIPseq peaks (after MACS calling) in a bedgraph file. The file has been loaded and read by R:

A1 <-read.delim("CTCF_ChIPseq_MACSPeaks_mm9.bedgraph")
head(A1)
chrom chromstart chromend value
1 chr1 3002401 3003150 1.828043
2 chr1 3258151 3258900 1.546266
3 chr1 3315001 3315900 1.731632
4 chr1 3434551 3435300 1.581921
5 chr1 3438001 3438900 1.843272
6 chr1 3601051 3601950 10.541419
Plotting for mouse CTCF peaks (mm9) chrom12 as follows:
chrom = "12"
chromstart = 41402356
chromend = 41412356
Plotting command is below, where A1 = CTCF_ChIPSeq_MACsPeaks_mm9genome.bedgraph file which for all chromosomes:
plotBedgraph(A1, chrom, chromstart, chromend, range = NULL, color = SushiColors(2)(2)[1], lwd = 1, linecolor = NA, addscale = FALSE, overlay = FALSE, rescaleoverlay = FALSE, transparency = 0.5, flip = FALSE, xaxt = "none", yaxt = "none", xlab = "Chrom 12 mm9", ylab = "CTCF ChIPSeq Peaks", xaxs = "i", yaxs = "i", bty = "n", ymax = 1.04, colorbycol = NULL)

Error received is: [1] "not enough data within range to plot"

The plot is BLANK! Not even the x-axis or the y-axis, only the xlab as "Chrom 12 mm9" and the ylab as CTCF ChIPSeq Peaks". Using Rstudio Version 1.1.463 in Ubuntu 14.04 or 16.04 environments. What do I need to do now to get the plot I wish? Thanks in advance!