when I turn on GGSASHIMI_DEBUG and run the R script it produces, I get the following when I run the script in R:
library(ggplot2)
library(grid)
library(gridExtra)
library(data.table)
library(gtable)
scale_lwd = function(r) {
lmin = 0.1
lmax = 4
return( r*(lmax-lmin)lmin )
}
base_size = 14
height = ( 5.5 base_size*0.352777778/67 ) * 1.02
width = 10
theme_set(theme_bw(base_size=base_size))
theme_update(
plot.margin = unit(c(15,15,15,15), "pt"),
panel.grid = element_blank(),
panel.border = element_blank(),
axis.line = element_line(size=0.5),
axis.title.x = element_blank(),
axis.title.y = element_text(angle=0, vjust=0.5)
)
labels = list("#000000"="#000000","WT"="WT","1KO"="1KO")
density_list = list()
junction_list = list()
color_list = list("WT"="#000000","1KO"="#000000")
# data table with exons
ann_list = list(
"exons" = data.table(),
"introns" = data.table()
)
gtfp = ggplot()
if (length(ann_list[['introns']]) 0) {
gtfp = gtfp geom_segment(data=ann_list[['introns']], aes(x=start, xend=end, y=tx, yend=tx), size=0.3)
gtfp = gtfp geom_segment(data=txarrows, aes(x=V1,xend=V2,y=tx,yend=tx), arrow=arrow(length=unit(0.02,"npc")))
}
if (length(ann_list[['exons']]) 0) {
gtfp = gtfp geom_segment(data=ann_list[['exons']], aes(x=start, xend=end, y=tx, yend=tx), size=5, alpha=1)
}
gtfp = gtfp scale_y_discrete(expand=c(0,0.5))
gtfp = gtfp scale_x_continuous(expand=c(0,0.25))
gtfp = gtfp coord_cartesian(xlim = c(6064339,6070090))
gtfp = gtfp labs(y=NULL)
gtfp = gtfp theme(axis.line = element_blank(), axis.text.x = element_blank(), axis.ticks = element_blank())
[density lists and junction lists omitted bc too long]
junction_list[["1KO"]] = data.frame(x=c(6064422,6064422,6065846,6064422,6064422,6065846,6064422,6064422,6065846), xend=c(6065705,6069831,6069831,6065705,6069831,6069831,6065705,6069831,6069831), y=c(39,39,13,39,39,13,39,39,13), yend=c(14,34,34,14,34,34,14,34,34), count=c(13,22,8,13,22,8,13,22,8))
Error: unexpected symbol in:
"20,6064821,6064822,6064823,6064824,6064825,6064826,6064827,6064828,6064829,6064830,6064831,6064832,6064833,6064834,6064835,6064836,6064837,6064838,6064839,6064840,6
junction_list"
pdf(NULL) # just to remove the blank pdf produced by ggplotGrob
if(packageVersion('ggplot2') = '3.0.0'){ # fix problems with ggplot2 vs 3.0.0
vs = 1
} else {
vs = 0
}
if(TRUE) {
print(max(unlist(lapply(density_list, function(df){max(df$y)}))))
maxheight = max(unlist(lapply(density_list, function(df){max(df$y)})))
breaks_y = labeling::extended(0, maxheight, m = 4)
}
[1] -Inf
Error in seq.default(from = best$lmin, to = best$lmax, by = best$lstep) :
'from' must be of length 1
In addition: Warning messages:
1: In max(unlist(lapply(density_list, function(df) { :
no non-missing arguments to max; returning -Inf
2: In max(unlist(lapply(density_list, function(df) { :
no non-missing arguments to max; returning -Inf
if(exists('coord_dict')){
all_pos_shrinked = do.call(rbind, density_list)$x
s2r = merge(intersected_introns, coord_dict, by.x = 'real_xend', by.y = 'real')
s2r = merge(s2r, coord_dict, by.x = 'real_x', by.y = 'real', suffixes = c('_xend', '_x'))
breaks_x_shrinked = labeling::extended(min(all_pos_shrinked), max(all_pos_shrinked), m = 5)
breaks_x = c()
out_range = c()
for (b in breaks_x_shrinked){
iintron = FALSE
for (j in 1:nrow(s2r)){
l = s2r[j, ]
if(b = l$shrinked_x && b <= l$shrinked_xend){
# Intersected intron
p = (b-l$shrinked_x)/(l$shrinked_xend - l$shrinked_x)
realb = round(l$real_x p*(l$real_xend - l$real_x))
breaks_x = c(breaks_x, realb)
iintron = TRUE
break
}
}
if (!iintron){
# Exon, upstream/downstream intergenic region or intron (not intersected)
if(b <= min(s2r$shrinked_x)) {
l <- s2r[which.min(s2r$shrinked_x), ]
if(any(b == all_pos_shrinked)){
# Boundary (subtract)
s = l$shrinked_x - b
realb = l$real_x - s
breaks_x = c(breaks_x, realb)
} else {
out_range <- c(out_range, which(breaks_x_shrinked == b))
}
} else if (b = max(s2r$shrinked_xend)){
l <- s2r[which.max(s2r$shrinked_xend), ]
if(any(b == all_pos_shrinked)){
# Boundary (sum)
s = b - l$shrinked_xend
realb = l$real_xend s
breaks_x = c(breaks_x, realb)
} else {
out_range <- c(out_range, which(breaks_x_shrinked == b))
}
} else {
delta = b-s2r$shrinked_xend
delta[delta < 0] = Inf
l = s2r[which.min(delta), ]
# Internal (sum)
s = b - l$shrinked_xend
realb = l$real_xend s
breaks_x = c(breaks_x, realb)
}
}
}
if(length(out_range)) {
breaks_x_shrinked = breaks_x_shrinked[-out_range]
}
}
density_grobs = list();
for (bam_index in 1:length(density_list)) {
id = names(density_list)[bam_index]
d = data.table(density_list[[id]])
junctions = data.table(junction_list[[id]])
# Density plot
gp = ggplot(d) geom_bar(aes(x, y), width=1, position='identity', stat='identity', fill=color_list[[id]], alpha=0.05)
gp = gp labs(y=labels[[id]])
if(exists('coord_dict')) {
gp = gp scale_x_continuous(expand=c(0, 0.25), breaks = breaks_x_shrinked, labels = breaks_x)
} else {
gp = gp scale_x_continuous(expand=c(0, 0.25))
}
if(!TRUE){
maxheight = max(d[['y']])
breaks_y = labeling::extended(0, maxheight, m = 4)
gp = gp scale_y_continuous(breaks = breaks_y)
} else {
gp = gp scale_y_continuous(breaks = breaks_y, limits = c(NA, maxheight))
}
# Aggregate junction counts
row_i = c()
if (nrow(junctions) 0 ) {
junctions$jlabel = as.character(junctions$count)
junctions = setNames(junctions[,.(max(y), max(yend),round(mean(count)),paste(jlabel,collapse=",")), keyby=.(x,xend)], names(junctions))
if ("mean" != "") {
junctions = setNames(junctions[,.(max(y), max(yend),round(mean(count)),round(mean(count))), keyby=.(x,xend)], names(junctions))
}
# The number of rows (unique junctions per bam) has to be calculated after aggregation
row_i = 1:nrow(junctions)
}
for (i in row_i) {
j_tot_counts = sum(junctions[['count']])
j = as.numeric(junctions[i,1:5])
if ("mean" != "") {
j[3] = ifelse(length(d[x==j[1]-1,y])==0, 0, max(as.numeric(d[x==j[1]-1,y])))
j[4] = ifelse(length(d[x==j[2]1,y])==0, 0, max(as.numeric(d[x==j[2]1,y])))
}
# Find intron midpoint
xmid = round(mean(j[1:2]), 1)
ymid = max(j[3:4]) * 1.2
# Thickness of the arch
lwd = scale_lwd(j[5]/j_tot_counts)
curve_par = gpar(lwd=lwd, col=color_list[[id]])
# Arc grobs
# Choose position of the arch (top or bottom)
nss = i
if (nss%%2 == 0) { #bottom
ymid = -0.3 * maxheight
# Draw the arcs
# Left
curve = xsplineGrob(x=c(0, 0, 1, 1), y=c(1, 0, 0, 0), shape=1, gp=curve_par)
gp = gp annotation_custom(grob = curve, j[1], xmid, 0, ymid)
# Right
curve = xsplineGrob(x=c(1, 1, 0, 0), y=c(1, 0, 0, 0), shape=1, gp=curve_par)
gp = gp annotation_custom(grob = curve, xmid, j[2], 0, ymid)
}
if (nss%%2 != 0) { #top
# Draw the arcs
# Left
curve = xsplineGrob(x=c(0, 0, 1, 1), y=c(0, 1, 1, 1), shape=1, gp=curve_par)
gp = gp annotation_custom(grob = curve, j[1], xmid, j[3], ymid)
# Right
curve = xsplineGrob(x=c(1, 1, 0, 0), y=c(0, 1, 1, 1), shape=1, gp=curve_par)
gp = gp annotation_custom(grob = curve, xmid, j[2], j[4], ymid)
}
# Add junction labels
gp = gp annotate("label", x = xmid, y = ymid, label = as.character(junctions[i,6]),
vjust=0.5, hjust=0.5, label.padding=unit(0.01, "lines"),
label.size=NA, size=(base_size*0.352777778)*0.6
)
# gp = gp annotation_custom(grob = rectGrob(x=0, y=0, gp=gpar(col="red"), just=c("left","bottom")), xmid, j[2], j[4], ymid)
# gp = gp annotation_custom(grob = rectGrob(x=0, y=0, gp=gpar(col="green"), just=c("left","bottom")), j[1], xmid, j[3], ymid)
}
gpGrob = ggplotGrob(gp);
gpGrob$layout$clip[gpGrob$layout$name=="panel"] <- "off"
if (bam_index == 1) {
maxWidth = gpGrob$widths[2vs] gpGrob$widths[3vs]; # fix problems ggplot2 vs
maxYtextWidth = gpGrob$widths[3vs]; # fix problems ggplot2 vs
# Extract x axis grob (trim=F -- keep empty cells)
xaxisGrob <- gtable_filter(gpGrob, "axis-b", trim=F)
xaxisGrob$heights[8vs] = gpGrob$heights[1] # fix problems ggplot2 vs
x.axis.height = gpGrob$heights[7vs] gpGrob$heights[1] # fix problems ggplot2 vs
}
# Remove x axis from all density plots
kept_names = gpGrob$layout$name[gpGrob$layout$name != "axis-b"]
gpGrob <- gtable_filter(gpGrob, paste(kept_names, sep="", collapse="|"), trim=F)
# Find max width of y text and y label and max width of y text
maxWidth = grid::unit.pmax(maxWidth, gpGrob$widths[2vs] gpGrob$widths[3vs]); # fix problems ggplot2 vs
maxYtextWidth = grid::unit.pmax(maxYtextWidth, gpGrob$widths[3vs]); # fix problems ggplot2 vs
density_grobs[[id]] = gpGrob;
}
Error in density_list[[id]] :
attempt to select less than one element in get1index
# Add x axis grob after density grobs BEFORE annotation grob
density_grobs[["xaxis"]] = xaxisGrob
# Annotation grob
if (1.0 == 1) {
gtfGrob = ggplotGrob(gtfp);
maxWidth = grid::unit.pmax(maxWidth, gtfGrob$widths[2vs] gtfGrob$widths[3vs]); # fix problems ggplot2 vs
density_grobs[['gtf']] = gtfGrob;
}
Error in ans[ypos] <- rep(yes, length.out = len)[ypos] :
replacement has length zero
In addition: Warning message:
In rep(yes, length.out = len) : 'x' is NULL so the result will be NULL
# Reassign grob widths to align the plots
for (id in names(density_grobs)) {
density_grobs[[id]]$widths[1] <- density_grobs[[id]]$widths[1] maxWidth - (density_grobs[[id]]$widths[2vs] maxYtextWidth); # fix problems ggplot2 vs
density_grobs[[id]]$widths[3vs] <- maxYtextWidth # fix problems ggplot2 vs
}
Error: object of type 'closure' is not subsettable
# Heights for density, x axis and annotation
heights = unit.c(
unit(rep(2, length(density_list)), "in"),
x.axis.height,
unit(1.5*1.0, "in")
)
Error in unit(rep(2, length(density_list)), "in") :
'x' and 'units' must have length 0
# Arrange grobs
argrobs = arrangeGrob(
grobs=density_grobs,
ncol=1,
heights = heights,
);
Error in arrangeGrob(grobs = density_grobs, ncol = 1, heights = heights, :
object 'heights' not found
# Save plot to file in the requested format
if ("pdf" == "tiff"){
# TIFF images will be lzw-compressed
ggsave("/orange/swanson/carter.h/ChP_Splicing_Project/Mbnl_12KO_ChP/rMATS_output_1KO_ChP_Triplicate/testing_SE_sashimi_plots/Arhgap12_chr18_6064340-6070091_meanj_-.pdf", plot = argrobs, device = "tiff", width = width, height = height, units = "in", dpi = 300, compression = "lzw", limitsize = FALSE)
} else {
ggsave("/orange/swanson/carter.h/ChP_Splicing_Project/Mbnl_12KO_ChP/rMATS_output_1KO_ChP_Triplicate/testing_SE_sashimi_plots/Arhgap12_chr18_6064340-6070091_meanj_-.pdf", plot = argrobs, device = "pdf", width = width, height = height, units = "in", dpi = 300, limitsize = FALSE)
}
Error in grDevices::pdf(file = filename, ..., version = version) :
cannot open file '/orange/swanson/carter.h/ChP_Splicing_Project/Mbnl_12KO_ChP/rMATS_output_1KO_ChP_Triplicate/testing_SE_sashimi_plots/Arhgap12_chr18_6064340-6070091_meanj_-.pdf'
dev.log = dev.off()
React
Typescript
Django
Laravel
Facebook
Microsoft
Google
Alibaba
D3
Tencent