BOL

Log in
Jit

Jit

Enthusiast computational biologist with business dream !!!

Our Sponsors


Download BioinformaticsOnline(BOL) Apps in your chrome browser.


Plot custom gene density with R

Jit
  • Public
By Jit 2243 days ago
library(karyoploteR) pp <- getDefaultPlotParams(plot.type=2) pp$data1outmargin <- 100 pp$data2outmargin <- 100 pp$topmargin <- 450 gff.file <- "http://plasmodb.org/common/downloads/Current_Release/PvivaxP01/gff/data/PlasmoDB-35_PvivaxP01.gff" data.points.colors <- c("#FFBD07", "#00A6ED", "#FF3366", "#8EA604", "#C200FB") header.lines <- readLines(gff.file, n = 30) ## Error in file(con, "r"): cannot open the connection to 'http://plasmodb.org/common/downloads/Current_Release/PvivaxP01/gff/data/PlasmoDB-33_PvivaxP01.gff' #The lines with the standard chromosomes start with "##sequence-region PvP01". #Select them. ll <- header.lines[grepl(header.lines, pattern = "##sequence-region PvP01")] ## Error in eval(expr, envir, enclos): object 'header.lines' not found #split them by space, and create a data.frame gg <- data.frame(do.call(rbind, strsplit(ll, split = " "))) ## Error in strsplit(ll, split = " "): object 'll' not found gg[,3] <- as.numeric(as.character(gg[,3])) ## Error in eval(expr, envir, enclos): object 'gg' not found gg[,4] <- as.numeric(as.character(gg[,4])) ## Error in eval(expr, envir, enclos): object 'gg' not found #and create a GRanges with the information PvP01.genome <- toGRanges(gg[,c(2,3,4)]) ## Error in is(A, "GRanges"): object 'gg' not found PvP01.genome ## Error in eval(expr, envir, enclos): object 'PvP01.genome' not found #kp <- plotKaryotype(genome=PvP01.genome) #kp <- plotKaryotype(genome=PvP01.genome, ideogram.plotter = NULL, plot.type=2, plot.params = pp) kp <- plotKaryotype(genome=PvP01.genome, plot.type=2, plot.params = pp) #kp <- plotKaryotype(genome=PvP01.genome, ideogram.plotter = NULL) #kpAddCytobandsAsLine(kp) features <- import(gff.file) table(features$type) genes <- features[features$type=="gene"] #kp <- plotKaryotype(genome=PvP01.genome, ideogram.plotter = NULL) #kpAddCytobandsAsLine(kp) #kpPlotRegions(kp, data=genes) #kp <- plotKaryotype(genome=PvP01.genome, ideogram.plotter = NULL, plot.type=2) kpAddMainTitle(kp, "Plasmodium Vivax - PvP01 with genes", cex=2) kpPlotRegions(kp, data=genes[strand(genes)=="+"], avoid.overlapping = FALSE, col="deepskyblue") kpPlotRegions(kp, data=genes[strand(genes)=="-"], avoid.overlapping = FALSE, col="gold", data.panel=2) kpAddLabels(kp, "strand +", cex=0.8, col="#888888") kpAddLabels(kp, "strand -", data.panel=2, cex=0.8, col="#888888") #plot all pp <- getDefaultPlotParams(plot.type = 4) pp$data1inmargin <- 0 pp$bottommargin <- 20 kp <- plotKaryotype(genome=PvP01.genome, plot.type=4, ideogram.plotter = NULL, labels.plotter = NULL, plot.params = pp, main="Gene Density") kpAddCytobandsAsLine(kp) kpAddChromosomeNames(kp, srt=45) kpPlotDensity(kp, genes, window.size = 10e2, col="#ddaacc") #plot all pp <- getDefaultPlotParams(plot.type = 2) pp$data1inmargin <- 0 pp$bottommargin <- 20 kp <- plotKaryotype(genome=PvP01.genome, plot.params = pp) #kp <- kpPlotDensity(kp, genes) kpAddChromosomeNames(kp, srt=45) kpPlotDensity(kp, genes, col="#ddaacc", data.panel = 1) kpAbline(kp, h=0.4, data.panel = 2, r0=0.2, r1=0, col=data.points.colors[3])
Powered by Elgg