https:\/\/github.com\/suhrig\/arriba\/blob\/master\/draw_fusions.R<\/a><\/p>\n#!\/usr\/bin\/env Rscript\n\n# print warnings as they happen instead of collecting them for after a loop ends\noptions(warn=1)\n\n# define valid parameters\nparameters <- list(\n fusions=list("fusionsFile", "file", "fusions.tsv", T),\n annotation=list("exonsFile", "file", "annotation.gtf", T),\n output=list("outputFile", "string", "output.pdf", T),\n alignments=list("alignmentsFile", "file", "Aligned.sortedByCoord.out.bam"),\n cytobands=list("cytobandsFile", "file", "cytobands.tsv"),\n minConfidenceForCircosPlot=list("minConfidenceForCircosPlot", "string", "medium"),\n proteinDomains=list("proteinDomainsFile", "file", "protein_domains.gff3"),\n sampleName=list("sampleName", "string", ""),\n squishIntrons=list("squishIntrons", "bool", T),\n printExonLabels=list("printExonLabels", "bool", T),\n render3dEffect=list("render3dEffect", "bool", T),\n pdfWidth=list("pdfWidth", "numeric", 11.692),\n pdfHeight=list("pdfHeight", "numeric", 8.267),\n color1=list("color1", "string", "#e5a5a5"),\n color2=list("color2", "string", "#a7c4e5"),\n mergeDomainsOverlappingBy=list("mergeDomainsOverlappingBy", "numeric", 0.9),\n optimizeDomainColors=list("optimizeDomainColors", "bool", F),\n fontSize=list("fontSize", "numeric", 1),\n fontFamily=list("fontFamily", "string", "Helvetica"),\n showIntergenicVicinity=list("showIntergenicVicinity", "string", "0"),\n transcriptSelection=list("transcriptSelection", "string", "provided"),\n fixedScale=list("fixedScale", "numeric", 0),\n coverageRange=list("coverageRange", "string", "0")\n)\n\n# print help if necessary\nargs <- commandArgs(trailingOnly=T)\nif (any(grepl("^--help", args)) || length(args) == 0) {\n usage <- "Usage: draw_fusions.R"\n for (parameter in names(parameters)) {\n usage <- paste0(usage, " ")\n if (length(parameters[[parameter]]) <= 3 || !parameters[[parameter]][[4]])\n usage <- paste0(usage, "[")\n usage <- paste0(usage, "--", parameter, "=", parameters[[parameter]][[3]])\n if (length(parameters[[parameter]]) <= 3 || !parameters[[parameter]][[4]])\n usage <- paste0(usage, "]")\n }\n message(usage)\n quit("no", ifelse(length(args) == 0, 1, 0))\n}\n\n# make sure mandatory arguments are present\nfor (parameter in names(parameters))\n if (length(parameters[[parameter]]) > 3 && parameters[[parameter]][[4]])\n if (!any(grepl(paste0("^--", parameter, "="), args), perl=T))\n stop(paste0("Missing mandatory argument: --", parameter))\n\n# set default values\nfor (parameter in names(parameters))\n assign(parameters[[parameter]][[1]], ifelse(parameters[[parameter]][[2]] == "file", "", parameters[[parameter]][[3]]))\n\n# parse command-line parameters\nfor (arg in args) {\n argName <- sub("=.*", "", sub("^--", "", arg, perl=T), perl=T)\n argValue <- sub("^[^=]*=", "", arg, perl=T)\n if (!(argName %in% names(parameters)) || !grepl("^--", arg, perl=T))\n stop(paste("Unknown parameter:", arg))\n if (parameters[[argName]][[2]] == "bool") {\n if (argValue %in% c("TRUE", "T", "FALSE", "F")) {\n assign(parameters[[argName]][[1]], as.logical(argValue))\n } else {\n stop(paste0("Invalid argument to --", argName))\n }\n } else if (parameters[[argName]][[2]] == "string") {\n assign(parameters[[argName]][[1]], argValue)\n } else if (parameters[[argName]][[2]] == "numeric") {\n if (is.na(suppressWarnings(as.numeric(argValue))))\n stop(paste0("Invalid argument to --", argName))\n assign(parameters[[argName]][[1]], as.numeric(argValue))\n } else if (parameters[[argName]][[2]] == "file") {\n if (file.access(argValue) == -1)\n stop(paste("Cannot read file:", argValue))\n assign(parameters[[argName]][[1]], argValue)\n }\n}\n\n# validate values of parameters\nif (cytobandsFile == "")\n warning("Missing parameter '--cytobands'. No ideograms and circos plots will be drawn.")\nif (!(minConfidenceForCircosPlot %in% c("none", "low", "medium", "high")))\n stop("Invalid argument to --minConfidenceForCircosPlot")\nshowIntergenicVicinity <- as.list(unlist(strsplit(showIntergenicVicinity, ",", fixed=T)))\nif (!(length(showIntergenicVicinity) %in% c(1,4)))\n stop(paste0("Invalid argument to --showIntergenicVicinity"))\nshowIntergenicVicinity <- lapply(showIntergenicVicinity, function(x) {\n if (x == "closestGene") {\n return("exon")\n } else if (x == "closestProteinCodingGene") {\n return("CDS")\n } else if (is.na(suppressWarnings(as.numeric(x))) || as.numeric(x) < 0) {\n stop(paste0("Invalid argument to --showIntergenicVicinity"))\n } else {\n return(as.numeric(x))\n }\n})\nif (length(showIntergenicVicinity) == 1)\n showIntergenicVicinity <- rep(showIntergenicVicinity, 4)\nif (squishIntrons)\n if (any(!is.numeric(unlist(showIntergenicVicinity))) || any(showIntergenicVicinity > 0))\n stop("--squishIntrons must be disabled, when --showIntergenicVicinity is > 0")\nif (!(transcriptSelection %in% c("coverage", "provided", "canonical")))\n stop("Invalid argument to --transcriptSelection")\nif (fixedScale < 0)\n stop("Invalid argument to --fixedScale")\nif (!(fontFamily %in% names(pdfFonts())))\n stop(paste0("Unknown font: ", fontFamily, ". Available fonts: ", paste(names(pdfFonts()), collapse=", ")))\ncoverageRange <- suppressWarnings(as.numeric(unlist(strsplit(coverageRange, ",", fixed=T))))\nif (!(length(coverageRange) %in% 1:2) || any(is.na(coverageRange)) || any(coverageRange < 0))\n stop("Invalid argument to --coverageRange")\n\n# check if required packages are installed\nif (!suppressPackageStartupMessages(require(GenomicRanges)))\n warning("Package 'GenomicRanges' is not installed. No protein domains and circos plots will be drawn.")\nif (!suppressPackageStartupMessages(require(circlize)))\n warning("Package 'circlize' is not installed. No circos plots will be drawn.")\nif (alignmentsFile != "")\n if (!suppressPackageStartupMessages(require(GenomicAlignments)))\n stop("Package 'GenomicAlignments' must be installed when '--alignments' is used")\n\n# define colors\nchangeColorBrightness <- function(color, delta) {\n rgb(\n min(255,max(0,col2rgb(color)["red",]+delta)),\n min(255,max(0,col2rgb(color)["green",]+delta)),\n min(255,max(0,col2rgb(color)["blue",]+delta)),\n maxColorValue=255\n )\n}\ngetDarkColor <- function(color) { changeColorBrightness(color, -100) }\ngetBrightColor <- function(color) { changeColorBrightness(color, +190) }\ndarkColor1 <- getDarkColor(color1)\ndarkColor2 <- getDarkColor(color2)\ncircosColors <- c(translocation="#000000", duplication="#00bb00", deletion="#ff0000", inversion="#0000ff")\n\n# convenience functions to add\/remove "chr" prefix\naddChr <- function(contig) {\n ifelse(contig == "MT", "chrM", paste0("chr", contig))\n}\nremoveChr <- function(contig) {\n sub("^chr", "", sub("^chrM", "MT", contig, perl=T), perl=T)\n}\n\n# convenience function to check if a value is between two others\nbetween <- function(value, start, end) {\n value >= start & value <= end\n}\n\n# read fusions\nfusions <- read.table(fusionsFile, stringsAsFactors=F, sep="\\t", header=T, comment.char="", quote="")\nif (colnames(fusions)[1] == "X.gene1") { # Arriba output\n colnames(fusions)[colnames(fusions) %in% c("X.gene1", "strand1.gene.fusion.", "strand2.gene.fusion.")] <- c("gene1", "strand1", "strand2")\n fusions$display_contig1 <- sub(":[^:]*$", "", fusions$breakpoint1, perl=T)\n fusions$display_contig2 <- sub(":[^:]*$", "", fusions$breakpoint2, perl=T)\n fusions$contig1 <- removeChr(fusions$display_contig1)\n fusions$contig2 <- removeChr(fusions$display_contig2)\n fusions$breakpoint1 <- as.numeric(sub(".*:", "", fusions$breakpoint1, perl=T))\n fusions$breakpoint2 <- as.numeric(sub(".*:", "", fusions$breakpoint2, perl=T))\n fusions$split_reads1 <- fusions$split_reads1\n fusions$split_reads2 <- fusions$split_reads2\n fusions$type <- sub(".*(translocation|duplication|deletion|inversion).*", "\\\\1", fusions$type, perl=T)\n fusions$fusion_transcript <- gsub("[()^$]", "", fusions$fusion_transcript)\n} else if (colnames(fusions)[1] == "X.FusionName") { # STAR-Fusion\n fusions$gene1 <- sub("\\\\^.*", "", fusions$LeftGene, perl=T)\n fusions$gene2 <- sub("\\\\^.*", "", fusions$RightGene, perl=T)\n fusions$strand1 <- sub(".*:(.)$", "\\\\1\/\\\\1", fusions$LeftBreakpoint, perl=T)\n fusions$strand2 <- sub(".*:(.)$", "\\\\1\/\\\\1", fusions$RightBreakpoint, perl=T)\n fusions$display_contig1 <- sub(":[^:]*:[^:]*$", "", fusions$LeftBreakpoint, perl=T)\n fusions$display_contig2 <- sub(":[^:]*:[^:]*$", "", fusions$RightBreakpoint, perl=T)\n fusions$contig1 <- removeChr(fusions$display_contig1)\n fusions$contig2 <- removeChr(fusions$display_contig2)\n fusions$breakpoint1 <- as.numeric(sub(".*:([^:]*):[^:]*$", "\\\\1", fusions$LeftBreakpoint, perl=T))\n fusions$breakpoint2 <- as.numeric(sub(".*:([^:]*):[^:]*$", "\\\\1", fusions$RightBreakpoint, perl=T))\n fusions$direction1 <- ifelse(grepl(":\\\\+$", fusions$LeftBreakpoint, perl=T), "downstream", "upstream")\n fusions$direction2 <- ifelse(grepl(":\\\\+$", fusions$RightBreakpoint, perl=T), "upstream", "downstream")\n fusions$gene_id1 <- sub(".*\\\\^", "", fusions$LeftGene, perl=T)\n fusions$gene_id2 <- sub(".*\\\\^", "", fusions$RightGene, perl=T)\n fusions$transcript_id1 <- ifelse(rep(!("CDS_LEFT_ID" %in% colnames(fusions)), nrow(fusions)), ".", fusions$CDS_LEFT_ID)\n fusions$transcript_id2 <- ifelse(rep(!("CDS_RIGHT_ID" %in% colnames(fusions)), nrow(fusions)), ".", fusions$CDS_RIGHT_ID)\n fusions$fusion_transcript <- ifelse(rep(!("FUSION_CDS" %in% colnames(fusions)), nrow(fusions)), ".", toupper(sub("([a-z]*)", "\\\\1|", fusions$FUSION_CDS, perl=T)))\n fusions$reading_frame <- ifelse(rep(!("PROT_FUSION_TYPE" %in% colnames(fusions)), nrow(fusions)), ".", ifelse(fusions$PROT_FUSION_TYPE == "INFRAME", "in-frame", ifelse(fusions$PROT_FUSION_TYPE == "FRAMESHIFT", "out-of-frame", ".")))\n fusions$split_reads <- fusions$JunctionReadCount\n fusions$discordant_mates <- fusions$SpanningFragCount\n fusions$site1 <- rep("exon", nrow(fusions))\n fusions$site2 <- rep("exon", nrow(fusions))\n fusions$confidence <- rep("high", nrow(fusions))\n fusions$type <- ifelse(fusions$contig1 != fusions$contig2, "translocation", ifelse(fusions$direction1 == fusions$direction2, "inversion", ifelse((fusions$direction1 == "downstream") == (fusions$breakpoint1 < fusions$breakpoint2), "deletion", "duplication")))\n} else {\n stop("Unrecognized fusion file format")\n}\n\npdf(outputFile, onefile=T, width=pdfWidth, height=pdfHeight, title=ifelse(sampleName != "", sampleName, fusionsFile))\npar(family=fontFamily)\n\nif (nrow(fusions) == 0) {\n plot(0, 0, type="l", xaxt="n", yaxt="n", xlab="", ylab="")\n text(0, 0, "empty input file")\n warning("empty input file")\n dev.off()\n quit("no")\n}\n\n# read cytoband annotation\ncytobands <- NULL\nif (cytobandsFile != "") {\n cytobands <- read.table(cytobandsFile, header=T, colClasses=c("character", "numeric", "numeric", "character", "character"))\n cytobands <- cytobands[order(cytobands$contig, cytobands$start, cytobands$end),]\n}\n\n# read exon annotation\nmessage("Loading annotation")\nexons <- scan(exonsFile, what=list(contig="",src="",type="",start=0,end=0,score="",strand="",frame="",attributes=""), sep="\\t", comment.char="#", quote='"', multi.line=F)\nattr(exons, "row.names") <- .set_row_names(length(exons[[1]]))\nclass(exons) <- "data.frame"\nexons <- exons[exons$type %in% c("exon","CDS"),c("contig","type","start","end","strand","attributes")]\nexons$contig <- removeChr(exons$contig)\nparseGtfAttribute <- function(attribute, gtf) {\n parsed <- sub(paste0(".*", attribute, "[ =]([^;]+).*"), "\\\\1", gtf$attributes, perl=T)\n failedToParse <- parsed == gtf$attributes\n if (any(failedToParse)) {\n warning(paste0("Failed to parse '", attribute, "' attribute of ", sum(failedToParse), " record(s)."))\n parsed <- ifelse(failedToParse, "", parsed)\n }\n return(parsed)\n}\nexons$geneID <- parseGtfAttribute("gene_id", exons)\nexons$geneName <- parseGtfAttribute("gene_name", exons)\nexons$geneName <- ifelse(exons$geneName == "", exons$geneID, exons$geneName)\nexons$transcript <- parseGtfAttribute("transcript_id", exons)\nexons$exonNumber <- ifelse(rep(printExonLabels, nrow(exons)), parseGtfAttribute("exon_number", exons), "")\n\n# read protein domain annotation\nproteinDomains <- NULL\nif (proteinDomainsFile != "") {\n message("Loading protein domains")\n proteinDomains <- scan(proteinDomainsFile, what=list(contig="",src="",type="",start=0,end=0,score="",strand="",frame="",attributes=""), sep="\\t", comment.char="", quote="", multi.line=F)\n attr(proteinDomains, "row.names") <- .set_row_names(length(proteinDomains[[1]]))\n class(proteinDomains) <- "data.frame"\n proteinDomains$color <- parseGtfAttribute("color", proteinDomains)\n proteinDomains$proteinDomainName <- sapply(parseGtfAttribute("Name", proteinDomains), URLdecode)\n proteinDomains$proteinDomainID <- parseGtfAttribute("protein_domain_id", proteinDomains)\n}\n\n# insert dummy annotations for intergenic breakpoints\nif (any(fusions$site1 == "intergenic" | fusions$site2 == "intergenic")) {\n intergenicBreakpoints <- rbind(\n setNames(fusions[fusions$site1 == "intergenic",c("gene1", "strand1", "contig1", "breakpoint1")], c("gene", "strand", "contig", "breakpoint")),\n setNames(fusions[fusions$site2 == "intergenic",c("gene2", "strand2", "contig2", "breakpoint2")], c("gene", "strand", "contig", "breakpoint"))\n )\n exons <- rbind(exons, data.frame(\n contig=intergenicBreakpoints$contig,\n type="intergenic",\n start=sapply(intergenicBreakpoints$breakpoint-1000, max, 1),\n end=intergenicBreakpoints$breakpoint+1000,\n strand=".",\n attributes="",\n geneName=intergenicBreakpoints$gene,\n geneID=paste0(intergenicBreakpoints$contig, ":", intergenicBreakpoints$breakpoint),\n transcript=paste0(intergenicBreakpoints$contig, ":", intergenicBreakpoints$breakpoint),\n exonNumber="intergenic"\n ))\n fusions[fusions$site1 == "intergenic","gene_id1"] <- paste0(fusions[fusions$site1 == "intergenic","contig1"], ":", fusions[fusions$site1 == "intergenic","breakpoint1"])\n fusions[fusions$site2 == "intergenic","gene_id2"] <- paste0(fusions[fusions$site2 == "intergenic","contig2"], ":", fusions[fusions$site2 == "intergenic","breakpoint2"])\n}\n\ndrawVerticalGradient <- function(left, right, y, color, selection=NULL) {\n # check if gradient should only be drawn in part of the region\n if (!is.null(selection)) {\n y <- y[selection]\n left <- left[selection]\n right <- right[selection]\n }\n # draw gradient\n for (i in 1:length(y)) {\n polygon(\n c(left[1:i], right[1:i]),\n c(y[1:i], y[i:1]),\n border=NA,\n col=rgb(col2rgb(color)["red",], col2rgb(color)["green",], col2rgb(color)["blue",], col2rgb(color, alpha=T)["alpha",]*(1\/length(y)), max=255)\n )\n }\n}\nif (!render3dEffect) # nullify function, if no 3D effect should be drawn\n drawVerticalGradient <- function(left, right, y, color, selection=NULL) { }\n\ndrawCurlyBrace <- function(left, right, top, bottom, tip) {\n smoothness <- 20\n x <- cumsum(exp(-seq(-2.5, 2.5, len=smoothness)^2))\n x <- x\/max(x)\n y <- seq(top, bottom, len=smoothness)\n lines(left+(tip-left)+x*(left-tip), y)\n lines(tip+x*(right-tip), y)\n}\n\ndrawIdeogram <- function(adjust, left, right, y, cytobands, contig, breakpoint) {\n # define design of ideogram\n bandColors <- setNames(rgb(100:0, 100:0, 100:0, maxColorValue=100), paste0("gpos", 0:100))\n bandColors <- c(bandColors, gneg="#ffffff", acen="#ec4f4f", stalk="#0000ff")\n cytobands$color <- bandColors[cytobands$giemsa]\n arcSteps <- 30 # defines roundness of arc\n curlyBraceHeight <- 0.03\n ideogramHeight <- 0.04\n ideogramWidth <- 0.4\n # extract bands of given contig\n bands <- cytobands[cytobands$contig==contig,]\n if (nrow(bands) == 0) {\n warning(paste("Ideogram of contig", contig, "cannot be drawn, because no Giemsa staining information is available."))\n return(NULL)\n }\n # scale width of ideogram to fit inside given region\n bands$left <- bands$start \/ max(cytobands$end) * ideogramWidth\n bands$right <- bands$end \/ max(cytobands$end) * ideogramWidth\n # left\/right-align cytobands\n offset <- ifelse(adjust=="left", left, right - max(bands$right))\n bands$left <- bands$left + offset\n bands$right <- bands$right + offset\n # draw curly braces\n tip <- min(bands$left) + (max(bands$right)-min(bands$left)) \/ (max(bands$end)-min(bands$start)) * breakpoint\n drawCurlyBrace(left, right, y-0.05+curlyBraceHeight, y-0.05, tip)\n # draw title of chromosome\n text((max(bands$right)+min(bands$left))\/2, y+0.07, paste("chromosome", contig), font=2, cex=fontSize, adj=c(0.5,0))\n # draw name of band\n bandName <- bands[which(between(breakpoint, bands$start, bands$end)), "name"]\n text(tip, y+0.03, bandName, cex=fontSize, adj=c(0.5,0))\n # draw start of chromosome\n leftArcX <- bands[1,"left"] + (1+cos(seq(pi\/2,1.5*pi,len=arcSteps))) * (bands[1,"right"]-bands[1,"left"])\n leftArcY <- y + sin(seq(pi\/2,1.5*pi,len=arcSteps)) * (ideogramHeight\/2)\n polygon(leftArcX, leftArcY, col=bands[1,"color"])\n # draw bands\n centromereStart <- NULL\n centromereEnd <- NULL\n for (band in 2:(nrow(bands)-1)) {\n if (bands[band,"giemsa"] != "acen") {\n rect(bands[band,"left"], y-ideogramHeight\/2, bands[band,"right"], y+ideogramHeight\/2, col=bands[band,"color"])\n } else { # draw centromere\n if (is.null(centromereStart)) {\n polygon(c(bands[band,"left"], bands[band,"right"], bands[band,"left"]), c(y-ideogramHeight\/2, y, y+ideogramHeight\/2), col=bands[band,"color"])\n centromereStart <- bands[band,"left"]\n } else {\n polygon(c(bands[band,"right"], bands[band,"left"], bands[band,"right"]), c(y-ideogramHeight\/2, y, y+ideogramHeight\/2), col=bands[band,"color"])\n centromereEnd <- bands[band,"right"]\n }\n }\n }\n # draw end of chromosome\n band <- nrow(bands)\n rightArcX <- bands[band,"right"] - (1+cos(seq(1.5*pi,pi\/2,len=arcSteps))) * (bands[band,"right"]-bands[band,"left"])\n rightArcY <- y + sin(seq(pi\/2,1.5*pi,len=arcSteps)) * ideogramHeight\/2\n polygon(rightArcX, rightArcY, col=bands[band,"color"])\n # if there is no centromere, make an artificial one with length zero\n if (is.null(centromereStart) || is.null(centromereEnd)) {\n centromereStart <- bands[1,"right"]\n centromereEnd <- bands[1,"right"]\n }\n # draw gradients for 3D effect\n drawVerticalGradient(leftArcX, rep(centromereStart, arcSteps), leftArcY, rgb(0,0,0,0.8), 1:round(arcSteps*0.4)) # black from top on p-arm\n drawVerticalGradient(leftArcX, rep(centromereStart, arcSteps), leftArcY, rgb(1,1,1,0.7), round(arcSteps*0.4):round(arcSteps*0.1)) # white to top on p-arm\n drawVerticalGradient(leftArcX, rep(centromereStart, arcSteps), leftArcY, rgb(1,1,1,0.7), round(arcSteps*0.4):round(arcSteps*0.6)) # white to bottom on p-arm\n drawVerticalGradient(leftArcX, rep(centromereStart, arcSteps), leftArcY, rgb(0,0,0,0.9), arcSteps:round(arcSteps*0.5)) # black from bottom on p-arm\n drawVerticalGradient(rightArcX, rep(centromereEnd, arcSteps), rightArcY, rgb(0,0,0,0.8), 1:round(arcSteps*0.4)) # black from top on q-arm\n drawVerticalGradient(rightArcX, rep(centromereEnd, arcSteps), rightArcY, rgb(1,1,1,0.7), round(arcSteps*0.4):round(arcSteps*0.1)) # white to top on q-arm\n drawVerticalGradient(rightArcX, rep(centromereEnd, arcSteps), rightArcY, rgb(1,1,1,0.7), round(arcSteps*0.4):round(arcSteps*0.6)) # white to bottom on q-arm\n drawVerticalGradient(rightArcX, rep(centromereEnd, arcSteps), rightArcY, rgb(0,0,0,0.9), arcSteps:round(arcSteps*0.5)) # black from bottom on q-arm\n}\n\ndrawCoverage <- function(left, right, y, coverage, start, end, color) {\n maxResolution <- 5000 # max number of data points to draw coverage\n # draw coverage as bars\n if (!is.null(coverage)) {\n coverageData <- as.numeric(coverage[IRanges(sapply(start, max, min(start(coverage))), sapply(end, min, max(end(coverage))))])\n # downsample to maxResolution, if there are too many data points\n coverageData <- aggregate(coverageData, by=list(round(1:length(coverageData) * (right-left) * maxResolution\/length(coverageData))), mean)$x\n polygon(c(left, seq(left, right, length.out=length(coverageData)), right), c(y, y+coverageData*0.1, y), col=color, border=NA)\n }\n}\n\ndrawStrand <- function(left, right, y, color, strand) {\n if (strand %in% c("+", "-")) {\n # draw strand\n lines(c(left+0.001, right-0.001), c(y, y), col=color, lwd=2)\n lines(c(left+0.001, right-0.001), c(y, y), col=rgb(1,1,1,0.1), lwd=1)\n # indicate orientation\n if (right - left > 0.01)\n for (i in seq(left+0.005, right-0.005, by=sign(right-left-2*0.005)*0.01)) {\n arrows(i, y, i+0.001*ifelse(strand=="+", 1, -1), y, col=color, length=0.05, lwd=2, angle=60)\n arrows(i, y, i+0.001*ifelse(strand=="+", 1, -1), y, col=rgb(1,1,1,0.1), length=0.05, lwd=1, angle=60)\n }\n }\n}\n\ndrawExon <- function(left, right, y, color, title, type) {\n gradientSteps <- 10 # defines smoothness of gradient\n exonHeight <- 0.03\n if (type == "CDS") {\n # draw coding regions as thicker bars\n rect(left, y+exonHeight, right, y+exonHeight\/2-0.001, col=color, border=NA)\n rect(left, y-exonHeight, right, y-exonHeight\/2+0.001, col=color, border=NA)\n # draw border\n lines(c(left, left, right, right), c(y+exonHeight\/2, y+exonHeight, y+exonHeight, y+exonHeight\/2), col=getDarkColor(color), lend=2)\n lines(c(left, left, right, right), c(y-exonHeight\/2, y-exonHeight, y-exonHeight, y-exonHeight\/2), col=getDarkColor(color), lend=2)\n # draw gradients for 3D effect\n drawVerticalGradient(rep(left, gradientSteps), rep(right, gradientSteps), seq(y+0.03, y+0.015, len=gradientSteps), rgb(0,0,0,0.2))\n drawVerticalGradient(rep(left, gradientSteps), rep(right, gradientSteps), seq(y-0.03, y-0.015, len=gradientSteps), rgb(0,0,0,0.3))\n } else if (type == "exon") {\n rect(left, y+exonHeight\/2, right, y-exonHeight\/2, col=color, border=getDarkColor(color))\n # draw gradients for 3D effect\n drawVerticalGradient(rep(left, gradientSteps), rep(right, gradientSteps), seq(y, y+exonHeight\/2, len=gradientSteps), rgb(1,1,1,0.6))\n drawVerticalGradient(rep(left, gradientSteps), rep(right, gradientSteps), seq(y, y-exonHeight\/2, len=gradientSteps), rgb(1,1,1,0.6))\n # add exon label\n text((left+right)\/2, y, title, cex=0.9*fontSize)\n }\n}\n\ndrawCircos <- function(fusion, fusions, cytobands, minConfidenceForCircosPlot, circosColors) {\n # check if Giemsa staining information is available\n for (contig in unlist(fusions[fusion,c("contig1", "contig2")])) {\n if (!any(cytobands$contig==contig)) {\n warning(paste0("Circos plot cannot be drawn, because no Giemsa staining information is available for contig ", contig, "."))\n # draw empty plots as placeholder\n plot(0, 0, type="l", xlim=c(0, 1), ylim=c(0, 1), bty="n", xaxt="n", yaxt="n", xlab="", ylab="")\n plot(0, 0, type="l", xlim=c(0, 1), ylim=c(0, 1), bty="n", xaxt="n", yaxt="n", xlab="", ylab="")\n return(NULL)\n }\n }\n # initialize with empty circos plot\n circos.clear()\n circos.initializeWithIdeogram(cytoband=cytobands, plotType=NULL)\n # use gene names as labels or <contig>:<position> for intergenic breakpoints\n geneLabels <- data.frame(\n contig=c(fusions[fusion,"contig1"], fusions[fusion,"contig2"]),\n start=c(fusions[fusion,"breakpoint1"], fusions[fusion,"breakpoint2"])\n )\n geneLabels$end <- geneLabels$start + 1\n geneLabels$gene <- c(fusions[fusion,"gene1"], fusions[fusion,"gene2"])\n geneLabels$gene <- ifelse(c(fusions[fusion,"site1"], fusions[fusion,"site2"]) == "intergenic", paste0(c(fusions[fusion,"display_contig1"], fusions[fusion,"display_contig2"]), ":", geneLabels$start), geneLabels$gene)\n # draw gene labels\n circos.genomicLabels(geneLabels, labels.column=4, side="outside", cex=fontSize, labels_height=0.27)\n # draw chromosome labels in connector plot\n for (contig in unique(cytobands$contig)) {\n set.current.cell(track.index=2, sector.index=contig) # draw in gene label connector track (track.index=2)\n circos.text(CELL_META$xcenter, CELL_META$ycenter, contig, cex=0.85)\n }\n # draw ideograms\n circos.genomicIdeogram(cytoband=cytobands)\n # draw arcs\n confidenceRank <- c(low=0, medium=1, high=2)\n for (i in c(setdiff(1:nrow(fusions), fusion), fusion)) { # draw fusion of interest last, such that its arc is on top\n f <- fusions[i,]\n if (any(cytobands$contig==f$contig1) && any(cytobands$contig==f$contig2)) # ignore viral contigs, because we have no cytoband information for them\n if (minConfidenceForCircosPlot != "none" && confidenceRank[f$confidence] >= confidenceRank[minConfidenceForCircosPlot] || i==fusion)\n circos.link(\n f$contig1, f$breakpoint1,\n f$contig2, f$breakpoint2,\n lwd=2, col=ifelse(i==fusion, circosColors[f$type], getBrightColor(circosColors[f$type]))\n )\n }\n # draw legend\n plot(0, 0, type="l", xlim=c(0, 1), ylim=c(0, 1), bty="n", xaxt="n", yaxt="n", ylab="", xlab="")\n legend(x="top", legend=names(circosColors), col=sapply(circosColors, getBrightColor), lwd=3, ncol=2, box.lty=0)\n}\n\ndrawProteinDomains <- function(fusion, exons1, exons2, proteinDomains, color1, color2, mergeDomainsOverlappingBy, optimizeDomainColors) {\n\n exonHeight <- 0.2\n exonsY <- 0.5\n geneNamesY <- exonsY - exonHeight\/2 - 0.05\n\n # find coding exons\n codingExons1 <- exons1[exons1$type == "CDS" & fusion$site1 != "intergenic",]\n codingExons2 <- exons2[exons2$type == "CDS" & fusion$site2 != "intergenic",]\n\n # cut off coding regions beyond breakpoint\n if (fusion$direction1 == "upstream") {\n codingExons1 <- codingExons1[codingExons1$end >= fusion$breakpoint1,]\n codingExons1$start <- ifelse(codingExons1$start < fusion$breakpoint1, fusion$breakpoint1, codingExons1$start)\n } else {\n codingExons1 <- codingExons1[codingExons1$start <= fusion$breakpoint1,]\n codingExons1$end <- ifelse(codingExons1$end > fusion$breakpoint1, fusion$breakpoint1, codingExons1$end)\n }\n if (fusion$direction2 == "upstream") {\n codingExons2 <- codingExons2[codingExons2$end >= fusion$breakpoint2,]\n codingExons2$start <- ifelse(codingExons2$start < fusion$breakpoint2, fusion$breakpoint2, codingExons2$start)\n } else {\n codingExons2 <- codingExons2[codingExons2$start <= fusion$breakpoint2,]\n codingExons2$end <- ifelse(codingExons2$end > fusion$breakpoint2, fusion$breakpoint2, codingExons2$end)\n }\n\n # find overlapping domains\n exonsGRanges1 <- GRanges(codingExons1$contig, IRanges(codingExons1$start, codingExons1$end), strand=codingExons1$strand)\n exonsGRanges2 <- GRanges(codingExons2$contig, IRanges(codingExons2$start, codingExons2$end), strand=codingExons2$strand)\n domainsGRanges <- GRanges(proteinDomains$contig, IRanges(proteinDomains$start, proteinDomains$end), strand=proteinDomains$strand)\n domainsGRanges$proteinDomainName <- proteinDomains$proteinDomainName\n domainsGRanges$proteinDomainID <- proteinDomains$proteinDomainID\n domainsGRanges$color <- proteinDomains$color\n domainsGRanges <- domainsGRanges[suppressWarnings(unique(queryHits(findOverlaps(domainsGRanges, union(exonsGRanges1, exonsGRanges2)))))]\n\n # group overlapping domains by domain ID\n domainsGRangesList <- GRangesList(lapply(unique(domainsGRanges$proteinDomainID), function(x) { domainsGRanges[domainsGRanges$proteinDomainID == x] }))\n\n # trim protein domains to exon boundaries\n trimDomains <- function(domainsGRangesList, exonsGRanges) {\n do.call(\n "rbind",\n lapply(\n domainsGRangesList,\n function(x) {\n intersected <- as.data.frame(reduce(suppressWarnings(intersect(x, exonsGRanges))))\n if (nrow(intersected) > 0) {\n intersected$proteinDomainName <- head(x$proteinDomainName, 1)\n intersected$proteinDomainID <- head(x$proteinDomainID, 1)\n intersected$color <- head(x$color, 1)\n } else {\n intersected$proteinDomainName <- character()\n intersected$proteinDomainID <- character()\n intersected$color <- character()\n }\n return(intersected)\n }\n )\n )\n }\n retainedDomains1 <- trimDomains(domainsGRangesList, exonsGRanges1)\n retainedDomains2 <- trimDomains(domainsGRangesList, exonsGRanges2)\n\n # calculate length of coding exons\n codingExons1$length <- codingExons1$end - codingExons1$start + 1\n codingExons2$length <- codingExons2$end - codingExons2$start + 1\n\n # abort, if there are no coding regions\n if (sum(exons1$type == "CDS") + sum(exons2$type == "CDS") == 0) {\n text(0.5, 0.5, "Genes are not protein-coding.")\n return(NULL)\n }\n codingLength1 <- sum(codingExons1$length)\n codingLength2 <- sum(codingExons2$length)\n if (codingLength1 + codingLength2 == 0) {\n text(0.5, 0.5, "No coding regions retained in fusion transcript.")\n return(NULL)\n }\n if ((codingLength1 == 0 || grepl("\\\\.$", fusion$strand1)) && (codingLength2 == 0 || grepl("\\\\.$", fusion$strand2))) {\n text(0.5, 0.5, "Failed to determine retained protein domains due to lack of strand information.")\n return(NULL)\n }\n antisenseTranscription1 <- sub("\/.*", "", fusion$strand1) != sub(".*\/", "", fusion$strand1)\n antisenseTranscription2 <- sub("\/.*", "", fusion$strand2) != sub(".*\/", "", fusion$strand2)\n if ((codingLength1 == 0 || antisenseTranscription1) && (codingLength2 == 0 || antisenseTranscription2)) {\n text(0.5, 0.5, "No coding regions due to antisense transcription.")\n return(NULL)\n }\n\n # remove introns from protein domains\n removeIntronsFromProteinDomains <- function(codingExons, retainedDomains) {\n if (nrow(codingExons) == 0) return(NULL)\n cumulativeIntronLength <- 0\n previousExonEnd <- 0\n for (exon in 1:nrow(codingExons)) {\n if (codingExons[exon,"start"] > previousExonEnd)\n cumulativeIntronLength <- cumulativeIntronLength + codingExons[exon,"start"] - previousExonEnd\n domainsInExon <- which(between(retainedDomains$start, codingExons[exon,"start"], codingExons[exon,"end"]))\n retainedDomains[domainsInExon,"start"] <- retainedDomains[domainsInExon,"start"] - cumulativeIntronLength\n domainsInExon <- which(between(retainedDomains$end, codingExons[exon,"start"], codingExons[exon,"end"]))\n retainedDomains[domainsInExon,"end"] <- retainedDomains[domainsInExon,"end"] - cumulativeIntronLength\n previousExonEnd <- codingExons[exon,"end"]\n }\n # merge adjacent domains\n retainedDomains <- do.call(\n "rbind",\n lapply(\n unique(retainedDomains$proteinDomainID),\n function(x) {\n domain <- retainedDomains[retainedDomains$proteinDomainID == x,]\n merged <- reduce(GRanges(domain$seqnames, IRanges(domain$start, domain$end), strand=domain$strand))\n merged$proteinDomainName <- head(domain$proteinDomainName, 1)\n merged$proteinDomainID <- head(domain$proteinDomainID, 1)\n merged$color <- head(domain$color, 1)\n return(as.data.frame(merged))\n }\n )\n )\n return(retainedDomains)\n }\n retainedDomains1 <- removeIntronsFromProteinDomains(codingExons1, retainedDomains1)\n retainedDomains2 <- removeIntronsFromProteinDomains(codingExons2, retainedDomains2)\n\n # abort, if no domains are retained\n if (is.null(retainedDomains1) && is.null(retainedDomains2)) {\n text(0.5, 0.5, "No protein domains retained in fusion.")\n return(NULL)\n }\n\n # merge domains with similar coordinates\n mergeSimilarDomains <- function(domains, mergeDomainsOverlappingBy) {\n if (is.null(domains)) return(domains)\n merged <- domains[F,] # create empty data frame\n domains <- domains[order(domains$end - domains$start, decreasing=T),] # start with bigger domains => bigger domains are retained\n for (domain in rownames(domains)) {\n if (!any((abs(merged$start - domains[domain,"start"]) + abs(merged$end - domains[domain,"end"])) \/ (domains[domain,"end"] - domains[domain,"start"] + 1) <= 1-mergeDomainsOverlappingBy))\n merged <- rbind(merged, domains[domain,])\n }\n return(merged)\n }\n retainedDomains1 <- mergeSimilarDomains(retainedDomains1, mergeDomainsOverlappingBy)\n retainedDomains2 <- mergeSimilarDomains(retainedDomains2, mergeDomainsOverlappingBy)\n\n # if desired, reassign colors to protein domains to maximize contrast\n if (optimizeDomainColors) {\n uniqueDomains <- unique(c(retainedDomains1$proteinDomainID, retainedDomains2$proteinDomainID))\n # make rainbow of pretty pastell colors\n colors <- rainbow(length(uniqueDomains))\n colors <- apply(col2rgb(colors), 2, function(x) { 0.3 + y\/255 * 0.7 }) # make pastell colors\n colors <- apply(colors, 2, function(x) {rgb(x["red"], x["green"], x["blue"])}) # convert back to rgb\n # reassign colors\n names(colors) <- uniqueDomains\n retainedDomains1$color <- colors[retainedDomains1$proteinDomainID]\n retainedDomains2$color <- colors[retainedDomains2$proteinDomainID]\n }\n\n # reverse exons and protein domains, if on the reverse strand\n if (any(codingExons1$strand == "-")) {\n codingExons1$length <- rev(codingExons1$length)\n temp <- retainedDomains1$end\n retainedDomains1$end <- codingLength1 - retainedDomains1$start\n retainedDomains1$start <- codingLength1 - temp\n }\n if (any(codingExons2$strand == "-")) {\n codingExons2$length <- rev(codingExons2$length)\n temp <- retainedDomains2$end\n retainedDomains2$end <- codingLength2 - retainedDomains2$start\n retainedDomains2$start <- codingLength2 - temp\n }\n\n # normalize length to 1\n codingExons1$length <- codingExons1$length \/ (codingLength1 + codingLength2)\n codingExons2$length <- codingExons2$length \/ (codingLength1 + codingLength2)\n retainedDomains1$start <- retainedDomains1$start \/ (codingLength1 + codingLength2)\n retainedDomains1$end <- retainedDomains1$end \/ (codingLength1 + codingLength2)\n retainedDomains2$start <- retainedDomains2$start \/ (codingLength1 + codingLength2)\n retainedDomains2$end <- retainedDomains2$end \/ (codingLength1 + codingLength2)\n\n # draw coding regions\n rect(0, exonsY-exonHeight\/2, sum(codingExons1$length), exonsY+exonHeight\/2, col=color1, border=NA)\n rect(sum(codingExons1$length), exonsY-exonHeight\/2, sum(codingExons1$length) + sum(codingExons2$length), exonsY+exonHeight\/2, col=color2, border=NA)\n\n # indicate exon boundaries as dotted lines\n exonBoundaries <- cumsum(c(codingExons1$length, codingExons2$length))\n if (length(exonBoundaries) > 1) {\n exonBoundaries <- exonBoundaries[1:(length(exonBoundaries)-1)]\n for (exonBoundary in exonBoundaries)\n lines(c(exonBoundary, exonBoundary), c(exonsY-exonHeight, exonsY+exonHeight), col="white", lty=3)\n }\n\n # find overlapping domains\n # nest if one is contained in another\n # stack if they overlap partially\n nestDomains <- function(domains) {\n if (length(unlist(domains)) == 0) return(domains)\n domains <- domains[order(domains$end - domains$start, decreasing=T),]\n rownames(domains) <- 1:nrow(domains)\n # find nested domains and make tree structure\n domains$parent <- 0\n for (domain in rownames(domains))\n domains[domains$start >= domains[domain,"start"] & domains$end <= domains[domain,"end"] & rownames(domains) != domain,"parent"] <- domain\n # find partially overlapping domains\n maxOverlappingDomains <- max(1, as.integer(coverage(IRanges(domains$start*10e6, domains$end*10e6))))\n padding <- 1 \/ maxOverlappingDomains * 0.4\n domains$y <- 0\n domains$height <- 0\n adjustPositionAndHeight <- function(parentDomain, y, height, padding, e) {\n for (domain in which(e$domains$parent == parentDomain)) {\n overlappingDomains <- which((between(e$domains$start, e$domains[domain,"start"], e$domains[domain,"end"]) |\n between(e$domains$end , e$domains[domain,"start"], e$domains[domain,"end"])) &\n e$domains$parent == parentDomain)\n e$domains[domain,"height"] <- height\/length(overlappingDomains) - padding * (length(overlappingDomains)-1) \/ length(overlappingDomains)\n e$domains[domain,"y"] <- y + (which(domain==overlappingDomains)-1) * (e$domains[domain,"height"] + padding)\n adjustPositionAndHeight(domain, e$domains[domain,"y"]+padding, e$domains[domain,"height"]-2*padding, padding, e)\n }\n }\n adjustPositionAndHeight(0, 0, 1, padding, environment())\n domains <- domains[order(domains$height, decreasing=T),] # draw nested domains last\n return(domains)\n }\n retainedDomains1 <- nestDomains(retainedDomains1)\n retainedDomains2 <- nestDomains(retainedDomains2)\n retainedDomains1$y <- exonsY - exonHeight\/2 + 0.025 + (exonHeight-2*0.025) * retainedDomains1$y\n retainedDomains2$y <- exonsY - exonHeight\/2 + 0.025 + (exonHeight-2*0.025) * retainedDomains2$y\n retainedDomains1$height <- retainedDomains1$height * (exonHeight-2*0.025)\n retainedDomains2$height <- retainedDomains2$height * (exonHeight-2*0.025)\n\n # draw domains\n drawProteinDomainRect <- function(left, bottom, right, top, color) {\n rect(left, bottom, right, top, col=color, border=getDarkColor(color))\n # draw gradients for 3D effect\n gradientSteps <- 20\n drawVerticalGradient(rep(left, gradientSteps), rep(right, gradientSteps), seq(top, bottom, len=gradientSteps), rgb(1,1,1,0.7))\n drawVerticalGradient(rep(left, gradientSteps), rep(right, gradientSteps), seq(bottom, bottom+(top-bottom)*0.4, len=gradientSteps), rgb(0,0,0,0.1))\n }\n if (length(unlist(retainedDomains1)) > 0)\n for (domain in 1:nrow(retainedDomains1))\n drawProteinDomainRect(retainedDomains1[domain,"start"], retainedDomains1[domain,"y"], retainedDomains1[domain,"end"], retainedDomains1[domain,"y"]+retainedDomains1[domain,"height"], retainedDomains1[domain,"color"])\n if (length(unlist(retainedDomains2)) > 0)\n for (domain in 1:nrow(retainedDomains2))\n drawProteinDomainRect(sum(codingExons1$length)+retainedDomains2[domain,"start"], retainedDomains2[domain,"y"], sum(codingExons1$length)+retainedDomains2[domain,"end"], retainedDomains2[domain,"y"]+retainedDomains2[domain,"height"], retainedDomains2[domain,"color"])\n\n # draw gene names, if there are coding exons\n if (codingLength1 > 0)\n text(sum(codingExons1$length)\/2, geneNamesY, fusion$gene1, font=2, cex=fontSize)\n if (codingLength2 > 0)\n text(sum(codingExons1$length)+sum(codingExons2$length)\/2, geneNamesY, fusion$gene2, font=2, cex=fontSize)\n\n # calculate how many non-adjacent unique domains there are\n # we need this info to know where to place labels vertically\n countUniqueDomains <- function(domains) {\n uniqueDomains <- 0\n if (length(unlist(domains)) > 0) {\n uniqueDomains <- 1\n if (nrow(domains) > 1) {\n previousDomain <- domains[1,"proteinDomainID"]\n for (domain in 2:nrow(domains)) {\n if (previousDomain != domains[domain,"proteinDomainID"])\n uniqueDomains <- uniqueDomains + 1\n previousDomain <- domains[domain,"proteinDomainID"]\n }\n }\n }\n return(uniqueDomains)\n }\n if (length(unlist(retainedDomains1)) > 0)\n retainedDomains1 <- retainedDomains1[order(retainedDomains1$start),]\n uniqueDomains1 <- countUniqueDomains(retainedDomains1)\n if (length(unlist(retainedDomains2)) > 0)\n retainedDomains2 <- retainedDomains2[order(retainedDomains2$end, decreasing=T),]\n uniqueDomains2 <- countUniqueDomains(retainedDomains2)\n\n # draw title of plot\n titleY <- exonsY + exonHeight\/2 + (uniqueDomains1 + 2) * 0.05\n text(0.5, titleY+0.01, "RETAINED PROTEIN DOMAINS", adj=c(0.5, 0), font=2, cex=fontSize)\n text(0.5, titleY, ifelse(fusion$reading_frame %in% c("in-frame", "out-of-frame"), paste(fusion$reading_frame, "fusion"), ifelse(fusion$reading_frame == "stop-codon", "stop codon before fusion junction", "reading frame unclear")), adj=c(0.5, 1), cex=fontSize)\n\n # draw domain labels for gene1\n if (length(unlist(retainedDomains1)) > 0) {\n previousConnectorX <- -1\n previousLabelX <- -1\n labelY <- exonsY + exonHeight\/2 + uniqueDomains1 * 0.05\n for (domain in 1:nrow(retainedDomains1)) {\n # if possible avoid overlapping lines of labels\n connectorX <- min(retainedDomains1[domain,"start"] + 0.01, (retainedDomains1[domain,"start"] + retainedDomains1[domain,"end"])\/2)\n if (connectorX - previousConnectorX < 0.01 && retainedDomains1[domain,"end"] > previousConnectorX + 0.01)\n connectorX <- previousConnectorX + 0.01\n labelX <- max(connectorX, previousLabelX) + 0.02\n # use a signle label for adjacent domains of same type\n adjacentDomainsOfSameType <- domain + 1 <= nrow(retainedDomains1) && retainedDomains1[domain+1,"proteinDomainID"] == retainedDomains1[domain,"proteinDomainID"]\n if (adjacentDomainsOfSameType) {\n labelX <- retainedDomains1[domain+1,"start"] + 0.015\n } else {\n text(labelX, labelY, retainedDomains1[domain,"proteinDomainName"], adj=c(0,0.5), col=getDarkColor(retainedDomains1[domain,"color"]), cex=fontSize)\n }\n lines(c(labelX-0.005, connectorX, connectorX), c(labelY, labelY, retainedDomains1[domain,"y"]+retainedDomains1[domain,"height"]), col=getDarkColor(retainedDomains1[domain,"color"]))\n if (!adjacentDomainsOfSameType)\n labelY <- labelY - 0.05\n previousConnectorX <- connectorX\n previousLabelX <- labelX\n }\n }\n\n # draw domain labels for gene2\n if (length(unlist(retainedDomains2)) > 0) {\n previousConnectorX <- 100\n previousLabelX <- 100\n labelY <- exonsY - exonHeight\/2 - (uniqueDomains2+1) * 0.05\n for (domain in 1:nrow(retainedDomains2)) {\n # if possible avoid overlapping connector lines of labels\n connectorX <- sum(codingExons1$length) + max(retainedDomains2[domain,"end"] - 0.01, (retainedDomains2[domain,"start"] + retainedDomains2[domain,"end"])\/2)\n if (previousConnectorX - connectorX < 0.01 && sum(codingExons1$length) + retainedDomains2[domain,"start"] < previousConnectorX - 0.01)\n connectorX <- previousConnectorX - 0.01\n labelX <- min(connectorX, previousLabelX) - 0.02\n # use a signle label for adjacent domains of same type\n adjacentDomainsOfSameType <- domain + 1 <= nrow(retainedDomains2) && retainedDomains2[domain+1,"proteinDomainID"] == retainedDomains2[domain,"proteinDomainID"]\n if (adjacentDomainsOfSameType) {\n labelX <- sum(codingExons1$length) + retainedDomains2[domain+1,"end"] - 0.015\n } else {\n text(labelX, labelY, retainedDomains2[domain,"proteinDomainName"], adj=c(1,0.5), col=getDarkColor(retainedDomains2[domain,"color"]), cex=fontSize)\n }\n lines(c(labelX+0.005, connectorX, connectorX), c(labelY, labelY, retainedDomains2[domain,"y"]), col=getDarkColor(retainedDomains2[domain,"color"]))\n if (!adjacentDomainsOfSameType)\n labelY <- labelY + 0.05\n previousConnectorX <- connectorX\n previousLabelX <- labelX\n }\n }\n\n}\n\nfindExons <- function(exons, contig, geneID, direction, breakpoint, coverage, transcriptId, transcriptSelection) {\n # use the provided transcript if desired\n if (transcriptSelection == "provided" && transcriptId != "." && transcriptId != "") {\n candidateExons <- exons[exons$transcript == transcriptId,]\n if (nrow(candidateExons) == 0) {\n warning(paste0("Unknown transcript given in fusions file (", transcriptId, "), selecting a different one"))\n } else {\n return(candidateExons)\n }\n }\n\n if (transcriptSelection == "canonical") {\n candidateExons <- exons[exons$geneID == geneID & exons$contig == contig,]\n } else {\n # look for exon with breakpoint as splice site\n transcripts <- exons[exons$geneID == geneID & exons$contig == contig & exons$type == "exon" & (direction == "downstream" & abs(exons$end - breakpoint) <= 2 | direction == "upstream" & abs(exons$start - breakpoint) <= 2),"transcript"]\n candidateExons <- exons[exons$transcript %in% transcripts,]\n # if none was found, use all exons of the gene closest to the breakpoint\n if (nrow(candidateExons) == 0)\n candidateExons <- exons[exons$geneID == geneID & exons$contig == contig,]\n # if we have coverage information, use the transcript with the highest coverage if there are multiple hits\n if (!is.null(coverage)) {\n highestCoverage <- -1\n transcriptWithHighestCoverage <- NULL\n lengthOfTranscriptWithHighestCoverage <- 0\n for (transcript in unique(candidateExons$transcript)) {\n exonsOfTranscript <- candidateExons[candidateExons$transcript==transcript,]\n exonsOfTranscript$start <- sapply(exonsOfTranscript$start, max, min(start(coverage)))\n exonsOfTranscript$end <- sapply(exonsOfTranscript$end, min, max(end(coverage)))\n lengthOfTranscript <- sum(exonsOfTranscript$end - exonsOfTranscript$start + 1)\n coverageSum <- sum(as.numeric(coverage[IRanges(exonsOfTranscript$start, exonsOfTranscript$end)]))\n # we prefer shorter transcripts over longer ones, because otherwise there is a bias towards transcripts with long UTRs\n # => a longer transcript must have substantially higher coverage to replace a shorter one\n substantialDifference <- (1 - min(lengthOfTranscript, lengthOfTranscriptWithHighestCoverage) \/ max(lengthOfTranscript, lengthOfTranscriptWithHighestCoverage)) \/ 10\n if (lengthOfTranscript > lengthOfTranscriptWithHighestCoverage && coverageSum * (1-substantialDifference) > highestCoverage ||\n lengthOfTranscript < lengthOfTranscriptWithHighestCoverage && coverageSum > highestCoverage * (1-substantialDifference)) {\n highestCoverage <- coverageSum\n transcriptWithHighestCoverage <- transcript\n lengthOfTranscriptWithHighestCoverage <- lengthOfTranscript\n }\n }\n if (highestCoverage > 0)\n candidateExons <- candidateExons[candidateExons$transcript==transcriptWithHighestCoverage,]\n }\n # if the gene has multiple transcripts, search for transcripts which encompass the breakpoint\n if (length(unique(candidateExons$transcript)) > 1) {\n transcriptStart <- aggregate(candidateExons$start, by=list(candidateExons$transcript), min)\n rownames(transcriptStart) <- transcriptStart[,1]\n transcriptEnd <- aggregate(candidateExons$end, by=list(candidateExons$transcript), max)\n rownames(transcriptEnd) <- transcriptEnd[,1]\n encompassingExons <- between(breakpoint, transcriptStart[candidateExons$transcript,2], transcriptEnd[candidateExons$transcript,2])\n if (any(encompassingExons))\n candidateExons <- candidateExons[encompassingExons,]\n }\n }\n\n # find the consensus transcript, if there are multiple hits\n if (length(unique(candidateExons$transcript)) > 1) {\n consensusTranscript <-\n ifelse(grepl("appris_principal_1", candidateExons$attributes), 12,\n ifelse(grepl("appris_principal_2", candidateExons$attributes), 11,\n ifelse(grepl("appris_principal_3", candidateExons$attributes), 10,\n ifelse(grepl("appris_principal_4", candidateExons$attributes), 9,\n ifelse(grepl("appris_principal_5", candidateExons$attributes), 8,\n ifelse(grepl("appris_principal", candidateExons$attributes), 7,\n ifelse(grepl("appris_candidate_longest", candidateExons$attributes), 6,\n ifelse(grepl("appris_candidate", candidateExons$attributes), 5,\n ifelse(grepl("appris_alternative_1", candidateExons$attributes), 4,\n ifelse(grepl("appris_alternative_2", candidateExons$attributes), 3,\n ifelse(grepl("appris_alternative", candidateExons$attributes), 2,\n ifelse(grepl("CCDS", candidateExons$attributes), 1,\n 0\n ))))))))))))\n candidateExons <- candidateExons[consensusTranscript == max(consensusTranscript),]\n }\n # use the transcript with the longest coding sequence, if there are still multiple hits\n if (length(unique(candidateExons$transcript)) > 1) {\n codingSequenceLength <- ifelse(candidateExons$type == "CDS", candidateExons$end - candidateExons$start, 0)\n totalCodingSequenceLength <- aggregate(codingSequenceLength, by=list(candidateExons$transcript), sum)\n rownames(totalCodingSequenceLength) <- totalCodingSequenceLength[,1]\n candidateExons <- candidateExons[totalCodingSequenceLength[candidateExons$transcript,2] == max(totalCodingSequenceLength[,2]),]\n }\n # use the transcript with the longest overall sequence, if there are still multiple hits\n if (length(unique(candidateExons$transcript)) > 1) {\n exonLength <- candidateExons$end - candidateExons$start\n totalExonLength <- aggregate(exonLength, by=list(candidateExons$transcript), sum)\n rownames(totalExonLength) <- totalExonLength[,1]\n candidateExons <- candidateExons[totalExonLength[candidateExons$transcript,2] == max(totalExonLength[,2]),]\n }\n # if there are still multiple hits, select the first one\n candidateExons <- unique(candidateExons[candidateExons$transcript == head(unique(candidateExons$transcript), 1),])\n return(candidateExons)\n}\n\nfindClosestGene <- function(exons, contig, breakpoint, extraConditions) {\n\n # find exons near breakpoint (extraConditions must define what is considered "near")\n closestExons <- exons[exons$contig == contig & extraConditions,] # find closest exon\n closestExons <- exons[exons$contig == contig & exons$geneID %in% closestExons$geneID,] # select all exons of closest gene\n\n # when more than one gene found with the given name, use the closest one\n if (length(unique(closestExons$geneID)) > 1) { # more than one gene found with the given name => use the closest one\n distanceToBreakpoint <- aggregate(1:nrow(closestExons), by=list(closestExons$geneID), function(x) { min(abs(closestExons[x,"start"]-breakpoint), abs(closestExons[x,"end"]-breakpoint)) })\n closestGene <- head(distanceToBreakpoint[distanceToBreakpoint[,2] == min(distanceToBreakpoint[,2]),1], 1)\n closestExons <- closestExons[closestExons$geneID == closestGene,]\n }\n\n # when no gene was found, return default values\n if (nrow(closestExons) == 0) {\n return(IRanges(max(1, breakpoint-1000), breakpoint+1000))\n } else {\n return(IRanges(min(closestExons$start), max(closestExons$end)))\n }\n}\n\n# main loop starts here\nfor (fusion in 1:nrow(fusions)) {\n\n message(paste0("Drawing fusion #", fusion, ": ", fusions[fusion,"gene1"], ":", fusions[fusion,"gene2"]))\n\n # if showIntergenicVicinity is a number, take it as is\n # if it is a keyword (closestGene\/closestProteinCodingGene), determine the range dynamically\n showVicinity <- rep(0, 4)\n if (fusions[fusion,"site1"] == "intergenic") {\n showVicinity[1] <- ifelse(\n is.numeric(showIntergenicVicinity[[1]]),\n showIntergenicVicinity[[1]],\n fusions[fusion,"breakpoint1"] - start(findClosestGene(exons, fusions[fusion,"contig1"], fusions[fusion,"breakpoint1"], exons$end < fusions[fusion,"breakpoint1"] & exons$type == showIntergenicVicinity[[1]]))\n )\n showVicinity[2] <- ifelse(\n is.numeric(showIntergenicVicinity[[2]]),\n showIntergenicVicinity[[2]],\n end(findClosestGene(exons, fusions[fusion,"contig1"], fusions[fusion,"breakpoint1"], exons$start > fusions[fusion,"breakpoint1"] & exons$type == showIntergenicVicinity[[2]])) - fusions[fusion,"breakpoint1"]\n )\n }\n if (fusions[fusion,"site2"] == "intergenic") {\n showVicinity[3] <- ifelse(\n is.numeric(showIntergenicVicinity[[3]]),\n showIntergenicVicinity[[3]],\n fusions[fusion,"breakpoint2"] - start(findClosestGene(exons, fusions[fusion,"contig2"], fusions[fusion,"breakpoint2"], exons$end < fusions[fusion,"breakpoint2"] & exons$type == showIntergenicVicinity[[3]]))\n )\n showVicinity[4] <- ifelse(\n is.numeric(showIntergenicVicinity[[4]]),\n showIntergenicVicinity[[4]],\n end(findClosestGene(exons, fusions[fusion,"contig2"], fusions[fusion,"breakpoint2"], exons$start > fusions[fusion,"breakpoint2"] & exons$type == showIntergenicVicinity[[4]])) - fusions[fusion,"breakpoint2"]\n )\n }\n\n # compute coverage from alignments file\n coverage1 <- NULL\n coverage2 <- NULL\n if (alignmentsFile != "") {\n # determine range in which we need to compute the coverage\n determineCoverageRegion <- function(exons, geneID, contig, breakpoint, showVicinityLeft, showVicinityRight) {\n closestGene <- findClosestGene(exons, contig, breakpoint, exons$geneID == geneID)\n return(IRanges(min(start(closestGene), breakpoint-showVicinityLeft), max(end(closestGene), breakpoint+showVicinityRight)))\n }\n coverageRegion1 <- determineCoverageRegion(exons, fusions[fusion,"gene_id1"], fusions[fusion,"contig1"], fusions[fusion,"breakpoint1"], showVicinity[1], showVicinity[2])\n coverageRegion2 <- determineCoverageRegion(exons, fusions[fusion,"gene_id2"], fusions[fusion,"contig2"], fusions[fusion,"breakpoint2"], showVicinity[3], showVicinity[4])\n # function which reads alignments from BAM file with & without "chr" prefix\n readCoverage <- function(alignmentsFile, contig, coverageRegion) {\n coverageData <- tryCatch(\n {\n alignments <- readGAlignments(alignmentsFile, param=ScanBamParam(which=GRanges(contig, coverageRegion)))\n coverage(alignments)[[contig]]\n },\n error=function(e) {\n alignments <- readGAlignments(alignmentsFile, param=ScanBamParam(which=GRanges(addChr(contig), coverageRegion)))\n coverage(alignments)[[addChr(contig)]]\n }\n )\n if (exists("alignments")) rm(alignments)\n return(coverageData)\n }\n # get coverage track\n coverage1 <- readCoverage(alignmentsFile, fusions[fusion,"contig1"], coverageRegion1)\n coverage2 <- readCoverage(alignmentsFile, fusions[fusion,"contig2"], coverageRegion2)\n # shrink coverage range to chromosome boundaries to avoid subscript out of bounds errors\n coverageRegion1 <- IRanges(max(start(coverageRegion1), min(start(coverage1))), min(end(coverageRegion1), max(end(coverage1))))\n coverageRegion2 <- IRanges(max(start(coverageRegion2), min(start(coverage2))), min(end(coverageRegion2), max(end(coverage2))))\n }\n\n # find all exons belonging to the fused genes\n exons1 <- findExons(exons, fusions[fusion,"contig1"], fusions[fusion,"gene_id1"], fusions[fusion,"direction1"], fusions[fusion,"breakpoint1"], coverage1, fusions[fusion,"transcript_id1"], transcriptSelection)\n if (nrow(exons1) == 0) {\n par(mfrow=c(1,1))\n plot(0, 0, type="l", xaxt="n", yaxt="n", xlab="", ylab="")\n text(0, 0, paste0("exon coordinates of ", fusions[fusion,"gene1"], " not found in\\n", exonsFile))\n warning(paste("exon coordinates of", fusions[fusion,"gene1"], "not found"))\n next\n }\n exons2 <- findExons(exons, fusions[fusion,"contig2"], fusions[fusion,"gene_id2"], fusions[fusion,"direction2"], fusions[fusion,"breakpoint2"], coverage2, fusions[fusion,"transcript_id2"], transcriptSelection)\n if (nrow(exons2) == 0) {\n par(mfrow=c(1,1))\n plot(0, 0, type="l", xaxt="n", yaxt="n", xlab="", ylab="")\n text(0, 0, paste0("exon coordinates of ", fusions[fusion,"gene2"], " not found in\\n", exonsFile))\n warning(paste("exon coordinates of", fusions[fusion,"gene2"], "not found"))\n next\n }\n\n # in case of intergenic breakpoints, show the vicinity\n if (sum(showVicinity) > 0) {\n if (fusions[fusion,"site1"] == "intergenic") {\n for (geneID in unique(exons[exons$contig == fusions[fusion,"contig1"] & exons$exonNumber != "intergenic" &\n (between(exons$end, fusions[fusion,"breakpoint1"]-showVicinity[1], fusions[fusion,"breakpoint1"]+showVicinity[2]) |\n between(exons$start, fusions[fusion,"breakpoint1"]-showVicinity[1], fusions[fusion,"breakpoint1"]+showVicinity[2])),"geneID"]))\n exons1 <- rbind(exons1, findExons(exons, fusions[fusion,"contig1"], geneID, fusions[fusion,"direction1"], fusions[fusion,"breakpoint1"], coverage1, fusions[fusion,"transcript_id1"], transcriptSelection))\n # crop genes that are only partially within user-defined vicinity, because the coverage data is incomplete for those\n exons1 <- exons1[exons1$start >= fusions[fusion,"breakpoint1"]-showVicinity[1] & exons1$end <= fusions[fusion,"breakpoint1"]+showVicinity[2] | exons1$exonNumber == "intergenic",]\n }\n if (fusions[fusion,"site2"] == "intergenic") {\n for (geneID in unique(exons[exons$contig == fusions[fusion,"contig2"] & exons$exonNumber != "intergenic" &\n (between(exons$end, fusions[fusion,"breakpoint2"]-showVicinity[3], fusions[fusion,"breakpoint2"]+showVicinity[4]) |\n between(exons$start, fusions[fusion,"breakpoint2"]-showVicinity[3], fusions[fusion,"breakpoint2"]+showVicinity[4])),"geneID"]))\n exons2 <- rbind(exons2, findExons(exons, fusions[fusion,"contig2"], geneID, fusions[fusion,"direction2"], fusions[fusion,"breakpoint2"], coverage2, fusions[fusion,"transcript_id2"], transcriptSelection))\n # crop genes that are only partially within user-defined vicinity, because the coverage data is incomplete for those\n exons2 <- exons2[exons2$start >= fusions[fusion,"breakpoint2"]-showVicinity[3] & exons2$end <= fusions[fusion,"breakpoint2"]+showVicinity[4] | exons2$exonNumber == "intergenic",]\n }\n }\n\n # normalize coverage\n if (alignmentsFile != "") {\n coverageNormalization <- function(coverage, coverageRegion, exons) {\n max(1, ifelse(\n squishIntrons, # => ignore intronic coverage\n max(as.numeric(coverage[IRanges(sapply(exons$start,max,min(start(coverage))),sapply(exons$end,min,max(end(coverage))))])),\n round(quantile(coverage[coverageRegion], 0.9999)) # ignore coverage spikes from read-attracting regions\n ))\n }\n coverageNormalization1 <- ifelse(head(coverageRange,1) == 0, coverageNormalization(coverage1, coverageRegion1, exons1), head(coverageRange,1))\n coverageNormalization2 <- ifelse(tail(coverageRange,1) == 0, coverageNormalization(coverage2, coverageRegion2, exons2), tail(coverageRange,1))\n if (length(coverageRange) == 1 && coverageRange == 0) { # harmonize scales of gene1 and gene2\n coverageNormalization1 <- max(coverageNormalization1, coverageNormalization2)\n coverageNormalization2 <- max(coverageNormalization1, coverageNormalization2)\n }\n coverage1 <- coverage1\/coverageNormalization1\n coverage2 <- coverage2\/coverageNormalization2\n coverage1[coverage1 > 1] <- 1\n coverage2[coverage2 > 1] <- 1\n }\n\n # sort coding exons last, such that they are drawn over the border of non-coding exons\n exons1 <- exons1[order(exons1$start, -rank(exons1$type)),]\n exons2 <- exons2[order(exons2$start, -rank(exons2$type)),]\n\n # insert dummy exons, if breakpoints are outside the gene (e.g., in UTRs)\n # this avoids plotting artifacts\n breakpoint1 <- fusions[fusion,"breakpoint1"]\n breakpoint2 <- fusions[fusion,"breakpoint2"]\n if (breakpoint1 < min(exons1$start)) {\n exons1 <- rbind(c(exons1[1,"contig"], "dummy", max(1,breakpoint1-1000), max(1,breakpoint1-1000), exons1[1,"strand"], "", "dummy", exons1[1,"geneID"], exons1[1,"transcript"], ""), exons1)\n } else if (breakpoint1 > max(exons1$end)) {\n exons1 <- rbind(exons1, c(exons1[1,"contig"], "dummy", breakpoint1+1000, breakpoint1+1000, exons1[1,"strand"], "", "dummy", exons1[1,"geneID"], exons1[1,"transcript"], ""))\n }\n if (breakpoint2 < min(exons2$start)) {\n exons2 <- rbind(c(exons2[1,"contig"], "dummy", max(1,breakpoint2-1000), max(1,breakpoint2-1000), exons2[1,"strand"], "", "dummy", exons2[1,"geneID"], exons2[1,"transcript"], ""), exons2)\n } else if (breakpoint2 > max(exons2$end)) {\n exons2 <- rbind(exons2, c(exons2[1,"contig"], "dummy", breakpoint2+1000, breakpoint2+1000, exons2[1,"strand"], "", "dummy", exons2[1,"geneID"], exons2[1,"transcript"], ""))\n }\n exons1$start <- as.integer(exons1$start)\n exons1$end <- as.integer(exons1$end)\n exons2$start <- as.integer(exons2$start)\n exons2$end <- as.integer(exons2$end)\n\n exons1$left <- exons1$start\n exons1$right <- exons1$end\n exons2$left <- exons2$start\n exons2$right <- exons2$end\n\n squishedIntronSize <- 200\n if (squishIntrons) {\n # hide introns in gene1\n cumulativeIntronLength <- 0\n previousExonEnd <- -squishedIntronSize\n for (exon in 1:nrow(exons1)) {\n if (breakpoint1 > previousExonEnd+1 && breakpoint1 < exons1[exon,"left"])\n breakpoint1 <- (breakpoint1-previousExonEnd) \/ (exons1[exon,"left"]-previousExonEnd) * squishedIntronSize + previousExonEnd - cumulativeIntronLength\n if (exons1[exon,"left"] > previousExonEnd) {\n cumulativeIntronLength <- cumulativeIntronLength + exons1[exon,"left"] - previousExonEnd - squishedIntronSize\n previousExonEnd <- exons1[exon,"right"]\n }\n if (breakpoint1 >= exons1[exon,"left"] && breakpoint1 <= exons1[exon,"right"]+1)\n breakpoint1 <- breakpoint1 - cumulativeIntronLength\n exons1[exon,"left"] <- exons1[exon,"left"] - cumulativeIntronLength\n exons1[exon,"right"] <- exons1[exon,"right"] - cumulativeIntronLength\n }\n\n # hide introns in gene2\n cumulativeIntronLength <- 0\n previousExonEnd <- -squishedIntronSize\n for (exon in 1:nrow(exons2)) {\n if (breakpoint2 > previousExonEnd+1 && breakpoint2 < exons2[exon,"left"])\n breakpoint2 <- (breakpoint2-previousExonEnd) \/ (exons2[exon,"left"]-previousExonEnd) * squishedIntronSize + previousExonEnd - cumulativeIntronLength\n if (exons2[exon,"left"] > previousExonEnd) {\n cumulativeIntronLength <- cumulativeIntronLength + exons2[exon,"left"] - previousExonEnd - squishedIntronSize\n previousExonEnd <- exons2[exon,"right"]\n }\n if (breakpoint2 >= exons2[exon,"left"] && breakpoint2 <= exons2[exon,"right"]+1)\n breakpoint2 <- breakpoint2 - cumulativeIntronLength\n exons2[exon,"left"] <- exons2[exon,"left"] - cumulativeIntronLength\n exons2[exon,"right"] <- exons2[exon,"right"] - cumulativeIntronLength\n }\n } else { # don't squish introns\n # shift exon coordinates to align the gene to the left border of the plot\n exons1$right <- exons1$right - min(exons1$left)\n breakpoint1 <- breakpoint1 - min(exons1$left)\n exons1$left <- exons1$left - min(exons1$left)\n exons2$right <- exons2$right - min(exons2$left)\n breakpoint2 <- breakpoint2 - min(exons2$left)\n exons2$left <- exons2$left - min(exons2$left)\n }\n\n # scale exon sizes to fit on page\n scalingFactor <- max(exons1$right) + max(exons2$right)\n if (fixedScale > 0) {\n if (fixedScale >= scalingFactor) {\n scalingFactor <- fixedScale\n } else {\n warning(paste("fallback to automatic scaling, because value for --fixedScale is too small to fit transcripts on canvas (increase it to", scalingFactor, "to avoid this)"))\n }\n }\n exons1$left <- exons1$left \/ scalingFactor\n exons1$right <- exons1$right \/ scalingFactor\n exons2$left <- exons2$left \/ scalingFactor\n exons2$right <- exons2$right \/ scalingFactor\n breakpoint1 <- breakpoint1 \/ scalingFactor\n breakpoint2 <- breakpoint2 \/ scalingFactor\n\n # shift gene2 to the right border of the page\n gene2Offset <- 1 + 0.05 - max(exons2$right)\n\n # center fusion horizontally\n fusionOffset1 <- (max(exons1$right)+gene2Offset)\/2 - ifelse(fusions[fusion,"direction1"] == "downstream", breakpoint1, max(exons1$right)-breakpoint1)\n fusionOffset2 <- fusionOffset1 + ifelse(fusions[fusion,"direction1"] == "downstream", breakpoint1, max(exons1$right)-breakpoint1)\n\n # layout: fusion on top, circos plot on bottom left, protein domains on bottom center, statistics on bottom right\n layout(matrix(c(1,1,1,2,4,5,3,4,5), 3, 3, byrow=TRUE), widths=c(1.1, 1.2, 0.7), heights=c(1.55, 1.2, 0.25))\n par(mar=c(0, 0, 0, 0))\n plot(0, 0, type="l", xlim=c(-0.12, 1.12), ylim=c(0.4, 1.1), bty="n", xaxt="n", yaxt="n", xlab="", ylab="")\n\n # vertical coordinates of layers\n ySampleName <- 1.04\n yIdeograms <- ifelse(alignmentsFile != "", 0.94, 0.84)\n yBreakpointLabels <- ifelse(alignmentsFile != "", 0.86, 0.76)\n yCoverage <- 0.72\n yExons <- 0.67\n yGeneNames <- 0.58\n yFusion <- 0.5\n yTranscript <- 0.45\n yScale <- 0.407\n yTrajectoryBreakpointLabels <- yBreakpointLabels - 0.035\n yTrajectoryExonTop <- yExons + 0.03\n yTrajectoryExonBottom <- yExons - 0.055\n yTrajectoryFusion <- yFusion + 0.03\n\n # print sample name (title of page)\n text(0.5, ySampleName, sampleName, font=2, cex=fontSize*1.5, adj=c(0.5,0))\n\n # draw ideograms\n if (!is.null(cytobands)) {\n drawIdeogram("left", min(exons1$left), max(exons1$right), yIdeograms, cytobands, fusions[fusion,"contig1"], fusions[fusion,"breakpoint1"])\n drawIdeogram("right", gene2Offset, gene2Offset+max(exons2$right), yIdeograms, cytobands, fusions[fusion,"contig2"], fusions[fusion,"breakpoint2"])\n }\n\n # draw gene & transcript names\n if (fusions[fusion,"gene1"] != ".")\n text(max(exons1$right)\/2, yGeneNames, fusions[fusion,"gene1"], font=2, cex=fontSize, adj=c(0.5,0))\n if (fusions[fusion,"site1"] != "intergenic")\n text(max(exons1$right)\/2, yGeneNames-0.01, head(exons1$transcript,1), cex=0.9*fontSize, adj=c(0.5,1))\n if (fusions[fusion,"gene2"] != ".")\n text(gene2Offset+max(exons2$right)\/2, yGeneNames, fusions[fusion,"gene2"], font=2, cex=fontSize, adj=c(0.5,0))\n if (fusions[fusion,"site2"] != "intergenic")\n text(gene2Offset+max(exons2$right)\/2, yGeneNames-0.01, head(exons2$transcript,1), cex=0.9*fontSize, adj=c(0.5,1))\n\n # if multiple genes in the vicinity are shown, label them\n if (fusions[fusion,"site1"] == "intergenic")\n for (gene in unique(exons1$geneName)) {\n exonsOfGene <- exons1[exons1$geneName == gene & exons1$type != "dummy",]\n if (any(exonsOfGene$type == "exon"))\n text(mean(c(min(exonsOfGene$left), max(exonsOfGene$right))), yExons-0.04, gene, cex=0.9*fontSize, adj=c(0.5,1))\n }\n if (fusions[fusion,"site2"] == "intergenic")\n for (gene in unique(exons2$geneName)) {\n exonsOfGene <- exons2[exons2$geneName == gene & exons2$type != "dummy",]\n if (any(exonsOfGene$type == "exon"))\n text(gene2Offset+mean(c(min(exonsOfGene$left), max(exonsOfGene$right))), yExons-0.04, gene, cex=0.9*fontSize, adj=c(0.5,1))\n }\n\n # label breakpoints\n text(breakpoint1+0.01, yBreakpointLabels-0.03, paste0("breakpoint1\\n", fusions[fusion,"display_contig1"], ":", fusions[fusion,"breakpoint1"]), adj=c(1,0), cex=fontSize)\n text(gene2Offset+breakpoint2-0.01, yBreakpointLabels-0.03, paste0("breakpoint2\\n", fusions[fusion,"display_contig2"], ":", fusions[fusion,"breakpoint2"]), adj=c(0,0), cex=fontSize)\n\n # draw coverage axis\n if (alignmentsFile != "") {\n # left axis (gene1)\n lines(c(-0.02, -0.01, -0.01, -0.02), c(yCoverage, yCoverage, yCoverage+0.1, yCoverage+0.1))\n text(-0.025, yCoverage, "0", adj=c(1,0.5), cex=0.9*fontSize)\n text(-0.025, yCoverage+0.1, coverageNormalization1, adj=c(1,0.5), cex=0.9*fontSize)\n text(-0.05, yCoverage+0.08, "Coverage", srt=90, cex=0.9*fontSize, adj=c(1,0.5))\n\n # right axis (gene2)\n if (length(coverageRange) == 2) { # separate axes for gene1 and gene2\n rightCoverageAxisX <- gene2Offset+max(exons2$right)\n lines(c(rightCoverageAxisX+0.02, rightCoverageAxisX+0.01, rightCoverageAxisX+0.01, rightCoverageAxisX+0.02), c(yCoverage, yCoverage, yCoverage+0.1, yCoverage+0.1))\n text(rightCoverageAxisX+0.025, yCoverage, "0", adj=c(0,0.5), cex=0.9*fontSize)\n text(rightCoverageAxisX+0.025, yCoverage+0.1, coverageNormalization2, adj=c(0,0.5), cex=0.9*fontSize)\n text(rightCoverageAxisX+0.05, yCoverage+0.08, "Coverage", srt=90, cex=0.9*fontSize, adj=c(1,0.5))\n }\n\n # plot coverage 1\n rect(min(exons1$left), yCoverage, max(exons1$right), yCoverage+0.1, col="#eeeeee", border=NA)\n if (squishIntrons) {\n for (exon in 1:nrow(exons1))\n if (exons1[exon,"type"] != "CDS") # don't draw coverage twice for coding regions\n drawCoverage(exons1[exon,"left"], exons1[exon,"right"], yCoverage, coverage1, exons1[exon,"start"], exons1[exon,"end"], color1)\n } else {\n drawCoverage(min(exons1$left), max(exons1$right), yCoverage, coverage1, min(exons1$start), max(exons1$end), color1)\n }\n\n # plot coverage 2\n rect(gene2Offset+min(exons2$left), yCoverage, gene2Offset+max(exons2$right), yCoverage+0.1, col="#eeeeee", border=NA)\n if (squishIntrons) {\n for (exon in 1:nrow(exons2))\n if (exons2[exon,"type"] != "CDS") # don't draw coverage twice for coding regions\n drawCoverage(gene2Offset+exons2[exon,"left"], gene2Offset+exons2[exon,"right"], yCoverage, coverage2, exons2[exon,"start"], exons2[exon,"end"], color2)\n } else {\n drawCoverage(gene2Offset+min(exons2$left), gene2Offset+max(exons2$right), yCoverage, coverage2, min(exons2$start), max(exons2$end), color2)\n }\n }\n\n # plot gene 1\n lines(c(min(exons1$left), max(exons1$right)), c(yExons, yExons), col=darkColor1)\n for (gene in unique(exons1$geneName))\n drawStrand(min(exons1[exons1$geneName == gene,"left"]), max(exons1[exons1$geneName == gene,"right"]), yExons, darkColor1, head(exons1[exons1$geneName == gene,"strand"],1))\n for (exon in 1:nrow(exons1))\n drawExon(exons1[exon,"left"], exons1[exon,"right"], yExons, color1, exons1[exon,"exonNumber"], exons1[exon,"type"])\n\n # plot gene 2\n lines(c(gene2Offset, gene2Offset+max(exons2$right)), c(yExons, yExons), col=darkColor2)\n for (gene in unique(exons2$geneName))\n drawStrand(gene2Offset+min(exons2[exons2$geneName == gene,"left"]), gene2Offset+max(exons2[exons2$geneName == gene,"right"]), yExons, darkColor2, head(exons2[exons2$geneName == gene,"strand"],1))\n for (exon in 1:nrow(exons2))\n drawExon(gene2Offset+exons2[exon,"left"], gene2Offset+exons2[exon,"right"], yExons, color2, exons2[exon,"exonNumber"], exons2[exon,"type"])\n\n # plot gene1 of fusion\n if (fusions[fusion,"direction1"] == "downstream") {\n # plot strands\n lines(c(fusionOffset1, fusionOffset1+breakpoint1), c(yFusion, yFusion), col=darkColor1)\n for (gene in unique(exons1$geneName)) {\n exonsOfGene <- exons1[exons1$geneName == gene,]\n if (min(exonsOfGene$start) <= fusions[fusion,"breakpoint1"])\n drawStrand(fusionOffset1+min(exonsOfGene$left), fusionOffset1+min(breakpoint1, max(exonsOfGene$right)), yFusion, col=darkColor1, exonsOfGene$strand[1])\n }\n # plot exons\n for (exon in 1:nrow(exons1))\n if (exons1[exon,"start"] <= fusions[fusion,"breakpoint1"])\n drawExon(fusionOffset1+exons1[exon,"left"], fusionOffset1+min(breakpoint1, exons1[exon,"right"]), yFusion, color1, exons1[exon,"exonNumber"], exons1[exon,"type"])\n # plot trajectories\n lines(c(0, 0, fusionOffset1), c(yTrajectoryExonTop, yTrajectoryExonBottom, yTrajectoryFusion), col="red", lty=2)\n lines(c(breakpoint1, breakpoint1, fusionOffset1+breakpoint1), c(yTrajectoryBreakpointLabels, yTrajectoryExonBottom, yTrajectoryFusion), col="red", lty=2)\n } else if (fusions[fusion,"direction1"] == "upstream") {\n # plot strands\n lines(c(fusionOffset1, fusionOffset2), c(yFusion, yFusion), col=darkColor1)\n for (gene in unique(exons1$geneName)) {\n exonsOfGene <- exons1[exons1$geneName == gene,]\n if (max(exonsOfGene$end+1) >= fusions[fusion,"breakpoint1"])\n drawStrand(fusionOffset2-max(exonsOfGene$right)+breakpoint1, min(fusionOffset2, fusionOffset2-min(exonsOfGene$left)+breakpoint1), yFusion, col=darkColor1, chartr("+-", "-+", exonsOfGene$strand[1]))\n }\n # plot exons\n for (exon in 1:nrow(exons1))\n if (exons1[exon,"end"]+1 >= fusions[fusion,"breakpoint1"])\n drawExon(fusionOffset1+max(exons1$right)-exons1[exon,"right"], min(fusionOffset2, fusionOffset1+max(exons1$right)-exons1[exon,"left"]), yFusion, color1, exons1[exon,"exonNumber"], exons1[exon,"type"])\n # plot trajectories\n lines(c(max(exons1$right), max(exons1$right), fusionOffset1), c(yTrajectoryExonTop, yTrajectoryExonBottom, yTrajectoryFusion), col="red", lty=2)\n lines(c(breakpoint1, breakpoint1, fusionOffset1+max(exons1$right)-breakpoint1), c(yTrajectoryBreakpointLabels, yTrajectoryExonBottom, yTrajectoryFusion), col="red", lty=2)\n }\n\n # plot gene2 of fusion\n if (fusions[fusion,"direction2"] == "downstream") {\n # plot strands\n lines(c(fusionOffset2, fusionOffset2+breakpoint2), c(yFusion, yFusion), col=darkColor2)\n for (gene in unique(exons2$geneName)) {\n exonsOfGene <- exons2[exons2$geneName == gene,]\n if (min(exonsOfGene$start) <= fusions[fusion,"breakpoint2"])\n drawStrand(max(fusionOffset2, fusionOffset2+breakpoint2-max(exonsOfGene$right)), fusionOffset2+breakpoint2-min(exonsOfGene$left), yFusion, col=darkColor2, chartr("+-", "-+", exonsOfGene$strand[1]))\n }\n # plot exons\n for (exon in 1:nrow(exons2))\n if (exons2[exon,"start"] <= fusions[fusion,"breakpoint2"])\n drawExon(max(fusionOffset2, fusionOffset2+breakpoint2-exons2[exon,"right"]), fusionOffset2+breakpoint2-exons2[exon,"left"], yFusion, color2, exons2[exon,"exonNumber"], exons2[exon,"type"])\n # plot trajectories\n lines(c(gene2Offset, gene2Offset, fusionOffset2+breakpoint2), c(yTrajectoryExonTop, yTrajectoryExonBottom, yTrajectoryFusion), col="red", lty=2)\n lines(c(gene2Offset+breakpoint2, gene2Offset+breakpoint2, fusionOffset2), c(yTrajectoryBreakpointLabels, yTrajectoryExonBottom, yTrajectoryFusion), col="red", lty=2)\n } else if (fusions[fusion,"direction2"] == "upstream") {\n # plot strands\n lines(c(fusionOffset2, fusionOffset2+max(exons2$right)-breakpoint2), c(yFusion, yFusion), col=darkColor2)\n for (gene in unique(exons2$geneName)) {\n exonsOfGene <- exons2[exons2$geneName == gene,]\n if (max(exonsOfGene$end+1) >= fusions[fusion,"breakpoint2"])\n drawStrand(max(fusionOffset2, fusionOffset2+min(exonsOfGene$left)-breakpoint2), fusionOffset2+max(exonsOfGene$right)-breakpoint2, yFusion, col=darkColor2, exonsOfGene$strand[1])\n }\n # plot exons\n for (exon in 1:nrow(exons2))\n if (exons2[exon,"end"]+1 >= fusions[fusion,"breakpoint2"])\n drawExon(max(fusionOffset2, fusionOffset2+exons2[exon,"left"]-breakpoint2), fusionOffset2+exons2[exon,"right"]-breakpoint2, yFusion, color2, exons2[exon,"exonNumber"], exons2[exon,"type"])\n # plot trajectories\n lines(c(gene2Offset+max(exons2$right), gene2Offset+max(exons2$right), fusionOffset2+max(exons2$right)-breakpoint2), c(yTrajectoryExonTop, yTrajectoryExonBottom, yTrajectoryFusion), col="red", lty=2)\n lines(c(gene2Offset+breakpoint2, gene2Offset+breakpoint2, fusionOffset2), c(yTrajectoryBreakpointLabels, yTrajectoryExonBottom, yTrajectoryFusion), col="red", lty=2)\n }\n\n if (fusions[fusion,"fusion_transcript"] != ".") {\n # print fusion transcript colored by gene of origin\n fusion_transcript1 <- gsub("\\\\|.*", "", fusions[fusion,"fusion_transcript"], perl=T)\n fusion_transcript1 <- substr(fusion_transcript1, max(1, nchar(fusion_transcript1)-30), nchar(fusion_transcript1))\n fusion_transcript2 <- gsub(".*\\\\|", "", fusions[fusion,"fusion_transcript"], perl=T)\n fusion_transcript2 <- substr(fusion_transcript2, 1, min(nchar(fusion_transcript2), 30))\n # check for non-template bases\n non_template_bases <- gsub(".*\\\\|([^|]*)\\\\|.*", "\\\\1", fusions[fusion,"fusion_transcript"], perl=T)\n if (non_template_bases == fusions[fusion,"fusion_transcript"]) # no non-template bases found\n non_template_bases <- ""\n # divide non-template bases half-and-half for centered alignment\n non_template_bases1 <- substr(non_template_bases, 1, floor(nchar(non_template_bases)\/2))\n non_template_bases2 <- substr(non_template_bases, ceiling(nchar(non_template_bases)\/2+0.5), nchar(non_template_bases))\n # transcript 1\n text(fusionOffset2, yTranscript, bquote(.(fusion_transcript1) * phantom(.(non_template_bases1))), col=darkColor1, adj=c(1,0.5), cex=fontSize)\n # transcript 2\n text(fusionOffset2, yTranscript, bquote(phantom(.(non_template_bases2)) * .(fusion_transcript2)), col=darkColor2, adj=c(0,0.5), cex=fontSize)\n # non-template bases\n text(fusionOffset2, yTranscript, non_template_bases1, adj=c(1,0.5), cex=fontSize)\n text(fusionOffset2, yTranscript, non_template_bases2, adj=c(0,0.5), cex=fontSize)\n }\n\n # draw scale\n realScale <- max(exons1$end - exons1$start, exons2$end - exons2$start)\n mapScale <- max(exons1$right - exons1$left, exons2$right - exons2$left)\n # choose scale which is closest to desired scale length\n desiredScaleSize <- 0.2\n realScale <- desiredScaleSize \/ mapScale * realScale\n mapScale <- desiredScaleSize\n realScaleOptimalFit <- signif(realScale, 1) # round to most significant digit\n mapScaleOptimalFit <- realScaleOptimalFit \/ realScale * mapScale\n # draw scale line\n lines(c(1-mapScaleOptimalFit, 1), c(yScale, yScale)) # scale line\n lines(c(1-mapScaleOptimalFit, 1-mapScaleOptimalFit), c(yScale-0.007, yScale+0.007)) # left whisker\n lines(c(1, 1), c(yScale-0.007, yScale+0.007)) # right whisker\n # draw units above scale line\n realScaleThousands <- max(0, min(3, floor(log10(realScaleOptimalFit)\/3)))\n scaleUnits <- c("bp", "kbp", "Mbp", "Gbp")\n scaleLabel <- paste(realScaleOptimalFit\/max(1,1000^realScaleThousands), scaleUnits[realScaleThousands+1])\n text(1-mapScaleOptimalFit\/2, yScale+0.005, scaleLabel, adj=c(0.5,0), cex=fontSize*0.9)\n if (squishIntrons)\n text(1-mapScaleOptimalFit\/2, yScale-0.005, "introns not to scale", adj=c(0.5,1), cex=fontSize*0.9, font=3)\n\n # draw circos plot\n if (is.null(cytobands) || !("circlize" %in% names(sessionInfo()$otherPkgs)) || !("GenomicRanges" %in% names(sessionInfo()$otherPkgs))) {\n plot(0, 0, type="l", xlim=c(0, 1), ylim=c(0, 1), bty="n", xaxt="n", yaxt="n", xlab="", ylab="")\n plot(0, 0, type="l", xlim=c(0, 1), ylim=c(0, 1), bty="n", xaxt="n", yaxt="n", xlab="", ylab="")\n } else {\n par(mar=c(2,4,0,0), xpd=NA)\n drawCircos(fusion, fusions, cytobands, minConfidenceForCircosPlot, circosColors)\n par(mar=c(0,0,0,0), xpd=F)\n }\n\n # draw protein domains\n plot(0, 0, type="l", xlim=c(-0.1, 1.1), ylim=c(0, 1), bty="n", xaxt="n", yaxt="n", xlab="", ylab="")\n par(xpd=NA)\n if (!is.null(proteinDomains) && "GenomicRanges" %in% names(sessionInfo()$otherPkgs))\n drawProteinDomains(fusions[fusion,], exons1, exons2, proteinDomains, color1, color2, mergeDomainsOverlappingBy, optimizeDomainColors)\n par(xpd=F)\n\n # print statistics about supporting alignments\n plot(0, 0, type="l", xlim=c(0, 1), ylim=c(0, 1), bty="n", xaxt="n", yaxt="n", xlab="", ylab="")\n text(0, 0.575, "SUPPORTING READ COUNT", font=2, adj=c(0,0), cex=fontSize)\n if ("split_reads" %in% colnames(fusions)) { # STAR-Fusion reports split reads from both breakpoints combined\n text(0, 0.525, paste0("Split reads = ", fusions[fusion,"split_reads"], "\\n", "Discordant mates = ", fusions[fusion,"discordant_mates"]), adj=c(0,1), cex=fontSize)\n } else { # Arriba reports split reads separately for the two breakpoints\n text(\n 0, 0.525,\n paste0(\n "Split reads at breakpoint1 = ", fusions[fusion,"split_reads1"], "\\n",\n "Split reads at breakpoint2 = ", fusions[fusion,"split_reads2"], "\\n",\n "Discordant mates = ", fusions[fusion,"discordant_mates"]\n ),\n adj=c(0,1), cex=fontSize\n )\n }\n\n}\n\ndevNull <- dev.off()\nmessage("Done")\n<\/code><\/pre>\n\u6784\u5efabash\u811a\u672c\u7528\u4e8e\u6279\u91cf\u8fd0\u884cArriba<\/h3>\ncd arriba_fusion_vis\n\n# \u6ce8\u610fcytobands\u7248\u672c\u8981\u4e0estar-fusion\u7248\u672c\u4e00\u81f4\uff0cannotation\u8981\u4e0estar-fusion\u7684annotation\u7248\u672c\u4e00\u81f4\n# \u6bd4\u5982\u6b64\u5904starfusion\u7528\u7684\u662fGRCh38_gencode_v37_CTAT_lib_Mar012021.plug-n-play\n# \u90a3\u4e48\uff0c\u8fd9\u91cc\u7684annotation\u7528\u7684\u5c31\u8981\u662fgencode.v37.annotation.gtf\n\nfor n in `ls ..\/fusion_hg38\/*star-fusion.fusion_predictions.tsv`;do echo Rscript .\/draw_fusions.R --fusions=$n --output=$(basename "$n" star-fusion.fusion_predictions.tsv)_fusion.pdf --cytobands=\/home\/zhoukaiwen\/software\/anaconda3\/envs\/arriba\/var\/lib\/arriba\/cytobands_hg38_GRCh38_v2.4.0.tsv --annotation=\/home\/zhoukaiwen\/database\/GENCODE\/gencode.v37.annotation.gtf --proteinDomains=\/home\/zhoukaiwen\/software\/anaconda3\/envs\/arriba\/var\/lib\/arriba\/protein_domains_hg38_GRCh38_v2.4.0.gff3>>.\/arriba_vis.sh;done\n\nnohup .\/arriba_vis.sh>arriba_vis.log 2>&1 &<\/code><\/pre>\n","protected":false},"excerpt":{"rendered":"Arriba\u63d0\u4f9b\u7684draw_fusions.R https:\/\/github.com\/suhrig\/arrib…<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"closed","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[1],"tags":[],"class_list":["post-135","post","type-post","status-publish","format-standard","hentry","category-uncategorized"],"_links":{"self":[{"href":"https:\/\/www.kz-hub.tech\/index.php\/wp-json\/wp\/v2\/posts\/135","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.kz-hub.tech\/index.php\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.kz-hub.tech\/index.php\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.kz-hub.tech\/index.php\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/www.kz-hub.tech\/index.php\/wp-json\/wp\/v2\/comments?post=135"}],"version-history":[{"count":2,"href":"https:\/\/www.kz-hub.tech\/index.php\/wp-json\/wp\/v2\/posts\/135\/revisions"}],"predecessor-version":[{"id":137,"href":"https:\/\/www.kz-hub.tech\/index.php\/wp-json\/wp\/v2\/posts\/135\/revisions\/137"}],"wp:attachment":[{"href":"https:\/\/www.kz-hub.tech\/index.php\/wp-json\/wp\/v2\/media?parent=135"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.kz-hub.tech\/index.php\/wp-json\/wp\/v2\/categories?post=135"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.kz-hub.tech\/index.php\/wp-json\/wp\/v2\/tags?post=135"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}